We review the statistical hadronization picture for charmonium production in ultra-relativistic nuclear collisions. Our starting point is a brief reminder of the status of the thermal model description of hadron production at high energy. Within this framework an excellent account is achieved of all data for hadrons built of (u,d,s) valence quarks using temperature, baryo-chemical potential and volume as thermal parameters. The large charm quark mass brings in a new (non-thermal) scale which is explicitely taken into account by fixing the total number of charm quarks produced in the collision. Emphasis is placed on the description of the physical basis for the resulting statistical hadronization model. We discuss the evidence for statistical hadronization of charmonia by analysis of recent data from the SPS and RHIC accelerators. Furthermore we discuss an extension of this model towards lower beam energies and develop arguments about the prospects to observe medium modifications of open and hidden charm hadrons. With the imminent start of the LHC accelerator at CERN, exciting prospects for charmonium production studies at the very high energy frontier come into reach. We present arguments that, at such energies, charmonium production becomes a fingerprint of deconfinement: even if no charmonia survive in the quark-gluon plasma, statistical hadronization at the QCD phase boundary of the many tens of charm quarks expected in a single central Pb-Pb collision could lead to an enhanced, rather than suppressed production probability when compared to results for nucleon-nucleon reactions scaled by the number of hard collisions in the Pb-Pb system.