Materials and methods

First, the shapefile of Mexico’s potential vegetation map (^{Rzedowski 1990}) was downloaded from CONABIO’s geoportal (http://www.conabio.gob.mx/informacion/gis/?vns=gis_root/usv/otras/vpr4mgw). The vegetation types and the surface area they occupy are shown in Table 1. The shapefile of the vegetation and land use map, series VII (^{INEGI 2021}), was downloaded from the same geoportal, at a scale of 1: 250,000. The 183 classes included in this latter map were assigned to the nine Rzedowski’s vegetation types (Supplementary material 1). Classes that could not be assigned this way -for example, polygons without vegetation or with agricultural or urban use- were assigned to the “No vegetation” class. Polygons that were recognized as belonging to one of the nine Rzedowski’s vegetation types were separated into a new layer. Subsequently, the centroid (coordinate calculated as the arithmetic mean of its minimum and maximum latitude and longitude coordinates) was calculated for each polygon and used to calculate Thiessen polygons (Figure 2A).

Figure 2 Example of the creation of Thiessen polygons and the updated vegetation map. A. The centroids of each polygon with an assigned vegetation type are shown as a point, and the polygons without recognized vegetation (in gray). B. The centroids were used to create Delaunay triangles and their sides to construct the bisector lines that constitute the margins of the Thiessen polygons. C. Thiessen polygons are used to interpolate the polygons without vegetation. D. The polygons with interpolated vegetation are merged to the polygons with vegetation in A. 

Thiessen polygons, also known as Voronoi diagrams or Dirichlet tessellations, can be used as a spatial interpolation method that predicts the attributes of unsampled areas based on the known attributes of neighboring sampled areas (^{Brassel & Reif 1979}, ^{Yang 2009}). Each trio of neighboring centroids was used to form a triangle (Delauney triangles). A perpendicular line was drawn in the middle of each side of the triangle, forming the edges of the Thiessen polygons and their intersections became the vertices of the polygons (Figure 2B). The technique has been used in various areas where it is necessary to know the spatial influence of a given phenomenon, including previous applications for vegetation mapping (^{Sanders et al. 2004}, ^{Xiao et al. 2005}, ^{Biró et al. 2006}).

Once the Thiessen polygons were calculated, covering the entire territory of Mexico, they were intersected with the “No vegetation” polygons, to interpolate the corresponding vegetation type (Figure 2C). When a “No vegetation” polygon fell into two or more Thiessen polygons of different vegetation, the vegetation type with the largest area was assigned to form larger, more continuous vegetation patches. Finally, the original polygons with vegetation and the polygons with interpolated vegetation were integrated into a single layer (Figure 2D). Spatial operations were performed in ArcGIS 10.5 (^{ESRI 2013}). The updated potential vegetation map was exported as a shapefile at a scale of 1:250,000 (Supplementary material 2).

The accuracy of a map describes how well the mapped features match the features that occur on the ground. The accuracy can be assessed quantitatively or qualitatively. Quantitative methods use independent observation points to estimate how often the mapped vegetation attribute matches that observed in the field (^{Story & Congalton 1986}). To assess the reliability of the generated map, it was compared with the vegetation types reported for 10,000 georeferenced sites from a database that records herbarium collections of Asteraceae species from Mexico, managed by one of the authors (JLV). This database uses the nomenclature proposed by ^{Rzedowski (1978)} in the vegetation file (Table 1). The Asteraceae family was selected because it has a wide distribution in the country and is well-represented in all vegetation types (^{Villaseñor 2018}).

The correspondence or discrepancy between the validation points and the map was quantified and used to construct confusion matrices to summarize the number of true positive (TP), true negative (TN), false negative (FN), and false positive (FP) cases (^{Story & Congalton 1986}, ^{Sánchez-Muñoz 2016}) for all vegetation types. Using the mesophyllous montane forest (bosque mesófilo de montaña) as an example, a true positive (TP) occurred when a validation point of this vegetation type coincided with a polygon of the same vegetation type. A true negative (TN) occurred when a point where the mesophyllous montane forest was not recorded coincided with a polygon not assigned to this vegetation type. A false negative (FN) occurred when a validation point of this forest type coincided with a polygon not associated with mesophyllous montane forest. A false positive (FP) occurred when a validation point not recorded as mesophyllous montane forest coincided with a polygon of this vegetation type. These quantifications were performed for each vegetation type. The confusion matrices allowed the calculation of three attributes of quantitative precision: the Kappa statistic, sensitivity, and specificity.

Cohen's Kappa statistic (^{Cohen 1960}) compares the observed agreement between a pair of categorical data sets to that which would occur by chance alone. Kappa ranges from -1 to 1, with a value closer to -1 reflecting a greater degree of disagreement, 0 indicating agreement due to chance, and 1 indicating complete agreement between the data sets. ^{Landis & Koch (1977)} and ^{Monserud & Leemans (1992)} used the following qualitative descriptors to characterize the degree of agreement suggested by the Kappa statistic: very poor to poor (< 0.4), fair (0.4-0.55), good (0.55-0.7), very good (0.7-0.85), and excellent (> 0.85). If the number of agreements observed is smaller than expected by chance, Kappa will be negative, indicating disagreement. Sensitivity (TP / (TP + FN)) is the proportion of cases correctly classified to a vegetation type (TP) relative to the total number of cases classified to that vegetation type, either correctly (TP) or incorrectly (FN). In contrast, specificity (TN / (FP + TN)) is the proportion of cases correctly classified as not belonging to a vegetation type (TN) relative to the total number of cases classified as not belonging to a vegetation type, either correctly (TN) or incorrectly (FP). The Kappa statistic, specificity, and sensitivity were also calculated for Rzedoswki´s (1990) potential vegetation map. These analyses were implemented in R 4.3.1 (^{R Core Team 2024}), using the “yardstick” package (^{Kuhn et al. 2024}).

Determining which of the categories evaluated has more or fewer classification errors is also informative. ^{Rosenfield & Fitzpatrick-Lins (1986)} proposed the Conditional Coefficient of Agreement (conditional Kappa) to measure each category’s precision. To perform this calculation, we used the same 10,000 Asteraceae sites described above (according to their geographic coordinates), with their vegetation type observed in the field and the vegetation type of the map interpolated with the Thiessen polygons. The conditional Kappa value was also calculated for ^{Rzedowski’s (1990)} potential vegetation map. These calculations were performed in Excel, using the formula proposed by ^{Wang & Hu (2017)}.

For both maps, the percentages of omission and commission errors were also calculated for each vegetation type. Returning to the example of the mesophyllous montane forest, an omission error occurs when a validation point of this vegetation type coincides with a polygon different from the focal vegetation type (FN). A commission error occurs when a validation point not coming from this forest type coincides with a mesophyllous montane forest polygon (FP). In other words, omission and commission errors measure the fraction of validation points incorrectly classified by the map. The desirable goal is a proportion of error as low as possible (^{Fielding & Bell 1997}).

Results

Figure 1B shows the updated potential vegetation map, and Table 1 shows the area corresponding to each vegetation type. Xerophilous scrub (matorral xerófilo) is the most widely distributed vegetation type, covering 37 % of the country's surface area, whereas mesophyllous montane forest occupies the smallest area, only 1.4 % of the national territory.

Table 2 shows the confusion matrices of the 10,000 Asteraceae sites used as verification points for Rzedowski’s map and the Thiessen polygons map produced here. The Cohen's Kappa values showed a substantial consensus increase in the map produced with the Thiessen polygons (0.927, P < 0.001) relative to Rzedowski’s map (0.588, P < 0.001). The map obtained here had high sensitivity (0.912) and specificity (0.992) values, while Rzedowski’s map had medium sensitivity (0.583) yet high specificity (0.956). The conditional Kappa value, which shows the precision values for each vegetation type (Table 3), was also higher for our map than for Rzedowski’s. In the latter, the vegetation type with the best rating was xerophilous scrub (0.82), while the Thiessen polygons increased the precision for this vegetation type to 0.99. The vegetation type with the lowest accuracy in Rzedowski’s map was evergreen tropical forest (bosque tropical perennifolio) (0.3), while the map produced here had a consensus of 0.85. Regarding the proportions of omission and commission errors, the Thiessen polygons map generally had lower values compared to Rzedowski’s map (Table 3).

Discussion

Vegetation maps are an effective way to process and illustrate spatial information describing vegetation cover. They play a significant role in real-world applications, such as their relationship with climate, soil, or other ecological factors; they are also relevant for understanding biodiversity at the plant community level, and in the design and planning of protected areas, land use, and environmental assessment (^{Pedrotti 2013}).

Vegetation maps can be constructed using very different philosophies and methods, subject to debate. One of these controversies is whether vegetation maps should be represented as vegetation mosaics with defined boundaries (like the maps analyzed in this work) or as vegetation gradients with fuzzy boundaries (^{Triepke 2017}). The advantages or disadvantages of each of these approaches are beyond the scope of this study since only maps with sharp edges were analyzed here, although the interested reader is referred to ^{Feilhauer et al. (2021)}.

More often than not, users must rely on existing maps, which may be based on data collected at different times, drawn at different scales, and with different accuracy levels (^{Küchler & Zonnelveld 1988}, ^{Demers 1991}, ^{Franklin 1995}, ^{Bredenkamp et al. 1998}). Before using one map for a given analysis, users should evaluate the strengths and weaknesses of various maps, and then select one that meets their criteria for accuracy, precision, and currency that are compatible with their objectives. Selecting the most suitable map can ensure a certain degree of reliability in the results.

Although vegetation maps are useful in various applications, according to ^{Demers (1991)}, they have been focused more on communication than analysis; in other words, the predominant purpose of vegetation maps has been to produce thematic visual patterns that correspond to the vegetation classification used. This communicative orientation underlying the production of some maps limits their analytical utility because attributes such as the map's accuracy, precision, and validity are left aside. These attributes are essential in the analytical paradigm of cartography, which focuses on the spatial operations that can be performed with a map, mainly with a geographic information system (^{Moellering 2000}).

The process of creating vegetation maps of sufficient quality to be used for analytical operations requires abundant human and material resources. For example, ^{Biró et al. (2006)} worked for four years to create a vegetation map for a region in Hungary; their team consisted of 56 people, who used paper and digital cartography and satellite images processed in geographic information systems and field validation. Similarly, ^{Su et al. (2020)} updated China's vegetation map, involving 250 researchers over 20 years of work in the office and the field. In Mexico, INEGI regularly publishes (¹⁹⁹⁷, ²⁰⁰¹, ²⁰⁰⁵, ²⁰⁰⁹, ²⁰¹³, ²⁰¹⁶, ²⁰²¹) updates to the country's land use and vegetation map, applying the institute's considerable human and material resources. Some alternatives consume fewer resources; for example, ^{Eva et al. (2002)} created a vegetation map for South America, using satellite images and an unsupervised automatic classification. Similarly, ^{Vaca et al. (2011)} produced a vegetation map for southeastern Mexico, using sites obtained from Rzedowski’s map, to which they associated climatic variables, and with this, they trained a Random Forest model to predict the distribution of vegetation types. These applications can mean savings in the resources used to create vegetation maps.

Map accuracy describes the degree to which mapped attributes match attributes that occur on the ground. This can be assessed quantitatively or qualitatively. Quantitative methods use independent point observations to estimate how often the mapped vegetation state matches the observed state. Confusion matrices, which summarize the frequencies of matches and mismatches, can be used to calculate indices, such as the Kappa statistic, to estimate the agreement between mapped and observed attributes relative to that expected by chance (^{Legendre & Legendre 1998}). To our knowledge, this study is the first to fully evaluate the accuracy of ^{Rzedowski’s (1990)} digital vegetation map.

When comparing the predictive success using Cohen's Kappa coefficient, Rzedowski’s potential vegetation map, with a Kappa value of 0.588, can be rated as good according to the criteria of ^{Landis & Koch (1977)} and ^{Monserud & Leemans (1992)}, while the map based on Thiessen polygons (0.927) can be rated as excellent. Rzedowski’s map also showed lower sensitivity (0.583) compared to the interpolated map with Thiessen polygons (0.912); this result indicates a lower capacity of the former to correctly classify the validation sites with the polygon, including the corresponding vegetation type. On the other hand, both maps have a high specificity, successfully identifying true negatives (TN).

When evaluating each vegetation type individually using conditional Kappa, Rzedowski’s map shows comparatively lower agreement between the validation points and the vegetation polygons than the map obtained with the Thiessen polygons. This lower accuracy is also noted in the percentages of omission and commission errors, which are again higher in Rzedowski’s map. A notable example of these differences is the tropical evergreen forest (TEF), with a conditional Kappa value of 0.3 (poor) in Rzedowski’s map compared to the value of 0.85 (very good) in our map. This means that the former map performs poorly at correctly assigning the validation points of this vegetation type to the polygon corresponding to the TEF. ^{Márcia-Barbosa et al. (2013)} related omission and commission errors with the degree of underestimation or overestimation of a map. A high rate of omission errors is related to an underestimation of the map and a high rate of commission errors, with a map overestimation. The TEF polygon in Rzedowski’s map has a commission error rate of 68 %; this high rate is mainly due to the high number of validation records from tropical sub-deciduous forests (TSF, 361 records, Table 2) that fall within it, which is substantially higher than the number of true positives (278 records). The high percentage of commission errors suggests that Rzedowski’s TEF polygon is possibly overestimated. This overestimation is noticeable in the area of this vegetation type since it is more than twice as large as that of the polygon interpolated using Thiessen (Table 1), which has a lower percentage of commission errors (14 %). The confusion between these two types of vegetation may also stem from the fact that they are not readily distinguished, given their strong floristic and structural similarities. Such similarities often cause errors in identifying the local vegetation type by botanical collectors, who are not always sufficiently well-trained to discriminate against similar vegetation types. Another possible explanation for the difference between Rzedowski’s map and ours is the change in vegetation due to successional or anthropogenic processes in the interval between both works. Similar changes have already been reported in other regions of the world (^{Pancer-Koteja et al. 2009}, ^{Kong et al. 2018}, ^{Zhang et al. 2022}).

Another example of the contrasting accuracy of these maps is the mesophyllous montane forest (MMF); although Rzedowski’s map has a very good conditional Kappa value (0.77), its percentage of omission errors is high (66 %). The confusion matrix (Table 2) shows that the field records to validate the MMF classified by Rzedowski’s map fall in TEF polygons (93 records) and coniferous and oak forests (65 records), which omit MMF polygons. The high percentage of omission errors suggests underestimating the polygon’s area. The area of Rzedowski’s MMF was only 66 % of that projected by our polygon (Table 1), which has better agreement with the validation points, as shown by an excellent conditional Kappa value (0.92) and low rates of omission (6 %) and commission errors (7 %). ^{Vaca et al. (2011)} and ^{Cruz-Cárdenas et al. (2012)} pointed out the underestimation of the surface occupied by the MMF in Mexico.

The vegetation map created by interpolation with Thiessen polygons generally shows better accuracy than the Rzedowski map. The tests indicate greater correspondence between the validation points relative to the polygons of the vegetation types in the interpolated map and lower accuracy in Rzedowski’s map. Similar results were found by ^{Vaca et al. (2011)}, who made a vegetation map of southeastern Mexico and compared it against that of Rzedowski using field sampling sites.

The map presented here is not error-free; Tables 2 and 3 list several of them. These flaws may arise from the misassignment of the vegetation types of the Land Use and Vegetation map (^{INEGI 2021}) to one of the nine vegetation types of ^{Rzedowski (1990)} (Supplementary material 1). We may also have erroneously assigned “No vegetation” polygons that fell between two Thiessen polygons with interpolated vegetation. A further source of error could be the field validation records of Asteraceae, where the collector may not have recognized the vegetation type correctly or when the point fell in an ecotone (the fuzzy limits mentioned above); such is the case of tropical evergreen forest and tropical sub-deciduous forest discussed above.

A map created mainly with the aim of communication, such as Rzedowski’s vegetation map, cannot be used under the analytical paradigm just because it can be incorporated into a geographic information system. Rzedowski's intention in 1978 was likely to provide a broad overview of the distribution of vegetation types in the country rather than providing a guide to deciding what vegetation type was present in a given location. It is important to recognize that maps mainly oriented toward communication may not prioritize accuracy and that the ability to upload the map into a GIS does not guarantee that it is appropriate for analytical applications; unfortunately, despite this reality, using such maps for purposes different from their original intention is becoming a common practice. In this paper, we calculated the same precision values for Rzedowski’s map and the map interpolated with Thiessen polygons; these values are important because they will allow users to estimate the precision in future analyses that use them. This paper presents a method for creating a vegetation map, and the result (Supplementary material 2) is accessible as the product of this method. The methodological procedure can be adjusted to produce a vegetation map under any other classification scheme of vegetation types. This effort thus constitutes a heuristic exercise that can serve as a hypothesis to support or refute future vegetation analyses at different scales in the country. The updated map is available to users to evaluate whether it is appropriate for their analysis. To date, there is no other proposal for a vegetation map (potential or natural) for Mexico that would replace the great effort made by Rzedowski almost five decades ago. We are only aware of regional or state efforts (for example, ^{González-Elizondo et al. 2007}, ^{Vaca et al. 2011}, ^{González-Elizondo et al. 2012}, ^{Siqueiros-Delgado et al. 2016}). However, Rzedowski's map is the main source used to discuss vegetation types in the country, despite its coarse scale and the difficulty of specifying many of the types he used in his classification locally. Proof of this is the number of omission and commission errors found. A better alternative is undoubtedly the land use and vegetation maps generated by ^{INEGI (2022)}, although these, besides natural vegetation, include areas that humans have transformed. Unfortunately, those maps come with their own difficulties for interpretation and use (too many categories, many of them slight variations of one or another, or poor knowledge of interpretation using a geographic information system), which limits its usefulness for many users interested in the subject.

We hope that the map presented here, and the accompanying shapefile format, will be useful in future studies that require interpretation of the vegetation types in Mexico. But above all, we hope it will be a basis for critique, review, and update as new information is obtained. We hope that the vegetation synthesis presented here, which updates Dr. J. Rzedowski's key contribution to the understanding of vegetation types in Mexico, will help to achieve, in a short time, a more generalized proposal that serves many other objectives. As ^{Ibarra-Manríquez et al. (2022}) argue, this would greatly impact ecology and many other areas of biology and related sciences.

Supplementary material

The supplementary material for this article can be consulted here: https://doi.org/10.17129/botsci.3598.

Supplementary material 1 and 2.