Study area

Figure 1 shows the BG, VM, and BB regions, which have contrasting hydrological characteristics and represent different geographic and climatic zones in Mexico (^{Rascón & Rivera, 2005}). The BG basin is part of the northwest zone of Hydrological Region Grijalva-Usumacinta at the south of Mexican territory (^{CONAGUA, 2007}). This region has a humid tropical climate (^{García, 1998}), a low degree of alteration, and a high-pressure level on water resources (^{Cotler, Garrido, Bunge, & Cuevas, 2010}). The salaried employees in this region are mainly engaged in industrial activities but about 30 % of the total receive less than a minimum wage (^{INEGI, 2009}). The most important population center is the city of Villahermosa with more than 300,000 habitants (^{INEGI, 2015}). The hydrological parameters of this exoreic region are MP 1 500-4 500 mm (^{Cuervo-Robayo et al., 2014}), mean annual runoff of 1 000-2 000 mm (^{Jiménez & Maderey, 1992}), evapotranspiration of 1 401-1 500 (^{Maderey, 1990}) and RC of 0-30 % (^{INEGI, 2010}). VM is in the central high plains at the homonymous Hydrological Region (^{CONAGUA, 2007}). This region has a temperate climate (^{García, 1998}), an extreme degree of alteration, and medium pressure on water resources pressure level (^{Cotler et al., 2010}). It is highly urbanized and has the highest population density among the compared basins in this study. In this area, approximately 70 % of the population is employed in shops and services activities (^{INEGI, 2009}). The most important population center in Mexico City and its metropolitan area with more than 19 million inhabitants (^{INEGI, 2015}). This exorheic basin has the hydrological parameters: MP of 400-1 500 mm (^{Cuervo-Robayo et al., 2014}), mean annual runoff of 10-500 mm (^{Jiménez & Maderey, 1992}), evapotranspiration of 400-800 (^{Maderey, 1990}), and RC of 0-30 % (^{INEGI, 2010}). BB is in the north at the Bravo Conchos y Rio San Juan Hydrological Regions (^{CONAGUA, 2007}). It has a semi-arid climate (García, 1998), a high degree of alteration, and a high-pressure level on water resources (^{Cotler et al., 2010}). In this part of the territory, the population density varies in an interval of 23 and 80 inhabitants/km². The population is mostly employed in industrial activities and other salaried employments (^{INEGI, 2009}). The most important population center is the city of Monterrey, with 1 109 171 inhabitants (^{INEGI, 2015}). The hydrological parameters of this exorheic basin are MP of 500-800 mm (^{Cuervo-Robayo et al., 2014}), mean annual runoff of 10-50 mm (^{Jiménez & Maderey, 1992}), evapotranspiration of 400-500 (^{Maderey, 1990}), and RC of 0-20 % (^{INEGI, 2010}).

Figure 1 Location of study regions BB, BG, and VM. Source: Own design with data from ^{Conagua (2007)}. 

Information resources

Table 1 shows information from the L8, S1, and S2 images that were selected for each region, which was downloaded from https://earthexplorer.usgs.gov/ and https://scihub.copernicus.eu/. These images correspond to the end of the rainy season to include LC during a period that has a high level of runoff. For comparison purposes, the closest possible dates in September were used. The areas in the images having the largest spatial overlap were selected for each region. They were also selected based on the criterion of low cloud cover. Table 1 shows the minimum and maximum coordinates, and the areas with image overlap of the three sensors for each of the regions. Also, the identifier and the date of each of the downloaded images are included.

Base criteria for image selection

To identify if each sensor can help to determine the LC of the Mexican territory and what the advantages are given by each of them, a comparative analysis was for SR, TR, and SpR was done. The analysis of SpR of optical images included the sensibility of the sensor for the spectral signatures involved in the RC: tap water, concrete, brownish-gray sand, and vegetation. For this, a graphic (fig. 2 included in the results section) was generated with the mentioned signatures and the ranges of the electromagnetic spectrum of each sensor (^{Baldridge, Hook, Grove, & Rivera, 2009}; ^{ESA, 2013a}; ^{ESA, 2013b}; ^{USGS, 2016}). For S1, the polarizations Vertical-Horizontal (VH), Vertical-Vertical (VV), and its subtraction (VV-VH) were used, following ^{Abdikan et al. (2016)}.

Processing and evaluation of the LC classification

For L8 and S2, a radiometric and atmospheric correction was performed with the DOS method, using the QGIS 2.18 “Semi-Automatic Classification Plugin” (SCP) (^{Song, Woodcock, Seto, Lenney, & Macomber 2001}; ^{Congedo, 2016}). For S1, noise elimination, calibration, multi-sweep, and terrain corrections were performed using the “SNAP Tool” (^{ESA, 2013a}). The LC classes were defined based on the characteristics of each region, using as a reference official data INEGI on land use and “series V” vegetation, as well as the LC-K correspondence table (^{INEGI, 2013}). The SCP was used to create the regions of interest (^{Solis et al., 2005}), which were defined as: water bodies, human settlements, bare soil, cropland-pasture, clouds, shadows, xerophytic mesquite, thorny shrubs, submerged vegetation, flood-prone zones, forest, and shrubs.

While various supervised classification algorithms exist, our study is limited to Minimum Distance (MD), Maximum Likelihood (ML), Spectral Angle Mapper (SAM), and SVM. MD is mathematically simple. It calculates the Euclidian distance from the spectral signatures, assigning the class corresponding to the signature that is closest to the pixel (^{Wacker & Landgrebe, 1972}; ^{Richards & Xiuping, 2013}; ^{Congedo, 2016}). ML calculates probability distributions in accordance with the Bayes theorem (^{Sudhakar, Sridevi, Ramana, Rao, & Raha, 1999}; ^{Bruzzone & Prieto, 2001}; Richards & Xiuping, 2013; ^{Zhao & Zhao, 2013}; ^{Jia et al., 2014}). SAM calculates the angle between the spectral signatures of image pixels and training spectral signatures (^{Kruse et al., 1993}). SVM defines the separation surface between two classes based on the definition of support vectors (^{Manning, Raghavan, & Schuetze, 2009}; ^{Theodoridis & Koutroumbas, 2008}; ^{Mountrakis, Im, & Ogole, 2011}). The four algorithms were applied to each region using 25 to 30 training areas per class, which corresponded to a sampling of 2 to 5 % of the total area. The training areas were selected based on a stratified random sampling of each type of cover. MD and SAM were applied with SCP-QGIS, ML with GRASS GIS, and SVM with Orfeo-QGIS.

The accuracy that a classification presents with the actual value in the terrain depends on the SR, TR, and radiometric resolution (RR) of each sensor (^{Rogan & Chen, 2004}). The validation of the obtained LC has been performed taking reference data from the National Institute of Statistical Geography (INEGI, acronym in Spanish), as well as with high-resolution images from Google Earth and Bing Aerial (^{Abdikan et al., 2016}; ^{Yu & Gong, 2012}). At least 30 validation points per LC class were also randomly generated, and the individual results of each algorithm-sensor-region combination were validated with those points according to ^{Abdikan et al. (2016)}, and ^{Yu and Gong (2012)}. After that, the unbiased population matrix and the accuracy calculations have been performed (^{Pontius Jr. & Santacruz, 2014}) to validate the classifications correcting sampling bias (^{Pontius Jr. & Millones, 2011}; ^{Pontius Jr. & Santacruz, 2014}; ^{Estoque et al., 2018}; ^{Huang, Huang, Pontius, & Tu, 2018}). The equations proposed by ^{Pontius Jr. and Santacruz (2014)} were applied to each combination of algorithm-sensor for each region to obtain the following indicators: unbiased population matrix, overall accuracy (OA), overall omission (OO), overall commission (OComm), the accuracy of cover (AC), omission of cover (OC), and commission of cover (CC). To synthesize the results, the average OA was obtained for each hydrological region, as well as for each sensor and each classification algorithm.

Calculation of RC

The first step to calculate RC was to select the LC image. To this end, the criteria resulting from section “Information resources” are used, as well as the results of the classifications with the highest OA. After that, the calculation was done accordingly to the following equations.

RC is expressed in function of MP and the K parameter, as shown by (^{DOF, 2015}):

K is obtained according to tables that combine LC and ST, such as those reported by ^{Solis et al. (2005)}.

Equation (1) and Equation (2) were applied to the selected images. They were transformed into a decision algorithm to calculate RC that used a K factor raster obtained from the correspondence tables by ^{Solis et al. (2005)}, a raster of annual MP obtained with regression-kriging interpolation obtained accordingly to the reported by ^{Tapia-Silva (2019)}, and an ST raster (^{INIFAP, 1995}). The use of correspondence tables is one method that is suggested in countries such as Mexico for determining the RC (^{DOF, 2015}). These tables show the relationship between the K parameter and ST (^{INIFAP, 1995}; ^{Solís et al., 2005}).

Base criteria for image selection

Figure 2 presents a comparison of the spectral bands selected for optical sensors L8 and S2. The top section shows the behavior of the spectral signatures of tap water, concrete, brownish-gray sand, and vegetation (^{Baldridge et al., 2009}). The wavelength, in micrometers, is shown on the x-axis and reflectance on the y-axis. Dotted lines run through and intersect the wavelengths that are perceived by each sensor. This graphic enabled us to select the optical sensor bands. The bands selected were those that are sensitive to the electromagnetic spectrum range, where reflectance varies in function of the type of LC. We selected bands 2, 3, 4, 5, 6, and 7 to classify LC using L8. Since band 1 is used to study shallow water, we included that in the BG region only. We selected bands 2, 3, 4, 5, 6, 7, 8, 8A, 11, and 12 to classify LC using S2. We also included 1 with S2 for the BG region to improve the definition of water bodies.

Figure 2 Spectral comparison of sensors L8 and S2 to determine LC. Developed by the authors with information from ^{Baldridge et al. (2009)}, ^{ESA (2013a)}, ^{ESA (2013b)}, and ^{USGS (2016)}. 

Processing and Evaluation of the LC Classification

Table 2 shows the OA of all the algorithm-sensor combinations for each of the regions, in order of high to low. Region BG had the highest average OA for all the classifications (70.3 %), while BB resulted in 61.5 %, and VM in 58.2 %. The size of the area selected for each region may be a factor of influence, given that the highest OA corresponds to the area of analysis that is smallest in size. In BG, the water bodies cover a large area, which could have influenced the high OA obtained since an extra band helped to identify them using optic sensors and radar images are useful to detect water. In the case of the VM, the dynamics of change and the complexity of the area may have influenced the OA obtained for that region, which had the lowest average of the study regions. The bottom of Table 2 shows the results for each sensor. L8 presented the highest average OA, with 76.1 %, even though it had the lowest SR. That gives it the distinction of being a very useful sensor for defining LC types as well as other biophysical parameters. The performance of S2 was similar, with an OA of 75 %. S1 resulted in an OA of 38.8 %, which is much lower than the sensors mentioned above. The sensors that had the highest OA were those with the highest SpR. The low accuracy found with S1 may be due to the use of only three bands for the classification, which strongly differs from a report by ^{Abdikan et al. (2016)}, who obtained the best results with the combination of bands VV, VH, and the subtraction of VV-VH. The MD was the classifier with the best performance, with an OA of 71.2 %, followed by SVM with 70.5 %, and SAM with 60.9%. The poorest performance was ML, with 50.7 %.

For the BB region, eight combinations had an OA over 60 %. The higher OA in this region was obtained with S2 combined with MD, SVM, and SAM (in that order), and the lowest OA was obtained with ML combined with L8 and S1 sensors. Based on these results, S2 is recommended for regions that are similar to BB (little precipitation and a large presence of irrigated cropland). The SR of S2 facilitates the identification of small objects of interest and was therefore helpful for this region, which contains large territorial contrasts.

For the BG region, eight combinations had an OA over 80 %: L8 and S2 combined with the four algorithms. The highest OA was 90.5 %, which was obtained with L8-SVM which presented the best performance of all combinations in this study. S2 with MD ended up in fourth place overall, even though its average OA was higher than those of the other algorithms. In contrast, the ML for this region had the second-highest OA, although it performed poorly in other combinations. The BG images have 1 to 3 % cloud cover, and images without clouds could not be found for that date. In addition, cloud cover and their shadows presented the highest accuracy in this region, which positively affected overall performance. Normally, S1 is advantageous for regions such as BG, where there is high precipitation and cloud cover since its technology is insensitive to those factors. However, in our case, this advantage did not mean improvements in the OA and this may be related to the fact that the selected images presented the above-mentioned percentages of cloudiness.

For the VM region, six combinations had an OA of over 60 %. L8 combined with the four algorithms resulted in the highest OA, even though more bands were used with S2. This indicates that this combination is suitable for regions with spatial configurations that are similar to VM (predominantly urban). This combination presented higher OA with SAM and MD. The OA was less than 60 % for S2 combined with SVM and ML. In this region, the results were poorest for S1-ML, S1-SAM, and S2-ML.

Figure 3 presents (in % of the domain) the presence of cover (PC), represented by a black square. AC (green bars), OC (red bars), and CC (blue bars) for the classifications with the best performance in each region are shown in the same figure. For BB, cropland-pasture presented the highest PrC (55.7) and AC (51.4) with the S2-MD combination. AC for human settlements was 3.0. For this category, CC was greater than OC, which indicates that this category was overestimated. Xerophytic mesquite and thorny shrubs have very similar spectral signatures, which resulted in high OO and CC values. The errors for bare soil may be due to confusion with areas where crops had been recently planted. Thus, a high number of pixels was mistakenly assigned to bare soil (1.8 of CC), and the category was overestimated. For BG, bare soil was the land cover with the lowest AC and PC, with 1.1 % in both values. The CC for cropland-pasture (0.9) and flood prone zones (2.5) may be due to their having very similar spectral signatures when crops begin to grow in those areas. Mangroves were distinguishable in the image and had a high AC (13 % for 15.5 % of PC). Water bodies had a low CC (0.1 %) since there was virtually no spectral confusion. Human settlements were underestimated, since CC for this class was lower than OC, with a PC of 2.8 %.

Figure 3 PC, AC, CC, and OC classifications of best performance in the analysis (as % of domain). 

In the VM, urban zones were dominant, and the AC was high (46.3 % with a PC of 49.9 %), and they were correctly classified by the L8-SAM combination. Bare soil, forest, and water bodies had lower AC and PC. Cropland-pasture and shrubs had a low PC. For these cases, AC exceeded half of the PC. As in the other regions, OC was higher for bare soil and cropland-pasture, and the latter was confused with the shrub. Shrub and urban areas were overestimated, with a CC of 4.9 and 3.6 %, respectively. The remaining classes were underestimated, as indicated by an OC less than CC.

In all of the study regions, the dominant class made it possible to identify training areas with a high number of pixels and a homogenous spectral signature, resulting in high AC values. This proportional relation between PC and AC was reflected by bare soil, which had a low PC and AC. While shadows and clouds resulted in high AC levels, they do not contribute to the RC value, and furthermore, they make it difficult to adequately determine RC by making it difficult to identify the LC on the site.

Table 3 presents the PC results from the sensor-classifier combinations for the three regions. The combinations with the highest OA are highlighted in bold (S2-MD-BB, L8-SVM-BG, and L8-SAM-VM). For BB, L8-SVM and S1-SVM were the combinations with PC that were most similar to the one with the highest AC. The combinations with PC that were least similar to those with the highest AC were L8-ML and S1-ML. S1-SAM and S1-MD are not recommended because they overestimated or underestimated one of the classes studied.

For the BG region, the PC obtained with S2-SVM and S2-MD were similar to those obtained with L8-SVM (which had the highest OA). The combinations that overestimated or underestimated PC were all the combinations with S1. As previously observed, using S1 for regions similar to BG (humid tropical) with cloudy conditions less than 3 % would result in considerable errors that would need to be corrected. Frequent cloud cover significantly reduces the ability to use optical images and according to our findings, these should be occupied with cloud cover less than 3 %. Additionally, we do not recommend the L8-ML combination to study human settlements in regions similar to BG, since that would overestimate PC.

For the VM region, the PC obtained with L8-SAM was comparable to S1-SVM and S2-SAM. The combinations that differed the most from what was obtained with L8-SAM were S1-SAM, S1-MD, and S2-ML, and therefore, they are not recommended. The PC obtained with L8 was similar to all the algorithms tested, which is a strong argument for recommending the use of this sensor with any of the classifiers for regions where urban land cover is predominant.

Calculation of RC

Table 4 shows the factors that this study took into account to select the most useful image for calculating RC: start of the study period, TR, SR, SpR, and Sensitivity to weather conditions. At the bottom of the table, the obtained average OA by the sensor is included as a reference. The TR of S2 (five days) was better than that of L8 (16 days) and S1 (12 days). S2 began to acquire images in 2015 and L8 in 2012. L8 can be complemented with the Landsat series, which has been available since 1972. S1 images have been available for all of Mexico since 2016. While S1 and S2 have the disadvantage of being only recently launched their SR is higher than that of L8 which makes it possible to distinguish LC having a small presence. For example, in BB, the SR of S2 made it possible to classify areas that had high contrast and covered a small area. Nevertheless, the size of the studied hydrological region should be considered, since a high SR requires a longer processing time and higher computing capacity.

All three sensors were able to distinguish the four main LC categories in the study: water, bare soil, human settlements, and vegetation. The greater number of bands in S2 did not result in a higher OA. Also, the three bands in S1 did not provide good classification results. Although this sensor is associated with the good classification of urban areas with a high SR (^{Goldewijk, 2001}; ^{Liu & Tian, 2010}; ^{Loperfido, Noe, Jarnagin, & Hogan, 2014}), its use in the present study did not improve the results in the VM region (highly urbanized). Given the TR of the sensors used in this study, any one of them could be used to perform the analysis presented by this investigation. S1 can be used with high cloud cover, in which case it will be helpful to use the combinations that the previous section indicated as having the highest accuracy. Certain actions are recommended to improve the classification results, such as post-processing techniques, increasing the number of bands, and selecting training areas that have a larger number of pixels and spectral homogeneity (^{Congedo, 2016}; ^{Richards & Xiuping, 2013}; ^{Chuvieco, 1991}). Table 4 can be used as a reference for selecting the image. Data from S1 and S2 can be used if the project being conducted involves the study of current LC, but if the study includes historical changes it will be limited to the data from the Landsat series, which has generated images since 1972. All three sensors can be used for regional or local analyses, and enable identifying water bodies, vegetation, bare soil, and human settlements with a good degree of accuracy. After consulting Table 4 to select the most suitable sensor, we recommend that it be combined with the algorithm having the highest degree of accuracy, as reported by this study’s findings (Table 2).

Figure 4 shows the RC results for each region, including the LC classification obtained with the sensor-algorithm combination having the highest OA, and four geographic layers: MP, ST, the resulting RC, and the RC from ^{INEGI (2010)} (official information) for comparison purposes. The following observations are applicable to all the study areas. The spatial distribution of the RC obtained by the present study is drastically different than those from INEGI’s hydrological chart, 2010. The information presented in the INEGI layers is much more general than the layers that resulted from our study. The information resulting from our study has a higher SR (10 m). In the case of INEGI, the data are at a scale of 1:50 000. INEGI’s RC values were not determined in a precise manner, but rather, there are five discrete categories, the highest of which is between 0.10 and 0.20 for BB and VM, and ≥ 0.3 for BG.

Figure 4 LC classifications were obtained with the sensor-algorithm combination having the highest OA and the layers MP, TS, resulting RC, and the RC from INEGI (official information), for the three regions studied. 

In the case of BB, the combination that was selected was S2-MD. MP ranged from 470 to 620 mm and increased towards the southern portion of the hydrological region. The dominant ST was medium, with small areas of fine texture and very little thick texture in the north. The RC ranged from 0.17 to 0.28, with the latter exceeding the upper limit reported by INEGI (0.20). The spatial distribution was quite different from INEGI’s hydrological chart. The RC in the chart ranged from 0 to 0.5 % in the southeast and north, while our study resulted in a smaller RC for a larger area of the northern portion of the region. In our study, the highest values were found in the south, whereas the INEGI chart identified the northwest as having the highest values.

For the BG region, the LC was obtained with L8-SVM. The lowest MP was 1 780 mm, which increased to nearly 2 000 mm towards the west. The dominant ST was fine, medium in the south, and thick in the north. The RC ranged from 0 to 1, with a predominance of values over 0.85. The highest RC values were in the western portion of the hydrological region, whereas the INEGI chart identified the central region as having the highest values.

For the VM region, LC was obtained with the L8-SAM combination. The minimum MP was 570 mm, which increased to 1 190 mm towards the west. The dominant ST was fine, with the medium in the central region and thick in the northeast. The RC ranged from 0.17 to 0.6, with a predominance of values over 0.24. The highest values were found in the central portion of the region, which may be because the calculation used by the present study included MP as a variable that is important to the RC value, considering that an immediate relationship between the two can be seen when visually comparing the RC results with the MP map.