This section is divided into four main subsections: First, the analysis of the land use and its spatial distribution at macro-class level is presented. Second, an analysis at class level is performed, similar to the one conducted in the first section. Third, an assessment of the land use change from 2015 to 2021 is carried out as a study case, and, lastly, a discussion of the possible drivers forcing land use change is established.
3.1. LULC Classifier Assessment and Selection
The three classification methods were run for the seven-year period proposed in this analysis. Results show a similar behavior for extracting LULC information using the same level of ancillary data. This condition allowed a fair comparison between the classification algorithms.
Figure 3 shows the macro-class classification for the year 2017, which was the year which presented higher discrepancy among LULC classification methods. Graphically, one can observe that algorithms performed comparably. Significant variations were found within the clusters sharing bare soil and and built-up zones, whereas a more accurate definition was found for those composed of built-up vegetation and bare soil vegetation. These variations were attributed mostly to the spectral signature dispersion and departure within the pixels defining each class and the nearby areas. Bare soil and vegetation areas, being those mostly distributed over the west and east parts of the city, respectively, presented similar spatial distributions.
Numerically, all the classification methods presented good behavior for both overall precision and Kappa values, all being above 85% and 0.85, respectively (
Table 6). The spectral angle mapping (SAM) was the method that showed the smallest overall precision and Kappa values, 89% and 0.86, respectively. The highest accuracy was obtained using the maximum likelihood classifier (MLC), which presented overall precision above 90% and Kappa coefficients above 0.90 for all the years in the study period. It was observed that the MLC required less field data to obtain accuracy compared to the other classification methods.
The three classification methods showed the advantages of a perfect decision boundary to distinguish the macro-classes and a consistent mathematical expression based on the decision boundary for further classifications. However, these methods might overtrain the decision tree, as the training set needs several examples to cover all the possible cases within a specific class. Lastly, training tends to require significant computational time to be effective. As the MLC presented the highest level of accuracy under the same training set, it was selected as the main classification method for this study.
3.2. Macro-Class-Based LULC Classification
The classification based on macro-classes represented a course scheme to generate four of the main land uses found in Coatzacoalcos. These macro-classes are water, built-up, vegetation and bare soil. The spatial distribution of them can be observed in
Figure 4. Graphically, one can see that the terrestrial surface is proportionally divided into three land uses without considering water. Most of the urbanization or built-up areas are in the north and south of the city near the river, as they were sites of the first settlers since the foundation of the city. However, the south potential flood areas and their development have been delayed, and small changes can be observed. It can also be seen that both north and south areas contain significant patches of vegetation, which, including the creation of green areas, is a good practice in urbanization development. However, urbanization in the west area of the city seemed to increment over the study period, but no vegetation buffers were observed in those zones. This indicates the absence of good practices in urbanization and government decisions in terms of territorial planning as well as the economic, social, cultural and environmental unconsciousness of new developers [
70]. These zones may provoke a higher temperature sensation or heat islands, which might translate into health problems for the inhabitants [
71]. Additionally, bare soil areas seemed to increase over time and were more prominent in the west area of the study. As seen in
Figure 4, bare soil replaced areas of vegetation and also increased close to developed areas. These changes can be attributed to water availability during the growing season and afforestation during the development of urbanization.
In terms of vegetation, the east part of the city seemed to remain unaltered, except for zones near Allende Village in the northeastern and southeastern parts of the city, where the industries and petrochemical complexes are located and trigger the soil degradation. Numerically,
Table 7 supports the claim observed on the above maps. Pixels representing water show small variations attributed to possible changes in the surface water balance, but a more refined classification is needed to understand how these variations occurred. The maximum percentage of water land use was
in 2017, whereas its minimum was
in 2019. Urbanization or man-made construction increased from
in 2015 to
in 2021. This increment can be explained by the trend observed in
Figure 5, which shows the last six 5-year censuses. Coatzacoalcos had a mean population growth rate from 1995 to 2015 equal to 2.68%. However, we can observe that a rate equal −0.54% was shown in 2020. This means that population growth decreased in the last census. This variation can be illustrated with the reported built-up information, where man-made construction showed similar increments on a year-to-year basis. Nonetheless, the growth from 2019 to 2020 was its minimum, and no change was seen in 2021. This decrease in population and the variation within estimated in this study have socio-economic drivers. The city of Coatzalcoalcos was the second most dangerous city in 2019 [
72]. This fact provoked serious economic and social issues due to extortion, racketeering and kidnapping. People moved to nearby cities to feel safe, and big companies closed their services; as a result, unemployment rate grew and investments in the city slowed down, which, in turn, decelerated urbanization and industry development.
Bare soil also seems to increase visually in
Figure 4 and is confirmed in
Table 7. Bare soil variations are linked to a decrease in vegetation. These two variables are highly dependent on water availability, droughts and weather fluctuations. We can see that, in 2015, the area of study was 25.67% of bare soil, while, in 2021, it increased up to 26.37%. On the other hand, vegetation presented the inverse behavior, being 26% and 24.96% in 2015 and 2021, respectively. Analyzing the water availability, which is proportional to the precipitation, one can see in
Figure 6 that rainfall tended to fluctuate from a minimum of 1107.6 mm in 2015 to a maximum of 1730.90 mm in 2017. This fluctuation, in general, tells us that, of the last three years of the period of analysis studied here, precipitation reduced drastically in 2019, which can explain the cause of the bare soil increment.
Additionally,
Figure 7 shows that although minimum temperature tended to be steady from 2013 to 2021, mean temperature and maximum temperature increased during the study period. For instance, mean temperature was 26.60
C in 2018, whereas the maximum temperature increased from 28.20 to 30.20
C in 2018 and 2020, respectively. Increments in temperature increases soil evaporation and plant transpiration. If the latter is limited due to lack of water, water stress occurs in plants and they do not reach maturity, which might explain changes in pixel color and the move from vegetation to bare soil. It can be seen that the use of remote sensing for LULC evolution can be related and explained with inland information. This not only promotes the use of remote sensing but also validates its implementation. All the macro-class-based classification showed average user accuracy (UA) and producer accuracy (PA) higher than 90% (
Table 8).
This means that both producers and users of these maps can rely on this classification with at least 90% confidence. These values are comparable to the ones found in [
47,
48]. The mean overall precision for the study period is 94.81%. The Kappa coefficients for all the studied years are located in an excellent range, as they all are higher than 0.81, according to Viera and Garret [
66]. Although this classification was appropriate for the proposed macro-classes, a more refined algorithm was implemented to understand more about the evolution of the LULC of Coatzacoalcos, as it is necessary for decision making purposes.
3.3. Class-Based LULC Classification
A class-based LULC classification was conducted to identify a more discrete evolution of the land use within the study period. This identification will help prioritize the location of areas that require ecological attention, the inclusion of best management practices (BMP) for water and soil conservation as well as a more appropriate location of the upcoming urbanization due to the Interoceanic Corridor project execution. At this level, we identified that pixels representing the sea did not change significantly. This behavior can be observed in
Figure 8 and
Table 9, in which the percentage area representing the sea is about 20.1% for the entire period. However, water courses and water bodies seemed to change over time. Although these changes might not be significant, these changes can be confusing. We observed in
Figure 6 that the annual precipitation throughout the study period decreased. As a result, bare soil area increased due to lack of water. One could assume that water bodies and water courses cannot increase in area. Nonetheless, this claim might not be entirely true, as increments of bare soil decrease infiltration and enhance the overland flow of the areas draining to the closest water courses. For that reason, water courses and water bodies in western areas increased due to increment of overland runoff in areas classified into bare soil with different vegetation densities (i.e., bare surfaces, sparse vegetation or grassland) and some over-flooded wetlands (i.e., swamps).
Man-made construction was divided into industry, residential, commercial and roads classes. Industry shows a constant percentage from 2016 to 2021. The only change occurred in 2015. The industrial area called Etileno XXI, the largest petrochemical complex in Latin America, concluded 99.2% of its construction in 2015 [
75]. This, once again, confirms the accuracy of the generated map. Residence land use increased from 11.85% in 2015 to 12.24% in 2019 and maintained a steady value in 2021 (
Table 9). This, once again, matches the period of violence described above that eventually provoked a decrease in population, as seen in
Figure 5, which means that residence land use remained as it was in the previous years. Commercial areas showed fluctuations due to the city’s economical variation. It can be observed that, in 2018, commercial areas reached their maximum of 0.55%. However, the value decreased at the end of the study period (0.43%). A small increment was observed in 2021, when commercial areas increased to 0.50%. The city experienced transitions from high commercial zones to abandoned places due to the socio-economical problems previously described. These areas transformed either to low density vegetation areas or abandoned buildings. Lastly, in this class, roads were the most difficult to identify due to the Landsat satellite image resolution equaling 30 m. Roads close to the coast line and main avenues crossing from north to south and west to east were detected easily. However, narrow streets were not characterized as roads and they were either characterized as residence, commercial or industry classes. For that reason, although roads showed increments due to increase in population and urbanization development, no claim is expressed in order to avoid confusion in this class.
Vegetation plays an important role in our natural ecosystem and also holds up the biosphere in various ways. Vegetation helps to regulate the flow of numerous biogeochemical cycles, most importantly those of nitrogen, carbon and water. It also contributes to the local and global energy balances. In this study, the classes within vegetation included mixed forest, shrubland, wetland and dune.
Figure 8 shows that mixed forest is the largest extension of the study area. One can see that the eastern and southwestern areas are dominated by this class, especially for those areas where no urbanization is present. Mixed forest occupied about 18% of city and its metropolitan area for the analyzed period, reaching its minimum in 2019 when the precipitation reached its lowest value (
Figure 6). As previously mentioned, at this point, one can confirm that vegetation established in the urban areas were mostly mixed forest, having more pixels where the city was initially developed than where the more recent residential areas in the west of the city have been developed. Shrublands, in general, showed an increment over time. They reached their maximum percentage area, equal to 5.31%, in 2020 with some fluctuations in 2017, a year preceded by annual precipitation below 1500 mm, which represents the driest year of the analyzed period.
One of the most representative classes in this classification is wetland (i.e, swamps) because it represents multiple biological, economic and social values. Swamps provide services to ecological well-being, such as groundwater recharge, water purification, micro-climate regulation, food resources, biodiversity and carbon storage [
76,
77,
78]. In the last decades, in Coatzacoalcos, urbanization has degraded swamps indiscriminately through industrial development or residential areas. These developments followed the erroneous idea that swamps were areas with dangerous species, such as snakes, alligators, mosquitoes, etc., which represented risks to the nearby areas [
79]. This study evidences how wetland or swamp degradation continues. Visually, swamps located south and southeast of the city have mostly transitioned to some degree of vegetation density, either from water stress or landfill. In 2015, swamps represented about 2% of the study area, whereas, in 2021, this fraction reached 0.92%. This fact makes this study a good indicator to prevent and preserve swamps and wetlands within Coatzacoalcos and its metropolitan area due to the important ecological services that they represent.
Dunes, in
Figure 8, are located along the coastline of the study area. In addition, there are some banks utilized by the industry (sandblasting) in the east area across the Coatzacoalcos river. One can see that dunes were replaced mainly by residential areas in the northwest part of the city. In 2015, dunes occupied almost 3% of the area in question. However, this area declined to 1.96% in 2021. The most significant evolution of urbanization occurred from 2017 to 2019. Additionally, the western coastline presented a transition from dunes to bare surfaces because that area is naturally preserved and the absence of tourism has improved the growth of native vegetation. Lastly, bare surfaces are the dominant land use along with mixed forest and residential areas (
Figure 8). This land use increased from 19.31 to 23.19% through the study period in
Table 9. Visually, most of the transitions were from sparse vegetation and high grassland to bare surfaces. This might have occurred because of the decrease in infiltration and precipitation previously mentioned and explored in
Figure 6 and
Figure 7. Sparse vegetation decreased from 2.79% to 0.79% and increased to 0.84% in 2021, whereas grasslands reduced from 1.51% to 0.42%.
Table A1 shows the UA, PA and
K values from each class throughout the analyzed period. User accuracy and producer accuracy show remarkable performance, with values within 90–100% and 80.48–100%, respectively. These individual accuracies, PA and UA, represent how well referenced pixels of the ground cover class are classified and the probability that a pixel classified into a given category actually represents that category on the ground, respectively. As expected, water and its classes showed the best performance, as it is the macro-class that presented the least variation while the other classes presented more significant variability. The Kappa coefficient shown for all classes was located in the range of excellent, according to the suggested values by Viera and Garret [
66]. One can expect from these results that all the classes provide at least 90% confidence to the users of this information. The mean UA and PA values are presented in
Table 10. Both UA and PA also are greater than 90%, as shown in each of the selected classes. The overall precision of the maps is also above 90%, which indicates that more than 90% of the reference pixels were correctly classified, and the
K values validate the quality of this study, as they are all above 0.95.
3.4. Land Use and Land Cover Change through LULC Transition Matrix
The last part of the analysis of these results includes the generation of the LULC transition matrix (
Table A2), which indicates the transitions of each of the classes with respect to each other, for instance, how much area of the dunes converted to low vegetation density areas [
69]. This matrix summarizes the changes already validated by the precision matrix but in terms of actual surface area. One disadvantage of this analysis is that it does not consider what happened throughout the two isolated years. However, it brings up an excellent tool for studies where LULC changes want to be predicted for future scenarios along with the possible drivers forcing those changes. The main diagonal contains the reference areas of the classes between the analyzed years. Columns named loss and gain represent the area lost or gained, respectively, for a particular class. Total gain or loss are calculated by adding up columns or rows of the reference area for each class. Particularly, this study selected the ends of the period, the years 2015 and 2021. The analysis of this matrix is straightforward. For instance, the water bodies gained 0.24 km
but lost 0.0297 km
. Gains were associated with transitions of wetland (0.22 km
), shrubland (0.0063 km
), mixed forest (0.0009 km
) and bare surfaces (0.0081 km
) to water bodies. On the other hand, losses are due to a mixture of transitions, including sparse vegetation (0.009 km
), bare surfaces (0.009 km
) and wetlands (0.00117 km
). All these transitions occurred between the two selected classifications.
For the sake of simplicity, only the most significant changes are discussed here. Residential land use grew only 0.88 km
during the last six years, which reflects a very small growth in comparison with that of industrial cities in Mexico, which has been characterized as about 5% per year [
80], while Coatzacoalcos only grew 1% per year. Commercial areas reduced by 0.16 km
, resulting in abandoned areas due to the previous mentioned socio-economic issues. One of the prominent changes was in areas of mixed forest, which gained 0.52 km
and lost 0.43 km
, indicating reforestation practices and more sustainable urbanization development in some new residential, industrial and commercial areas but some activities of deforestation. Swamps, on the other hand, lost 2.09 km
. This loss warns us of the possible future loss of ecological services provided by wetlands and their habitat, which are essential for biochemical processes and water purification because some of the non-point pollution in the city is contained by them. Lastly, vegetation densities fluctuated the most among them, namely, grassland and sparse vegetation density areas transitioned to bare surfaces. The latter gained 9.12 km
, whereas the former two lost 2.48 and 4.66 km
, respectively, which represents 80% of the bare surfaces gain.