Global navigation satellite systems (GNSS) provide a great data source about the ionosphere state. These data can be used for testing ionosphere models. We studied the performance of nine ionospheric models (Klobuchar, NeQuickG, BDGIM, GLONASS, IRI-2016, IRI-2012, IRI-Plas, NeQuick2, and GEMTEC) both in the total electron content (TEC) domain—i.e., how precise the models calculate TEC—and in the positioning error domain—i.e., how the models improve single frequency positioning. The whole data set covers 20 years (2000–2020) from 13 GNSS stations, but the main analysis involves data during 2014–2020 when calculations are available from all the models. We used single-frequency positioning without ionospheric correction and with correction via global ionospheric maps (IGSG) data as expected limits for errors. Improvements against noncorrected solution were as follows: GIM IGSG—22.0%, BDGIM—15.3%, NeQuick2—13.8%, GEMTEC—13.3%, NeQuickG and IRI-2016—13.3%, Klobuchar—13.2%, IRI-2012—11.6%, IRI-Plas—8.0%, GLONASS—7.3%. TEC bias and mean absolute TEC errors for the models are as follows: GEMTEC—−0.3 and 2.4 TECU, BDGIM—−0.7 and 2.9 TECU, NeQuick2—−1.2 and 3.5 TECU, IRI-2012—−1.5 and 3.2 TECU, NeQuickG—−1.5 and 3.5 TECU, IRI-2016—−1.8 and 3.2 TECU, Klobuchar—1.2 and 4.9 TECU, GLONASS—−1.9 and 4.8 TECU, and IRI-Plas—3.1 and 4.2 TECU. While TEC and positioning domains differ, new-generation operational models (BDGIM and NeQuickG) could overperform or at least be at the same level as classical empirical models.