Application of Remote Sensing Tools to Assess the Land Use and Land Cover Change in Coatzacoalcos, Veracruz, Mexico

Land use and land cover (LULC) change has become an important research topic for global environmental change and sustainable development. As an important part of worldwide land conservation, sustainable development and management of water resources, developing countries must ensure the use of innovative technology and tools that support their various decision making systems. This study provides the most recent LULC change analysis for the last six years (2015–2021) of Coatzacoalcos, Veracruz, Mexico, one of the most important petrochemical cities in the world and host of the ongoing Interoceanic Corridor project. The analysis was carried out using Landsat 8 Operational Land Imager (OLI) satellite images, ancillary data and ground-based surveys and the Normalized Difference Vegetation Index (NDVI) to identify and to ameliorate the discrimination between four main macro-classes and fourteen classes. The LULC classification was performed using the maximum likelihood classifier (MLC) to produce maps for each year, as it was found to be the best approach when compared to minimum distance (MDM) and spectral angle mapping (SAM) methods. The macro-classes were water, built-up, vegetation and bare soil, whereas the classes were an improved classification within those. Our study achieved both user accuracy (UA) and producer accuracy (PA) above 90% for the proposed macro-classes and classes. The average Kappa coefficient for macro-classes was 0.93, while for classes it was 0.96, both comparable to previous studies. The results from the LULC analysis show that residential, industry and commercial areas slowed down their growth throughout the study period. These changes were associated with socioeconomical drivers such as insecurity and lack of economic investments. Groves and trees presented steady behaviors, with small increments during the five-year period. Swamps, on the other hand, significantly degraded, being about 2% of the study area in 2015 and 0.93% in 2021. Dunes and medium and high vegetation densities (∼ 80%) transitioned mostly to low vegetation densities. This behavior is associated with rainfall below the annual reference and increments of surface runoff due to the loss of vegetation cover. Lastly, the present study seeks to highlight the importance of remote sensing for a better understanding of the dynamics between human–nature interactions and to provide information to assist planners and decision-makers for more sustainable land development.


Introduction
Land use and land cover (LULC) changes are among the research topics more recurrent in remote sensing [1][2][3][4]. Remote sensing provides comprehensive and extensive information to understand the interaction between terrestrial ecosystems and their responses to environmental factors [5][6][7][8][9][10][11]. This technology is considered a powerful source to obtain information from terrestrial surface characteristics at different temporal and spatial scales [12][13][14][15][16]. For more than half a century, scientists have used several types of remote sensing data. From those, the most commonly used source is the Landsat satellite based-data [17]. One of main uses of the Landsat-based satellite images is the identification or classification of LULC changes [18][19][20][21][22]. Some interesting conclusions have been driven through the use of remote sensing in LULC changes. It has been observed that changes in land uses are dynamic, non-linear interactions between humans and nature and are driven by complex stochastic processes [23]; worldwide changes in LULC for the last 300 years have shown gains in agriculture and losses in forest [24]; rapid population growth in Africa has been identified as the main driving force for the expansion of agricultural land, whereas in developing countries, urbanization dynamics are attributed to demographic factors [25] and transitions from cropland to urban and water bodies areas have been identified as an important trend in Europe [26].
Some of the more remarkable research found in literature is the work conducted by Yuan et al. [27]. Their research consists of the use of multi-temporal Landsat Thematic Mapper data to map and quantify the LULC changes in seven counties of the Twin Cities located in the metropolitan area of the state of Minnesota, United States. They used satellite images corresponding to the years 1986, 1991, 1998 and 2002. Yuan et al. demonstrated an urban zone increment from 23.7 to 32.8% in the study area. Surfaces such as rural areas, croplands, forests and wetlands decreased from 69.6 to 60.5%. Another unique example is the work carried out by Demissie et al. [28], where Landsat satellite data from 1973 to 2015 were classified to study land use change and its possible causes in Gonder, Ethiopia.
The study showed that about 60.1% of the area experienced land use changes within the study period.
In Mexico, little research related to the quantification of LULC changes has been found within the literary review. One of the few works found was carried out by Colditz et al. [29]. Their research presents a methodology to develop a land use map in Mexico for the year 2005. The scheme was based on time series from a Moderate Resolution Imaging Spectroradiometer (MODIS) with a 250 m resolution and an extensive sampling of data for the different geographic zones of the Mexican surface. The results showed a map with an overall precision of 82.5% (Kappa statistic = 0.79). A further evaluation with 780 randomly generated samples within the classes with referenced field data indicated a precision equal to 83.4% (Kappa statistic = 0.80). Another study found in the literature includes the generation of a land use map of the Latin American and Caribbean zone for the year 2008. This project was developed within the framework of the project of the Latina American Network for Monitoring and Studying of Natural Resources (SERENA) [30]. Similar to Colditz et al. [29], time series and decision trees with MODIS data (250 m) were used for the classification of land use changes. The discrete SERENA model showed an overall accuracy of 84%. Other uses of remote sensing in Mexico include the detection of chinampas in the Xaltocan area in the Northern Basin of Mexico [31], the evolution of sea temperature on the west coast of Baja California [32] and/or hydropower assessment [33].
Due to the limitations observed in the use of remote sensing in Mexico, such as limited research using remote sensing and map course resolution, this study has integrated remote sensing, Geographic Information Systems (GIS) analysis, field sampling and image revision to assess the potential of the use of Landsat 8 satellite data to assess and monitor land use changes and coverage changes in the southeastern region of the state of Veracruz, Mexico, located in the southern coastal area of the Gulf of Mexico. Three widely used classification methods are explored to select the more appropriate algorithm for LULC classification under the same level of information previously extracted from field data, digital maps and ancillary data. This comparison allows for conclusions on the advantages and drawbacks of current classification algorithms. The analysis encompasses the definitions of macro-classes of land use and their derived classes. The Normalized Difference Vegetation Index, also known as NDVI, was used to identify and improve discrimination between macro-classes and classes within the study area. A supervised classification method, the Maximum Likelihood Classification (MLC), was applied in this study due to its availability within the GIS applications in addition to it not requiring an extensive training process [34,35]. The most important advantage of the MLC as a parametric classifier is that it takes into account the variance within the selected classes and that, for normally distributed data, the MLC method performs a better qualification than other methods, namely, decision trees [34,36].
Additionally, this study proposes a methodology to obtain updated temporal information of the change in land use and land cover for the Coatzacoalcos area using Landsat 8 satellite imagery, which has a higher resolution with the available MODIS maps. The mapping of land use and its changes are critical and important factors for sustainable development as well as the monitoring of environmental impacts. The final product of this research will help government agencies to make decisions on urban development and preservation of available natural resources, since one of the emblematic projects in Mexico known as the Interoceanic Corridor is being executed in the present year. Lastly, this piece of research is going to serve as the baseline for additional contributions on LULC change prediction with more complex analyses of the drivers forcing these landscape modifications.

Study Area
The study area is the southeastern region of the state of Veracruz, Mexico, in which one of the most important pretrochemical ports of the country is located at the city of Coatzacoalcos. This city is surrounded by the Coatzacoalcos and Tonalá rivers and extends between latitudes 18.06 • and 18.21 • and longitudes −94.22 • and −94.64 • (Figure 1). The area of study extends over an area about 220 km 2 . The port of Coatzacoalcos is dominated by the petrochemical sector [37]. Four of the petrochemical complexes located near the city make the city one of the most important oil areas in the world [38]. Coatzacoalcos has been a transportation hub for hundreds of years; it is connected by air, land, sea and railroad to the rest of the world. Currently, the Mexican government is working on a project called the Interoceanic Corridor of the Isthmus of Tehuantepec, which represents the modernization and installation of industrial parks inside the industrial areas of the ports of Coatzacoalcos and Salina Cruz, Mexico [39]. This project will represent economic and population growth within the city. As a result, land use assessments and projections have become priorities to help both state and federal government make decisions for the projected economical growth.
Physiographically, the study area has an average elevation of 10 meters above sea level. The soils found in the region are clayey lateritic Acrisol and Luvisol in high elevations and Gleysol, Cambisol, Vertisol and Nitosol in plains [40]. The vegetation consists of riparian vegetation, swamps and wetlands in floodable areas, mesophilic mountain pine oak forests in high areas and high evergreen forests in hills and acahual areas in abandoned and cultivated pastures. The climate of the region is a warm humid climate with abundant rains in summer. The average annual temperature varies between 24 and 25°C, while the average annual precipitation varies between 1500 and 2500 mm [41].

Data Acquisition and Processing
The images used in this study were obtained from the Landsat 8 satellite, which captures multispectral images of 30 m resolution sampled with the OLI (Operational Land Imager) sensor. The Landsat 8 satellite also provides a 15 m panchromatic band, allowing for a greater spatial resolution. All images were geometrically corrected and acquired in level 1 (L1). The images can be directly downloaded from the EarthExplorer gateway of the United States Geological Survey (USGS) [17]. Table 1 displays the details of the images used in this study.  [42]. All images were processed using the free software Quantum GIS (QGIS) version 3.14.16. Atmospheric correction was performed through the Semi-Automatic Classification Plugin (SCP) tool based on the Dark Object Subtraction (DOS1) algorithm [43]. The panchromatic band was used to provide higher spatial resolution. This is achieved using Gram-Schmidt Pan-Sharpening to produce a high-resolution color image that improves the training set and the classification process. Google Earth images were acquired according to acquisition dates (Table 1) and geoferenced using QGIS. The closest dates to the acquisition date were selected. Field data included collection of georeferenced points and their corresponding land use for the last two years of the analyzed periods. Figure 2 shows the methodology used for the study.  The DOS1 atmospheric correction method was selected, as it is the most widely used method for LULC change detection [42]. This method uses the properties of the images. For instance, pixels conforming elements such as shadows, water or forests are considered dark objects when the magnitude of the reflectance is close to zero, namely, when the reflectance is less than or equal to 1%, and all the analyzed pixels are assumed dark elements. This assumption indicates that pixels receiving low values of solar radiation (∼100% shade) registered by the satellite correspond to atmospheric dispersion, which generally is caused by effects of the topography [44]. Once a dark object is found in the image, the minimum reflectance value of the total Digital Number (DN) histogram is assigned to the object. This minimum is then modified by the effects of atmospheric correction [45]. The surface reflectance is calculated using the following expression: where ρ s is the surface reflectance, d is the distance to the sun, ESUN λ is the mean solar exo-atmospheric irradiance, L λ is the solar spectral radiance to the satellite, θ s is the zenith angle and λ p is the path radiance given by [46] L p = L min − L DO1% where L min is the radiance corresponding to a digital count value for which the sum of all the pixels with digital counts lower than or equal to this value is equal to 0.01% of all the pixels from the image considered [46] and L DO1% is the radiance of Dark Object assumed to have reflectance value of 0.01.

NDVI Classification
Based on the interpretation of the Normalized Difference Vegetation Index (NDVI), the analysis of Google-Earth-extracted images, RGB (red, green, blue) images, panchromatic bands and ground-based georeferenced information, NDVI thresholds were proposed to classify the area of study into four macro-classes during the analysis period (2015-2021), as shown in Table 2. The NDVI thresholds implemented in the proposed classification were selected after an extensive literature review of previous studies where the NDVI was used as the spectral index of diverse land covers [47][48][49][50][51][52]. The NDVI values were calculated by the equation proposed in [53,54]: where φ nir is the near-infrared reflectance and φ red is the red-band reflectance. Band 4 (0.64-0.67 µm) and band 5 (0.85-0.88 µm) represent φ nir and φ red , respectively, in the Landsat products. The NDVI values vary from −1 to 1. Values close to 1 represent vegetation in optimal environmental conditions, whereas low values of NDVI indicate low vegetation density or a different land use. After the macro-class definition, a second classification was performed to refine the LULC characteristics through a more extensive discrimination of pixels. These new classes were defined according to the CORINE land cover inventory proposed by the European Environmental Agency [55]. This classification process required the use of the ancillary data and field data to assign new information to the training set. This assignation was repeated until a suitable spatial distribution was obtained. This new classification consisted of 14 classes, as shown in Table 3.

Classifier Comparison
In addition to its social and environmental objectives, this study seeks to evaluate the potential and the drawbacks of the minimum distance classifier [56], the spectral angle mapping classifier [57] and the maximum likelihood classifier [34,35,58] in deriving information on LULC. Although there are more classifications, such as artificial neural networks [34] and parallelepiped classification [56], the comparison was limited to the three selected methods because they are often used to classify LULC. However, limited information exists when the same training input and ancillary data are used to perform the classification process. Here, the comparison was set under the assumption that the same amount of information was available. The minimum distance method (MDM) calculates the Euclidean distance d(x, y) between the spectral signatures given in the training data set and the spectral signatures of the image pixels. The spectral distance is calculated through the following expression: where x is the spectral signature vector of an image, y is the spectral signature vector of the training area and n is the number of bands of the image. Once the spectral distance is computed for every pixel, the class with the closest spectral signature to the training set is assigned according to the following discrimination function [56]: where C k is the land cover macro-class or class, y k is the spectral signature of class k and y j is the spectral signature of class j. This equation is valid when k = j.

Spectral Angle Mapping Classifier
The spectral angle mapping (SAM) algorithm computes the spectral angle between the spectral signatures of the image pixels and the training spectral signatures. The spectral angle θ is given by Thus, a pixel resides in the macro-class or class that has the lowest spectral angle, as provided by where k = j.

Maximum Likelihood Classifier
The maximum likelihood classifier (MLC) is based on the probability that a pixel belongs to or is within macro-classes or a particular class [59,60].
The MLC algorithm calculates the weighted distance or probability D that an unknown value in the vector M p belongs to one of the macro-classes or classes M c . This likelihood is based on the Bayesian equation [61]: where D is the weighted distance or probability, c is a particular macro-class or class, M p is the measurement vector of the candidate pixel, M c is the mean vector of the sample of macro-class or class c, a c is the percentage probability of any pixel belonging to the macro-class or class c, Cov c is the covariance matrix of the pixels in the sampled macro-class or class c, |Cov c | is the determinant of the matrix Cov c , Cov c−1 is the inverse of Cov c and T is the transposition function.

Precision Assessment
The analysis of the accuracy is an important step for the evaluation of the resulting classification because the users of the information, once the classification is performed, need to know how accurate the result is in order to use the data in their decision making. A minimum level of overall precision for a selected macro-class or class of at least 85%, according to the recommendations found in similar works [49,60,62,63], was proposed in this study. The selected ratio between a training set and validation set was 70/30. This means that 70% of the available information was used for training purposes, whereas 30% was used for validation. Acceptable precision measurements used in this work included producer precision (PA > 85%), user precision (UA > 85%), overall precision (OA > 85%) and the Kappa coefficient (K) [49,64,65]. The reference values of the Kappa coefficient proposed by Viera and Garrett [66] are shown in Table 4. For precision evaluation, a set of verification points are required. The sample should be designed to achieve low standard errors in precision estimates, and this is generally achieved by random selection of points. The number of samples should be calculated by where N is the sample size, W i is the mapped area proportion of the class i, S i is the standard deviation of the stratum i and So is the expected standard deviation of overall accuracy, often valued at 0.01 [67][68][69].

Validation Point Estimation
The number of samples calculated with Equation (9) is shown in Table 5. The assumption of this expression claims that randomly generated validation points are proportional to the size of the selected class. In other words, the larger the class is, the more sample points are needed to verify the correct LULC assignation. For instance, as vegetation represented the larger area of the city and its surroundings, more points were needed to validate the land use classification. This behavior can be seen in Table 5 for the macro-classes in question. Similarly, the number of samples required to validate the classes presents the same behavior; the classes representing mixed forest (vegetation) and bare surfaces (bare soil) needed 79 and 88 samples, respectively, as they are the largest land uses in the study area. It is important to mention that pixels representing the sea area were not counted, as they did not change significantly during this period of time and we found no factor for any small increment of decrease that suggested any environmental effect on ties. However, it was taken into consideration to illustrate the applicability of the classification algorithm. On the other hand, rivers and water bodies might be subjected to changes due to water availability or variations in hydrologic processes.

Results
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.

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.

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 28.58% in 2017, whereas its minimum was 28.43% in 2019. Urbanization or man-made construction increased from 19.88% in 2015 to 20.21% 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.

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, microclimate 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 Figures 6 and 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.

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 2 but lost 0.0297 km 2 . Gains were associated with transitions of wetland (0.22 km 2 ), shrubland (0.0063 km 2 ), mixed forest (0.0009 km 2 ) and bare surfaces (0.0081 km 2 ) to water bodies. On the other hand, losses are due to a mixture of transitions, including sparse vegetation (0.009 km 2 ), bare surfaces (0.009 km 2 ) and wetlands (0.00117 km 2 ). 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 2 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 2 , 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 2 and lost 0.43 km 2 , 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 2 . 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 2 , whereas the former two lost 2.48 and 4.66 km 2 , respectively, which represents 80% of the bare surfaces gain.

Discussion
LULC changes measure the transitions of different land uses in complex interactions between humans and the physical environment [81,82]. Analysing LULC changes helps facilitate sustainable land use planning to protect and conserve the natural habitat and resources. The present study applied available remote sensing technology to classified course land uses, similarly to what has been presented in several studies [23,81,[83][84][85][86]. The macro-classes identified were water, built-up, vegetation and bare soil. This first division allowed us to discriminate NDVI spectral indexes to improve a latter classification. Fractions of the area of study divided into macro-classes help us to observe the evolution of the city and its metropolitan area development and to explore the possible drivers forcing the changes in this initial classification. Broadly, some socio-economic factors, climate and topography were proposed as possible drivers. These drivers have also been identified as important factors for land use evolution [23,[87][88][89]. This paper identified that security, precipitation, climate and topography were the main drivers causing LULC changes. Measurements such as producer accuracy (PA) and user accuracy (UA) showed high confidence, since they were higher than 90% for most of the macro-classes. These values are comparable and even superior to previous studies found in literature [49,84,89,90]. One can observe that the maximum likelihood classifier (MLC) is a suitable classification method to obtain acceptable results, as cited by several scientific contributions [34,35,58]. After the analysis of the macro-classes, 14 classes were exhaustively found to analyze and identify the land use in detail. Among those classes, some discussions are driven. Industry has experienced no change since the last petrochemical complex opened in 2015. Commercial areas have declined and residential zones present a slow increment due to the abovementioned socio-economic drivers, which served to confirm the proposed hypothesis. Mixed forest and bare surfaces with low vegetation density tended to increase in surface area because reforestation programs and vegetation implemented in the coastline by the government have been a priority, as wind erosion has provoked much damage to the city's infrastructure [91]. Swamps (i.e., wetlands) also showed an important decrease in spatial contribution, losing about 2 km 2 . This, in turn, might become an important environmental issue, since swamps are seen as important regulators of water pollution, carbon storage and habitat of species [92][93][94][95]. For that reason, actions need to be considered to avoid the continuous degradation of local wetlands and swamps. Additionally, it was observed that about 80% of the bare surface areas came from sparse vegetation and grassland areas. Associated changes in precipitation and high temperatures throughout the assessment period were responsible for these changes, as shown in results section. As a result of increments of bare surfaces, infiltration decreases and overland flow is exacerbated because soil properties such density, porosity and hydraulic conductivity depend on the level of established vegetation [96,97]. Lastly, a LULC transition matrix reflected the changes in surface area for each of the classes identified in this study. This is the baseline for future research that involves local problems in the city of Coatzacoalcos. A subsequent prediction of the land use using more extended spatial information of the drivers forcing the changes observed here will be proposed as an extension of this work. Additionally, the city is facing an issue with the relocation of solid residues, which naturally correlated with the actual land use. As a result, this study and methodology using higher-resolution imagery (i.e., sentinel satellite) help to study possible and suitable landfill sites.

Conclusions
This paper determined land use and land cover (LULC) over the period 2015-2021 for the city of Coatzacoalcos and its metropolitan area, located in the state of Veracruz, Mexico. Based on images from Landsat 8, MLC, Geographic Information Systems (GIS), NVDI spectral index and field data, the annual land use variability was produced up-todate as reliable information for decision making and ecosystem preservation during the execution of the ongoing project called the Interoceanic Corridor. This project represents Mexico's narrowest stretch between the Pacific and Atlantic oceans and the expansion and modernization of Coatzacoalcos and Salina Cruz ports.
The objective of this study was to provide the most recent information on land use spatial distribution that Coatzacoalcos experienced in the last six years and to improve the current available course resolution maps. The satisfactory results can be summarized in four main aspects: (1) built-up areas, including industry, residential and commercial areas, have slowed down their growth due to socio-economical drivers such as security and null monetary investments; (2) vegetation such as mixed forest and low density vegetation (bare surfaces, sparse vegetation and grassland) has been sustained and increased over time due to reforestation or migration from other classes; (3) swamps experienced considerable degradation over the past five years and (4) high and medium vegetation densities have transformed mostly to low vegetation densities due to climate drives such as low precipitation and possible high soil evaporation, which might also increase the overland flow for those areas.
Lastly, this study demonstrated the use of free available Landsat data and their processing by open source tools. It provided an accurate approach to mapping and assessing LULC changes over time. This methodology can be applied similarly for longer periods of time and other satellite products and contributes to improving the number of applications of remote sensing and research in Mexico and other developing countries.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments:
The authors would like to thank the anonymous reviewers for their valuable comments and feedback on this article.

Conflicts of Interest:
The authors declare no conflict of interest.