We study the possibility to probe the unpolarized and linearly polarized transverse momentum- dependent gluon distributions in unpolarized pp collision in charmonium and bottomonium produc- tion, employing non-relativistic QCD (NRQCD) based color octet model within transverse momentum dependent (TMD) factorization framework. The transverse momentum (p T ) and rapidity distributions of J/ψ and Υ(1S) at LHCb, RHIC and AFTER energies are estimated. Significant modulations in transverse momentum spectrum of quarkonium in the low p T region is obtained when contributions of linearly polarized gluons inside an unpolarized proton are taken into account. The results of quarko- nium production in color octet model and color evaporation model are compared.