Extraction of Old Towns in Hangzhou (2000–2018) from Landsat Time Series Image Stacks

: With rapid urbanization in recent decades, more and more urban renewal has taken place in China. Meanwhile, the early developed areas without change have become old towns, which need special attention in future city planning. However, other than ﬁeld surveys, there is no speciﬁc method to identify old towns. To ﬁll this gap, we used time-series image stacks established from Landsat Surface Reﬂectance Tier 1 data on the Google Earth Engine (GEE) platform, facilitated by Global Urban Boundary (GUB), Essential Urban Land Use Categories (EULUC) and Global Artiﬁcial Impervious Area (GAIA) data. The LandTrendr change detection algorithm was applied to extract detailed information from 14 band/index trajectories. These features were then used as inputs to two methods of old town identiﬁcation: statistical thresholding and random forest classiﬁcation. We assessed these two methods in a rapidly developing large city, Hangzhou, and subsequently obtained overall accuracies of 81.33% and 90.67%, respectively. Red band, NIR band and related indices show higher importance in random forest classiﬁcation, and the magnitude feature plays an outstanding role. The ﬁnal map of Hangzhou during the 2000–2018 period shows that the old towns were concentrated in the downtown region near West Lake within the urban boundaries in 2000, and far fewer than the renewed areas. The results could serve as references in the provincial and national planning of future urban developments.


Introduction
Old towns (OT), namely old construction areas, refer to constant urban impervious areas where structural features have not changed within a certain time frame. In contrast, renewed areas (RA) refer to urban surface structures that have been through subsequent changes at least once. The conversion between impervious and non-impervious surfaces, or the demolition and reconstruction of old buildings, both belong to urban renewal.
During the past 40 years, China's urban land has expanded more than 12 times [1,2]. In the meantime, some built-up areas older than 35 years have also undergone re-construction. The average building age in China is only 23 years [3]. Old towns as a whole or in part may have preservation values for cultural and historical reasons. However, in rapidly developing cities, old built-up areas are disappearing quickly due to urban-rebuild projects. Excessively fragmented land blocks and landscapes resulting from disordered construction seriously affect the appearance of a city, causing a loss of the historical style in cities. What is even worse, disordered re-constructions over old town areas are particularly disturbing socially and culturally, destroying the original characteristics and unique glamour of cities [4] in people's precious memories of their hometowns. The problem has become so serious that the Chinese government has decided to monitor the level of destruction of old towns.
Since 2012, Ecological Civilization has been officially adopted as a means to achieve sustainable development and build Beautiful China. It has now been widely recognized that many old but historic parts of cities deserve preservation in the trend of urban renewal. The Chinese government has formulated the Arable Land Minimum, Ecological Redline Policy and other urban development boundaries [5] to protect basic farmland and the ecological environment. In other words, local governments need to activate the stock and strictly control the increased use of new land resources that are currently non-urban. Hence there is a requirement to go for re-construction in built-up areas to improve the quality of urban development. Thus, urban planners and developers face the issue of which part of an old town can be replaced and which part should be preserved. However, many cities do not have maps of old construction areas or the age of built-up areas, nor do they have adequate methods for the mapping of old towns.
Although there are existing annual maps of urban boundaries [2], they do not provide information about urban renewal. The purpose of this study is to come up with a suitable method to distinguish re-built areas in original built-up areas of a certain age. In China, constructions older than 20 years, are considered old towns. To discover whether an area has experienced renewal in old built-up areas, change detection is indispensable. Relevant algorithms include VCT [6], LandTrendr [7,8], BFAST [9], CCDC [10], COLD [11], TSCCD [12], etc., mostly for forests [13][14][15], grasslands [16], croplands [17] and impervious surfaces [18]. When fitting time-series spectral trajectories and depicting interannual dynamics, LandTrendr can capture gradual long-term subtle changes and short-term mutations to satisfy various needs and purposes, thus achieving effective segmentation [7,8]. This method allows arbitrary cutting of band trajectories to obtain intuitive line segments and vertices with practical meaning. The line segments represent different land surface stages and their turning vertices represent actual ground change events. The segments and vertices can provide rich information about the events, such as start and end values, start and end years, magnitude, duration, change rate and so on. In this way, we can extract abundant land surface change signals from the spectral trajectories, which are critical for subsequent identification. Moreover, this algorithm controls noise by empirical parameter setting. Therefore, we used LandTrendr for modeling and segmentation, focusing on the interannual variability of spectral features. Afterward, statistical thresholding and random forest classification methods were adopted to utilize the extracted signals and features and judge from them.
Since 2000, the rapid growth of China's economy and gradual adjustment of land policy have continued to influence the urbanscape. Taking Hangzhou as an example, we developed a mapping framework for extracting old towns and renewal areas after 2000. We hope the method developed here can be applied to other cities for old town mapping.

Study Area and Samples
Hangzhou, the capital city of Zhejiang Province (Figure 1a), has a subtropical monsoon climate with plains and hills as the main terrain characteristics. By the end of 2018, the land area of the city had reached 16,850 km 2 , with a total population of 9.806 million and a population density of 582 person/km 2 (from Hangzhou Statistical Yearbook 2019, http://tjj.hangzhou.gov.cn/art/2019/10/23/art_1229453592_3819412.html, accessed on 18 December 2020). The developed area in the urban district increased rapidly from 177.18 km 2 in 2000 to 615.22 km 2 in 2018. There are 13 jurisdictions in Hangzhou, including 10 districts, 2 counties and 1 county-level city (Figure 1b). The region of interest (ROI) in the research was determined by Global Urban Boundaries (GUB) [19] and Essential Urban Land Use Categories (EULUC) [20] datasets. GUB data provides clear urban-rural boundaries, enabling us to focus on human activities within the urban built-up outlines. The overlaying of urban boundaries in 2000 and 2018 can highlight the difference between the extent of old built-up areas and the newly expanded regions. Rural areas included in the rapid urban expansion may constitute the phenomenon of villages-in-city [21]. On this basis, EULUC parcels remove main roads, mountains and waters, specifically orienting to human settlements.
According to multi-year high-resolution remote sensing images available in Google Earth, two categories of samples, namely, old towns and renewed areas, were collected separately through manual interpretation. Plots were first randomly selected on the images inside urban boundaries in 2018. Only those with a complete image series from 2000 to 2018 and clear land-use status were then collected as sample plots. Here, a complete image series means land surface structures and changes of plots can be checked from continuous images without lengthy time gaps. As a result, a total of 450 sample units were divided into a training set and a validation set with a 2:1 ratio ( Figure 2). The land surface features in old towns were consistent in dense images, and did not change with time, while those in renewed areas had changed at least once between 2000 and 2018.
The flow of data processing is shown in Figure 3. Google Earth Engine (GEE) was used as the platform to construct our mapping framework. It almost covered all the procedures, including Landsat data preparation, image stacks generation, LandTrendr fitting, segmentation, band trajectory inspection and random forest classification.

Data Preparation
Considering the integrity of satellite imagery and the balance of temporal frequency and spatial resolution, Landsat images are a relatively optimal source of remote sensing. For urban studies, finer resolution data is desirable, but it may be hard to include most data sources from early observations due to the later launch time of satellites, which is critical to thoroughly understand a pixel's history. The 30-m Landsat Enhanced Thematic Mapper Plus (ETM+) and Operational Land Imager (OLI) data are sufficient to cover the whole period 2000-2018. We used atmospherically corrected Surface Reflectance Tier 1 products, which include the mask of cloud, shadow, water and snow with CFMASK [22,23] and the spectral transformation from OLI to ETM+ [24]. A medoid approach [25] assigns the spectral value closest to the median among all candidate images with different timestamps to a given pixel, which can reduce data volume and minimize atmospheric impact [26]. Annual compositing in this way was applied through the green season in a certain year, to generate image collections of 14 available bands/indices in the functions library of the LandTrendr JavaScript module on GEE, including B1, B2, B3, B4, B5, B7, NDMI, NBR, NDVI, NDSI, TCB, TCG, TCW, TCA (Table 1). Therefore, for a given band, each year corresponds to one image layer. Arranging and stacking all the image layers within the study period in chronological sequence, every pixel will also correspondingly have a spectral value sequence. Afterward, its spectral value sequence will be transferred to LandTrendr for fitting and segmentation. The same stacking process was applied for all the bands, that is, each pixel can obtain 14 spectral value sequences. We included 14 bands because we expected to integrate, extract and use spectral information as much as possible to improve classification accuracy and avoid omissions.

LandTrendr
Since the spectral variation within built-up areas may not be significant due to the similarity of building materials, we expected to capture disturbance to the maximum extent [7] and improve the sensitivity of emergency response, rather than smoothing the overall trend. After sufficient testing and adjustment, the main parameters in the function were set as shown in Table 2 based on foregoing considerations. It is hard to define the change direction explicitly due to complicated human intervention, so both gain and loss were taken into account. Our target was the greatest one with the largest change magnitude (mag, multiplied by 1000) among all segments of a trajectory (just like successive parts of red broken lines in Figure 4). Then the change rate of this segment can be calculated from its duration (dur): rate = mag dur (1)  Figure 4a shows the time series trajectories of a typical renewed area. The bold crimson lines in the curves refer to the greatest segments, the slope of which is equivalent to the rate mentioned before. In this region, the original residence was demolished in 2007 and rebuilt afterward, so the indices fluctuated greatly after 2007, with high mag. Mag of NDMI was usually less than that of NBR and NDVI, and there may be time differences between indices. On the contrary, we expected the ideal trajectory of an old town to converge to a horizontal straight line. Therefore, in Figure 4b, the index trajectories of a typical old town were much smoother with low mag. It is worth noting that if the curve is stable enough, the greatest segment will cover the entire period, like NDVI in Figure 4b.

Statistical Thresholding (ST)
Normalized Difference Built-up Index (NDBI) [37] is an excellent indicator for built-up areas. From Equation (2) and Table 1, we can see that NDMI is equivalent to reverse NDBI. We do not mind the change direction of segments in the research, so NDMI can replace NDBI to obtain relevant information about built-up areas. NDVI contains much information about vegetation. NBR shows a high Signal to Noise Ratio (SNR) and prediction ability when detecting disturbance [14]. Therefore, we took statistical analysis for the three representative indices to determine the thresholds between the old and the renewed.
Each index has mag in two directions, gain and loss, and we take the higher value as ∆Index in Equation (3), where Index refers to NDMI, NBR or NDVI. Figure 5a shows the density curves of ∆Index for old town and renewed area sets separately. It is obvious that notable differences exist between the two sample categories. The ∆Index values of old towns are quite concentrated and generally low, while those of renewed areas are relatively more scattered and higher. The greatest absolute difference between the two categories is in NDVI, and the least is in NDMI, which can be seen from their distinct averages. However, there are also partial overlaps between them leading to incomplete separation. The box plots in Figure 5b show similar differences. Additionally, by Kolmogorov-Smirnov test based on cumulative distribution function, there are significant differences between ∆Index of old and renewed samples, with no significant differences between the training and validation set in each category. Therefore, it is feasible to determine general thresholds for classification according to the differences between old and renewed training samples. For NDMI or NBR, the threshold is determined by the maximum in the upper quartile (Q3) of training samples of old towns and the lower quartile (Q1) of training samples of renewed areas (Equation (4)). Because NDVI is affected by some external factors (such as vegetation growth and landscape greening), the limit for it should be relaxed. Therefore, its threshold is determined by the upper edge (Q3 + 1.5IQR) of training samples of old towns (Equation (4)). The final thresholds are 142 for NDMI, 260 for NBR and 275 for NDVI. Only when the mags of three indices are all no more than the corresponding thresholds, can a pixel be judged as an old town. Index = MAX gain mag Index , loss mag Index Threshold Index = MAX{ Index ot : Q3, Index rt : Q1}, Index = NDMI or NBR Index ot : (Q3 + 1.5IQR) , Index = NDVI where, ot: training samples of old towns rt: training samples of renewed areas.

Random Forest (RF)
Random forest was first introduced by Leo Breiman and Adele Cutler [38,39], and it belongs to the range of supervised machine learning. According to the idea of ensemble learning, the random forest classifier contains numerous decision trees, and meanwhile trains the samples and makes predictions based on various features. It shows stable, robust and excellent classification performance [40] and is able to deal with huge and high-dimensional datasets, equipped with rapid learning ability and high flexibility.
We expected to make full use of multi-dimensional spectral space in mapping. Therefore, a total of 84 features composed of the mag/dur/rate in gain/loss of the 14 bands/indices were used as inputs to the random forest classifier on GEE. The same training set was used to train the model, then generating the classification result for the whole ROI.

Post-Processing
Since some non-built-up areas that should not be included were still mixed in the results, which would affect the statistics of old towns and renewed areas, we used Global Artificial Impervious Area (GAIA) [1,2] data to mask out those areas. In GAIA, local maps of China cover 1978-2017 with more precise details, while global maps cover 1985-2018 with more macro patterns. To satisfy the study period 2000-2018 and meanwhile improve the precision, a non-impervious (non-urban) pixel can be confirmed only when it exists in both maps. The obtained non-impervious areas were composited with classification results to remove undesired non-built-up areas.

Validation and Comparison
The randomly distributed validation set includes 50 old sample units and 100 renewed sample units. As shown in Table 3, the overall accuracy (OA) of ST is 81.33%, and producer's accuracy (PA) and user's accuracy (UA) of renewed areas are much higher than those of old towns. The overall accuracy of RF is 90.67%, which is nearly 10% higher than the former, and PA/UA of the two categories are also largely improved. Moreover, there is more balance between omission and commission errors, indicating that the RF result is more robust. In contrast, ST appears relatively plain and limited. Because the sample units unevenly fell into districts of Hangzhou, we screened out districts with more than 6 sample units to calculate local overall accuracies and to observe the spatial distribution pattern of mapping accuracies ( Figure 6). The overall accuracies in the districts with fewer validation sample units are not representative enough and biased, so they were not taken into account. Despite the general similarity, certain differences still exist in the two mapping results. On the whole, all the available accuracies are higher than 70%. However, the local accuracies of the random forest result perform better than those of the statistical thresholding result, just like the situation for the entire study area. In terms of the detailed common ground between the two results, both Yuhang and Binjiang show remarkable accuracies (all 100%), which are also the main expanded regions. Additionally, Shangcheng and Gongshu in the RF result also have similarly ideal accuracies (94.87% and 100%, respectively). However, Shangcheng and Jianggan in ST and Xiacheng in RF were mapped relatively poorly with higher errors than other districts, where the land use is more fragmented, complex and mixed.
In addition, we compared the actual classification quality of the two methods through local viewing. The selected part is close to West Lake, with dense commercial and residential land (Figure 7a). Old towns in the ST result are relatively fragmented, presenting poor identification for the renewed core commercial blocks with small spectral variation, such as the lakeside zone in the lower-left corner. This may be related to the spectral similarity of soil, cement and other building materials before, during and after reconstruction, as well as the resolution limitation and spectral mixing. In the RF result, the connectivity of homogeneous areas is higher, forming obviously concentrated blocks, which can provide more intuitive and practical information. From three description perspectives of mag/dur/rate simultaneously, the full use of change signals in 14 indices enables more accurate judgment for commercial blocks. The red/black boxes indicate the West Lake Cultural Plaza. During the 2000-2018 period, its original residential buildings were demolished and then rebuilt into commercial buildings and a square ahead. In the RF result, the whole region was effectively identified as renewed areas. However, in the ST result, in spite of successful identification for the commercial buildings, the method failed to identify the square with mainly cement surfaces and greenbelts. On the whole, the RF method is more accurate, robust and reliable. We also magnified the northwestern, northeastern, southwestern and southeastern parts of this extent and covered the 2000 orthophoto with semi-transparent ST result and RF result, to inspect the details (Figure 7b). It can be seen that in the magnified perspective, the RF result similarly presents better accuracy and connectivity in many homogeneous neighborhoods, which is generally consistent with the ground truth. On the other hand, the RF method is relatively more aggressive, and consequently, some actual old towns were misjudged as renewed areas.

Supplement of Random Forest Result
From the importance of 84 features in random forest classification (Figure 8), mag is the decisive factor in classification and contains a lot of information related to changes. Contrarily, dur appears less important. Summing the importance for each index, B3, TCG, NBR, B4 and NDVI perform better in classification, so the Red and NIR bands are more effective and powerful in identifying old towns. Figure 9 shows partial random forest results in Hangzhou. It can be intuitively felt that old towns are far fewer than renewed areas in main municipal districts, and are mostly distributed in Shangcheng, Xiacheng and Xihu Districts. These are also the earliest developed areas in Hangzhou, where old and new buildings often coexist adjacently. Old towns are concentrated in the urban center near the West Lake within the urban boundaries in 2000, and gradually extend and diffuse outwards with density decreasing. The traditional central market position of the lakeside zone can be traced back to the 1910s, which also laid the foundation for the subsequent spreading urbanization and modernization [41]. Between the urban boundaries in 2000 and 2018 are the expanded zones, most of which used to be countryside. These areas usually experienced urban renewal and became builtup areas. Many agricultural lands still remain on the east side of Xiaoshan, that is, the lower-right corner in the figure showing the main districts. Old towns, mainly in the form of cottages, are arranged systematically and distributed in the gaps among fields. In Lin'an, most old towns are also distributed inside the urban boundaries in 2000 to the west of Qingshan Lake. Old towns in Tonglu mainly gather along the Fuchun River, especially where it joins the tributary Fenshui River. In Chun'an, the residents also live near water. Old towns are compactly located on the northeastern lakeshore of Qiandao Lake, at the foot of the mountains, and the human habitation has gradually extended eastward into gentle valleys, developing livable built-up areas through artificial transformation.   Most regions in Hangzhou have experienced urban renewal (1411.77 km 2 ), and the total area of old towns is 134.24 km 2 , accounting for 8.68% of its urban area (OT + RA) ( Table 4). Within the urban boundaries in 2000, the area of old towns is 59.06 km 2 , accounting for 19.79% of the urban area, which is mainly residential. Within the expanded zones between the urban boundaries in 2000 and 2018, although the actual area of old towns (75.18 km 2 ) is more than the former, its percentage is only 6.03%. This means that much earlier urban centers are more likely to remain the same as before, and renewal occurs more frequently in the suburbs. Urban expansion involves not only planning adjustment and population migration, but also surface reconstruction. In 13 jurisdictions, Shangcheng District located in the core of Hangzhou, has the highest percentage of old towns, accounting for 39.47% (Table 4), in spite of its actual area being small. The next is the adjacent Xiacheng District, accounting for more than one-fifth (22.27%). The percentage of Xiaoshan is not high (7.80%), but it possesses the largest old town area (35.88 km 2 ). The smallest area is in Chun'an, which is also the seat of Qiandao Lake, with many surrounding hills in the territory. The lowest percentage is in Yuhang. As we can see in Figure 10, from downtown outwards, the area of old towns in districts presents a low-high-low change process in general, while the corresponding percentage presents an opposite high-low-high change tendency.

Discussion
In order to capture more significant interannual differences, we did not include the non-green season images in annual compositing, which may meanwhile omit the surface changes before June or after September in a year. Observation conditions and noise will also influence the image quality. The fitting of LandTrendr cannot guarantee that it is consistent with the actual change, such as smoothing small fluctuations when detecting a disturbance. Therefore, the setting of control parameters is particularly critical, which should be adjusted cautiously according to our specific needs.
Currently, the inputs of the random forest method are features of the greatest segments in band/index trajectories. When considering introducing more information as a classification basis, we also tried to add the fastest segments to double the feature quantity. However, the classification effect was not improved as well as we expected, and the accuracy was even lower. The reason may be that some of the added features did not contain more valuable information, and even belonged to noise for the classifier, or an excessive quantity of features led to the over-fitting problem. Therefore, we only kept the greatest segments.
The spatial scope of multi-year high-resolution images is limited in Google Earth, and is mostly concentrated in densely populated areas. High-quality satellite imagery is often missing in the urban fringe, so there is a certain spatial bias in the obtained samples. Additionally, we tried to expand the period to 1990, but the high-resolution images from 1990 to 2000 in most areas were hard to acquire. Further, 30-m resolution Landsat imagery is not fine enough to observe the internal changes inside built-up areas, leading to insufficient samples. Thus, it is impossible to conduct effective training and validation and ensure the credibility of the results. Next, we hope to obtain other useful data sources, such as official approval documents on real estate development, enabling us to judge accurately when collecting samples without fine images. Therefore, maps of old towns with different duration can be composited to increase information capacity in the data.
Uncertainty also exists on the map, especially in the city center. House repair, roof material replacement, nearby shed construction, surrounding plant growth over years and tree canopy sheltering are all possible error sources, resulting in false positives, while sometimes the spectral curves cannot accurately indicate surface changes, which will result in false negatives. The mixture inside a pixel will also affect the judgment of the classifier.
Auxiliary data may influence the results. GUB does not necessarily coincide with the actual urban boundaries. The compositing of the GAIA product may lead to a certain degree of error accumulation. Therefore, more efforts need to be made to ensure that urban boundaries are precisely mapped.

Conclusions
In the process of rapid urban expansion and renewal, some old towns have been preserved. We used the LandTrendr algorithm to extract knowledge from Landsat timeseries image stacks during the period 2000-2018, and adopted statistical thresholding and random forest classification methods to map old towns in Hangzhou respectively. Random forest is clearly superior, with better overall accuracy of 90.67%. In the classification, B3, TCG, NBR, B4 and NDVI had higher importance, indicating that the Red and NIR bands contain more significant information about built-up areas. In the map of Hangzhou, old towns are concentrated in the city center near West Lake, which are mainly residential and only account for about one-tenth of the renewed areas. The percentage of old towns in urban areas within the 2000 urban boundaries was about three times that within the expanded zones during the period 2000-2018. Shangcheng District has the highest percentage, and Xiaoshan District has the largest area, relative to the lowest percentage in Yuhang and the smallest area in Chun'an. Although this research takes Hangzhou as an example, the mapping framework based on the GEE platform can be easily migrated to other regions and different periods with targeted regional samples. Through Cloud Computing, national and even global maps could be produced at little cost.
Because the classification completely relies on spectral information, image quality and disturbance detection appear particularly critical, and to some extent, the errors can be reduced through parameter optimization. Moreover, the sample quality depends on the artificial visual interpretation of observers and affects the capacity of the classifier. In future, the rules of the statistical thresholding method will be further revised and expanded to meet higher accuracy and interpretability. We will also test alternative change detection, mapping algorithms and multi-algorithm integration to extract old towns. The distribution of old towns in this research can provide some guidance for urban planning and highlight the remaining areas lagging behind urban development.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to China's geographic data restriction policy.