Urban Land Mapping Based on Remote Sensing Time Series in the Google Earth Engine Platform: A Case Study of the Teresina-Timon Conurbation Area in Brazil

: Teresina-Timon conurbation (TTC) area is an example of urban agglomeration, situated in the semiarid environment of the northeast region of Brazil, which has shown an accelerated process of urban development over the last four decades (1985–2019). In this study, we developed a semi-automatic urban land mapping framework at the Google Earth Engine (GEE) platform to (a) evaluate spatiotemporal sprawl of the TTC area (1985–2018); and (b) quantify current urban fabric structures of TTC area (2019). The main empirical results demonstrate that the use of the Landsat historical dataset is a suitable option for generating consistent urban land maps across the years in semiarid environments. Teresina and Timon expanded, respectively, from 70.34 km 2 and 12.20 km 2 in 1985 to 159.02 km 2 and 30.68 km 2 in 2018, increasing annually at 3.05% and 3.69% averaged rate, showing an underlying tendency of continuous growth, and magnitude similar to Asian cities. The results of the urban fabric (UF) structures mapping demonstrates a high complexity of the urbanized surfaces, characterized by irregular shapes and variability of urban coverage. In 2019, the TTC metropolitan area was covered by urban land use classes as ceramic roofs, other types of roofs, and impervious surface, in the proportions of 28.02%, 11.97%, and 5.67%, respectively. roofs, and Impervious surface.


Introduction
Urbanization has permanently transformed ecosystems, causing massive amounts of carbon dioxide (CO 2 ) emissions, land degradation, air pollution, social inequalities, and diminishing human health. The process of deforestation and native vegetation clearing in urban settlements had accelerated mechanisms of net carbon losses from the natural vegetation, biodiversity disturbance, and soil rarefaction [1][2][3][4]. In 2018, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) identified urban expansion as one of the significant drivers of land-use change [5], contributing to the loss of ecosystem services and fall of forest-based livelihoods in peri-urban areas. Combined with population growth, rapid processes of urban sprawl (US) are pointed worldwide as a crucial challenge for human well-being on Earth.
Between 2000 and 2014, urban areas sprawled globally faster than the world population, resulting in a fall of inhabitant densities in cities, which poses an additional threat to environmental sustainability. In this context, the Sustainable Development Goal (SDG) associated with sustainable cities prescribes that to guarantee sustainable urbanization, we need to improve our understanding of managing urban expansion [6,7]. Besides, the United Nations New Urban Agenda (UN-Habitat) established the necessity to optimize the spatial dimension of the urban form and extension, prioritizing urban renewal [8]. On the other hand, the discussion regarding the optimal urban design, its sustainable population

Study Area
The two adjacent cites of Teresina, in the Brazilian state of Piauí (PI), and Timon, in Maranhão (MA), have been growing at an accelerated rate over the last few decades. Teresina-Timon conurbation (TTC) area is one of the three Brazilian Integrated Economic Development Regions (IEDR) of bi-state urban areas and more than one federative unit. The Brazilian IEDR's are similar to metropolitan areas, though, the individual cities at the IEDR's are under distinct political-institutional arrangements. Teresina and Timon are connected by three bridges, across the 300 m wet border of the Parnaíba River. These cities share flows of commodities, capital, and inhabitants every day, and their urban dynamics over the years is an interdependent system.
The TTC urban agglomeration is located in the semiarid region of the northeast of Brazil, placed among the Cerrado and Caatinga biomes ( Figure 1). This ecotone zone comprises natural vegetation characterized by a mix of semi-humid forest and savannas, with most of them degraded by mosaics of agriculture or pasture. Teresina-Timon conurbation area has a semiarid climate with recurrent droughts, which historically had motivated major waves of rural-urban migration [12,27]. The total population estimated in 2019 was equal to one million inhabitants, with averaged population densities varying from 89. 18 persons per km 2 in Timon municipality to 584.94 persons per km 2 in Teresina (Table 1). In terms of the total population, Timon is larger than the second most populated city in the state of Piauí, after Teresina. Both cities have most of the urban population concentrated into their continuous fabric, circumcised by an expressive rural area.  Each one of these two cities has implemented different urban strategies over the past several years, although both cities have been following the main dimensions of the 2001 City Statute, the legal framework governing urban development and management in Brazil [28]. Teresina, as the capital of Piauí, and a major center for commerce, industry, and Remote Sens. 2021, 13, 1338 4 of 16 services, concentrates the majority of efforts on urban municipal legislation, regarding the delimitation of urban districts and guidelines for urban occupation, recently increasing its attention to environmental sustainability. The building areas in Teresina are much more concentrated than in Timon, and the arborization in the streets in the Teresina downtown are, as a result, also less evident ( Figure 2) [27,29,30]. However, despite the interdependency among these two cities, there are still little efforts to understand the TTC urban dynamics as a unit. We lack information on the timing and magnitude of their changes, and although the Human Development Index (HDI) among these two cities is relatively close (Table 1: 0.751 in Teresina, and 0.649 in Timon), the percentage of wastewater sanitation network, for example, varies mostly ( Each one of these two cities has implemented different urban strategies over the past several years, although both cities have been following the main dimensions of the 2001 City Statute, the legal framework governing urban development and management in Brazil [28]. Teresina, as the capital of Piauí, and a major center for commerce, industry, and services, concentrates the majority of efforts on urban municipal legislation, regarding the delimitation of urban districts and guidelines for urban occupation, recently increasing its attention to environmental sustainability. The building areas in Teresina are much more concentrated than in Timon, and the arborization in the streets in the Teresina downtown are, as a result, also less evident ( Figure 2) [27,29,30]. However, despite the interdependency among these two cities, there are still little efforts to understand the TTC urban dynamics as a unit. We lack information on the timing and magnitude of their changes, and although the Human Development Index (HDI) among these two cities is relatively close (    Figure 3 presents the proposed framework to generate urban land maps. Our framework is virtually the same for our two major procedures: (I) classification of yearly mosaics of Landsat (LS) collection from 1985 to 2018; and (II) classification of the yearly mosaic of Sentinel-2 (S2) collection in 2019. Regarding the territorial coverage, we classified the entire two municipalities in (I); and the two adjacent urban lands in (II). In addition to the selection of different datasets, we also made adaptations on each procedure regarding the thematic land use and land cover (LULC) classes and the sample points used ( Table 2). Our workflow procedure was based in Rosa [31] and Alencar et al. [32].
We produced yearly best-pixel mosaics by merging pixels of distinct images collected from July 1st to September 30th [32], which covers the dry season. The date filtering returned a range of five to ten LS scenes for each path/row (219/063 and 219/064) in each year. We applied cloud and cloud shadow mask functions with a quality control band for each image, which returned pixels with minimum cloud coverage and better atmospheric conditions ( Figure 3).
For each band, we used the GEE median function in each selected pixel-location into

Shuttle Radar Topography Mission (SRTM) Slope (in degrees)
We produced yearly best-pixel mosaics by merging pixels of distinct images collected from 1st July to 30th September [32], which covers the dry season. The date filtering returned a range of five to ten LS scenes for each path/row (219/063 and 219/064) in each year. We applied cloud and cloud shadow mask functions with a quality control band for each image, which returned pixels with minimum cloud coverage and better atmospheric conditions ( Figure 3).
For each band, we used the GEE median function in each selected pixel-location into the period window. The resulting value was allocated for the final respective pixel-location in the annual mosaic. The median value is aimed at data dimensionality reduction and eliminating pixels corrupted by noise caused by clouds or cloud shadows not removed in the initial mask. Table 3 shows the years considered for filtering each data collection.

Landsat and Sentinel-2 Data Processing
Additional input variables were included in the procedure workflow, as shown in Figure 3 and Table 4. From the Landsat and Sentinel-2 datasets, we added values from the: normalized difference vegetation index (NDVI), the median of enhanced vegetation index 2 (median EVI2) [36], range of variation of EVI2 relative to the difference between the minimum and maximum values (amplitude EVI2), normalized difference built-up index (NDBI), and spatial texture computed by the entropy GEE function for NDBI values using a 5-pixel kernel centered (entropy NDBI). Furthermore, we added the slope value from the shuttle radar topography mission (SRTM). We resampled and regularized the range validity of values before inputting these new bands.
Rosa [31] and Alencar et al. [32] have shown that the addition of variables from remote sensing indexes and topographic data has contributed to the increment of the overall accuracy of classification results. For example, the amplitude EVI2 allows to better differentiate targets that remain stable throughout the year (buildings, roads, or industrial areas) from targets that vary during the year (shadows, vegetation, agriculture, or water). The textured band (entropy NDBI) allows a differentiation by the surrounding pixels, adding context information, and it helped to classify more homogeneous regions, reflecting a good generalization that is expected for classifications of urban areas.
For the classification processing, we selected the optimal parameters for the random forest (RF) algorithm from the literature (Figure 3) [37], expecting to get the highest overall classification accuracy possible. Several studies have used the number of trees equal to 100, and the number of variables in each split equal to the square root of the number of total variables [31,32,38].
Our training and validation sample data was collected based on the remote sensing visual interpretation of high-resolution imagery available in Google Earth. Based on this visual interpretation, we created polygons for each thematic land use category in areas that had not changed along the years. 80% of the total sampling in each thematic class was selected randomly from these polygons and used for training the RF algorithm, and 20% of the sample data was used for validation, considering the values described in Table 2. To assess the accuracy of classification, we used the two most popular metrics in the literature: overall accuracy (OA) and Kappa coefficient (KC), derived from the confusion matrix of the random sampling. The OA and KC coefficients were obtained for each year.

Urban Dynamics Metrics
As suggested by Espindola et al. [12] and Wu et al. [39], we calculated two indexes to quantify the magnitude of urban expansion: the annual increase (AI) and the normalized annual growth rate (AGR). AI (pixels per year) and AGR (%) is the annual growth extent and the annual growth rate of the urban class, respectively. P end and P start are the numbers of pixels labeled as urban at the end and start periods, respectively, and d is the time interval in years.

Spatial Dynamic of Urban Sprawl
The TTC area experienced notable urban expansion between 1985 and 2018. To systematically evaluate the spatial dynamics of urban sprawl over these last four decades, we present the urban land maps generated by the processing of the Landsat historical dataset. the consistency of the LULC classification among these sensors and across the years. We see that the proportion of the urban area has increased in both municipalities, ranging from 5.00% (1985) to 11.31% (2018)  correspond to carbon stocks of 68.6 t C/ha [40][41][42]. As an outcome, urban sprawl directly impacted the carbon and water cycles in the TTC area.   Here, all classification results showed high overall accuracies (OA) and Kappa coefficients (KC), ranging from 90% to 95%. The combined LULC map annually demonstrates that the urban land continued to increase and expand to suburban areas. Besides, over these 33 years, Teresina and Timon have undergone rapid processes of land consumption in almost every direction from the original urban core. Our results demonstrate the feature and the trends of the urban landscape in the TTC area, as these two cities presented a monotonic increasing trajectory. In Figure 6, we analyze the main directions of the sprawling in the TTC area in more detail. Although we have a significant increase in all trends, the expansion movement to the south (S) and east (E) directions was prominent. The center of these directions is the point labeled as 3 in Figure 2. Our results also demonstrate that the increment of the urban area was based on the replacement of the natural vegetation (Savanna and Forest), which has decreased considerably since 1985. Between 1985 and 2018, 161.31 km 2 of Forest and 39.75 km 2 of Savanna were lost. In Teresina and Timon municipalities, the predominant savanna vegetation is characterized by tree-shrub stratum, with a canopy cover of around 60%, and an averaged plant biomass and carbon stocks of 39.9 t C/ha [40]. This vegetation type differs from the forest canopy, characterized by more dense vegetation with relatively large trees, which correspond to carbon stocks of 68.6 t C/ha [40][41][42]. As an outcome, urban sprawl directly impacted the carbon and water cycles in the TTC area.
To further understand the annual patterns of urban sprawl (US) in Teresina, Timon, and in the TTC area as a unit, we show in Figure 5 the annual and continuous movement of expansion of their urban land, presenting the urban land coverage for 1985 (L5), 2012 (L7), and 2018 (L8), and each annual urban land map combined. Here, all classification results showed high overall accuracies (OA) and Kappa coefficients (KC), ranging from 90% to 95%. The combined LULC map annually demonstrates that the urban land continued to increase and expand to suburban areas. Besides, over these 33 years, Teresina and Timon have undergone rapid processes of land consumption in almost every direction from the original urban core. Our results demonstrate the feature and the trends of the urban landscape in the TTC area, as these two cities presented a monotonic increasing trajectory. In Figure 6, we analyze the main directions of the sprawling in the TTC area in more detail. Although we have a significant increase in all trends, the expansion movement to the south (S) and east (E) directions was prominent. The center of these directions is the point labeled as 3 in Figure 2.     Moreover, from Figure 7A-C, we can draw that, although the increase was constant in the TTC area, the speed of growth varied over the years. We present the magnitude of urban sprawl revealed by the annual increase (AI, in pixels) metric and normalized annual urban growth rate (AGR, %) for Teresina, Timon, and the TTC area as a unit. As the proportion of the urban area of Teresina is higher in the TTC area, the pattern of the AI in the graph ( Figure 7A) is similar to the pattern presented for Teresina. In 2018, 84% of the urban land of the TTC area was in Teresina, and 16% in Timon. Then, the most massive AI occurred during the 2017-2018 period for all three areas, followed by 1998-1999 and 1993-1994. Furthermore, Figure 7A shows a positive increasing trend in the 2006-2016 decade.
Moreover, from Figure 7A-C, we can draw that, although the increase was constant in the TTC area, the speed of growth varied over the years. We present the magnitude of urban sprawl revealed by the annual increase (AI, in pixels) metric and normalized annual urban growth rate (AGR, %) for Teresina, Timon, and the TTC area as a unit. As the proportion of the urban area of Teresina is higher in the TTC area, the pattern of the AI in the graph ( Figure 7A) is similar to the pattern presented for Teresina. In 2018, 84% of the urban land of the TTC area was in Teresina, and 16% in Timon. Then, the most massive AI occurred during the 2017-2018 period for all three areas, followed by 1998-1999 and 1993-1994. Furthermore, Figure 7A shows a positive increasing trend in the 2006-2016 decade.
After removing the effect of the city size, the TTC area, Teresina and Timon had an average AGR of 3.12%, 3.05%, and 3.69% , ranging from 0.30% to 12.31%, 0.51% to 11.74%, and 0.17% to 15.38% in 2018-1985 period, respectively ( Figure 7B). From Figure 7C, we note that during the three decades selected, the AGR's were almost constant across the two cities and in the conurbation area. These numbers are similar in magnitude to the ones found by Wu et al. [39] in Asian cities, such as Beijing, which had an average AGR of 3.70% in the 1980-2010 period. In another study conducted by Fu et al. [43], the annual expansion rate in another cluster of cities in China reported was between 4% and 10% during the 2000-2010 period.

Characterization of the Urban Fabric
In Figure 8, we illustrate the 2019 LULC classification generated from the Sentinel-2 dataset. For this result, we got an OA and KC of 0.94 and 0.91, respectively. To characterize the urban fabric (UF) in the TTC area, we present the proportion of the LULC classes directly associated with intra-urban characteristics: Residential-Ceramic roofs, Residential-Other roofs, and Impervious surfaces. In the TTC area, 28.02% of the urban land was covered by ceramic roof tiles, and 11.97% by other roofs, including roofs of concrete, metallic, or asbestos. Usually, red ceramic roof tiles are associated with households, and our results show a predominance of this material with a high built-up density across Teresina (30.84%) and Timon (21.38%). On the other hand, other types of roofs are also frequent, and in this area, they are mostly associated with commercial buildings or apartments. This coverage shows a considerable difference among Teresina (15.42%) and Timon (3.85%). Indeed, part of the recent urban development of Teresina was based on vertical growth, especially in its high-income district.
In the literature, impervious surface cover is usually used as an indicator of the builtup urban land density [44]. We defined impervious surfaces as the surfaces that water After removing the effect of the city size, the TTC area, Teresina and Timon had an average AGR of 3.12%, 3.05%, and 3.69% (2018-1985), ranging from 0.30% to 12.31%, 0.51% to 11.74%, and 0.17% to 15.38% in 2018-1985 period, respectively ( Figure 7B). From Figure 7C, we note that during the three decades selected, the AGR's were almost constant across the two cities and in the conurbation area. These numbers are similar in magnitude to the ones found by Wu et al. [39] in Asian cities, such as Beijing, which had an average AGR of 3.70% in the 1980-2010 period. In another study conducted by Fu et al. [43], the annual expansion rate in another cluster of cities in China reported was between 4% and 10% during the 2000-2010 period.

Characterization of the Urban Fabric
In Figure 8, we illustrate the 2019 LULC classification generated from the Sentinel-2 dataset. For this result, we got an OA and KC of 0.94 and 0.91, respectively. To characterize the urban fabric (UF) in the TTC area, we present the proportion of the LULC classes directly associated with intra-urban characteristics: Residential-Ceramic roofs, Residential-Other roofs, and Impervious surfaces. In the TTC area, 28.02% of the urban land was covered by ceramic roof tiles, and 11.97% by other roofs, including roofs of concrete, metallic, or asbestos. Usually, red ceramic roof tiles are associated with households, and our results show a predominance of this material with a high built-up density across Teresina (30.84%) and Timon (21.38%). On the other hand, other types of roofs are also frequent, and in this area, they are mostly associated with commercial buildings or apartments. This coverage shows a considerable difference among Teresina (15.42%) and Timon (3.85%). Indeed, part of the recent urban development of Teresina was based on vertical growth, especially in its high-income district.
In the literature, impervious surface cover is usually used as an indicator of the builtup urban land density [44]. We defined impervious surfaces as the surfaces that water cannot infiltrate. They are usually made up of anthropogenic materials, e.g., streets and roads, rooftops, parking lots, and outdoor facilities [45]. Although roof tiles may also be included in this category, we differentiate roofs in other thematic classes, aiming to distinguish places without population settlement. Even showing a similar proportion of the impervious surface in Teresina (6.45%) and Timon (3.83), the visual interpretation of our results ( Figure 8B,C) demonstrates that in Timon, besides being less recurrent, this land use category is arranged with natural vegetation class. This category also corresponds to 5.67% of the total urban land in the TTC area.
The impervious surface category is relevant when associated with the urban land surface temperature. Impervious surface areas are reported to yield a temperature of around 10 degrees higher when compared with forest coverage [46]. Then, mapping this category is important for urban planning and environmental and resources management. As a result, the pattern of the land use found in Timon may contribute to an urban design that enhances the thermal comfort of the population [47].
The impervious surface category is relevant when associated with the urban land surface temperature. Impervious surface areas are reported to yield a temperature of around 10 degrees higher when compared with forest coverage [46]. Then, mapping this category is important for urban planning and environmental and resources management. As a result, the pattern of the land use found in Timon may contribute to an urban design that enhances the thermal comfort of the population [47].

Residential-Ceramic roofs
Residential-Other roofs Impervious surface

Discussion
Understanding the spatial dynamics of urban development is relevant in conurbation cities under distinct political-institutional arrangements, as in the case of the Teresina-Timon area. Most remote sensing datasets are currently freely available, giving us a vast panorama of how twins cities evolved and allowing more effective urban planning. Additionally, cloud computing platforms are revolutionizing the ways we deal with image processing features. However, frameworks that integrate freely available remote sensing data into cloud computing platforms designed for urban land monitoring are still scarce. Our results showed the robustness of Landsat and Sentinel-2 imagery for urban planning and management. The use of Landsat historical imagery is well documented in the literature for monitoring urban sprawl [48][49][50][51]. On the other hand, Sentinel-2 satellites are still being evaluated [52][53][54], and we demonstrated that the urban land classes selected for this study were well classified. The ability of Sentinel-2 imagery to distinguish different types of roofs opens new possibilities for modeling populational densities in the urban land. The landscape metrics used also promoted a deeper understanding of the environmental sustainability of TTC growth.
Besides, we showed that the urban land in the TTC area presented an underlying tendency of continuous growth. Analysis results from this study suggest that, compared with other cities in Latin America, this growth has been horizontal rather than vertical. A basic finding from this study is that the TTC area's sprawling pattern is still more relevant than the urban fabric's densification. During the last four decades, a large amount of vegetation removal and land consumption in the TTC area imposed a critical threat to urban sustainability, as this type of land-use change is often associated with a rapid decline in local ecosystem services. In this scenario, the opening of new spaces to urban development in Teresina and Timon cities had central environmental implications, especially knowing that the main drivers of the urban fringe sprawl are usually more related to economic and agrarian factors, rather than related to the well-being of the population. The landscape analysis provided with our results summarizes the size, shape, and spatial distribution of the urban sprawl dynamics in the TTC area, reinforcing that the observed patterns are a result of complex environmental and socio-political interactions. As a conurbation area located along rivers, thus the necessity of balancing urban development and water quality should also be stressed.
In the TTC area, urban dynamics were influenced by the diversification and growth of national and local economies and by population increase. This area presented cycles of heavy urban expansion, followed by periods of sparse growth. These cycles were directly related to the importance of the federal government's role in financing urban infrastructure. In our region, the federal grants typically followed the form of political and economic constraints.
During the 1990s, the population grew at a breakneck pace, accompanying a movement that was reported in other medium-sized cities in Brazil and Latin America [55][56][57][58][59]. After the 1990s, the TTC area sprawl was partially motivated by federal government programs that implemented national housing plans aimed to eradicate the large housing deficit in the country and focusing on low-income families.
After the 2008 economic crisis, massive investments in infrastructure displaced the population to peripheral zones without adequate urban planning, which reflected in alteration in land use and land cover (LULC), followed by environmental impacts and public health issues caused by thermal discomfort, particularly in semiarid regions. The implementation of the national housing policy in TTC contributed to the current shape of the urban land, as most of the housing settlements were in peripherical zones of Teresina and Timon. Nevertheless, the investments in urban infrastructure were never compatible with the amount of the population relocated to the new peripheral areas.
Despite empirical works that address the sprawl of Latin America medium-cities have been scarcer, commonly emphasizing only the major metropolitan areas such as São Paulo and Rio de Janeiro in Brazil, our results corroborate with Inostroza et al. [60]. Our findings show that, as a typical urban area of Latin America, TTC expressed spatial dynamics characterized by high rates of land consumption and fragmentation degrees. The TTC area manifested the same pattern of the radial growth of other Latin American cities, characterized by a low-density peripherical expansion with a lack of infrastructure, motivated by poorly regulated land use [12,[61][62][63][64]. Usually, this overall tendency of sprawling growth also contributes to the reduction of residential activities in well-serviced central areas, willing to the socio-spatial segregation [65,66].

Conclusions
Urbanization processes are one of the major drivers affecting ecosystem functioning and services locally. In semiarid environments, urban expansion contributes to the recurrence of extreme climate events that may further impact the landscape. This study aimed to improve the understanding of urbanization processes in one conurbation area in Latin America, using cloud computing platforms and free remote sensing data. Currently, the Landsat historical dataset is the most used source to track continuous regional land-use changes back to the 1980s. However, Sentinel-2 is arising as an alternative for more detailed analysis. The high-performance of the GEE computing platform makes the manipulation of a large amount of remote sensing data more user-friendly and suitable for urban assessments and management.
Then, we proposed a semi-automatic framework to map urban land in the TTC area. The produces maps resulting in spatially explicit inventories of the urban extents and conditions. The results demonstrated that the urban agglomeration in the TTC area experienced accelerated growth from 1985 to 2019. Despite our promising results, it is still necessary to incorporate these approaches into the mainframe of urban governance that accommodates each city's political-institutional arrangement. The joint governance of Teresian and Timon based on the use of the proposed framework would allow a more dynamic urban management by the creation of integrated urban public policies. As a future investigation, we suggest additional studies on the new dynamics of medium-sized cities in the context of urban sprawl and planning policies. We need to investigate deeper the main drivers of urban development in Latin America in contrast with distinct local, national, and international panoramas.
Author Contributions: E.C. analyzed the data and performed the experiments and computed the remote sensing analysis. W.L. and G.E. supervised the conception and design of the analysis; and worked on the final manuscript. All authors developed and discussed the manuscript together and finally wrote the paper. All authors have read and agreed to the published version of the manuscript.
Funding: UFPI and IFPI partially funded this research paper.