Quarkonium production offers good possibilities to study gluon TMDs. In this proceedings contribution this topic is explored for the linearly polarized gluons inside unpolarized hadrons and unpolarized gluons inside transversely polarized hadrons. It is argued that χ b0/2 and η b production at LHC are best to study the effects of linearly polarized gluons in hadronic collisions, by means of angular independent ratios of ratios of cross sections. This can be directly compared to cos 2φ asymmetries in heavy quark pair and dijet production in DIS at a future high-energy Electron-Ion Collider (EIC), which probe the same TMDs. In the small-x limit this corresponds to the Weizsäcker-Williams (WW) gluon distributions, which should show a change in behavior for transverse momenta around the saturation scale. Together with investigations of the dipole (DP) gluon distributions, this can provide valuable information about the polarization of the Color Glass Condensate if sufficiently small x-values are reached. Quarkonia can also be useful in the study of single transverse spin asymmetries. For transversely polarized hadrons the gluon distribution can be asymmetric, which is referred to as the Sivers effect. It leads to single spin asymmetries in for instance J/ψ (pair) production at AFTER@LHC, which probe the WW or f -type gluon Sivers TMD. It allows for a test of a sign-change relation w.r.t. the gluon Sivers TMD probed at an EIC in open heavy quark pair production. Single spin asymmetries in backward inclusive C-odd quarkonium production, such as J/ψ production, may offer probes of the DP or d-type gluon Sivers TMD at small x-values in the polarized proton, which in that limit corresponds to a correlator of a single Wilson loop, describing the spin-dependent odderon.