Recording Urban Land Dynamic and Its Effects during 2000–2019 at 15-m Resolution by Cloud Computing with Landsat Series

Cities, the core of the global climate change and economic development, are high impact land cover land use change (LCLUC) hotspots. Comprehensive records of land cover land use dynamics in urban regions are essential for strategic climate change adaption and mitigation and sustainable urban development. This study aims to develop a Google Earth Engine (GEE) application for high-resolution (15-m) urban LCLUC mapping with a novel classification scheme using pan-sharpened Landsat images. With this approach, we quantified the annual LCLUC in Changchun, China, from 2000 to 2019, and detected the abrupt changes (turning points of LCLUC). Ancillary data on social-economic status were used to provide insights on potential drivers of LCLUC by examining their correlation with change rate. We also examined the impacts of LCLUC on environment, specifically air pollution. Using this approach, we can classify annual LCLUC in Changchun with high accuracy (all above 0.91). The change detection based on the high-resolution wall-to-wall maps show intensive urban expansion with the compromise of cropland from 2000 to 2019. We also found the growth of green space in urban regions as the result of green space development and management in recent years. The changing rate of different land types were the largest in the early years of the observation period. Turning points of land types were primarily observed in 2009 and 2010. Further analysis showed that economic and industry development and population migration collectively drove the urban expansion in Changchun. Increasing built-up areas could slow wind velocity and air exchange, and ultimately led to the accumulation of PM2.5. Our implement of pan-sharpened Landsat images facilitates the wall-to-wall mapping of temporal land dynamics at high spatial resolution. The primary use of GEE for mapping urban land makes it replicable and transferable by other users. This approach is a first crucial step towards understanding the drivers of change and supporting better decision-making for sustainable urban development and climate change mitigation.

Driven by the demands of comprehensive records of urban land use land cover dynamics, this study aims to generate annual maps of urban LCLU at 15-m resolution by a multilevel decision rule method in GEE. This method would rely on a new Landsat pan-sharpening algorithm using GEE. The approach was applied to capture the annual change and turning point of LCLU in Changchun from 2000 to 2019. The LCLU maps were then used to further investigate the socio-economic drivers and environmental effects of LCLUC through statistical analysis. The implementation of our methods would help to acquire urban knowledge at a finer spatiotemporal resolution.

Study Area
Changchun, the major city in Northeast China, is the capital of Jilin province. The humid continental climate with the long cold winter governs the regional biomes, representing by temperate conifers. As lying in the continental plain of northeast Asia, Changchun became the modern city proto in China since the early 20th century due to railway construction by Russian and Japanese. As ordered by the central government, Changchun priorly developed heavy industrial and automobile manufacture from the 1950s. However, the economic depression in Changchun was impressed trail after the rapid development in coastal cities in recent decades. Winning the award of National Garden City in 2001 led Changchun to pay more attention to urban greenness, improving its livability and reputation. The built-up land and green space have filled Changchun ( Figure 1). The long-term city development made Changchun an important laboratory to study urbanization and environmental change [17,31].

Landsat Series and Pan-Sharpening in GEE
Landsat sensors have been providing multispectral images since 1972, serving earth observation at 30-60 m resolution for decades. Because lapped acquisition from sensors made dense observations available (Table 1), TM images were applied to fill the gap of ETM+ images' multispectral bands. The panchromatic band of the ETM+ and OLI sensors can apply to promote original bands to 15-m resolution, enhancing Landsat availability for urban study [26]. In GEE, to observe land dynamics, we employed Landsat top-of-atmosphere images (TOA) with the cloud cover <80% and the path/row of 128029, 128030, and 129030. The cloud-removing algorithm iterated all images, which were prepared to compose a yearly image by median stacking. A new algorithm of GEE was developed for Landsat pan-sharpening in this study. The JavaScript-based HSV method is the built-in function of GEE, which can sharpen three bands through one transformation that includes two-step algorithms of rgbToHsv() and hsvToRgb() [29]. Instead of using the HSV, we directly reassigned multispectral information to up-scaled pixels based on a sliding pixel-window, including five steps (1) stacking all bands in one image, (2) forming a mean array by sliding the window, (3) recording the variation from multispectral mean to panchromatic mean, (4) recording the variation of multispectral information from original to mean, and (5) resigning multispectral information into 15-m pixels by multiplying pan band and two variations in each pixel. By testing, the 5 × 5 (or 150-m length) pixel-window could reserve the original multispectral information, which is evidenced by a Pearson coefficient of more than 0.9 in each band ( Figure 2). As the pixel heterogeneity has been enhanced spatially, 15-m spectral indices could be obtained. This study pan-sharpened Landsat images by a 5 × 5 pixel-window that reserves the spectral information and enhances the spatial heterogeneity of pixels.
Spectral indices that include Normalized Difference Vegetation Index (NDVI) and Modified Normalized Difference Vegetation Index (MNDWI) were applied to extract LCLU pixels [32]. We further employed NDVI for phenology analysis by calibrating ETM+ and OLI images [33]. In addition, to avoid saturation of the NDVI, Fractional Vegetation Cover (FVC) calculated from NDVI was used for extracting vegetation.
where NDVI soil and NDVI veg denote the mean NDVI in bare ground and vegetation, respectively.

Integration of Automatic Thresholding and Multilevel Decision Rule Process
The urban LCLU category was defined by referencing peer-reviewed works [14,34,35] and visually investigating Google Earth images ( Table 2). The difference is a category of green space that contains three vegetation types of forest, grass, and shrub. Although individual vegetation classes were not classified, our scheme was meeting the research need in Changchun [17,31,35] and avoided commission error between vegetation types due to lack of structure definitions. A multilevel decision rule (MDR) algorithm was developed to generate urban LCLU maps by thresholding spectral index in this study. Despite machine learning classifiers are performing good LCLU classification [14,25], we aimed to timely obtain urban land information by saving the time taken in classifier debugging, sample collection, feature adjusting, etc. Previous studies indicated the effectiveness of decision rules, such as in mapping paddy rice and urban settlement [24,36]. Nonetheless, the multiple rules integration, namely MDR that orderly extracts LCLU [37,38], has not been effectively applied in GEE. This study applied and packaged MDR as an effective classifier in GEE to observe land dynamics, namely annual LCLUC [14].
MDR sequence maps LCLU through analyzing spectra, phenology, and ecology ( Figure 3). The vegetation extent was primarily extracted from the growing season FVC. Within the vegetation extent, green space and cropland could be differentiated by the NDVI of the pre-growing season, in which the growth of green space stays ahead of cropland ( Figure 4). On the other hand, band SWIR2 was applied to segment water-like pixels (including water and wetland) with high-albedo pixels (including built-up land and bare land) in non-vegetation areas [36]. Within water-like extents, we set a rule of MNDWI ≥ NDVI to map clean water space, and the other pixels accorded with characteristics of wetlands. Because building shadows have similar spectra with water bodies, we eliminated pseudo water by the cooling-effect from water space in urban. The effect indicates that urban water could lower its ambient temperature [39]. Built-up land and bare land were extracted separately by adjusting the threshold of brightness data that computed from the Tasseled-Cap Transformation. Thresholding was automatic in our MDR process. Thresholds used to be manually tested due to poor continuances of spectral information among images that differ in time and sensors [37,38]. GEE's computing stacked different images into one clean image, and thus the OTSU method can be applied for automatic thresholding [40]. The global scale study has evidenced the effectiveness of OTSU [23]. Within one decision rule, the OTSU applies to determine a value that segments two classes of pixels. Specifically, the algorithm could statistically search a number that maximizes the variances of inter-class pixels. Therefore, we could highlight the difference between the two target pixels by searching the spectral distinctness between the built-up area, bare area, water body, and vegetation [41]. OTSU thresholds could extract target land cover well (Figure 3), owing to well-analyzed urban phenology ( Figure 4) [42].  [42] for remotely sensed urban phenology.

Data and Method for Accuracy Assessment
The statistically robust method proposed by [43] was employed to assess the map accuracy. In brief, the accuracy is measured by adjusting with the areal extent of LCLU. Our previous work accumulated 311 unchanged samples that could be used to assess the accuracy [16,17]. In addition, through stratified sampling, we supplemented 677 samples by inspecting Google Earth images with very-high-resolution [24,37]. The accuracy of maps for 2004,2006,2009,2011,2012,2013,2014,2016,2017,2018, and 2019 was assessed because of an insufficient amount of very-high-resolution images. Figure A1 shows the distribution of samples. To test the availability of our methodology, we generated urban LCLU maps of nine big cities in the world for 2019. Correspondingly, we collected 5756 samples from Google Earth for accuracy assessment ( Figure A3).

Land Dynamic Analysis
The abrupt change (or turning-point) detection and change rate computation of LCLUC were conducted. The turning-point of LCLUC during the period considered was obtained by Pettitt's test, which is a nonparametric method [44]. Specifically, within a continuous observation (e.g., X 1 , X 2 , · · · , X N ), it is supposed that the r i is the rank of data. At each observing place j, the sum of ranks W j could be computed, where the j is ranged from 1 to N − 1. Hereafter, the statistic T j was introduced to infer the turning-point.
If the absolute value of the statistic T j reaches its maximum at j = J, then the turning-point occurs in time J [45]. For the turning time J, it is necessary to evident whether the estimation is significant using the sampling distribution of W j . The approximate statistic P was introduced, and the turning estimation has been performed with a significant level of 0.01.
Hereafter, the annual change rate (ACR) before and after the turning-point could be computed despite the temporal length might be different [37].
where the t 1 and t 2 denote the early year and later year, respectively. The Area t 1 and Area t 2 denote the LCLU area in the early year and later year, respectively. The ACR is in %·yr −1 .

Quantifying Cause and Effect of Land Dynamics
Time-series observations were applied to quantify the drivers and effects of land dynamic by computing correlations between LCLUC, socio-economic statistics, and climate records. The annual socio-economic statistics were collected from the statistic yearbook of Changchun, including urban population (UP), gross domestic product (GDP), GDP of agriculture sector (AS), GDP of industry sector (IS), and GDP of service sector (SS). The climate records of average temperature (ATEM) and average wind velocity (AWV) were collected from the urban station No.54161 of national meteorology (125.2 E, 43.9 N). Further, we quantified the relationship between land dynamic, finer particles (PM 2.5 ), and net primary production (NPP) to learn the effects of LCLUC. Because of lacking continuous PM 2.5 data for 2000-2019, we merged the raster dataset from 2000 to 2013 published by [46] and the official annual number from 2014 to 2019. The PM 2.5 raster dataset (0.01 • resolution) was used to compute the annual mean value within the study area. The product MOD13A3HGF (500-m resolution) was used to compute the annual mean value of NPP within the study area. Thus, the annual change of PM 2.5 in 2000-2013 and NPP in 2000-2019 was the mean value of Changchun that was computed based on published products.
We computed the built-up ratio (BR), cropland ratio (CR), green ratio (GR), and green space/built-up ratio (GBR) to echo the previous findings that suggested that green space, built-up area, and cropland changed dominantly in Changchun [8,30,31,35], GR = pixel green space /pixel total (10) where pixel built−up , pixel cropland , pixel green space , and pixel total denote the pixel area of built-up land, cropland, green space, and Changchun, respectively. Pearson analysis was applied to quantify the correlation coefficients between the observing variables mentioned above for analyzing the driver and effects of LCLUC [9], where r is the Pearson correlation coefficient. A higher absolute value of r could infer a closer relationship between n pairs of variables x i and y i .

Accuray Report of Annual Maps
The annual LCLU maps of Changchun from 2000 to 2019 were generated with a 15-m resolution. The accuracy assessment reported high overall accuracy (from 0.92 ± 0.03 to 0.98 ± 0.02) of maps (Table 3). Fewer omission errors were reflected by the producer's accuracy, indicating that our maps agreed with the ground truth. Especially, the user's accuracy of green space (≥0.84 ± 0.06) and built-up land (≥0.96 ± 0.03) had a slight commission error, laying a good foundation for observing green space and gray space. The bare land showed low accuracy, owing to the spectral similarity of bare ground, built-up area, and sparsely vegetated land. In addition, due to differences in cloud coverage, the yearly brightness image was constructed based on raw images of different numbers. Thus, the default threshold for extracting bare land might change, affecting mapping results ( Figure 3). Nonetheless, we temporally neglected that issue to ensure the entire automatic process because the bare land area was few [31,35]. Moreover, good overall accuracy was encouraging the further application of these maps.

Land Dynamics in Changchun during 2000-2019
The annual LCLUC of Changchun during 2000-2019 was recorded ( Figure 5). The built-up land expansion with the disappearance of cropland was obvious during 2000-2010. The green space exhibited a trend of descent-after-increase, turning at roughly 2011-2016. By contrast, the changes of bare land, water, and wetland were slight, altogether occupying a slight proportion.  Figure A2) [35]. The turning of cropland loss occurred in 2009 as well, and the ACR halved from −9.19%·yr −1 to −5.29%·yr −1 . Such an observation is in line with the previous study, which suggested the trend of urbanization encroaching on agriculture continually slowed since 2009 [19]. The growth of green space continued from 2000 to 2019 and shifted from rapid (2.66%·yr −1 ) to stable (0.43%·yr −1 ) in 2009. Though both the built-up land and green space were expanding, their obvious difference in ACR changes deserves in-depth study ( Figure A2). Other LCLU exhibited a slight proportion and areal changes. The bare land transitioned from gain to vanishing concurrently with the shifting of rapid-to-slow in built-up expansion. In addition, the water area kept expanding while the wetland kept losing.

Quantified Correlation of Land Dynamics, Socio-Economic Development, and Environment Change
Correlation computation indicated that though the built-up land was filling the urban space; the green space/built-up ratio stayed at~0.4 (Figure 6a), which related to the areal growth of green space (0.67) (Figure 6f). Regarding socio-economic development, Changchun's population generally kept increasing (especially after 2014), and GDP grew linearly (Figure 6b). Correspondingly, the land dynamic was significantly driven by immigration and economic process (Figure 6f). The increments of GDP (0.96) primarily encouraged built-up growth, which affirms previous knowledge [9,30], followed by migration (0.89) and industrial progress (0.8). On the other hand, the disappearance of cropland was negatively correlated with the economic and demographic changes, leading to the diving of agricultural production (0.97). The GBR was significantly related to IS (0.73), AS (−0.56), GDP (0.49), and UP (0.48), suggesting GBR was stabilized by industrial development, agricultural contraction, economic progress, and habitat growth. The urban environment changed over time. The NPP and ATEM changed complicatedly with general trends parallel mutually (Figure 6d). By contrast, the declining trend and increasing trend appeared in the AWV and PM 2.5 , respectively (Figure 6e). Pearson correlation indicated a close relationship between the AWV, PM 2.5 , and land dynamics (Figure 6f). Specifically, the AWV negatively related to BR (−0.92) and GR (−0.88), indicating built-up expansion and green space increment decelerated the urban wind. In addition, the PM 2.5 was positively related to GR, BR, and GBR, indicating that urban greening and constructing built-up land aggravated air pollution. Such a result deserves in-depth discussion because urban green space could offset emissions and improve air quality, which has been scientifically navigating the urban planning and greening in China during decades [8,17].

Availability of Recording Urban Land Dynamics by Integrating Landsat, MDR, and GEE
This study, to our knowledge, has accomplished the first record of urban land dynamics at 15-m resolution, attributing to GEE, Landsat, and MDR algorithm. Currently, GEE is the prevalent tool to observe LCLUC for global to local scale studies [6,11,23]. Large-scale LCLU data, as an important input for modeling and simulating, are furthering the understanding of global change and urban ecology. On the other hand, continuous LCLU data are helping to discover the real process of urban LCLU to correct past concepts derived from time-slice observations (e.g., impervious would not vanish). Especially, high spatiotemporal LCLU data could find the key causes of urban issues that were previously hidden by time-slice mapping [15].
Landsat-based urban mapping was promoted by GEE and MDR. First, GEE's supercomputing saved time taken in downloading, preprocessing, mosaic, etc. [17]. Second, multispectral bands were fully applied to construct spectral indices or rule nodes that could extract urban LCLU [13,24,36]. This study exactly utilized the urban phenology and ecology to generate annual maps with good accuracy (Table 2), indicating MDR is an effective classifier for mapping urban LCLU. What is more, compared with the built-in algorithm of GEE (i.e., HSV method), our pan-sharpening method saved the original multispectral information of Landsat while also enhancing the image resolution ( Figure 7). Although HSV could enhance pixel heterogeneity, it changed the original spectral information. Therefore, our method could be an alternative to HSV for the dual goals of saving spectral information and enhancing image resolution. Our pan-sharpening method implies a question: To what extent can 15-m LCLU maps enhance the overall accuracy from 30-m maps? Given that 30-m pixel concealed real land cover, it is expected that 15-m maps could promote mapping accuracy [21,27]. Moreover, we applied the 5 × 5 pixel-window to avoid the loss of spectral information (Figure 2), whether other window-size could help pixel unmixing is an interesting question.  [29]. The Landsat ID indicates that the compared image is a TOA image of the OLI sensor derived from 25 September 2017, and its path/row is 118/29. It is available and replicable to map LCLU of other cities by our methodology. For testing the availability of our methodology, we mapped the urban LCLU for nine big cities that share the same biomes with Changchun and have a population of more than five million (Figure 8). Within a 15 km radius, the urban LCLU can be mapped with high accuracy (0.85 ± 0.01 to 0.93 ± 0.01) (Figure 8). Thus, under the condition of temperate biomes, our method could be applied to map urban LCLU at 15-m resolution. Differentiating vegetation could be a difficult step because forest, grass, and crops exhibit similar spectral information during the growing season. However, non-vegetation land could be directly extracted by adjusting shortwave infrared bands or other spectral indices [36][37][38]. In this study, the non-vegetation mapping part that considered spectral variation and cooling-island could be directly used to map similar LCLU categories (Figure 3). Previous studies employed machine learning classifiers to map different vegetation through "black box" that learn knowledge from spectral indices [14,47]. The difference of this study is that the phenology of green space and cropland was well analyzed (Figure 4), implying the necessity of our method is the difference of growth time between different vegetation land covers. However, by our test, our method is currently not suitable for southern cities because of a lack of clean images and fewer phenology differences between the green space and crops. For example, in the past two decades, the total number of ETM+ and OLI images in Bangkok is 959, while the number of Beijing is 1593 (Figure 9a). For Beijing, the cloud cover in 60% of the images is less than or equal to 20%, and the phenology of green space and cropland in Beijing is distinct, laying a good foundation for mapping by the MDR method (Figure 9a). By contrast, for Bangkok, the ratio of low cloud cover images is smaller (Figure 9b). Especially, warmer climates make double-cropping possible in southern cities (Figure 9c), which reduced the chance of mapping forest and cropland directly based on the NDVI threshold (Figure 3). In addition, the similar structure of tall crops (e.g., oil palm) and forest in Southern countries makes their differentiation difficult [48]. Future work can start with building composite indices, thresholding by phenology fluctuation, adding texture information, or fitting image gaps to improve Southern urban LCLU mapping [6,41].

Time-Series Land Dynamic for Understanding Urbanization in Changchun
Changchun's 21st century exhibited characteristics of industrialization and de-agrarianization with the expansion of built-up land and the disappearance of cropland ( Figure 5). Previous studies reported a similar urbanization trend of Changchun by time-slice observation of the 2000s [30,31,35], but this study qualified the turning point of LCLU through continuous data of the entire 2000s-2010s. The year 2009 was the transition of Changchun as the ACR of major LCLU halved (Table 3), which could be the epitome of China's land use transformation from agriculture to city [19]. The temporal cost of China's urbanization reduced in the 2000s [49]. Therefore, it is necessary to ascertain what land transformed into new urban area and its effects after 2009 when the "red line" national policy proposed to maintain a basic area of cropland.
The overall economic progress promoted the growth of the BR and GR and the loss of CR, and the urban population growth and industrial production played a major role ( Figure 6). However, it is doubtful that the sustainable development in Changchun could achieve under those relationships. First, the industry-driven economy has faded in Northeast China and Jilin province in recent years [34], slowing down the expansion of built-up areas. This is a reasonable observation because the economic depression discourages investment, the emergence of high-tech industry, and the attraction of talent. The housing bubble in the early years has gradually disappeared, which could further deteriorate the economy. In addition, the administrative policy has become the major driver of urban development as the gradually decoupling between population, economy, and impervious in China [9]. In 2008, the National Eco-functional Zoning announced that the primary task of Northeast China was grain production and environment protection, which could further limit the industrial and urban development in Changchun. Though we neglected the urban volume increment as Landsat observations only reflect two-dimensional land changes, the second-rate positioning has destined the depressing economy and misty urban future in the short-term.

Implication from Building the City of Garden in Changchun
As the only northern one of four Garden Cities, Changchun has been working hard to make the urban environment more aesthetic, habitable, and famous. In 20 years, the GBR of~0.4 indicated that Changchun did well in greening under the urbanization (Figure 6a). This was partly because industrial production caused air pollution and urged more green space to promote air quality, which was inferred by a significant relationship between the IS, PM 2.5 , and GR ( Figure 6f). Nonetheless, though green space could mitigate emissions [17], the increasing GR significantly related to anabatic PM 2.5 ( Figure 6). This relationship deserves serious consideration because the national policy New-type Urbanization Plan published in 2014 has been asking cities to construct more green space to enhance the urban environment.
The increase in PM 2.5 was significantly related to the AWV descent (Figure 6f). Though our study showed that impervious expansion was slowing down in recent years ( Figure 5), the BR was significantly related to AWV descent (Figure 6f). The situation might be because of building density increasing by transforming inefficient land (e.g., urban villages) to economically efficient land (e.g., office building) under the old city renovation project in China. However, like what was observed for BR, the increasing GR significantly declined wind speed ( Figure 6). In other words, filling green space could likely hinder the airflow and thus centralize finer particles. This finding may be limited by the resolution of the PM 2.5 raster and the effectiveness of the official public data, which needs verification by integrating LCLU maps and finer PM 2.5 raster. By this study, the purposive thinning of forest density might be an effective way to strengthen the air exchange and air quality suggested by the concurrent descent of PM 2.5 and GBR after 2015 (Figure 6a,e).

Conclusions
We presented an approach to map annual urban LCLUC at high resolution (15-m) with Landsat images using GEE. With this method, we could capture the LCLUC in Changchun, and analyze the social-economic drivers and environmental effects of LCLUC. The continued urban expansion, green space growth, and cropland loss characterized LCLUC in Changchun. The LCLUC was primarily driven by economic development, immigration, and industrialization. However, inevitable urbanization under the environment of national economic development and ecological protection would weaken the sustainability of the urban development. With China's emphasis on economic development and ecosystem conservation such as the Belt and Road initiative, urban planners should explore strategies for urban sustainable development.
New knowledge about urban LCLUC in Changchun was revealed through our study. The year 2009 was the turning point of land dynamics from high-frequency to low-frequency in Changchun. The growth of built-up land and green space hindered airflow and thus increased PM 2.5 . Our research enhances the understanding of the impacts of urbanization on the environment. As air pollution has been one of the main concerns of the Chinese government, controlling air pollution and urban sustainable management requires research across different scales from local to national.  For demonstration purposes, the spatial statistical process was conducted in a pixel of 1.5 km × 1.5 km because the original pixel length is 15-m. ACR: annual change rate. Figure A3. Distribution of validation samples used for assessing the accuracy of big city maps. These samples were collected by stratified sampling from Google Earth. See Figure 8 for the maps.