We calculate the yields of quarkonia in heavy ion collisions at RHIC and the LHC as a function of their transverse momentum. Based upon non-relativistic quantum chromodynamics, our results include both color-singlet and color-octet contributions and feed-down effects from excited states. In reactions with ultra-relativistic nuclei, we focus on the consistent implementation of dynamically calculated nuclear matter effects, such as coherent power corrections, cold nuclear matter energy loss, and the Cronin effect in the initial state. In the final state, we consider radiative energy loss for the color-octet state and collisional dissociation of quarkonia as they traverse through the QGP. Theoretical results are presented for J/ψ and Υ and compared to experimental data where applicable. At RHIC, a good description of the high-pT J/ψ modification observed in central Cu+Cu and Au+Au collisions can be achieved within the model uncertainties. We find that measurements of J/ψ yields in proton-nucleus reactions are needed to constrain the magnitude of cold nuclear matter effects. At the LHC, a good description of the experimental data can be achieved only in mid-central and peripheral Pb+Pb collisions. The large five-fold suppression of prompt J/ψ in the most central nuclear reactions may indicate for the first time possible thermal effects at the level of the quarkonium wavefunction at large transverse momenta.