Land Use and Land Cover Changes , and Environment and Risk Evaluation of Dujiangyan City ( SW China ) Using Remote Sensing and GIS Techniques

Understanding of the Land Use and Land Cover (LULC) change, its transitions and Landscape risk (LR) evaluation in earthquake-affected areas is important for planning and urban sustainability. In the present study, we have considered Dujiangyan City and its Environs (DCEN), a seismic-prone area close to the 2008 Wenchuan earthquake (8.0 Mw) during 2007–2018. Five different multi-temporal data sets for the years 2007, 2008, 2010, 2015, and 2018 were considered for LULC mapping, followed by the maximum likelihood supervised classification technique. The individual LULC maps were further used in four time periods, i.e., 2007–2018, 2008–2018, 2010–2018, and 2015–2018, to evaluate the Land Use and Land Cover Transitions (LULCT) using combined remote sensing and GIS (Geographical Information System). Furthermore, multi-criteria evaluation (MCE) techniques were applied for LR mapping. The results of the LULC change data indicate that built-up, agricultural area, and forest cover are the prime categories that had been changed by the natural and anthropogenic activities. LULCT, along with multi-parameters, are suggested to avoid development in fault-existing areas that are seismically vulnerable for future landscape planning in a sustainable manner.


Introduction
The knowledge of spatio-temporal Land Use and Land Cover (LULC) is important for understanding the dynamics of the environments of urban areas that are associated with landscape transitions [1][2][3][4][5][6][7].The local, regional, and global LULC are greatly impacted by natural hazards, such as forest fires, and catastrophic events, such as earthquakes, floods and landslides.Studies of LULC are important for understanding many of the observed phenomena that are responsible for changes [8][9][10][11][12][13][14][15][16][17][18][19][20], landscape changes [21][22][23][24][25][26], landscape fragmentation [27,28], changes in ecosystems [29,30], climate changes [31], the growth of urban areas [32,33] and LULC changes, based on sustainable development [34,35].The term 'LULC change' refers to the human modification of the terrestrial surface of the earth as well as the study of land surface change [36].On the other hand, the term 'Transition' is defined by the process of changing, or the change of something from one form or state to another [37].The LULC pattern of a region is an outcome of natural and socio-economic factors, and their utilization are considered a central component for managing natural resources and monitoring environmental changes [38].Furthermore, reviews of LULC change, the development of land science, change patterns, and causes of changes have been widely studied [15,18,[39][40][41][42][43][44].In the United States of America (USA), the NW-SW are particularly facing very high rates of LULC change, with large declines in agriculture and forest areas.Remote sensing data have been widely used to estimate continuous LULC changes in USA for the period 1973-2000 [45], and during 1973 and 2010 [46].In Canada, LULC data is available at a 30 m spatial resolution nationwide, with an update frequency of five years, and has been considered as the fundamental earth surface attributes associated with geology, hydrology, climatic, and atmospheric parameters [47].According to NOAA, land cover (LC) indicates the physical land type, such as forest or open water, or land-use-landscape for development, conservation, or mixed uses [48].However, an account of 25 years of land cover change was performed in Europe, with the theme of Landscapes in Transition [49].Based on the European Environment Agency (EEA) Report (no: 10/2017), the landscape is undergoing fast and large-scale transformations according to land use changes, and a detailed analysis of land cover provides information about the drivers of the transitions [49].Moreover, the study framework of LULC based on sustainable development in China as discussed by Sui and Ming [35], has attracted the attention of scientists and planners worldwide.Recently, the spatio-temporal dynamics of LULC changes in metropolitan areas over the world have been carried out [50][51][52][53][54][55][56][57] in watershed [58,59], coastal [60], and wildlife sanctuary areas [61] using multi-temporal remote sensing data.
The study of LULC change is very common, and is based on a comparison of present and historical data that are available through ground, airborne, and satellite means.The integrated techniques of remote sensing and GIS provide some of the most accurate and important information over time, which enables us to know the changes and their extent of area, and the pattern of changes.Using optical and microwave remote sensing data, Singh et al. [61][62][63][64][65][66], Dey and Singh [67], Okada et al. [68], and Dey et al. [69] have carried out studies related to surface, ocean, and atmospheric parameters that are associated with the Gujarat earthquake of 26 January 2001.Singh and Singh [70] have used remote sensing data for lineament mapping, and its relation with stress changes in the epicentral areas of the Gujarat earthquake.A detailed analysis of damages, and changes in numerous parameters that are associated with the deadly Wenchuan earthquakes were carried out by many groups [71][72][73][74][75][76][77][78].A number of papers were published by Singh [79] in the special issue of International Journal of Remote Sensing.Recently, a seismic hazard assessment was also carried out in Syria using seismicity, Digital Elevation Model (DEM), slope, active faults, and GIS techniques [80].Many studies have been carried out around the globe using remote sensing and GIS techniques to assess LULC changes, especially in mountainous regions where accessibility is limited [81][82][83][84].LULC change detection and its variations in types can be easily identified in a given location through the spectral signature characteristics over the time period.Recently, changes in the lineament structure, caused by earthquakes, were investigated in South America based on lineament analysis using ASTER (Terra) satellite data [85].In the present study, lineament extraction and analysis have been carried out using Landsat (5-TM and 8-OLI) satellite data.
Sustainable assessment of LULC changes provide environmental, social, and economic dimensions.The role of earth observation datasets (Landsat data) is now well known in land management, monitoring, and degradation; however, such studies have not been carried out in Sichuan province.In the present study, this is the first time that we have applied remote sensing and GIS techniques in the DCEN earthquake-affected areas.Hence, it is vital to draw attention to earthquake-affected cities more closely, to observe historical LULC changes, along with the stress patterns.
In this paper, we have studied LULC change, LULC transitions (LULCT), and landscape risk (LR) evaluations in Dujiangyan City and Environs (DCEN) located in Sichuan province, SW China.
LULC in DCEN has been impacted by the deadly Wenchuan earthquake of 12 May 2008, which caused major landscape changes over the past decade.A comparison and analysis of LULC prior and after the earthquake is important in understanding the changes in spatio-temporal LULC [86][87][88][89][90].After the 2008 Wenchuan earthquake, the Chinese government rebuilt the city during 2008-2010.Some of the surviving people, together with new people from outside, moved in the city.
In the present study, multi-temporal Landsat data for the period 2007-2018 are considered for the detailed analysis.The LULC of the year 2007 is considered before the earthquake (BEQ), and compared with LULC after the earthquake (AEQ) (8.0 M w ) in May 2008 for the analysis of LULC changes.We have also considered the 2010 and 2015 LULC for changes with rebuilding, whereas the current year (2018) is considered for studying the present status of LULC, LULCT, and LR of the study area.Based on LULC changes and LULCT maps, a change matrix (CM) was calculated for the period 2007-2018.These LULCT have been mapped in four different time frames, subtracting earlier images from the recent image (e.g., 2007-2018; 2008-2018, 2010-2018, and 2015-2018) for the change evaluation and the detailed analysis.The disaster chain effect is the main properties of landscape risk.Furthermore, along with these, LR areas have been investigated and evaluated with multiple parameters, such as earthquake epicenter locations, local geology, fault lines, geomorphological features, lineament density (LD), and LULC data.For instance, three different time periods (e.g., 2007 as BEQ, 2008 as AEQ, and 2018 as present scenario) were considered for risk evaluation using remote sensing and GIS tools.
This paper focuses on the main objective for the quantitative evaluation of spatio-temporal changes of LULC, LULCT since 2007-2018, and recent 2018 data was considered as an ancillary data to identify LR areas across the DCEN.However, the broad objectives of this study were to examine three issues: (1) the quantitative evaluation of LULC changes in the DCEN, considering multi-temporal datasets, (2) the quantitative evaluation of LULCT in different time periods, and (3) to identify and evaluate LR areas at the present context in consideration with important multiple environmental parameters.The present study will be useful in understanding the LULC changes, LULCT and LR in the highly earthquake-sensitive areas.The spatial pattern of LULC changes, and present landscape level risk results will help planners and decision makers with landscape planning, environment, and future land management of any area effectively, and in a sustainable manner.Policy responses are needed to support sustainable land management, thus contributing to achieving the United Nations (UN) 2030 Sustainable Development Goals (SDGs).We have considered Goal 11 of the SDGs in the present study which is "Sustainable Cities and Communities" to make cities inclusive, safe, resilient, and sustainable [91].

Study Area
The present study focuses on the LULC patterns of the heterogeneous DCEN is a county-level city of an area of 1208 km 2 , a sub-division of Chengdu located in Sichuan province, which is now a UNESCO world heritage site and has undergone major renovations in the year 2013.The region forms a border, with southern Nagawa Tibet in the north, attracting international tourists.Dujiangyan city (DC) was formerly a county named Guanxian or Guan County, known as "irrigation" County.In the year, 1988 this county was renamed after the Dujiangyan Irrigation System (DIS), as since around 250 BC, this famous irrigation system has provided water for Chengdu.DC is settled back against mountains, with five rivers flowing through the city, and mountain Yulei being located in the middle.Overall, the city is a combination of mountains, rivers, forests, weirs, bridges, and city buildings.The study area DCEN (30 • 53 0" N to 31 • 4 30" N and 103 • 27 5" E to 103 • 43 0" E), is an earthquake-affected city in SW China, 45 km north-west of the provincial capital of Chengdu of Sichuan province (Figure 1).
DCEN has a humid subtropical climate (Koppen Cwa) with cool, dry winters and hot, very wet summers.There are 17 towns in the study area (e.g., Guankou, Xingfu, Puyang, Yuyuan, Chongyi, Tianma, Shiyang, Liujie, Yutang, Zhongxing, Qingchengshan, Longchi, Xujia, Anlong, Daguan, Zipingpu, and Cuiyuehu including two townships (e.g., Xiange, and Hongkou).The study area extends over an area of 318.30 km 2 .The elevation and slope in the study area varies from 668 to 1861 m, with the average altitude being approximately 864.52 m above mean sea level, and 89.66 to 89.99 degrees with an average slope of 89.55 degrees, respectively.The mean high and mean low temperature in the complex heterogeneous landscape is 28 • C in the month of August, and 20 • C in the month of January, based on available records from 1981-2010 (Figure 2a).The highest recorded temperature was observed ≥ 30 °C during April to September.In the same time period, the mean precipitation varies in the range 230-260 mm during summer (July-August), and it further decreased to 170 mm in the month of September, and the mean highest relative humidity being observed to be approximately 85% (Figure 2b).The highest recorded temperature was observed ≥ 30 • C during April to September.In the same time period, the mean precipitation varies in the range 230-260 mm during summer (July-August), and it further decreased to 170 mm in the month of September, and the mean highest relative humidity being observed to be approximately 85% (Figure 2b).

Data
The remote sensing data used in this study was based on the Landsat 5 Thematic Mapper (TM), and Landsat 8 Operation Land Imager (OLI).The spectral resolution and bands of the data used in this study are shown in Table 1.The five temporal data sets were acquired in different time periods, and were obtained from the United States Geological Survey (USGS) Landsat archive: 2007-2018 [92].The changes in LULC, LULCT, and LR were estimated using LANDSAT data.In addition, Shuttle Radar Topographic Mission (SRTM) (30 m resolution) was used for the Digital Elevation Model (DEM) (Figure 1) and retrieved topographic features (i.e., slope).All maps were prepared at a 1: 130,000 scale, combining geological maps, fault lines, epicenter points, LULC, and lineament density (LD) maps to identify the LR areas.The area of interest (AOI) is demarcated in the satellite images with the help of the polygon feature, using the ArcGIS 10.5 software environment, and represents boundary area of the study area.These images were geometrically corrected and projected onto a World Geodetic System (WGS) 1984, into Universal Transverse Mercator (UTM), Zone 48N coordinate system.The pre-processing of Landsat 5 TM and 8 OLI images were conducted according to the following steps.First, the images were calibrated to at-sensor radiance images.Second, radiometric correction was performed for all six satellite images, using an atmospheric correction model to surface reflectance images with the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) tools, based on a MODTRAN radiative transfer module of ENVI 5.3 software.

LULC Mapping Training Sample Setup and Classification Accuracy Method
The heterogeneous landscape in the study area consists of a mixture of built-up area, forest area, agricultural area, reconstructed area, water bodies, and barren land, respectively.Spectral signatures of six different features were drawn on a Landsat image in an integrated manner with distinct line colors (Figure 4b).In the image signature and spectral profile, both used the same color code for clear identification, where: Pink = Built-up area (BU), Dark green = Forest area (F), Orange = Agricultural area (AG); Purple = Reconstructed area (RA), and Blue = Water bodies (WB); and Maroon = Barren land (BL), respectively.The spectral profiles of six different classes were further verified and validated with the USGS spectral library that is available in ENVI spectral library tool box with the ENVI 5.3 version software.To identify the unique pattern of each LULC, 20 signatures class are drawn over image for each class identification.A total of 120 signatures were collected that represent six LULC class used for image classification.Table 2 shows LULC classes and descriptions [93][94][95].A commonly pixel-based supervised image classification with a maximum likelihood classification (MLC) algorithm was used for the analysis [96].Based on this MLC classification approach, five temporal LULC maps were generated using ENVI 5.3 software.The signatures files (polygon file) have been further used for creating random points.For that, polygon-to-point conversion has been performed using the spatial analyst tool box of ArcGIS 10.5 software.For each (b)  A commonly pixel-based supervised image classification with a maximum likelihood classification (MLC) algorithm was used for the analysis [96].Based on this MLC classification approach, five temporal LULC maps were generated using ENVI 5.3 software.The signatures files (polygon file) have been further used for creating random points.For that, polygon-to-point conversion has been performed using the spatial analyst tool box of ArcGIS 10.5 software.For each image, total 120 random points were generated and validated with the corresponding classified image; as well, these points were converted into a Key Hole Mark Up Language (KML) file that is supported by Google Earth (GE).Each chosen point centered at a 30 × 30 m polygon which is one-pixel image size.Both points and corresponding polygons were overlayed in GE for visualization as well as validation.These points were identified and assigned classes based on LULC class types.GE has established itself as a powerful visualizer, and it has been found to be suitable for data validation because of its higher spatial resolution.A similar approach is applied for the rest of the images.The training samples represent the spatial distribution of training points geographically (Figure 4a).Classification accuracy assessments were performed, and a contingency errors matrix table were automatically created using ENVI 5.3 software.

LULC Dynamic Degree Estimation and Transition Matrices Computation Method
To estimate the dynamic degree changes between the various types of LULC to evaluate loss or gain in LULC-type areas, comparing different periods of the dynamic degree model is used to represent the spatio-temporal characteristics.The dynamic estimation is calculated by using the approach of Liu et al. [12][13][14] where, D represents the dynamic degree model, which refers to rate of change; A a is the area in the initial year; A b is the area in the terminal year; T is the temporal scale, which in our case are the time comparisons at 11, 10, 8 and 3 years, respectively.The dynamic degree model served to generate statistical data in a table with D (%), and Gain and Loss (%).Furthermore, transition matrices (TM) tables were computed using LULC data for different years (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).These data (vector format) were analyzed using ArcGIS 10.5 software, to identify areas with LULC changes.The results were exported to a text file and a database using GIS, for statistical analysis.

LULC Transitions Mapping and Analysis Method
In this section, the methodology used in the LULCT mapping and analysis for different periods using dynamic mapping and its corresponding transition matrixes tables [97,98] is described.LULCT is typically performed using LULC themes transformed into raster data.The five LULC maps were used for LULCT mapping and analysis.The ENVI 5.3 thematic change workflow tool was used to depict the dynamics of LULC change that had taken place in the DCEN during the periods 2007-2018.The tool measures the transition dynamics of the LULC class to another class at a certain extent.Firstly, an initial state image was considered as the time 1 (t1) image, and the final stage image as the time 2 (t 2 ) image.Only the changed areas were considered for visualizing the dynamic pattern of LULCT in different time periods.Unclassified areas were eliminated, as they have no valid change or '0' change value to visualize.After the execution, a cleanup operation was performed to eliminate speckling noise from the classified image, with the smooth Kernel size being specified as 3 × 3 pixels; next, the square kernel's center pixel was replaced with the majority class value of the kernel.
In the present study, we used the LULC data for the periods 2007, 2008, 2010, and 2015 as the initial image when compared with the recent 2018 LULC map.By allowing those time frames, LULCT maps for the periods 2007-2018, 2008-2018, 2010-2018, and 2015-2018 were prepared.Export thematic change vectors save the vectors that are created during classification to a shapefile.The supported vector output formats are a shapefile and a geodatabase.The default output area units are in square meters which was converted to square kilometers.

Landscape Risk Evaluation Method
In this study, the LR evaluation methodology [90,96,98] was used to determine the LR exposure, based on the past and present conditions of multiple features of the earth surface that could pose a potential threat to people, property, and the corresponding environment in which they live.Geology is one of the important parameters used for earthquake risk evaluation in DCEN for urban development.Tudes and Yigiter [98] have used this technique for the suitability analysis of Adana City in Turkey, and have suggested that development should take place where there is a low risk of earthquake.However, several steps were involved in our present LR mapping.The LULC and LULCT maps were prepared in the first stage, and lineament features and LD maps were necessary to construct at the 2 stage for LR evaluation.The following steps were followed for lineament construction and LD map preparation.

•
Lineaments were retrieved digitally by using a sufficient values of LINE algorithm parameters of PCI Geomatica 16.0 version software.The default parameters of the LINE algorithm were modified as an image, considering a small area rather than the whole satellite scene.Various computer-aided methods, such as edge detection, thresholding, and curve extraction steps, were carried out over derived PCA images (i.e., PC1 best band chosen) of Landsat 5 TM and Landsat 8 OLI data sets of three different times such as, 2007, 2008, and 2018, respectively.

•
In the next stage, LD mapping, the shapefile of the lineament data was exported and saved in the CADdxf file format, which was then used for lineament direction measurement by creating a rose diagram in the RockWorks 17 software environment.
Recently, the multi-criteria evaluation (MCE) method has been used for assessing and aggregating weighted maps (criterion) based on the expert's knowledge of factor influence and its interactions with LULC [96,97].An LR map is derived using MCE by combining the information from multiple criterions to form a single index of evaluation [90,96,98].
Multiple parameters, such as, LULC maps, lineament features, LD maps, fault lines, and epicenters were overlaid, along with local geology maps, slope, aspect, and DEM, in the integrated mapping.Meanwhile, the geology map was modified and zones I-IV were attributed according to their existing categories [99].At this stage, based on the above important layers, integrated overlay maps were prepared for three different time frames, 2007, 2008, and 2018, respectively.In the next stage, multiple buffer creations (i.e., 5, 10, 15, 20, 25 km distances) and weightings (5,4,3,2,1) were given for each layer, such as important places, past earthquake epicenters, and fault lines, along with LD-2018, and LULC-2018 data sets that were considered for LR evaluations.The weighted score value was assigned, followed by the conditions set as "higher risk value" when 'buffer distance is lower', and "lower risk value" when 'buffer distance is higher'.The recent LULC-2018 map was used for risk analysis of the study area.The LULC class was designated with risk priorities from very high to very low, followed by the weighted score, such as, BU = 1, RA = 1, F = 2, WB = 3, BL = 4, and AG = 5, respectively.
Details of the meaning and the individual weighted assigned score of each layer are given in Supplementary Table S7a-h, and the corresponding individual risk maps are shown in Figure S1a-f.Table S7a-h is based on certain conditions and should be interpreted as follows: if the category of LR is 'very high risk' and 'high risk', it is close to the earthquake epicenter.On the other hand, 'medium risk' represents proximity to the high-risk area.Similarly, low values are indicated by 'low risk' to 'very low risk', i.e., the risk is low because the area is far away from the historical epicenters and geological fault lines.Re-classification was performed for all of the layers, based on the individual assigned weighted values, with the aid of the field calculator of ArcGIS 10.5 software.A new total risk-weighted value field was created, and the weighted score was assigned automatically in the integrated database table, where the geo-processing wizard of the 'union' operator function of ArcGIS 10.5 was applied for the combined data values of all layers.The total LR was defined by the integration of all weighted values using Equation (2): Based on the layer integration approach, the total weighted score was calculated, and in this stage, for quality visual representation of the risk-weighted score, further classification procedures of the Natural Breaks (Jenks) method with five classes were implied.Thereafter, draw quantities using a graduated color ramp were applied with blue-to-red intensification, thus helping to prepare an LR map of the study area.In the final stage, the prepared LR map was masked out based on the study area boundary limits.The overall present study analyses were performed by the integration of multiple software platforms, such as ENVI 5.3, PCI Geomatica 16.0, ArcGIS 10.5, and RockWorks 16, respectively (Figure 5).multiple software platforms, such as ENVI 5.3, PCI Geomatica 16.0, ArcGIS 10.5, and RockWorks 16, respectively (Figure 5).

LULC Mapping, Accuracy Assessment, and Quantity Change Characteristics
Multi-temporal Landsat satellite data were used for image classification, and the reliability of the classification results depends on the overall accuracies of the classified images [25,79,80].In the present study, the accuracy assessment of the classified data of 2007, 2008, 2010, 2015, and 2018 confirmed that the results are reasonable and satisfactory for the present investigation (see Table S1).The overall accuracy is higher compared to the minimum acceptable standard of ≥ 85% , according to the USGS classification scheme.The results of the present analysis indicate that the LULC changes have been accurately identified and extracted.The results of the overall accuracy range from 86-97% for the LULC-classified maps of 2007-2018 (see Table S1), and over Kappa coefficient values of 0.83-0.96;user and producer accuracy for individual classess ranged from 65-100%, and from 38-100%, respectively.
The visual comparison method among the original satellite images, LULC maps, and Google earth high resolution images with a 1 km eye altitude for LULC individual class validation showed robust classification accuracy.Each LULC type area (in km 2 and % form) showed pronounced   S1).The overall accuracy is higher compared to the minimum acceptable standard of ≥ 85% , according to the USGS classification scheme.The results of the present analysis indicate that the LULC changes have been accurately identified and extracted.The results of the overall accuracy range from 86-97% for the LULC-classified maps of 2007-2018 (see Table S1), and over Kappa coefficient values of 0.83-0.96;user and producer accuracy for individual classess ranged from 65-100%, and from 38-100%, respectively.

Dynamic Characteristics of LULC Change Monitoring
The visual comparison method among the original satellite images, LULC maps, and Google earth high resolution images with a 1 km eye altitude for LULC individual class validation showed robust classification accuracy.Each LULC type area (in km 2 and % form) showed pronounced changes over the 11 years (2007-2018) as shown in Table S2.Results indicate that BU increased remarkably from 2007 to 2018, from 40.54 km 2 area in 2007, to 59.13 km 2 in 2018.Large-scale damage was observed in the BU area, due to the deadly Wenchuan earthquake.
Historical population statistics data suggests that the study area had a lower population in 2004 (598 person thousand), and from 2007 (BEQ) in 2008, the population showed an increasing trend.AEQ, the population loss was huge, and it remained low until 2010.From 2011-2014, the population in this area increased from 609 to 619 person per thousand, and the annual population growth rate increased at the rate of 0.6% (Figure 6).
remarkably from 2007 to 2018, from 40.54 km area in 2007, to 59.13 km in 2018.Large-scale damage was observed in the BU area, due to the deadly Wenchuan earthquake.
Historical population statistics data suggests that the study area had a lower population in 2004 (598 person thousand), and from 2007 (BEQ) in 2008, the population showed an increasing trend.AEQ, the population loss was huge, and it remained low until 2010.From 2011-2014, the population in this area increased from 609 to 619 person per thousand, and the annual population growth rate increased at the rate of 0.6% (Figure 6).
AG areas declined from 63.01 km 2 in 2007 to 47.22 km 2 in 2008, mainly due to reconstruction of city and its surroundings; rainfall was also responsible for these changes, and overall, the area declined to 30.66 km 2 in 2018.WB had an increase/decline from 17.60 km 2 in 2007, with overland flows of rivers, and lake water flows creating new channels to a new area of 18.17 km 2 in 2008 (after the 2008 earthquake); this declined to 13.54 km 2 in 2010.One of the reasons for the decline in water area was the accumulation of sediments in the river and channel.WB was found to be increasing in trend at present, and it was 17.14 km 2 in 2018 (see Table S2).From the western, southern, and northeast parts of the mountain to the north part of the study area, F region showed an increase from an area of 124.04 km 2 in 2007, to 175.45 km 2 in 2008 (Figure 7a,b), due to massive plantation activities after the earthquake, to maintain the ecological balance of the area.The areas declined to 148.63 km 2 in 2010 (Figure 7c) to 130.45 km 2 in 2015 (Figure 7d).The F region gradually declined due to reconstruction phase work, whereas the current scenario suggests an increasing state, at 139.25 km 2 in 2018 (Table S2).In contrast, BU area expansion continued until January 2018 (Figure 7e).One of the major geological fault lines, the Hanwang fault, which crosses the northeast region, and another that crosses in the SW-NE direction, flanking towards the NE (e.g., the Pengguan fault), form the BU region, which is seismically active.One of the aftershocks (Mw 5.2) (black dot small red circle) was located in the near half of these two fault lines, close to the BU area (see Figure 3a).AG areas declined from 63.01 km 2 in 2007 to 47.22 km 2 in 2008, mainly due to reconstruction of city and its surroundings; rainfall was also responsible for these changes, and overall, the area declined to 30.66 km 2 in 2018.WB had an increase/decline from 17.60 km 2 in 2007, with overland flows of rivers, and lake water flows creating new channels to a new area of 18.17 km 2 in 2008 (after the 2008 earthquake); this declined to 13.54 km 2 in 2010.One of the reasons for the decline in water area was the accumulation of sediments in the river and channel.WB was found to be increasing in trend at present, and it was 17.14 km 2 in 2018 (see Table S2).
From the western, southern, and northeast parts of the mountain to the north part of the study area, F region showed an increase from an area of 124.04 km 2 in 2007, to 175.45 km 2 in 2008 (Figure 7a,b), due to massive plantation activities after the earthquake, to maintain the ecological balance of the area.The areas declined to 148.63 km 2 in 2010 (Figure 7c) to 130.45 km 2 in 2015 (Figure 7d).The F region gradually declined due to reconstruction phase work, whereas the current scenario suggests an increasing state, at 139.25 km 2 in 2018 (Table S2).In contrast, BU area expansion continued until January 2018 (Figure 7e).One of the major geological fault lines, the Hanwang fault, which crosses the northeast region, and another that crosses in the SW-NE direction, flanking towards the NE (e.g., the Pengguan fault), form the BU region, which is seismically active.One of the aftershocks (M w 5.2) (black dot small red circle) was located in the near half of these two fault lines, close to the BU area (see Figure 3a).From Table 3, we observed that from 2007-2018, the DCEN showed higher dynamic degrees (DD) in RL compared to other periods, suggesting fast changes in LULC in the study area.The BU area showed a positive DD, with a 4.17% per year increase from 2007-2018, whereas, it showed an increase of 5.81% per year and a rapid change of LULC from 2008-2018 (AEQ until the present time), due to city development, reconstruction works, and the growth of the urban area.However, the change rate declined by 1.39% per year over the period 2010-2018, whereas during 2015-2018, BU showed very slow development.On the other hand, the F region had undergone quick changes.
After the earthquake, negative changes of DD, with a −2.03% decrease per year from 2008 to 2018, were observed in the study areas, because of landslides that were associated with the earthquake, two recorded aftershocks, and further expansion of the city in northwards, resulting in forest loss in the DCEN.The scenario of F further slows down, with a negative DD trend, which declined by −0.79% per year over the 2010-2018 time period (Table 3).The LULC change was faster, and it negatively impacted upon the AG of the DCEN area over different time intervals.From 2007 to 2018, it showed a −4.67% decline per year, whereas in the 2008-  After the earthquake, negative changes of DD, with a −2.03% decrease per year from 2008 to 2018, were observed in the study areas, because of landslides that were associated with the earthquake, two recorded aftershocks, and further expansion of the city in northwards, resulting in forest loss in the DCEN.The scenario of F further slows down, with a negative DD trend, which declined by −0.79% per year over the 2010-2018 time period (Table 3).
The LULC change was faster, and it negatively impacted upon the AG of the DCEN area over different time intervals.From 2007 to 2018, it showed a −4.67% decline per year, whereas in the 2008-2018 and 2010-2018 time periods, it decreased by −3.51% and −5.79% per year, respectively.On the contrary, AG severely declined by −12.93% per year during 2015 to 2018.
Table 4 shows temporal LULC area coverage (percentage %), with the gain/loss (%) estimation for the period of 2007-2018.Pronounced LULC changes were observed in different classes.Before earthquake in 2007, the major percentage of LULC in the DCEN were F (39.49%), BL (21.81%),AG (19.80%), and BU (12.74%), whereas in 2008 (AEQ), the data suggested that in particular, the F region showed (54.90%) an increase, due to massive plantation activities by the Chinese and local county government to maintain the ecological balance of the area.Figure 8 shows the two different outputs based on the Table 4 gain and loss (%) of different LULC classes.The rate of change (%) per year in different time periods are calculated.Figure 8a shows the gain and loss (%), and Figure 8b represents the rate of change (%) in LULC classes in different time periods.The overall gain accounted for 4.26% in the F region (Table 4, Figure 8a), and increased by 0.39% per year (Figure 8b) during 2007 to 2018, and the negative impact (−11.15%)(Table 4 and Figure 8a) decreased by 1.11% per year (Figure 8b) from AEQ until now (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).This change occurred due to the pace of city development work and massive urbanization through the expansion of built-up, residential areas, and the industrialization process.Therefore, the faster pace of changes observed with 6.83% positive changes (Table 4, Figure 8a), and it increased by 0.68% per year of BU, and 4.06% (0.41% per year) of RA from 2008-2018, compared to 2010-2018.In this time period (2010-2018), BU and RA showed positive changes of 1.86% (0.23% increase per year), and 2.51% (0.31% increase per year), respectively (Figure 8a,b).On the other hand, BU increased at a slower pace during 2015-2018, compared to earlier periods with a gain of 0.09% (Table 4, Figure 8a) and increased by 0.03% per year (Figure 8b).2018 and 2010-2018 time periods, it decreased by −3.51% and −5.79% per year, respectively.On the contrary, AG severely declined by −12.93% per year during 2015 to 2018.
Table 4 shows temporal LULC area coverage (percentage %), with the gain/loss (%) estimation for the period of 2007-2018.Pronounced LULC changes were observed in different classes.Before earthquake in 2007, the major percentage of LULC in the DCEN were F (39.49%), BL (21.81%),AG (19.80%), and BU (12.74%), whereas in 2008 (AEQ), the data suggested that in particular, the F region showed (54.90%) an increase, due to massive plantation activities by the Chinese and local county government to maintain the ecological balance of the area.Figure 8 shows the two different outputs based on the Table 4 gain and loss (%) of different LULC classes.The rate of change (%) per year in different time periods are calculated.Figure 8a shows the gain and loss (%), and Figure 8b represents the rate of change (%) in LULC classes in different time periods.The overall gain accounted for 4.26% in the F region (Table 4, Figure 8a), and increased by 0.39% per year (Figure 8b) during 2007 to 2018, and the negative impact (−11.15%)(Table 4 and Figure 8a) decreased by 1.11% per year (Figure 8b) from AEQ until now (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).This change occurred due to the pace of city development work and massive urbanization through the expansion of built-up, residential areas, and the industrialization process.Therefore, the faster pace of changes observed with 6.83% positive changes (Table 4, Figure 8a), and it increased by 0.68% per year of BU, and 4.06% (0.41% per year) of RA from 2008-2018, compared to 2010-2018.In this time period (2010-2018), BU and RA showed positive changes of 1.86% (0.23% increase per year), and 2.51% (0.31% increase per year), respectively (Figure 8a,b).On the other hand, BU increased at a slower pace during 2015-2018, compared to earlier periods with a gain of 0.09% (Table 4, Figure 8a) and increased by 0.03% per year (Figure 8b).

LULC Transition Mapping and Matrices Calculation
It was necessary to make a LULCT mapping based in different time periods to identify the transformation of each LULC through a visual, as well as a quantitative manner.LULCT mapping was performed with change in the workflow module of the ENVI 5.3 software, where each earlier LULC map was overlaid with recent one, in order to produce a transition matrix for each time period.Based on the six LULC categories, LULCT maps of four time periods (From→To approach) such as 2007-2018, 2008-2018, 2010-2018, and 2015-2018 have been prepared (Figure 9a-d), and corresponding change matrices of each time scales (where each conversion record was presented) were generated, with the change matrix calculation toolbar of the ENVI 5.3 software (Tables S3-S6).The bold values along the transition diagonals of each matrix table indicate the LULC types from time t 1 to time t 2 (t 2 > t 1 ), with the amount of LULC categories (area in km 2 ) that remained unchanged through each time period scale; while the off-diagonal data represent a transition from one LULC class category to another.Based on the LULC change class of each particular time period, the rate of change (km 2 ) per year of each LULC class was calculated in Microsoft Excel.
The LULCT took place between F to BU, AG to BU, and AL to BL, in particular.In the northeast, central, and southeast area of DCEN, the F transition, including BU, RA, AG, and BL are prominent, whereas the RL transition included BU, AG, and BL, and the AL transition included BL, F, and BU, respectively.Besides these, other minor transitions were also present in the study area (see Table S3).The LULCT during 2007-2018 is shown in Figure 9a.We observed increased area (204.35km 2 ), and changing by 18.58 km 2 per year from 2007-2018 in BU.Moreover, AL and BL also showed positive changes, which held second (80.57km 2 , changed 7.32 km 2 per year) and third (29.72 km 2 , changed 7.32 km 2 per year) positions, respectively.Forest area (F) also increased (10.33 km 2 ), with a change of 2.70 km 2 per year) during this period, along with both RA and WB, which were also observed to have minor changes.Similarly, based on the change matrices in Tables S3-S6, the rate of change per year was calculated for the corresponding time periods.Figure 9b-d

LULC Transition Mapping and Matrices Calculation
It was necessary to make a LULCT mapping based in different time periods to identify the transformation of each LULC through a visual, as well as a quantitative manner.LULCT mapping was performed with change in the workflow module of the ENVI 5.3 software, where each earlier LULC map was overlaid with recent one, in order to produce a transition matrix for each time period.Based on the six LULC categories, LULCT maps of four time periods (From→To approach) such as 2007-2018, 2008-2018, 2010-2018, and 2015-2018 have been prepared (Figure 9a-d), and corresponding change matrices of each time scales (where each conversion record was presented) were generated, with the change matrix calculation toolbar of the ENVI 5.3 software (Tables S3-S6).The bold values along the transition diagonals of each matrix table indicate the LULC types from time t1 to time t2 (t2 > t1), with the amount of LULC categories (area in km 2 ) that remained unchanged through each time period scale; while the off-diagonal data represent a transition from one LULC class category to another.Based on the LULC change class of each particular time period, the rate of change (km 2 ) per year of each LULC class was calculated in Microsoft Excel.
The LULCT took place between F to BU, AG to BU, and AL to BL, in particular.In the northeast, central, and southeast area of DCEN, the F transition, including BU, RA, AG, and BL are prominent, whereas the RL transition included BU, AG, and BL, and the AL transition included BL, F, and BU, respectively.Besides these, other minor transitions were also present in the study area (see Table S3).The LULCT during 2007-2018 is shown in Figure 9a.We observed increased area (204.35km 2 ), and changing by 18.58 km 2 per year from 2007-2018 in BU.Moreover, AL and BL also showed positive changes, which held second (80.57km 2 , changed 7.32 km 2 per year) and third (29.72 km 2 , changed 7.32 km 2 per year) positions, respectively.Forest area (F) also increased (10.33 km 2 ), with a change of 2.70 km 2 per year) during this period, along with both RA and WB, which were also observed to have minor changes.Similarly, based on the change matrices in Tables S3-S6, the rate of change per year was calculated for the corresponding time periods.Figure 9b-d

Landscape Risk Area Indentification
In this section, the present study further utilized an integrated approach based on the method discussed in Section 2.3.4,and we evaluated the LR zone.The study area is characterized by various limitations that appear in the form of faults and lineaments.The heterogeneity of the geological structure, and lithological types in the study areas during the Jurassic-Triassic, Jurassic, Karst, and Quaternary periods, respectively are shown in Figure 10a.These complex structures formed due to complex drainage networks that cut their way into the rock surface, and formed steep slopes (Figure 10b).The aspect (degrees) map was derived from DEM (30 m resolution) to visualize the terrain conditions of the study area (Figure 10c) with active fault lines in the Sichuan basin, including the study area [100].
In contrast, the structure of the study area and its surroundings are dominated by the Beichuan fault [93], Pengguan fault [101], Hanwang fault [101,102], and an unnamed fault [103], respectively.The major strike of these faults is in the NE-SW direction, and the corresponding fault layer is overlaid on the GE for visual representation (see Figure 10d).
For landscape planning and management in the earthquake prone areas, it is important to identify the LR areas, which must be prioritized for future sustainable development.We have further focused on the DCEN area for LR area mapping.In the past, the DCEN was marked as a tectonically active area.Therefore, the area needs vital attention to assess the present risk pattern at the landscape level.Lineament data extraction was necessary along with combinations of other ancillary data sets.

Landscape Risk Area Indentification
In this section, the present study further utilized an integrated approach based on the method discussed in Section 2.3.4,and we evaluated the LR zone.The study area is characterized by various limitations that appear in the form of faults and lineaments.The heterogeneity of the geological structure, and lithological types in the study areas during the Jurassic-Triassic, Jurassic, Karst, and Quaternary periods, respectively are shown in Figure 10a.These complex structures formed due to complex drainage networks that cut their way into the rock surface, and formed steep slopes (Figure 10b).The aspect (degrees) map was derived from DEM (30 m resolution) to visualize the terrain conditions of the study area (Figure 10c) with active fault lines in the Sichuan basin, including the study area [100].
In contrast, the structure of the study area and its surroundings are dominated by the Beichuan fault [93], Pengguan fault [101], Hanwang fault [101,102], and an unnamed fault [103], respectively.The major strike of these faults is in the NE-SW direction, and the corresponding fault layer is overlaid on the GE for visual representation (see Figure 10d).
For landscape planning and management in the earthquake prone areas, it is important to identify the LR areas, which must be prioritized for future sustainable development.We have further focused on the DCEN area for LR area mapping.In the past, the DCEN was marked as a tectonically active area.Therefore, the area needs vital attention to assess the present risk pattern at the landscape level.Lineament data extraction was necessary along with combinations of other ancillary data sets.Lineament data were further used for lineament density (LD) maps over the period 2007-2018.The LD maps were prepared based on five classes, from Very low to Very high-density distribution in three time periods, 2007 as BEQ, 2008 as AEQ, and 2018 as the present scenario, respectively (Figure 10e-g).The LD values changed and slightly increased in AEQ, categorized as Very High Density, which was 3.95 km/km 2 in 2007 (BEQ), and 4.04 km/km 2 in 2008 (AEQ), respectively; the current stress pattern shows almost double the BEQ values (8.27 km/km 2 : Very High Density) which is reflected in the recent LD map of 2018 (Figure 10g).

Landscape Risk Area Indentification
In this section, the present study further utilized an integrated approach based on the method discussed in Section 2.3.4,and we evaluated the LR zone.The study area is characterized by various limitations that appear in the form of faults and lineaments.The heterogeneity of the geological structure, and lithological types in the study areas during the Jurassic-Triassic, Jurassic, Karst, and Quaternary periods, respectively are shown in Figure 10a.These complex structures formed due to complex drainage networks that cut their way into the rock surface, and formed steep slopes (Figure 10b).The aspect (degrees) map was derived from DEM (30 m resolution) to visualize the terrain conditions of the study area (Figure 10c) with active fault lines in the Sichuan basin, including the study area [100].
In contrast, the structure of the study area and its surroundings are dominated by the Beichuan fault [93], Pengguan fault [101], Hanwang fault [101,102], and an unnamed fault [103], respectively.The major strike of these faults is in the NE-SW direction, and the corresponding fault layer is overlaid on the GE for visual representation (see Figure 10d).
For landscape planning and management in the earthquake prone areas, it is important to identify the LR areas, which must be prioritized for future sustainable development.We have further focused on the DCEN area for LR area mapping.In the past, the DCEN was marked as a tectonically active area.Therefore, the area needs vital attention to assess the present risk pattern at the landscape level.Lineament data extraction was necessary along with combinations of other ancillary data sets.Lineament data were further used for lineament density (LD) maps over the period 2007-2018.The LD maps were prepared based on five classes, from Very low to Very high-density distribution in three time periods, 2007 as BEQ, 2008 as AEQ, and 2018 as the present scenario, respectively (Figure 10e-g).The LD values changed and slightly increased in AEQ, categorized as Very High Density, which was 3.95 km/km 2 in 2007 (BEQ), and 4.04 km/km 2 in 2008 (AEQ), respectively; the current stress pattern shows almost double the BEQ values (8.27 km/km 2 : Very High Density) which is reflected in the recent LD map of 2018 (Figure 10g).The integrated layer overlay was applied in three different time periods, which provided information about the change in the landscape in the DCEN (see Figure 10h-j), thus indicating a fast forward move towards the necessity of LR map preparation.For the LR evaluations, the LULC maps and Landsat image source were used for lineament extraction.Thereafter, lineament direction maps of the concerned time periods, along with epicenter points (main event, aftershock event) and fault lines, were prepared with GIS mapping tools.All of the required individual layers have been initially evaluated and assigned according to the degree of risk weightage (see Table S7a-h).
Table 5 represents the LR area statistics details and its categorical risk level information prepared on the developed LR map.The final LR map is created (Figure 11a) based on the integrated layers weighted score, as discussed in Section 2.3.4.Based on the overall weighting, 14 is assigned as 'lowest', and 25 is marked as the 'highest' score.The quantile natural break (Jenks) method with a five-class range has been applied in the ArcGIS 10.Table 5 shows a maximum 149.23 km 2 area, marked as a high-risk area (see Figure 11a), which 46.88% of the total area of DCEN.The second highest area, 119.54 km 2 (37.56%), was identified as the medium LR area.On the other hand, the 25.37 km 2 (7.97%) area was marked as the very highrisk category (Figure 11a).The of the rest areas are shared by very low to low risk categories, which accounted for 4.59% and 3.00%, respectively, of the total risk area of DCEN (Table 5, Figure 11a).The medium risk zone is located in the NE, E, W, Central, and SE and SW parts of the DCEN (Figure 11a), close to adjoining fault lines, nearby WB, and BU where population density is extremely high.However, the associated LR is subject to change according to LD change and different LULC pattern changes; for example, the BU expansion towards fault lines and F clearance, and the shrinkage of AG, etc. Finally, the LR map was converted into the GE kml format, which is used for visual representation of the risk area in a real-earth perspective at different scales (Figure 11b).The integrated layer overlay was applied in three different time periods, which provided information about the change in the landscape in the DCEN (see Figure 10h-j), thus indicating a fast forward move towards the necessity of LR map preparation.For the LR evaluations, the LULC maps and Landsat image source were used for lineament extraction.Thereafter, lineament direction maps of the concerned time periods, along with epicenter points (main event, aftershock event) and fault lines, were prepared with GIS mapping tools.All of the required individual layers have been initially evaluated and assigned according to the degree of risk weightage (see Table S7a-h).
Table 5 represents the LR area statistics details and its categorical risk level information prepared based on the developed LR map.The final LR map is created (Figure 11a) based on the integrated layers weighted score, as discussed in Section 2.3.4.Based on the overall weighting, 14 is assigned as 'lowest', and 25 is marked as the 'highest' score.The quantile natural break (Jenks) method with a five-class range has been applied in the ArcGIS 10.Table 5 shows a maximum 149.23 km 2 area, marked as a high-risk area (see Figure 11a), which covers 46.88% of the total area of DCEN.The second highest area, 119.54 km 2 (37.56%), was identified as the medium LR area.On the other hand, the 25.37 km 2 (7.97%) area was marked as the very high-risk category (Figure 11a).The of the rest areas are shared by very low to low risk categories, which accounted for 4.59% and 3.00%, respectively, of the total risk area of DCEN (Table 5, Figure 11a).The medium risk zone is located in the NE, E, W, Central, and SE and SW parts of the DCEN (Figure 11a), close to adjoining fault lines, nearby WB, and BU where population density is extremely high.However, the associated LR is subject to change according to LD change and different LULC pattern changes; for example, the BU expansion towards fault lines and F clearance, and the shrinkage of AG, etc. Finally, the LR map was converted into the GE kml format, which is used for visual representation of the risk area in a real-earth perspective at different scales (Figure 11b).

Discussion
In developing and developed countries, the LULC changes studied in earthquake-affected regions have received much less attention and documentation in the past.Our study area, of Sichuan province in SW China, particularly the DCEN areas, have been frequently hit by earthquakes, landslides, etc.The Wenchuan earthquake that hit the study area extensively damaged buildings/structures and the landscape.This area has received little attention in the past, due to its relatively slight economic and infrastructure development, and even after the earthquake.After the 2008 earthquake, the Chinese Government focused on development in the area, with restructuring, redevelopment and expansion works occurring due to the inflow migration of people in city areas from neighboring rural areas.Only limited studies have been carried out to evaluate the landscape risk (LR).Therefore, the present assessments of LULC changes, LULCT, and LR evaluation are urgently needed to maintain the sustainability of SW China.
This study first time analyzed the LULC changes during 2007-2018 in the DCEN areas of SW China.The study shows an overall accuracy in the range range 86-97% for the LULC classification during 2007-2018 (see Table S1).The scenarios of changes reflect the LULC pattern of the earthquakeaffected areas.Our results show building reconstruction and BU area expansion with large-scale restructuring and rebuilding of the city, initiative untaken by the People's Republic of China government AEQ, which continued until the end of 2010.People who survived were forced to reenter the earthquake-affected areas during and after the reconstruction work was completed.This development was further accelerated and implemented in the surrounds of the urban area.The LULC changes impacted are observed in the N, NE, and SW directions of the study area, where important places such as Puyangzhen, Xujiazhen, Juyuanzhen, and Zhongxingzhen, respectively, are located.These areas will be further used for developmental work, including BU area expansion, to support more populations in the DCEN.This has been evident in the population statistics data from 2004-2014, showing an increasing trend compared to earlier times (see Figure 6).
The encroachment in AG areas was prominent, with positive DD changes, as observed in the periods 2008-2018 and 2010-2018, respectively while presently it remains as BL.However, the remaining LULC class, such as WB, has significantly improved during the periods 2010-2018 and 2015-2018, with positive trends, compared to its earlier period, which showed negative trends in DD (see Table 3).The changes of this LULC class is regulated by most of the river channels and upper dams located in the high mountain areas, which were badly damaged by the earthquake, thus increasing the water flow with high sedimentation rate, ultimately blocking the pathways of regular

Discussion
In developing and developed countries, the LULC changes studied in earthquake-affected regions have received much less attention and documentation in the past.Our study area, of Sichuan province in SW China, particularly the DCEN areas, have been frequently hit by earthquakes, landslides, etc.The Wenchuan earthquake that hit the study area extensively damaged buildings/structures and the landscape.This area has received little attention in the past, due to its relatively slight economic and infrastructure development, and even after the earthquake.After the 2008 earthquake, the Chinese Government focused on development in the area, with restructuring, redevelopment and expansion works occurring due to the inflow migration of people in city areas from neighboring rural areas.Only limited studies have been carried out to evaluate the landscape risk (LR).Therefore, the present assessments of LULC changes, LULCT, and LR evaluation are urgently needed to maintain the sustainability of SW China.
This study first time analyzed the LULC changes during 2007-2018 in the DCEN areas of SW China.The study shows an overall accuracy in the range range 86-97% for the LULC classification during 2007-2018 (see Table S1).The scenarios of changes reflect the LULC pattern of the earthquake-affected areas.Our results show building reconstruction and BU area expansion with large-scale restructuring and rebuilding of the city, initiative untaken by the People's Republic of China government AEQ, which continued until the end of 2010.People who survived were forced to re-enter the earthquake-affected areas during and after the reconstruction work was completed.This development was further accelerated and implemented in the surrounds of the urban area.The LULC changes impacted are observed in the N, NE, and SW directions of the study area, where important places such as Puyangzhen, Xujiazhen, Juyuanzhen, and Zhongxingzhen, respectively, are located.These areas will be further used for developmental work, including BU area expansion, to support more populations in the DCEN.This has been evident in the population statistics data from 2004-2014, showing an increasing trend compared to earlier times (see Figure 6).
The encroachment in AG areas was prominent, with positive DD changes, as observed in the periods 2008-2018 and 2010-2018, respectively while presently it remains as BL.However, the remaining LULC class, such as WB, has significantly improved during the periods 2010-2018 and 2015-2018, with positive trends, compared to its earlier period, which showed negative trends in DD (see Table 3).The changes of this LULC class is regulated by most of the river channels and upper dams located in the high mountain areas, which were badly damaged by the earthquake, thus increasing the water flow with high sedimentation rate, ultimately blocking the pathways of regular channels; thereafter, water flow was diverted into the new channels.We have also calculated the gain/loss (%) as well as the rate of change per year in the corresponding time periods.The results show that the LULC changes were driven by natural activities (such as earthquake-induced infrastructure loss), and also by human activities in the affected areas (through massive reconstruction work, BU area expansion, and forest loss, etc.).
Moreover, based on the LULC change maps, we have further examined and generated LULCT maps during 2007-2018 for the DCEN areas.The results showed the largest overall increase in LULC in BU 18.58 km 2 per year based, on 2007-2018.AG and BL also show positive changes, with increases of 7.32 km 2 and 2.70 km 2 per year, respectively.Forest area (F) also enhanced by 0.94 km 2 per year during this period, along with both RA and WB, which showed minor changes (see Table S3).After the deadly earthquake in 2008, the government and other agencies built new constructions, and also renovated damaged buildings in DCEN areas, which is clearly evident from the LULCT (2008-2018) map, and this increased by 24.26 km 2 per year (see Table S4, Figure 9b).During this period, the population increased, which increased the demand for BU areas.The spatio-temporal pattern of LULCT shows the conversion scenario among the six major LULC types in the next two time periods.In 2010-2018, BU increased by 2.81 km 2 per year (see Table S5, Figure 9c), and in recent times (2015-2018) it has increased by 4.86 km 2 per year (see Table S6, Figure 9d), respectively.
However, during 2010-2018, the LULC changes were associated with government policy, infrastructure development, built-up area expansion, industrial set ups, and so on.RA and AG were the two top major categories in the LULC change, at 191.87 km 2 and 30.41 km 2 , respectively, and they changed by 23.98 km 2 and 3.80 km 2 per year, respectively (see Table S5).On the contrary, BU and F changes were associated with 22.46 km 2 (increased 2.81 km 2 per year) and 21.15 km 2 (2.64 km 2 per year) land area, respectively (see Table S5).These changes were associated with the post-disaster improvement plan in the DCEN areas.
On the other hand, the LULCT change was slow during 2015-2018 compared to the period of 2010-2018 with respect to BU (changing 4.86 km 2 per year), RA (changing 8.71 km 2 per year), and AG (0.96 km 2 per year); except that F, WB, and BL showed an increasing trend and changed by 68.57km 2 , 0.88 km 2 and 7.97 km 2 per year, respectively (see Table S6).In the WB case, more channels formed after the earthquakes, and through these channels water diversion activities may have occurred that may have helped local farmers to enhance their agricultural activities.Moreover, the data suggests that the local government took major initiatives for the AG expansion within the close vicinity of the BU periphery zones, such as Puyangzhen, Xujiazhen, the E-SE part of Dujiangyan, and the E-W part of Zhongxingzhen, respectively.These areas are marked as promising areas that are considered for AG expansion, which supports the needs of growing populations (see Figure 9d).In the next, landscape factors, including structural distribution and integrated overlay mapping of the study area has been performed successfully (see Figure 10a-j).
Lineament density mapping in the DCEN area provided information about change in the stress pattern during 2007-2018.The main basic assumption was that, the higher the LD values, the more rock that was fractured in the surface area, and vice versa with regard to the stress.Furthermore, higher rock fracturing indicates high risk and impending disasters such as earthquakes, earthquake induced landslides, regular landslides, and liquefactions in the near future, in the DCEN areas.
In addition, this study finally prepared and analyzed an LR map (see Figure 11a), and geological zones I, II and IV are identified as risk categories (see Table 5).The landscape types such as BU, F, BL, AG, and WB, respectively, are found to be under stress.The geological zones I-IV are marked, and Landscape types such as WB, BU, F, AG, and BL are under risk, based on their priority levels, buffer distance from fault lines, and past earthquake epicenter buffer zone fusion.The pace of development occurred closely towards the direction of numerous fault lines, which makes the city vulnerable, and this has been given little attention in the past for risk assessments.In our present LR mapping, the BU is identified as being in a high LR risk category, and it needs great attention, to avoid any further damage in this area.For the data visualization, the final LR map is displayed in the GE viewer (see Figure 11b).
Sustainable assessment of LULC changes are often analyzed in four dimensions, such as economic, environmental, institutional, and social.Our present study deals with the environmental sustainability of DCEN areas, and it can be achieved through LULC change trajectories and LR analysis.Moreover, the integrated database and geo-integrated techniques have helped us to explore the disaster-prone area in a sustainable manner.Therefore, by following the aims and objectives of the present work, we successfully evaluated the DCEN areas and make suggestions to make city development more safe, resilient, and sustainable in the future, which may follow the UN sustainability development Goal No. 11.

Conclusions
The present study utilized spatial geo-integrated analysis techniques, which provide information about the pronounced changes in the epicentral and surrounding areas that are associated with the earthquake at a regular time interval.The study area is located in the highly seismic-prone areas in SW China, which is presently undergoing high-level stress, as revealed from the present study through LR analysis.We have studied the affected DCEN areas and quantitatively analyzed the period 2007-2018 using multi-temporal Landsat data.The results discussed in the present paper will of great help to city planners, local administrations, and earthquake engineers for future planning and development and technical-know-how about the landscape changes in the DCEN areas.
The results discussed in the present study show pronounced LULC changes in the 2008 earthquake-affected area (DCEN).The area has undergone major changes during the course of the earthquake time period and the AEQ.However, after the development planning was adopted by the Chinese government to rejuvenate the city again to its earlier status, with the inclusion of more people to live in the newly provided infrastructure, and better employment facilities in the regions.The city expansion and LULCT were not done properly, but the present results may provide some guidelines for the development of risk-prone areas.This change is evident from the recent satellite data of 2018, showing the present stress level towards the development of LR maps of the DCEN areas.
In conclusion, we have drawn from multi-temporal Landsat data for the period 2007-2018, and importantly, we have generated LULC maps of the DCEN based on MLC techniques with six LULC classes.The LULC maps are further successfully utilized to figure out the LULCT that had taken place in the DCEN environment.Moreover, based on the integrated layers, a final LR map was successfully developed for the DCEN area.As the DCEN is identified as risk-prone area, integrated LULC plans should be made by city planners and local administrators with a special consideration towards the presently developed LR map.A study of LULC changes with integrated LR mapping techniques have found ideal and significant relevant impact for sustainable developments in the DCEN areas, which might be effective for sustainable developments of any world regions.In the future, construction and developmental activities must be avoided in the existing fault lines.Furthermore, planners should focus on this area within the existing laws of LULC management, to prevent further landscape damage.Thus, planners and decision makers must work together for future LULC management, development activities, and allocation of land.These three things must be planned, and new policy frameworks should be adjusted with minimizing the risk factors in a sustainable manner.
Sustainability 2018, 10, x FOR PEER REVIEW 4 of 32 extends over an area of 318.30 km 2 .The elevation and slope in the study area varies from 668 to 1861 m, with the average altitude being approximately 864.52 m above mean sea level, and 89.66 to 89.99 degrees with an average slope of 89.55 degrees, respectively.The mean high and mean low temperature in the complex heterogeneous landscape is 28 °C in the month of August, and 20 °C in the month of January, based on available records from 1981-2010 (Figure2a).

Figure 1 .
Figure 1.Location of the study area (Dujiangyan city and its Environs); Location of Sichuan province and the map of the study area are shown in the upper panel.The lower panel shows the Digital Elevation Model (in meters) with important locations points and geological faults lines crossing over the study area.
The study area has suffered flood and landslide disasters in the past, and during 2008-2009, a major earthquake and aftershocks occurred.The frequency of earthquakes of Mw 5.0 and above in the epicentral and surrounding regions are shown for the periods 12 May 2008 to 15 June 2009 in Figure 3.

Figure 1 .
Figure 1.Location of the study area (Dujiangyan city and its Environs); Location of Sichuan province and the map of the study area are shown in the upper panel.The lower panel shows the Digital Elevation Model (in meters) with important locations points and geological faults lines crossing over the study area.
The study area has suffered flood and landslide disasters in the past, and during 2008-2009, a major earthquake and aftershocks occurred.The frequency of earthquakes of M w 5.0 and above in the epicentral and surrounding regions are shown for the periods 12 May 2008 to 15 June 2009 in Figure 3.

Figure 2 .
Figure 2. Historical climatic information of the study area from 1981-2010 represents different temperature information: (a) Temperature data (record high, mean high, and mean low); and (b) mean precipitation and mean relative humidity.Source: China Meteorological Data Service Center, and Weather China; http://en.m.wikipedia.org/wiki/Dujiangyan.

Figure 3 .Figure 2 . 32 Figure 2 .
Figure 3. Frequency of earthquakes in the study area during 12 May 2008 to 15 June 2009; (a) top right indicates the location of earthquakes points only within white bounding box; left panel shows the location in Google earth, and bottom right panel shows the locations of earthquakes of two different magnitudes, overlaid over the 2008 Landsat 5 TM false color image (bands used: R, G, B: 5,4,3), and

Figure 3 .Figure 3 .
Figure 3. Frequency of earthquakes in the study area during 12 May 2008 to 15 June 2009; (a) top right indicates the location of earthquakes points only within white bounding box; left panel shows the location in Google earth, and bottom right panel shows the locations of earthquakes of two different magnitudes, overlaid over the 2008 Landsat 5 TM false color image (bands used: R, G, B: 5,4,3), and
where, LR[T_Weight] refers to the Landscape Risk[Total Weight]; [IP_B_Weight] refers to the Important Places buffer weight; [EQ_B_Weight] refers to the Earthquake Point buffer weight; [F_B_Weight] refers to Fault Based buffer weight; [Risk_Geol] refers risk factor of geological zone; [LD2018_Risk] refers Lineament density based risk factor; and [LULC2018_Risk] indicates the Land use and Land cover 2018 risk factor value.

Figure 5 .
Figure 5.A schematic framework method used in the present study.

Figure 5 .
Figure 5.A schematic framework method used in the present study.
3.1.1.LULC Mapping, Accuracy Assessment, and Quantity Change Characteristics Multi-temporal Landsat satellite data were used for image classification, and the reliability of the classification results depends on the overall accuracies of the classified images [25,79,80].In the present study, the accuracy assessment of the classified data of 2007, 2008, 2010, 2015, and 2018 confirmed that the results are reasonable and satisfactory for the present investigation (see Table

Figure 6 .
Figure 6.Population distribution and annual population growth rate of the study area (2004-2014).Note: Annual population growth rate is calculated based on the standard formula.(Source: National Bureau of Statistics-Population County Level Region; Data updated on 10 March 2016, www.ceicdata.com).

Figure 6 .
Figure 6.Population distribution and annual population growth rate of the study area (2004-2014).Note: Annual population growth rate is calculated based on the standard formula.(Source: National Bureau of Statistics-Population County Level Region; Data updated on 10 March 2016, www.ceicdata.com).

Figure 7 .
Figure 7. LULC maps of the study area from 2007-2018; (a) LULC-2007, (b) LULC-2008, (c) LULC-2010, (d) LULC-2015, and (e) LULC-2018.3.1.2.LULC Dynamic Degree, and Gain and Loss Estimations from 2007-2018From Table3, we observed that from 2007-2018, the DCEN showed higher dynamic degrees (DD) in RL compared to other periods, suggesting fast changes in LULC in the study area.The BU area showed a positive DD, with a 4.17% per year increase from 2007-2018, whereas, it showed an increase of 5.81% per year and a rapid change of LULC from 2008-2018 (AEQ until the present time), due to city development, reconstruction works, and the growth of the urban area.However, the change rate declined by 1.39% per year over the period 2010-2018, whereas during 2015-2018, BU showed very slow development.On the other hand, the F region had undergone quick changes.

Figure 8 .
Figure 8. LULC change distribution in four different time periods; (a) Gain/Loss (%) of each LULC classes, (b) Per-year rate of change (%) of LULC classes in different time periods.Note: BU-Built-Up area; F-Forest area; AG-Agricultural area; RA-Reconstructed area; WB-Water bodies; and BL-Barren Land.

32 Figure 8 .
Figure 8. LULC change distribution in four different time periods; (a) Gain/Loss (%) of each LULC classes, (b) Per-year rate of change (%) of LULC classes in different time periods.Note: BU-Built-Up area; F-Forest area; AG-Agricultural area; RA-Reconstructed area; WB-Water bodies; and BL-Barren Land.
Lineament data were further used for lineament density (LD) maps over the period 2007-2018.The LD maps were prepared based on five classes, from Very low to Very high-density distribution in three time periods, 2007 as BEQ, 2008 as AEQ, and 2018 as the present scenario, respectively (Figure 10e-g).The LD values changed and slightly increased in AEQ, categorized as Very High Density, which was 3.95 km/km 2 in 2007 (BEQ), and 4.04 km/km 2 in 2008 (AEQ), respectively; the current stress pattern shows almost double the BEQ values (8.27 km/km 2 : Very High Density) which is reflected in the recent LD map of 2018 (Figure 10g).

Figure 11 .
Figure 11.LR mapping and visual display of the DCEN on Google Earth; (a) LR category map at the present context, (b) LR map overlaying on Google earth as a kml file for visualization and sharing.

Figure 11 .
Figure 11.LR mapping and visual display of the DCEN on Google Earth; (a) LR category map at the present context, (b) LR map overlaying on Google earth as a kml file for visualization and sharing.

Table 1 .
Details of data used in this study.

Table 3 .
Dynamic degree (%) estimation between different times based on the LULC total area coverage (km 2 ) from 2007-2018.

Table 3 .
Dynamic degree (%) estimation between different times based on the LULC total area coverage (km 2 ) from 2007-2018.
Note: BU-Built-Up area; F-Forest area; AG-Agricultural area; RA-Reconstructed area; WB-Water bodies; and BL-Barren Land.Source: Image statistical results retrieved from ENVI 5.3 and ArcGIS 10.5 software, based on LULC images from 2007-2018 and DD (%) values calculated by authors.
Note: BU-Built-Up area; F-Forest area; AL-Agricultural area; RL-Reconstructed area; WB-Water bodies; and BL-Barren Land.The per-year change is not displayed in this table, and the calculated results are referred to in the relevant text.

Table 4 .
Gain/Loss (%) estimation between different times, based on the LULC total area coverage (%) from 2007-2018.Barren Land.The per-year change is not displayed in this table, and the calculated results are referred to in the relevant text.

Table 5 .
Landscape risk (LR) identification and patterns of change based on integrated datasets in different zones.

Table 5 .
Landscape risk (LR) identification and patterns of change based on integrated datasets in different zones.