We consider the transverse momentum spectrum and the cos 2ϕ azimuthal distribution of J/ψ mesons produced in semi-inclusive, deep-inelastic electron-proton scattering, where the electron and the proton are unpolarized. At low transverse momentum, we propose factorized expressions in terms of transverse momentum dependent gluon distributions and shape functions. We show that our formulae, at the order αs, correctly match with the collinear factorization results at high transverse momentum. The latter are computed at the order αs2 in the framework of nonrelativistic QCD (NRQCD), with the inclusion of the intermediate 3S1 color-singlet Fock state, as well as the subleading color-octet ones that are relatively suppressed by a factor v4 in the NRQCD velocity parameter v. We show that the 1S0 and 3PJ (J = 0, 1, 2) contributions diverge in the small transverse momentum region and allow us to determine the perturbative tails of the shape functions, which carry the same quantum numbers. These turn out to be identical, except for the overall magnitude given by the appropriate NRQCD long distance matrix element.