Remote Sensing Derived Indices for Tracking Urban Land Surface Change in Case of Earthquake Recovery

: The study of post-disaster recovery requires an understanding of the reconstruction process and growth trend of the impacted regions. In case of earthquakes, while remote sensing has been applied for response and damage assessment, its application has not been investigated thoroughly for monitoring the recovery dynamics in spatially and temporally explicit dimensions. The need and necessity for tracking the change in the built-environment through time is essential for post-disaster recovery modeling, and remote sensing is particularly useful for obtaining this information when other sources of data are scarce or unavailable. Additionally, the longitudinal study of repeated observations over time in the built-up areas has its own complexities and limitations. Hence, a model is needed to overcome these barriers to extract the temporal variations from before to after the disaster event. In this study, a method is introduced by using three spectral indices of UI (urban index), NDVI (normalized di ﬀ erence vegetation index) and MNDWI (modiﬁed normalized di ﬀ erence water index) in a conditional algebra, to build a knowledge-based classiﬁer for extracting the urban / built-up features. This method enables more precise distinction of features based on environmental and socioeconomic variability, by providing ﬂexibility in deﬁning the indices’ thresholds with the conditional algebra statements according to local characteristics. The proposed method is applied and implemented in three earthquake cases: New Zealand in 2010, Italy in 2009, and Iran in 2003. The overall accuracies of all built-up / non-urban classiﬁcations range between 92% to 96.29%; and the Kappa values vary from 0.79 to 0.91. The annual analysis of each case, spanning from 10 years pre-event, immediate post-event, and until present time (2019), demonstrates the inter-annual change in urban / built-up land surface of the three cases. Results in this study allow a deeper understanding of how the earthquake has impacted the region and how the urban growth is altered after the disaster.


Introduction
The recovery process after a disaster has several components that are impacted by the nature of the hazardous event, the type and status of structures and infrastructures, and vulnerability, resilience, and socio-economic characteristics of the community [1][2][3][4][5][6][7]. One of the main indicators of how a community is recovering is depicted in the reconstruction and growth of city after the event. Earthquakes have unique characteristics compared with other natural hazards as the damaging impacts due to shaking depends on the building types and soil conditions in the affected region, and they have a short on-set time. The application of remote sensing for studying the impacted region is highly valuable and especially advantageous when other sources of data are unavailable or incomplete. account for the impact from barren lands [45]. Waqar et al. [46] proposed the BRBA (band ration for built-up area) and NBAI (normalized built-up area index) that were verified with a case study from Pakistan. The BRBA is the ratio of RED to SWIR, and NBAI is mathematical construct from Green, SWIR1, and SWIR2 of Landsat imagery.
Other urban indices to mention include BCI (biophysical composition index) based on tasseled cap transform; MBI (modified built-up index) using SWIR2, Red, and NIR; BAEI (built-up area extraction index); NDSV (normalized difference spectral vector) using Red, Green, and SWIR; and CBI (combinational build-up index) using the principal component analysis as well as SAVI and NDWI (normalized different wetness index). Valdiviezo et al. [39] conducted a comparative study of built-up index methods and found the highest accuracies for IBI, BAEI, and NDSV; with lowest accuracy for MBI and BRBA. Among these urban indices, one common issue is the difficulty in differentiating barelands and the built-up areas in one image because of their similar spectral responses [47]. The literature proves the complexity in barelands and built-up areas and the limitation of urban indices.
To accommodate for this complexity in urban/built-up land mapping, we propose to converge the aforementioned urban indices to distinguish between built-up areas and bare-land before and after an earthquake. A disaster-specific conditional algebra approach is developed, which has flexibility in determining thresholds due to local environmental variations. This study effectively delineates urban/built-up land surface change for the three heterogeneous case studies across the globe with dissimilar geographical and socio-economical contexts. Application of machine-learning techniques and support vector machine (SVM) classifiers in the future can further complement the results of our study and proposed methodology [48]. The extracted annual change trend in urban areas is an essential element in understanding the recovery and reconstruction process after each earthquake.

Study Areas and Datasets
Three heterogeneous study areas were selected based on three earthquakes in different geographical locations across the globe: the Christchurch earthquake 2010 in New Zealand; the L'Aquila earthquake 2009 in Italy; and the Bam earthquake 2003 in Iran. The three selected events have distinctive recovery characteristics. On September 3, 2010 a M7.0 earthquake struck approximately 50 km to the west-northwest of Christchurch, New Zealand [49]. This mainshock (Canterbury earthquake) was followed by several aftershocks, including a M6.3 shock in 2011 (Christchurch earthquake). Consequently, 87% of homes in Christchurch were damaged (30% of them had major damage), more than 180 people died based on New Zealand's police reports, and estimated damage from both shocks were about NZ$40 billion (2015$ value) according to the Reserve Bank of New Zealand [50,51]. The recovery encountered three to five year wait times for insurance claims payments, rezoning, and foundation standards [5]. On April 6, 2009, a M6.3 earthquake struck the city of L'Aquila, Italy [49]. As a result, 309 lost their lives, 1500 were injured, 3,000-10,000 buildings were damaged, and there was a total economic loss of EUR 17.4 billion (USD$ 20 billion, 2018$ value) according to Swiss Re estimates [52]. Within six months, the government built base isolated housing for 15,000 families. However, due to lack of funding, others lived in hotels and other towns for two to three years, and some relocated [5,53]. On December 26, 2003, an M6.6 earthquake struck city of Bam in southeastern Iran [49] that resulted in more than 26,000 loss of life, about 30,000 injuries, an estimated 70% complete loss of housing, and about USD$1.5 billion loss (2014$ value) from World Bank estimates [54,55].
All three earthquakes occurred in the time frame of 2000-2010 with Intensity > VIII and Magnitude > M6.0. The time frame of each case study starts from 10 years before the earthquake to present (2019). In this time frame, the pre-event and post-event imagery is collected. The 10 years of imagery before the earthquake is necessary for establishing a baseline estimate of the urban growth trend from pre-event conditions. Only urban areas are studied to monitor the urban/built-up land surface change before and after the event. The extent of area of interest (AOI) for each case is determined by Remote Sens. 2020, 12, 895 4 of 25 the boundary of intensity of six (VI) or higher from USGS ShakeMap [49] (i.e., 'slight damage' or more based on Modified Mercalli intensity scale) for the earthquake events, since it delineates the possible region of physical damage (Figure 1). Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 25 surface change before and after the event. The extent of area of interest (AOI) for each case is determined by the boundary of intensity of six (VI) or higher from USGS ShakeMap [49] (i.e., 'slight damage' or more based on Modified Mercalli intensity scale) for the earthquake events, since it delineates the possible region of physical damage (Figure 1).  The primary data in this study are Landsat imageries downloaded from USGS EarthExplorer (Landsat Collection Level 1) [56]. Table 1 below shows a summary of Landsat imagery used for  the case studies (more details in Table 2). In order to find best time of image acquisition in certain year, the atmospheric conditions and phenology were examined. The environmental consideration of phenological cycle characteristics can be either for vegetation phenology or urban-suburban phenological cycles [57].   The images for each case are collected based on temporal availability, cloud coverage, and type of sensor availability. The 2003 Bam earthquake occurred at 5:26 am local time on 26 December and the months of January and February are chosen for anniversary date. The 2009 L'Aquila earthquake occurred on 6 April and the months of June and July are more suitable for anniversary dates due to cloud and snow coverage in other months. Lastly, the 2010 New Zealand occurred at 4:35 am local time on 4 September, 2010 with a major aftershock at 12:51 pm on 22 February 2011. The months of October and November are more suitable for anniversary dates but the actual dates of gathered imagery vary between Septembers to February due to cloud coverage. The time period between January-March is also in cyclone season. All these effects are taken into consideration when selecting the images. For initial analysis and data acquisition, there are no threshold assignments for cloud-coverage and the image that had the least cloud-coverage for each year is selected.

Methods
The proposed methodology uses three indices of UI, NDVI, and MNDWI, and builds a knowledge-based classifier to extract the urban/built-up features with a conditional algebra. These three indices are chosen after testing several combinations of indices that have shown reliable results based on literature and test sample points from three different years in each case (minimum 50 points/pixels randomly checked with high-resolution Google Earth images in AOI, for the same three years described in accuracy assessment) to delineate the built-up areas and remove the vegetation, water, and barelands. The spectral profiles for the sample points from the three cases is shown in Figure 2 that depicts the local variations and application of the three chosen indices as an average of spectral responses for each of the four classes. For example, the UI is preferred over other urban indices such as NDBI or EBBI since the distinction between urban/built-up areas versus barelands is better represented from the difference between SWIR2 and NIR, than SWIR1 and NIR (Figure 2b). Some complex urban indices, such as IBI, combine several bands but do not provide flexibility to include the local variations observed in our study cases and therefore, are not further investigated in this study. Hence, the UI is used to extract built-up surfaces, while NDVI is applied for removing vegetation and MNDWI is employed to eliminate both water-bodies and bareland. The formulas for the three indices are   For each image, the three indices are calculated and composited into a three-layer RGB image: UI as Red, NDVI as Green, and MNDWI as Blue. This is the image used to classify urban/built-up and non-urban in this study. To extract built-up areas, a conditional algebraic test is applied to extract higher UI values and remove high NDVI and MNDWI values (Appendix A, Table A1). The Gaussian-based (Bayesian) automatic selection of thresholds for indices [58] was considered but due to the complexities of the urban landscape and earthquake cases, the pixel values do not follow a normal distribution. Also, as demonstrated with 100 randomly selected sample pixels in each case (Figure 3), urban index of built-up pixels reveals an extremely long range, reflecting heterogeneous urban structures. Brighter built-up pixels (higher UI values) are mixed with barelands, while the darker built-up pixels (lower UI values) could be mixed with some vegetation and water pixels. In the L'Aquila case (Figure 3.b), UI values of built-up pixels are more densely clustered with barelands, even at higher UI values. In the Bam case ( Figure 3.c), UI values can better delineate built-up from barelands, but they show higher mixture with water. Figure 3 reveals that the UI itself could not optimally extract built-up pixels from barelands, vegetation, and water. On the other hand, vegetation is more distinguishable from NDVI, and water could be better identified from MNDWI. Therefore, a multi-index approach is established in this study. Figure 3 also shows the complexity of urban surface in different earthquake cases. Conditional algebra provided the means to assign optimal thresholds based on each case characteristics and variability of indices' values that was especially effective in distinguishing barelands from built-up features. In the conditional statement, the threshold for high NDVI values removes vegetation and MNDWI is used to both remove both water bodies and the barelands, which is the reason that the conditional statement is applied to both higher values and lower values for this band. Since the histogram distributions are variable temporally and spatially, the thresholds start with a baseline defined based on literature [59][60][61], and the mean values for each band in each image to accommodate for the variability ( Figure  3, Table 3). Finally, the classifier assigns each pixel from the three-layer image to a binary value of built-up/non-urban for each year. The general conditional statement for extracting built-up surfaces, prior to any modifications, is as below, (where 1 is indicating built-up and 0 is indicating non-urban): The thresholds are different for each of the case studies due to the environmental, seasonal, and building variations. The variable phenological thresholds have been applied in past studies in analysis of NDVI time series [62]. In this study, Christchurch, New Zealand is in a flat area east of the Canterbury Plains and bounded by Pacific Ocean coast, volcanic slopes in the South and Waimakariri River in the North (Figure 3.a); while the case of L'Aquila, Italy is mountainous ( Figure  3.b), and Bam is in a desert environment surrounded by mountains (Figure 3.c). The values for a few cases are modified based on the histogram variation for the case and time of the year. It provides For each image, the three indices are calculated and composited into a three-layer RGB image: UI as Red, NDVI as Green, and MNDWI as Blue. This is the image used to classify urban/built-up and non-urban in this study. To extract built-up areas, a conditional algebraic test is applied to extract higher UI values and remove high NDVI and MNDWI values (Appendix A, Table A1). The Gaussian-based (Bayesian) automatic selection of thresholds for indices [58] was considered but due to the complexities of the urban landscape and earthquake cases, the pixel values do not follow a normal distribution. Also, as demonstrated with 100 randomly selected sample pixels in each case (Figure 3), urban index of built-up pixels reveals an extremely long range, reflecting heterogeneous urban structures. Brighter built-up pixels (higher UI values) are mixed with barelands, while the darker built-up pixels (lower UI values) could be mixed with some vegetation and water pixels. In the L'Aquila case (Figure 3b), UI values of built-up pixels are more densely clustered with barelands, even at higher UI values. In the Bam case (Figure 3c), UI values can better delineate built-up from barelands, but they show higher mixture with water. Figure 3 reveals that the UI itself could not optimally extract built-up pixels from barelands, vegetation, and water. On the other hand, vegetation is more distinguishable from NDVI, and water could be better identified from MNDWI. Therefore, a multi-index approach is established in this study. Figure 3 also shows the complexity of urban surface in different earthquake cases. Conditional algebra provided the means to assign optimal thresholds based on each case characteristics and variability of indices' values that was especially effective in distinguishing barelands from built-up features. In the conditional statement, the threshold for high NDVI values removes vegetation and MNDWI is used to both remove both water bodies and the barelands, which is the reason that the conditional statement is applied to both higher values and lower values for this band. Since the histogram distributions are variable temporally and spatially, the thresholds start with a baseline defined based on literature [59][60][61], and the mean values for each band in each image to accommodate for the variability ( Figure 3, Table 3). Finally, the classifier assigns each pixel from the three-layer image to a binary value of built-up/non-urban for each year. The general conditional statement for extracting built-up surfaces, prior to any modifications, is as below, (where 1 is indicating built-up and 0 is indicating non-urban): Remote Sens. 2020, 12, 895 8 of 25 3.b), and Bam is in a desert environment surrounded by mountains (Figure 3.c). The values for a few cases are modified based on the histogram variation for the case and time of the year. It provides flexibility in defining the indices' thresholds according to local characteristics, which enables more precise distinction of features based on environmental variability. Once classification is completed, the annual urban/built-up land surface change is calculated by subtracting the resulted binary pixel values of each year from its previous one. A histogram match process is applied prior to index calculations. The imagery for a reference year is picked for histogram match based on having the least cloud coverage. For each case, all bands for all years are histogram matched by building a model in ERDAS Imagine. Hence, the inter-year variation from non-earthquake impacts such as atmospheric noises and weather are removed before extracting thresholds. The radiometric calibration and atmospheric correction is the first step in image processing, which is done by using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) model in ENVI 5.4 software. Atmospheric correction is needed in this change detection method, since it is based on linear transformations of the data (e.g., NDVI) [54]. The FLAASH model incorporates the MODTRAN (MODerate spectral resolution atmospheric TRANsmittance) radiation transfer code to correct wavelengths in the visible through near-infrared and shortwave infrared regions [63]. If cloud coverage over the study area (AOI) was more that 10% of area for a year, those years were eliminated. The thresholds are different for each of the case studies due to the environmental, seasonal, and building variations. The variable phenological thresholds have been applied in past studies in analysis of NDVI time series [62]. In this study, Christchurch, New Zealand is in a flat area east of the Canterbury Plains and bounded by Pacific Ocean coast, volcanic slopes in the South and Waimakariri River in the North ( Figure 3a); while the case of L'Aquila, Italy is mountainous (Figure 3b), and Bam is in a desert environment surrounded by mountains ( Figure 3c). The values for a few cases are modified based on the histogram variation for the case and time of the year. It provides flexibility in defining the indices' thresholds according to local characteristics, which enables more precise distinction of features based on environmental variability. Once classification is completed, the annual urban/built-up land surface change is calculated by subtracting the resulted binary pixel values of each year from its previous one.
A histogram match process is applied prior to index calculations. The imagery for a reference year is picked for histogram match based on having the least cloud coverage. For each case, all bands for all years are histogram matched by building a model in ERDAS Imagine. Hence, the inter-year variation from non-earthquake impacts such as atmospheric noises and weather are removed before extracting thresholds. The radiometric calibration and atmospheric correction is the first step in image processing, which is done by using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) model in ENVI 5.4 software. Atmospheric correction is needed in this change detection method, since it is based on linear transformations of the data (e.g., NDVI) [54]. The FLAASH model incorporates the MODTRAN (MODerate spectral resolution atmospheric TRANsmittance) radiation transfer code to correct wavelengths in the visible through near-infrared and shortwave infrared regions [63]. If cloud coverage over the study area (AOI) was more that 10% of area for a year, those years were eliminated.
The accuracy assessment is the other part of methodology to estimate the rigor of classification results. The accuracy of classification is often reported as the percentage of correct identification, which is presented as errors of commission (presented as user's accuracy referring to the percentage that are actually correct) and errors of omission (presented as the producer's accuracy, is the percentage of the total pixels that belong to that category in the image). The confusion matrix is used with validation samples. The overall accuracy is calculated from the ratio of the sum of samples along the diagonal to the number of validation samples. In order to get representative random points for accuracy assessment, the number of points is selected by using a validation sample size for each case, and then they are distributed in two steps: 1) half of points are selected by equalized stratified random distribution that collects an equal number of random coordinates from each class; 2) the other half is a complete random distribution without replacement or any predefinitions. When accounting for the total number of pixels, with a 95% confidence level and confidence interval of 5, the minimum validation sample size is calculated for each case, which is rounded to 350 points for each case [56]. The random points are distributed in the pre-processed Landsat image and then compared with Google Earth satellite imagery and our classification.

Conditional statement thresholds for L'Aquila, Italy
Year

Conditional Algebra and Classification
The conditional algebra accommodates for local characteristics of each case with flexibility in defining the indices' thresholds based on environmental and seasonal variability. Several thresholds were tested for each case and based on the statistics and standard deviation from the mean for each index values (Appendix A, Table A1), histogram distributions and sample testing points from both built-up and non-urban features (Figure 3), the best fitted thresholds were chosen that had higher accuracy percentages. The fitted thresholds are extracted from samples like Figure 3 for three years in each case and matched with the mean and standard deviation values of each index for the study area (appendix A), so that the conditional statement is only extracting urban features, and the thresholds are adjusted based on the histogram for that year. All thresholds are summarized in Table 3.
The example on (Figure 3a) indicates that lower MNDWI values show barelands and can provide a threshold to distinguish them from built-up lands. The threshold between these two at MNDWI in the sample points, is seen at −0.48, and after matching with MNDWI values for the whole study area with the mean of −0.456 and standard deviation of 0.135, the value of mean is a chosen threshold. This is tested for sample points from three other years for Christchurch, New Zealand and then applied for all years. The same process is applied for other indices as well (comparison of final results is presented in Figure 4). For each year, in case of New Zealand the thresholds were set based on  For each year, in case of Italy the thresholds are assigned as The thresholds are set according to samples like Figure 3.b, as it can be seen for example the lowest NDVI value seen in the sample is 0.48 for vegetation, and matching with general distribution of NDVI value for the study region in 2017, the mean is 0.588, and standard deviation of 0.133 (Appendix A, Table A1), thus threshold of one standard deviation below the mean matches a threshold for NDVI at value of 0.455 for 2017. The threshold is tested with similar samples from three years, and same process is applied for other indices. A comparison of final results is presented in Figure 5. The Bam, Iran area, is surrounded by both desert (barelands) and mountains; therefore, the thresholds are slightly different from other cases to accommodate for the geographical features (Figure 3.c) and only extract built-up areas based on the histograms. The lower threshold for NDVI in this case is helpful for eliminating barelands of surrounding desert (comparison of final results is presented in Figure 6). For each year, in the Iranian case the thresholds are assigned as For each year, in case of Italy the thresholds are assigned as The thresholds are set according to samples like Figure 3b, as it can be seen for example the lowest NDVI value seen in the sample is 0.48 for vegetation, and matching with general distribution of NDVI value for the study region in 2017, the mean is 0.588, and standard deviation of 0.133 (Appendix A, Table A1), thus threshold of one standard deviation below the mean matches a threshold for NDVI at value of 0.455 for 2017. The threshold is tested with similar samples from three years, and same process is applied for other indices. A comparison of final results is presented in Figure 5.  For each year, in case of Italy the thresholds are assigned as The thresholds are set according to samples like Figure 3.b, as it can be seen for example the lowest NDVI value seen in the sample is 0.48 for vegetation, and matching with general distribution of NDVI value for the study region in 2017, the mean is 0.588, and standard deviation of 0.133 (Appendix A, Table A1), thus threshold of one standard deviation below the mean matches a threshold for NDVI at value of 0.455 for 2017. The threshold is tested with similar samples from three years, and same process is applied for other indices. A comparison of final results is presented in Figure 5. The Bam, Iran area, is surrounded by both desert (barelands) and mountains; therefore, the thresholds are slightly different from other cases to accommodate for the geographical features (Figure 3.c) and only extract built-up areas based on the histograms. The lower threshold for NDVI in this case is helpful for eliminating barelands of surrounding desert (comparison of final results is presented in Figure 6). For each year, in the Iranian case the thresholds are assigned as The Bam, Iran area, is surrounded by both desert (barelands) and mountains; therefore, the thresholds are slightly different from other cases to accommodate for the geographical features ( Figure 3c) and only extract built-up areas based on the histograms. The lower threshold for NDVI in this case is helpful for eliminating barelands of surrounding desert (comparison of final results is presented in Figure 6). For each year, in the Iranian case the thresholds are assigned as

Accuracy Assessment
The Kappa coefficient of agreement is computed from the confusion matrix and reported for each case, in selected sample years. Due to large number of imagery and consistent methodology for each case, one representative image from every 10 years is chosen to simplify the assessment process. Hence, three sample years from the beginning, middle, and end of the time frame for each case study is chosen for accuracy assessment. The two classes that are extracted and presented in binary images are evaluated: built-up and non-urban.
According to the results (Table 4), the overall accuracy of urban/built-up ranges between 92% to 96.29%; and the Kappa value varies from 0.79 to 0.91. The lower user accuracy percentages for New Zealand 2017 and Italy 2017 are due to errors introduced from cloud presence, (even though the clouds are not over the AOI). Differences between sensors (TM versus OLI) do not seem to change the accuracy, as is seen in case of Iran. Thus, for longitudinal studies of large geographic areas, the proposed method provides the flexibility to delineate built-up land.

Accuracy Assessment
The Kappa coefficient of agreement is computed from the confusion matrix and reported for each case, in selected sample years. Due to large number of imagery and consistent methodology for each case, one representative image from every 10 years is chosen to simplify the assessment process. Hence, three sample years from the beginning, middle, and end of the time frame for each case study is chosen for accuracy assessment. The two classes that are extracted and presented in binary images are evaluated: built-up and non-urban.
According to the results (Table 4), the overall accuracy of urban/built-up ranges between 92% to 96.29%; and the Kappa value varies from 0.79 to 0.91. The lower user accuracy percentages for New Zealand 2017 and Italy 2017 are due to errors introduced from cloud presence, (even though the clouds are not over the AOI). Differences between sensors (TM versus OLI) do not seem to change the accuracy, as is seen in case of Iran. Thus, for longitudinal studies of large geographic areas, the proposed method provides the flexibility to delineate built-up land.

Change Detection
With the resulting class maps for each year, one application of the proposed methodology is change detection, which is tested for the three earthquake cases here. The change detection method is based on the differences in indices through years. In this method, the change from built-up to non-urban or the opposite is calculated in a yearly interval. The results show the annual spatial variations of urban/built-up surface change at pixel level, and summarized in Table 5. According to these results the estimated damage (% change from pre-event built-up area) is about 32% for Christchurch, 16% for L'Aquila, and about 63% for Bam ( Table 5) that is close to the reports of extensive damages for these earthquakes. Also, results show that current built-up area is more than pre-event status in the study areas for Christchurch and L'Aquila, while total built-up area is still lower than pre-event for Bam. However, the annual growth trend can show the temporal dynamics of recovery process in these cases more in detail, thus, the results of urban land surface change are mapped for each case by pixels. In addition, the curves of percent urban land surface change for the study area are presented demonstrating the temporal trend. The maximum error detected from accuracy assessment for each case is included in the graphs to show the maximum extent of variability in the estimated urban land change. The years with imagery from Landsat 7 ETM+ after 31 May, 2003, (when the Scan Line Corrector (SLC) failed) are removed, since even the overlapped gap-filling images were introducing errors in change detection.

Case 1-Christchurch, New Zealand
The map below show the areas detected as built-up/non-urban before and after the Christchurch earthquake and aftershock in New Zealand (Figure 7). The year of the earthquake and aftershock (2010 and 2011) is added in black, to visualize the change in urban area post-event (Figure 7a). The developments after the earthquake are displayed in purple after eliminating the built-up areas in 2010 and 2011 (i.e., not damaged) that show more construction and development in the western part of Christchurch and in Rolleston (Figure 7b). The orange areas show built-up lands from before the earthquakes, while some of the damaged area is reconstructed, some pockets have remained vacant, especially in the eastern parts, Central City and by the Avon River (marked as 'Red Zones' due to extensive damage from the liquefaction caused by the earthquake [64]). areas in 2010 and 2011 (i.e., not damaged) that show more construction and development in the western part of Christchurch and in Rolleston (Figure 7.b). The orange areas show built-up lands from before the earthquakes, while some of the damaged area is reconstructed, some pockets have remained vacant, especially in the eastern parts, Central City and by the Avon River (marked as 'Red Zones' due to extensive damage from the liquefaction caused by the earthquake [64]).
(a) (b) As is apparent from Figure 8, the urban/built-up land surface change rate in Christchurch is relatively constant before the 2010 earthquake, where there is a drop due to a negative change rate that extends to 2011 due to the aftershock and consequent damages. The change in trend appears to start after 2008, but it is due to the missing data point for 2009 (because of partial cloud coverage over urban lands in the 2009 image). The urban/built-up change trend continues positively after 2011 with a slower rate until 2015 and then shows a faster rate of urban growth.  As is apparent from Figure 8, the urban/built-up land surface change rate in Christchurch is relatively constant before the 2010 earthquake, where there is a drop due to a negative change rate that extends to 2011 due to the aftershock and consequent damages. The change in trend appears to start after 2008, but it is due to the missing data point for 2009 (because of partial cloud coverage over urban lands in the 2009 image). The urban/built-up change trend continues positively after 2011 with a slower rate until 2015 and then shows a faster rate of urban growth. areas in 2010 and 2011 (i.e., not damaged) that show more construction and development in the western part of Christchurch and in Rolleston (Figure 7.b). The orange areas show built-up lands from before the earthquakes, while some of the damaged area is reconstructed, some pockets have remained vacant, especially in the eastern parts, Central City and by the Avon River (marked as 'Red Zones' due to extensive damage from the liquefaction caused by the earthquake [64]). As is apparent from Figure 8, the urban/built-up land surface change rate in Christchurch is relatively constant before the 2010 earthquake, where there is a drop due to a negative change rate that extends to 2011 due to the aftershock and consequent damages. The change in trend appears to start after 2008, but it is due to the missing data point for 2009 (because of partial cloud coverage over urban lands in the 2009 image). The urban/built-up change trend continues positively after 2011 with a slower rate until 2015 and then shows a faster rate of urban growth.  Year %Change in Urban/built-up Area, Christchurch New Zealand

Case 2-L'Aquila, Italy
The maps below show the areas near city of L'Aquila detected as built-up/non-urban for before and after the earthquake (Figure 9). The year of earthquake is added to visualize the change in urban area after the earthquake in black color (Figure 9a). The built-up area since 2010 is marked in purple and shows further development in the eastern part and some reconstruction near downtown areas (Figure 9b). The areas that remained as built-up are removed to only show new development in purple and the areas that have not been built back in orange. This shows that some areas in downtown have not been reconstructed yet.

Case 2-L'Aquila, Italy
The maps below show the areas near city of L'Aquila detected as built-up/non-urban for before and after the earthquake (Figure 9). The year of earthquake is added to visualize the change in urban area after the earthquake in black color (Figure 9.a). The built-up area since 2010 is marked in purple and shows further development in the eastern part and some reconstruction near downtown areas (Figure 9.b). The areas that remained as built-up are removed to only show new development in purple and the areas that have not been built back in orange. This shows that some areas in downtown have not been reconstructed yet.

Case 2-L'Aquila, Italy
The maps below show the areas near city of L'Aquila detected as built-up/non-urban for before and after the earthquake (Figure 9). The year of earthquake is added to visualize the change in urban area after the earthquake in black color (Figure 9.a). The built-up area since 2010 is marked in purple and shows further development in the eastern part and some reconstruction near downtown areas (Figure 9.b). The areas that remained as built-up are removed to only show new development in purple and the areas that have not been built back in orange. This shows that some areas in downtown have not been reconstructed yet.

Case 3-Bam, Iran
The maps below show the areas detected as built-up/non-urban for before and after the Bam earthquake, near city of Bam, Iran. The year of earthquake is added to visualize the change in urban area after the earthquake (Figure 11). The black areas represent undamaged and existing built-up surfaces after the earthquake in 2003 (Figure 11a), and the green areas show the pre-existing urban features that were not detected as built-up in 2003 (i.e., damaged or demolished). The areas that were not demolished in 2003 and remained built-up are subtracted from all years after the earthquake and presented in Figure 11b  The maps below show the areas detected as built-up/non-urban for before and after the Bam earthquake, near city of Bam, Iran. The year of earthquake is added to visualize the change in urban area after the earthquake (Figure 11). The black areas represent undamaged and existing built-up surfaces after the earthquake in 2003 (Figure 11.a), and the green areas show the pre-existing urban features that were not detected as built-up in 2003 (i.e., damaged or demolished). The areas that were not demolished in 2003 and remained built-up are subtracted from all years after the earthquake and presented in Figure 11 The urban/built-up land surface change in Bam is depicted in Figure 12 that shows a growth trend until 2003, which changes after the earthquake. The post-earthquake trend is rather constant until 2009 and shows growth until 2014, followed by an instant increase in growth rate until 2016 but then seems fairly constant until 2019. The curves of percent urban/built-up area change depict the development and growth trend in the three cases. The pre-event imagery is helpful in defining a baseline for the development trend The urban/built-up land surface change in Bam is depicted in Figure 12 that shows a growth trend until 2003, which changes after the earthquake. The post-earthquake trend is rather constant until 2009 and shows growth until 2014, followed by an instant increase in growth rate until 2016 but then seems fairly constant until 2019.

Case 3-Bam, Iran
The maps below show the areas detected as built-up/non-urban for before and after the Bam earthquake, near city of Bam, Iran. The year of earthquake is added to visualize the change in urban area after the earthquake (Figure 11). The black areas represent undamaged and existing built-up surfaces after the earthquake in 2003 (Figure 11.a), and the green areas show the pre-existing urban features that were not detected as built-up in 2003 (i.e., damaged or demolished). The areas that were not demolished in 2003 and remained built-up are subtracted from all years after the earthquake and presented in Figure 11.b to show new development and remaining areas that have not been built back in orange. The purple areas are the newly developed built-up areas since 2004 (a year after the event), which show new constructions in the eastern and southern part of Bam. The urban/built-up land surface change in Bam is depicted in Figure 12 that shows a growth trend until 2003, which changes after the earthquake. The post-earthquake trend is rather constant until 2009 and shows growth until 2014, followed by an instant increase in growth rate until 2016 but then seems fairly constant until 2019. The curves of percent urban/built-up area change depict the development and growth trend in the three cases. The pre-event imagery is helpful in defining a baseline for the development trend before the earthquake, while the post-event imagery shows how much the pre-event trend has been altered after the incident. Figure 13 juxtaposes the land surface change trends together with additional linear forecast lines from pre-event data to demonstrate a comparative measure of how these communities are recovering according to the pre-existing development processes. The slopes of projected trends from pre-event data are 0.0044 for Christchurch, 0.0035 for L'Aquila, and 0.0016 for Bam; which can be compared with each segment of post-earthquake trends to see the relative recovery and development process for different time-frames. On the other hand, the projected trends from 10-years of post-event data, has a slope of 0.0039 for Christchurch, 0.0033 for L'Aquila, and 0.0017 for Bam. These pre-event versus post-event trends indicate that Bam has just surpassed the pre-event growth, while the other two cases have not reached their pre-event growth trends in a decadal comparison. However, the diversity in general growth trends for each case is observable in Figure 13, where Christchurch shows a higher growth speed (steeper slope) for both pre-event and post-event trends. Also, since recovery is a process, looking at specific time-windows can be an instrument for monitoring the recovery process for each case, based on the baseline that is specific to that case from their pre-event trends. For example, within the first 5 years after the earthquake in Christchurch (depicted by blue line in Figure 13) the slope is 0.0016, which is much less that the pre-event growth trend, but the time-window of 5-year to 10-years after the earthquake has a slope of 0.0084 that is more than its post-event decadal growth trend. In other words, the trend of built-up land growth in Bam is caught up with pre-event trend; however, the total built-up area is still lower than pre-event (Table 5), thus still recovering. In contrast, even though the decadal post-event growth trend is slacking in the other two cases, the growth trend in more recent years has surpassed the pre-event trends in Christchurch and just caught up with pre-event trends in L'Aquila, and they both have higher percentage of built-up area when compared to pre-event (Table 5). These comparable curves show the relative trend of urban change that will inform more details about the recovery process if combined with the socio-economic information. The observed differences in the trends need further investigation to find explanatory variables of such variations through time.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 25 before the earthquake, while the post-event imagery shows how much the pre-event trend has been altered after the incident. Figure 13 juxtaposes the land surface change trends together with additional linear forecast lines from pre-event data to demonstrate a comparative measure of how these communities are recovering according to the pre-existing development processes. The slopes of projected trends from pre-event data are 0.0044 for Christchurch, 0.0035 for L'Aquila, and 0.0016 for Bam; which can be compared with each segment of post-earthquake trends to see the relative recovery and development process for different time-frames. On the other hand, the projected trends from 10-years of post-event data, has a slope of 0.0039 for Christchurch, 0.0033 for L'Aquila, and 0.0017 for Bam. These pre-event versus post-event trends indicate that Bam has just surpassed the pre-event growth, while the other two cases have not reached their pre-event growth trends in a decadal comparison. However, the diversity in general growth trends for each case is observable in Figure  13, where Christchurch shows a higher growth speed (steeper slope) for both pre-event and post-event trends. Also, since recovery is a process, looking at specific time-windows can be an instrument for monitoring the recovery process for each case, based on the baseline that is specific to that case from their pre-event trends. For example, within the first 5 years after the earthquake in Christchurch (depicted by blue line in Figure 13) the slope is 0.0016, which is much less that the pre-event growth trend, but the time-window of 5-year to 10-years after the earthquake has a slope of 0.0084 that is more than its post-event decadal growth trend. In other words, the trend of built-up land growth in Bam is caught up with pre-event trend; however, the total built-up area is still lower than pre-event (Table 5), thus still recovering. In contrast, even though the decadal post-event growth trend is slacking in the other two cases, the growth trend in more recent years has surpassed the pre-event trends in Christchurch and just caught up with pre-event trends in L'Aquila, and they both have higher percentage of built-up area when compared to pre-event (Table 5). These comparable curves show the relative trend of urban change that will inform more details about the recovery process if combined with the socio-economic information. The observed differences in the trends need further investigation to find explanatory variables of such variations through time. Depending on the extent of study, these results can be aggregated based on administrative boundaries to have the percent area change in urban surfaces for each administrative unit or neighborhood. A division of change trends by smaller administrative boundaries is specifically helpful in determining and comparing the location of new developments and reconstructions, since the total change trend for a large study area might be misleading as some neighborhoods might not been reconstructed or other new neighborhoods are developing that can be added to the analysis of recovery depending on the scale of study. For example, the percent built-up area change is summarized and presented by districts for the case of Christchurch, New Zealand ( Figure 14). Since these curves depict the percentage of changes in the built-up area, the changes for Christchurch city district shows a higher level of change compared with the other regions in the study area. Even though the post-event growth trend for Christchurch city is slow until 2015 (slope of 0.0037), which is much lower than its pre-event trend with slope of 0.0111, it shows a steep increase after 2015 until 2018 with a slope of 0.0147 that has surpassed its pre-event values. Also noteworthy is the growth trend for the Selwyn District that shows a higher slope in the recent years (slope of 0.0046) when compared with pre-event trend for this district (slope of 0.0011). A ground truthing study can further validate or evaluate the quality of these estimations.  Depending on the extent of study, these results can be aggregated based on administrative boundaries to have the percent area change in urban surfaces for each administrative unit or neighborhood. A division of change trends by smaller administrative boundaries is specifically helpful in determining and comparing the location of new developments and reconstructions, since the total change trend for a large study area might be misleading as some neighborhoods might not been reconstructed or other new neighborhoods are developing that can be added to the analysis of recovery depending on the scale of study. For example, the percent built-up area change is summarized and presented by districts for the case of Christchurch, New Zealand ( Figure 14). Since these curves depict the percentage of changes in the built-up area, the changes for Christchurch city district shows a higher level of change compared with the other regions in the study area. Even though the post-event growth trend for Christchurch city is slow until 2015 (slope of 0.0037), which is much lower than its pre-event trend with slope of 0.0111, it shows a steep increase after 2015 until 2018 with a slope of 0.0147 that has surpassed its pre-event values. Also noteworthy is the growth trend for the Selwyn District that shows a higher slope in the recent years (slope of 0.0046) when compared with pre-event trend for this district (slope of 0.0011). A ground truthing study can further validate or evaluate the quality of these estimations.

Discussion
The application of remote sensing for tracking urban land surface change in order to understand the reconstruction phase of the recovery process after a major disaster can inform decision makers for both recovery planning and resilience measurements. As indicated by the case studies' results, the urban/built-up land surface change could be monitored with the proposed methodology, which can depict the development process and growth trend through time. Even though this method only applies to the changes in the physical built environment and not the socio-economic characteristics of the community, it has the potential to be combined with such data from other sources to provide a more holistic view of the recovery process. In our upcoming research, these measures are applied to a combination of both physical and socio-economic aspects for modeling a more comprehensive recovery process and gauging the relative success of the recovery.
The pre-event imagery is essential to decipher the disaster impact and recovery process, which requires establishing a database and image collection with designed temporal intervals for hazardous regions. Since official records of building and structures are not uniformly available in different regions and countries due to complexities of data collection and regulations, a repository of built environment imagery is a promising alternative to monitor the changes and allocate development resources appropriately. This would enable longitudinal studies of development

Discussion
The application of remote sensing for tracking urban land surface change in order to understand the reconstruction phase of the recovery process after a major disaster can inform decision makers for both recovery planning and resilience measurements. As indicated by the case studies' results, the urban/built-up land surface change could be monitored with the proposed methodology, which can depict the development process and growth trend through time. Even though this method only applies to the changes in the physical built environment and not the socio-economic characteristics of the community, it has the potential to be combined with such data from other sources to provide a more holistic view of the recovery process. In our upcoming research, these measures are applied to a combination of both physical and socio-economic aspects for modeling a more comprehensive recovery process and gauging the relative success of the recovery.
The pre-event imagery is essential to decipher the disaster impact and recovery process, which requires establishing a database and image collection with designed temporal intervals for hazardous regions. Since official records of building and structures are not uniformly available in different regions and countries due to complexities of data collection and regulations, a repository of built environment imagery is a promising alternative to monitor the changes and allocate development resources appropriately. This would enable longitudinal studies of development trends that can be applied to studies other than disaster impacts as well. The importance of this method lies in establishing a baseline based on pre-existing conditions and growth trends in the community for a more realistic expectation of the temporal and spatial recovery trend. As it was seen in the three examples investigated in this study, the differential growth trends from pre-event to post-event provides a relative measure of recovery progress. Hence, application of this methodology for more cases will inform a comparative view of how communities recover through time after a disaster.
Moreover, the method provides the potential to aggregate the built-up surface change results with other layers of information about different indicators-such as the community traits and economic status, building materials, ground motion characteristics, hazard levels or damage status, etc.-for a more holistic view of the recovery process. The quality of rebuilding and reconstruction is missing in this method but can be achieved with ground truthing and field study. Clearly, higher resolution imagery can improve the results and precision of detecting the percent urban/built-up area change, and complement the results of this method. Additionally, the proposed method has flexibility in defining the thresholds with the conditional algebra statements according to local characteristics, which warrants more precise differentiation of features based on environmental variability. However, this flexibility introduces subjectivity in modifying the thresholds that is inevitable and more accurate in comparison with similar approaches. The logic calculation method has showed higher accuracy in similar cases, compared with a band spectral signature analysis, maximum likelihood supervised classification, or a principal component analysis (PCA), while being less subjective or time-consuming [60].
Future research directions would include application of higher spatial, spectral, and temporal resolution data to test the methodology and compare results. In addition, the results of this study can be juxtaposed with the object-based change detection models to further investigate their relative advantages and disadvantages. Finally, with larger number of images the machine learning methodologies would help processing more scenes and test the results for different geographical locations and contexts.

Conclusions
This study employs a conditional algebra in a knowledge-based classifier to extracts the urban/built-up features in three earthquake cases. Three spectral indices, UI, NDVI, and MNDWI, are used in the classifier. The results have an overall accuracy of 92% to 96.29%; and the Kappa value of 0.79 to 0.91 in three study cases.
The method's flexibility in threshold definition with the conditional algebra statements accounts for local characteristics and distinction of features based on environmental variability, which also can be seen as the shortcoming of the method due to the subjectivity in modifying the thresholds. Regardless, the results indicated the potential of the proposed model in delineating the built-up surfaces that can provide essential information about the urban change trend when other sources of data are not available or accessible. Also, this method enables a relative measure of recovery that-once combined with variations in ground motion characteristics, building type distribution, and socio-economic data-allows a more holistic view of the recovery process. According to the findings from the three earthquake cases, total built-up area has increased in two cases of Christchurch and L'Aquila compared with pre-event area, which is also observed from the annual growth trend with a steeper rate in recent years. However, in case of Bam, when comparing recent years to pre-event years, the growth trend has the same rate as pre-event, but the total built-up area is still less than pre-event status. The spatial distribution of built-up area change in all three cases indicated that the post-earthquake developments are not necessarily on the damaged sites, and some areas have not been built back.
Lastly, the curves of percent urban/built-up area change provide a comparable overview of the development and growth trend post-earthquake in the affected regions. Hence, the utilization of annual imagery for detecting urban surface change can inform researchers and policy makers about the temporal and spatial nature of the recovery process and the impacts of the disaster. Application of this methodology for more cases will provide more insight on the observed disparities in the recovery process based on pre-existing trends before the disaster event, both temporally and spatially.