In this paper we analyze azimuthal asymmetries in the processes of unpolarized and polarized J/ψ(Υ) production at an Electron-Ion Collider. Apart from giving access to various unknown gluon transverse momentum distributions, we suggest to use them as a new method to extract specific color-octet NRQCD long-distance matrix elements, i.e. ⟨0|O8J/ψ(1S0)|0⟩ and ⟨0|O8J/ψ(3P0)|0⟩, whose values are still quite uncertain and for which lattice calculations are unavailable. The new method is based on combining measurements of analogous asymmetries in open heavy-quark pair production which can be performed at the same energy. We also study for the first time the effects of transverse-momentum smearing in the quarkonium formation process. To enhance the gluon contribution one can consider smaller values of x and, in order to assess the impact of small-x evolution, we perform a numerical study using the MV model as a starting input and evolve it with the JIMWLK equations.