The central region of the Virgo cluster of galaxies contains thousands of globular clusters (GCs), an order of magnitude more than the numbers found in the Local Group. Relics of early star formation epochs in the universe, these GCs also provide ideal targets to test our understanding of the Spectral Energy Distributions (SEDs) of old stellar populations. Based on photometric data from the Next Generation Virgo cluster Survey (NGVS) and its near-infrared counterpart NGVS-IR, we select a robust sample of ~2000 GCs with excellent photometry and spanning the full range of colors present in the Virgo core. The selection exploits the well defined locus of GCs in the uiK diagram and the fact that the globular clusters are marginally resolved in the images. We show that the GCs define a narrow sequence in 5-dimensional color space, with limited but real dispersion around the mean sequence. The comparison of these SEDs with the predictions of eleven widely used population synthesis models highlights differences between models, and also shows that no single model adequately matches the data in all colors. We discuss possible causes for some of these discrepancies. Forthcoming papers of this series will examine how best to estimate photometric metallicities in this context, and compare the Virgo globular cluster colors with those in other environments.