Monitoring Quarry Area with Landsat Long Time-Series for Socioeconomic Study

Quarry sites result from human activity, which includes the removal of original vegetation and the overlying soil to dig out stones for building use. Therefore, the dynamics of the quarry area provide a unique view of human mining activities. Actually, the topographic changes caused by mining activities are also a result of the development of the local economy. Thus, monitoring the quarry area can provide information about the policies of the economy and environmental protection. In this paper, we developed a combined method of machine learning classification and quarry region analysis to estimate the quarry area in a quarry region near Beijing. A temporal smoothing based on the classification results of all years was applied in post-processing to remove outliers and obtain gently changing sequences along the monitoring term. The method was applied to Landsat images to derive a quarry distribution map and quarry area time series from 1984 to 2017, revealing significant inter-annual variability. The time series revealed a five-stage development of the quarry area with different growth patterns. As the study region lies on two jurisdictions—Tianjin and Hebei—a comparison of the quarry area changes in the two jurisdictions was applied, which revealed that the different policies in the two regions could impose different impacts on the development of a quarry area. An analysis concerning the relationship between quarry area and gross regional product (GRP) was performed to explore the potential application on socioeconomic studies, and we found a strong positive correlation between quarry area and GRP in Langfang City, Hebei Province. These results demonstrate the potential benefit of annual monitoring over the long-term for socioeconomic studies, which can be used for mining decision making.


Introduction
After the 1980s in China, especially in Beijing, there was a rapid expansion in economic growth with cities and urban areas consuming a huge amount of building materials.Between 1982 and 2015, the population in Beijing increased from 9.23 million to 21.71 million with an annual increase of 2.6%.In the same period, the built-up area in Beijing grew from 303.91 km 2 to 1401.01 km 2  (http://data.stats.gov.cn/english/easyquery.htm?cn=E0103).The growing population, the urban sprawl, and the improvement of infrastructures during the past three decades in Beijing have meant that quarries have progressively exploited the surrounding mountains.More specifically, quarries in Langfang City of Hebei Province went through an enormous expansion because of the great development of the major cities around this region.Metropolises such as Beijing and Tianjin experienced a wide upsurge in building construction and the materials of these buildings came mainly from the nearby quarries considering the transportation fees.Overall, the continuous urban sprawl has given rise to long-lasting quarry activity in the nearby mountains.Quarries contribute a great part to local government revenue, which pushes the industry of building material to a leading position in Hebei's economy.Quarries developed in a highly unsupervised manner during its early age, and its change was mainly subject to supply and demand at that time.With rapidly increasing house prices and in the absence of sound supervision from local governments, quarry areas may undergo continuous growth.However, with the intervention from local government, the direction and magnitude of its dynamics may change dramatically due to changes in economic policy and environmental initiatives.Therefore, monitoring quarry-rich regions at multi-temporal and high spatial resolution is important for local governments to gain information that is needed for the policies of economic and environmental protection.
Surface mining activities range from large open-cast coal and base metal mines to much smaller aggregate (rock, gravel, and sand), industrial minerals (potash, clay), and building material (granite, stone and marble) quarries [1].Mining operations are done in a place called a quarry or quarry site.A successful monitoring approach for evaluating surface mining processes and their dynamics on a regional scale requires observations with frequent temporal coverage over a long time to differentiate natural changes from those associated with human activities.To meet such challenges, accurate and up-to-date information is needed.Since the 1970s, analog aerial photographs have been used to map spatial changes in mined areas [2].Today, remote sensing with its capability of establishing calibrated observations over decadal-scale time periods and at sufficiently high spatial resolution (≤30 m) can provide this information to achieve better insight into quarry dynamics.Several local and regional studies have focused on the dynamics of quarries using a variety of remotely-sensed datasets with various spatial and temporal resolutions [1,[3][4][5][6][7][8][9][10].The findings from these studies have demonstrated the ability of remote sensing to detect and analyze the magnitude and spatial changes of natural environments as a result of mining activities [3,4].However, most prior remote sensing studies focusing on quarry area changes have been temporally limited, relying on the comparison of imagery from two to three time slices (e.g., [1,5,6,9,10]).Studies [5,6] with high spatial resolution images (<5 m) usually focused on rather short time scales, ranging from several months to a couple of years due to the limited extent and availability of very high-resolution imagery, while other analyses typically span a couple of decades.On the other hand, several studies have turned to multi-temporal analyses either with multi-sourced data [8] or a single source of data (Landsat) [7].Moreover, the diversity of data sources and the acquisition timing has resulted in limited comparability among the different studies.In this study, we used Landsat data series over a long time and with sufficient observations to fully investigate quarry dynamics.Unlike [7,8], we aimed to obtain an annual result of quarry area, trying to further investigate the relationship between quarry dynamics and local economy using the acquired quarry area series.
Since the entire Landsat archive is freely provided by USGS, it becomes viable to exploit this valuable and consistent data source to access multi-scaled land surface dynamics in most regions.Recently, with increasing computation capacities and novel processing techniques, more and more studies focusing on monitoring high-resolution land cover changes have shifted focus toward a high frequency multi-temporal analysis using the entire Landsat archive, with over 40 years of continuous acquisitions.Examples include mostly forestry applications.For example, Kennedy, Cohen, and Schroeder [11] applied a trajectory-based change detection method to a stack of Landsat images from 1984 to 2004 to monitor forest disturbance dynamics; Hansen et al. [12] mapped global forest loss and gain from 2000 to 2012 based on Landsat data; Rosenau, Scheinert, and Dietrich [13] monitored glacial flow velocities both at decadal and seasonal time scales using Landsat imagery; and Macander et al. [14] mapped snow persistence for Northwest Alaska using Landsat images from 1985 to 2011.These studies are predominantly based on the analysis of multi-spectral indices (MSI) or the original spectral bands.In Fraser et al. [15], robust linear trend analysis of temporal Landsat Tasseled Cap (TC) index time-series was proposed and applied in different studies of land changes, such as post-fire forest recovery [16], lake dynamics [17], and the land cover change classification [18].Recently, Google Earth Engine [19] has greatly facilitated the temporal applications of the Landsat archive with its powerful cloud computing ability, such as the monitoring of urban boundaries [20] and woody vegetation clearing [21].
We noticed that quarry changes might have a potential relationship with the regional economy for it is a result of human activities and economic development.However, with the scarce observations so far, the potentially strong relationship between quarry area and economic parameters have not been sufficiently analyzed.Few studies have examined the correlation between quarry area and GDP or gross regional product (GRP), whereas a large number of statistical studies have demonstrated that there is a positive correlation between GDP and remotely sensed nighttime light data captured by Operational Linescan System (OLS) of the US Defense Meteorological Satellite Program (DMSP) [22][23][24][25].This research indicates that it is feasible to use nighttime lights to estimate GDP or GRP.Therefore, we were inspired to research whether quarry area captured by Landsat satellites has a positive correlation with GDP or GRP and investigate the possibility of using daytime remote sensing data for socioeconomic study.
In this study, we first examined the areal expansion of aggregate quarries over 34 years in the study region in Northern China, using remote sensing techniques.Due to the small extension of quarries in one year, we chose a single imagery captured by the Landsat 5 or Landsat 8 satellite for each year between 1984 and 2017.We then applied a modified random forest classification for quarry range extraction and analyzed quarry dynamics in the study region.The studied regions governed by the Hebei and Tianjin governments allow us to compare different quarry dynamics responding to different policies in these two places.Our analysis aims to detect decadal-scale changes in the quarry sites, as well as the relationship between quarry area and socioeconomic parameters.We try to demonstrate in this study that quarry area extracted from remote sensing data can be used to sense society changes, for many land cover changes should be a result of human activities and social development.

Study Area
The study area we chose is part of the middle range of the Yanshan Mountains, located in the northeast of Langfang City, Hebei Province, China, where it has been exploited for aggregate since the 1980s.The study area is situated on the border of Beijing, Tianjin, and Hebei, east of Beijing city (Figure 1).Its location is in the north of China between latitudes 39 • 59 and 40 • 5 N and longitudes 117 • 5 and 117 • 14 E and covers an area of about 145 km 2 (Figure 1).The highest altitude in this area is 527 m above sea level and the lowest is 4 m.
There is a large amount of high-quality dolomite ore, limestone, shale, and other minerals under the hills in the study area.In the early 1980s, some locals had exploited these minerals for profit.With the development of urbanization in the surrounding region, a large amount of stone building materials have been required for a growing number of construction projects in the cities every year.Considering the cost of transportation, most stone building materials such as aggregates and cement come from quarry sites around the cities.The growing rock market triggered the rapid extraction of rocks from these hills.With the rise in housing prices, the price of stone in mining enterprises soared, resulting in mining activity that can be seen everywhere in the hills inside the study area.We choose this region as the study area because it contains quarry sites in both Langfang and Tianjin.This facilitates the comparison of quarry dynamics in the two regions, and the study period can effectively represent the quarry changes both in the early period when local people dominated the exploitation activity and in the later period when the government has managed the production and recovery in the quarry sites.The changes in the study area also reflect the many socioeconomic parameters and some big events that happened around this region.

Data
While there are numerous satellites orbiting the globe with a multitude of instruments onboard that explore various bandwidths in the electromagnetic spectrum, Landsat-based images are ideally suited for surface feature analyses because of the long history of earth observation and rich bands with ample spectral information [26].Note that Landsat 5 Thematic Mapper (TM) operational imaging ended in November 2011 [27]; subsequently, Landsat 8 was launched in February 2013 [28].Therefore, only Landsat 7 imagery with missing data are available for 2012 due to the failure of the Enhanced Thematic Mapper Plus (ETM+) Scan Line Corrector (SLC) on Landsat 7 post-2003.In this study, an annual observation is enough for tracking quarry areal extents because the processes that drive quarry area growth or decline are slow enough to be captured well with an annual repeat frequency.Thanks to the high temporal frequency of capturing images, we can easily find ideal annual observations with low cloud coverage from 1984 to present.Eventually, 33 Landsat 5 (TM) and Landsat 8 (Operational Land Imager-OLI) images were picked from the US Geological Survey's (USGS) EarthExplorer platform (https://earthexplorer.usgs.gov),representing years from 1984 to 2017 (except 2012), and WRS Paths 122 to 123 and Rows 32.All scenes were ordered as surface reflectance (SR) data.Landsat 5 SR data provided by USGS are generated using the atmospheric correction chain called Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [29].In 2010, USGS Earth Resources Observation and Science Center (EROS) released on-demand LEDAPS processing for generating SR products for the entire Landsat 5 archives (http://espa.cr.usgs.gov).On the other hand, Landsat 8 SR data are generated from the Landsat Surface Reflectance Code (LaSRC) [30].Currently, both Landsat 5 and Landsat 8 SR data can be ordered using EarthExplorer.Because of the SLC failure of Landsat 7, we chose not to pick data for 2012 in order to avoid the influence of the missing scan line data in ETM+ images.The images from TM and OLI sensors have a common

Data
While there are numerous satellites orbiting the globe with a multitude of instruments onboard that explore various bandwidths in the electromagnetic spectrum, Landsat-based images are ideally suited for surface feature analyses because of the long history of earth observation and rich bands with ample spectral information [26].Note that Landsat 5 Thematic Mapper (TM) operational imaging ended in November 2011 [27]; subsequently, Landsat 8 was launched in February 2013 [28].Therefore, only Landsat 7 imagery with missing data are available for 2012 due to the failure of the Enhanced Thematic Mapper Plus (ETM+) Scan Line Corrector (SLC) on Landsat 7 post-2003.
In this study, an annual observation is enough for tracking quarry areal extents because the processes that drive quarry area growth or decline are slow enough to be captured well with an annual repeat frequency.Thanks to the high temporal frequency of capturing images, we can easily find ideal annual observations with low cloud coverage from 1984 to present.Eventually, 33 Landsat 5 (TM) and Landsat 8 (Operational Land Imager-OLI) images were picked from the US Geological Survey's (USGS) EarthExplorer platform (https://earthexplorer.usgs.gov),representing years from 1984 to 2017 (except 2012), and WRS Paths 122 to 123 and Rows 32.All scenes were ordered as surface reflectance (SR) data.Landsat 5 SR data provided by USGS are generated using the atmospheric correction chain called Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [29].In 2010, USGS Earth Resources Observation and Science Center (EROS) released on-demand LEDAPS processing for generating SR products for the entire Landsat 5 archives (http://espa.cr.usgs.gov).On the other hand, Landsat 8 SR data are generated from the Landsat Surface Reflectance Code (LaSRC) [30].Currently, both Landsat 5 and Landsat 8 SR data can be ordered using EarthExplorer.
Because of the SLC failure of Landsat 7, we chose not to pick data for 2012 in order to avoid the influence of the missing scan line data in ETM+ images.The images from TM and OLI sensors have a common spatial resolution (30 m) and (six) overlapping spectral bands: blue, green, red, near-infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2).All common spectral bands were used for analysis while the remaining bands were excluded from further processing.Images from the first Landsat sensor generation (Multispectral Scanner/MSS) were not taken into consideration at this point because of their coarser spatial resolution and lower spectral fidelity.The image selection was filtered to acquisition dates in May as we observed that Landsat images of the study area always had no or very little cloud in May among these years.Finally, we got 33 scenes with a cloud coverage near to zero in the study area.Detailed image ID and acquired dates are shown in Appendix A.
For our analysis, we used GRP data of Langfang City from 1990 to 2016.All data came from the Hebei Economic Yearbook [31], which contains the annual economic development of Hebei Province as recorded by the Statistics Bureau of Hebei Province.

Methods
The method developed in this study follows a workflow of automatic image clipping and indices calculating, and subsequent classification and calculation of the quarry area based on a modified classification and quarry region analysis.The methodological framework is presented in Figure 2 and described in the following four subsections: pre-processing (Section 3.2.1),classification (Section 3.2.2),post-processing (Section 3.2.3),and calculation of quarry area (Section 3.2.4).
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 19 spatial resolution (30 m) and (six) overlapping spectral bands: blue, green, red, near-infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2).All common spectral bands were used for analysis while the remaining bands were excluded from further processing.Images from the first Landsat sensor generation (Multispectral Scanner/MSS) were not taken into consideration at this point because of their coarser spatial resolution and lower spectral fidelity.The image selection was filtered to acquisition dates in May as we observed that Landsat images of the study area always had no or very little cloud in May among these years.Finally, we got 33 scenes with a cloud coverage near to zero in the study area.Detailed image ID and acquired dates are shown in Appendix A.
For our analysis, we used GRP data of Langfang City from 1990 to 2016.All data came from the Hebei Economic Yearbook [31], which contains the annual economic development of Hebei Province as recorded by the Statistics Bureau of Hebei Province.

Methods
The method developed in this study follows a workflow of automatic image clipping and indices calculating, and subsequent classification and calculation of the quarry area based on a modified classification and quarry region analysis.The methodological framework is presented in Figure 2 and described in the following four subsections: pre-processing (Section 3.2.1),classification (Section 3.2.2),post-processing (Section 3.2.3),and calculation of quarry area (Section 3.2.4).

Pre-Processing
After being downloaded from USGS, all images were clipped into the study area.Then, we calculated several multispectral indices for each scene, serving as proxies for different land surface properties.We included the Landsat specific Tasseled Cap indices of brightness (TCB), greenness (TCG), and wetness (TCW) [32][33][34], as well as other general indices, such as the Normalized Difference Vegetation Index (NDVI) [35], the Normalized Difference Water Index (NDWI) [36], and the Normalized Difference Building Index (NDBI) [37].They were chosen to reflect a variety of surface characteristics such as vegetation status, surface moisture, or buildings.The TC indices were calculated with the sensor-specific formulas for reflectance data.The calculation was performed using the following equations:

Pre-Processing
After being downloaded from USGS, all images were clipped into the study area.Then, we calculated several multispectral indices for each scene, serving as proxies for different land surface properties.We included the Landsat specific Tasseled Cap indices of brightness (TCB), greenness (TCG), and wetness (TCW) [32][33][34], as well as other general indices, such as the Normalized Difference Vegetation Index (NDVI) [35], the Normalized Difference Water Index (NDWI) [36], and the Normalized Difference Building Index (NDBI) [37].They were chosen to reflect a variety of surface characteristics such as vegetation status, surface moisture, or buildings.The TC indices were calculated with the sensor-specific formulas for reflectance data.The calculation was performed using the following equations: Then, the six calculated index layers were stacked together with the six original bands to get a 12-band scene for each year.

Classification
We first applied a supervised machine learning classification approach, which separated the image into two target classes: quarry class and non-quarry class.Random forest [38] (RF) has been established as one of the most accurate and widely used algorithms [39][40][41][42] for remote sensing and other classification applications.One of its advantages is that it produces an estimation of classification accuracy based on out-of-bag cross-validation (OOB-CV) [43][44][45].However, OOB-CV strongly overestimates model accuracy when classifying remote sensing imagery using training areas with several pixels, which was recently demonstrated in [46].The authors of [46] think the reason is that pixels in a training patch are not independent-they have a strong spectral similarity-whereas the OOB-CV takes all samples as being independent.To deal with the problem, a modification of the RF algorithm, splitting the training patches instead of the pixels that compose them, is proposed in [46] and implemented in the R package called SDRF (Spatial Dependence Random Forest).For the classification process, we used the modified random forest function in the SDRF package.In this function, there are two important parameters that should be controlled to acquire a better classification result: ntree (number of trees to grow), and mtry (number of variables randomly sampled as candidates at each split).Since the random forest classifier is computationally efficient and does not overfit, the number of trees can be as large as possible [47].However, the classification no longer improved as the number of trees increased above a threshold, as demonstrated in some studies investigating the sensitivity of the RF classifier to the number of trees [42,48].A default value of ntree for remote sensing imagery classification is 500.Higher mtry will result in stronger individual decision trees but with an enhanced correlation between trees the model accuracy is reduced [49].Usually, a square root of the total number of variables is used in the classification task.

Training Sampling
We selected 20~50 training sample areas within and in close proximity to our study sites with respect to a set of remote sensing imagery with higher resolution, which mainly came from Google Earth.The ground truth selection process is based on a stratified manual sampling.Since the non-quarry class can have a variety of different surface conditions, we increased its sample size to attribute for this variety.We assumed that quarry range changes slightly over several years and its change limits to its boundary.Therefore, once a sample is identified as an inside part of the quarry according to high-resolution images, we can assume this sample remains unchanged in several near years unless it shows an extensive spectral change on these years' images.The Landsat data and in situ data, as well as high-resolution images, were used for determining suitable locations of the manually selected and determined ground truth locations.

Classification Method Details
The six common spectral bands (blue, green, red, NIR, SWIR1, and SWIR2), as well as the six calculated indices (TCB, TCG, TCW, NDVI, NDBI, NDMI), were used as input features for the pixel-based classification process.The SDRF classification model was trained for each image in every year using the respective training samples.The two parameters ntree and mtry were determined with a test based on the modified out-of-bag (M-OOB) accuracy in this study and we found that the number of 500 for ntree and 4 for mtry worked well for the classification task.As the quarry area covered only a small fraction of the study area before 1995, we used the class weights for the models before 1995 according to the prior of the sample areas.
The classification output contains a hard classification output with the two pre-defined classes (Table 1).This is an ensemble result voted by all decision trees in the RF classifier.The random forest can also output a ratio result for each class, which is calculated as the proportion of trees voting for one class against all votes.Each ratio result is output as a confidence layer for each class, which contains the classification probability (confidence index) of each class in a range from 0 to 1.The SDRF-internal M-OOB accuracy, Cohen's Kappa, and area under the ROC curve (AUC) were used for quality assessment of the classification.

Post-Processing
The output from the random forest is a per-pixel classification with two classes.The per-pixel results tend to be noisy due primarily to the effects of pixel-scale misclassification.These effects include couples of pixels among vegetation area where naturally bare stones are misclassified as quarry.There are also some speckles in the classification results, which is common to many per-pixel classification results.Thus, a post-classification editing step is necessary.
The presence of bright topographic features dominated by cement and stones poses a challenge when trying to extract quarry location.For example, cement ground, cement roads, and the depositories for raw rock materials in the flat terrain can be mistaken by the analysis algorithm as quarry sites as they present a very similar ground feature because of their similar surface materials.Thus, we need to introduce some degree of filtering that would eliminate anomalous pixels from the quarry pixels.For example, quarry sites only exist near or on hill slopes where the altitude is higher than other flat terrain.Therefore, we decided to process the classification result images with an altitude mask.We obtained the altitude mask from the ALOS DSM data by taking an area higher than 50 m as a valid field and then dilating this field by 9 pixels in order to cover quarry hollows at the foot of the hills.Figure 3 illustrates the altitude mask we processed from DSM image.Only quarry pixels inside the valid area were considered as true quarry area.

Temporal Smoothing
After all scenes are processed using the above steps, we applied a temporal smoothing to exclude some non-quarry pixels in the quarry class, such as dusty land caused by wind around the quarry in certain years and also to obtain a consistent change result.Once other classes are changed to quarry (always with the removal of vegetation and surface soil), it will remain as quarry for deep exploration of stones for the following years.For the reverse situation, once a quarry changes to other classes, whether naturally or artificially, it would not change back to a quarry in a short time because it is obviously forgotten or protected by humans.Based on this, it is safe to assume that outliers in the temporal sequence of confidence index is very possible a misclassified one.For each class, we first composed a confidence sequence for each pixel by putting together the confidence index of each year in order.Then, we applied a median filter to the sequences in order to remove outliers and obtain gently changing sequences along the observation period.An example of this method is expressed in Figure 4 where the lines indicate the confidence index of the quarry class during the monitoring period.Obviously, an outlier in 2005 is successfully removed using the temporal smoothing strategy.

Calculation of Quarry Area
During the final step, we calculated the areal extent of quarries in the study area.We noted that unlike other naturally occurring changes, such as land becoming a lake, quarries are exploited by humans in a short time, which means we cannot distinguish a transitional state of the land between quarry and other land-cover types.Therefore, it is reasonable that we obtained a per-pixel

Temporal Smoothing
After all scenes are processed using the above steps, we applied a temporal smoothing to exclude some non-quarry pixels in the quarry class, such as dusty land caused by wind around the quarry in certain years and also to obtain a consistent change result.Once other classes are changed to quarry (always with the removal of vegetation and surface soil), it will remain as quarry for deep exploration of stones for the following years.For the reverse situation, once a quarry changes to other classes, whether naturally or artificially, it would not change back to a quarry in a short time because it is obviously forgotten or protected by humans.Based on this, it is safe to assume that outliers in the temporal sequence of confidence index is very possible a misclassified one.For each class, we first composed a confidence sequence for each pixel by putting together the confidence index of each year in order.Then, we applied a median filter to the sequences in order to remove outliers and obtain gently changing sequences along the observation period.An example of this method is expressed in Figure 4 where the lines indicate the confidence index of the quarry class during the monitoring period.Obviously, an outlier in 2005 is successfully removed using the temporal smoothing strategy.

Temporal Smoothing
After all scenes are processed using the above steps, we applied a temporal smoothing to exclude some non-quarry pixels in the quarry class, such as dusty land caused by wind around the quarry in certain years and also to obtain a consistent change result.Once other classes are changed to quarry (always with the removal of vegetation and surface soil), it will remain as quarry for deep exploration of stones for the following years.For the reverse situation, once a quarry changes to other classes, whether naturally or artificially, it would not change back to a quarry in a short time because it is obviously forgotten or protected by humans.Based on this, it is safe to assume that outliers in the temporal sequence of confidence index is very possible a misclassified one.For each class, we first composed a confidence sequence for each pixel by putting together the confidence index of each year in order.Then, we applied a median filter to the sequences in order to remove outliers and obtain gently changing sequences along the observation period.An example of this method is expressed in Figure 4 where the lines indicate the confidence index of the quarry class during the monitoring period.Obviously, an outlier in 2005 is successfully removed using the temporal smoothing strategy.

Calculation of Quarry Area
During the final step, we calculated the areal extent of quarries in the study area.We noted that unlike other naturally occurring changes, such as land becoming a lake, quarries are exploited by humans in a short time, which means we cannot distinguish a transitional state of the land between quarry and other land-cover types.Therefore, it is reasonable that we obtained a per-pixel

Calculation of Quarry Area
During the final step, we calculated the areal extent of quarries in the study area.We noted that unlike other naturally occurring changes, such as land becoming a lake, quarries are exploited by humans in a short time, which means we cannot distinguish a transitional state of the land between quarry and other land-cover types.Therefore, it is reasonable that we obtained a per-pixel classification with just two classes.However, in the quarry margins, quarry expansion tends to occur on a sub-pixel scale (<30 m), resulting in mixed pixels in these areas.In order to obtain an accurate area result, we calculated the area of the internal zone and margin zone in different ways.
Based on the hard classification results, edges of quarries were identified.Then, we defined two zones to obtain a sub-pixel level accuracy in the quarry area calculation.The quarry margins are defined by the edges being dilated by 1 pixel to capture the sub-pixel expansion of the quarry.To obtain the internal zone, we used the hard-classified quarry mask and applied a morphological erosion to reduce the quarry's radius by 1 pixel to avoid the quarry margin.
The two defined zones were treated differently.The internal zone was regarded as pure pixels and its area was therefore calculated as a full quarry area.For margins, we chose a more dynamic sub-pixel analysis approach to account for its mixed pixels.Its area was represented by the classification probability (confidence index), which was assumed to be the fraction of each endmember/class and they were tested for plausibility on high-resolution imagery, as well as field measurements and observations, which were available for selected locations.These confidence indices were used as a weighing factor of quarry class within each pixel of the quarry margins.For example, a pixel in the margin zone has a confidence index of 0.705 for a quarry, representing 70.5% of the 900 m per-pixel; hence, 634.5 m 2 is calculated as quarry within this pixel.Figure 5 is an example taken from the classification result of 2002 to illustrate the calculation process.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 19 classification with just two classes.However, in the quarry margins, quarry expansion tends to occur on a sub-pixel scale (<30 m), resulting in mixed pixels in these areas.In order to obtain an accurate area result, we calculated the area of the internal zone and margin zone in different ways.
Based on the hard classification results, edges of quarries were identified.Then, we defined two zones to obtain a sub-pixel level accuracy in the quarry area calculation.The quarry margins are defined by the edges being dilated by 1 pixel to capture the sub-pixel expansion of the quarry.To obtain the internal zone, we used the hard-classified quarry mask and applied a morphological erosion to reduce the quarry's radius by 1 pixel to avoid the quarry margin.
The two defined zones were treated differently.The internal zone was regarded as pure pixels and its area was therefore calculated as a full quarry area.For margins, we chose a more dynamic sub-pixel analysis approach to account for its mixed pixels.Its area was represented by the classification probability (confidence index), which was assumed to be the fraction of each endmember/class and they were tested for plausibility on high-resolution imagery, as well as field measurements and observations, which were available for selected locations.These confidence indices were used as a weighing factor of quarry class within each pixel of the quarry margins.For example, a pixel in the margin zone has a confidence index of 0.705 for a quarry, representing 70.5% of the 900 m 2 per-pixel; hence, 634.5 m 2 is calculated as quarry within this pixel.Figure 5

Classification Accuracy
The classification of two defined classes yielded a very good separation on the ground truth data of the defined classes.M-OOB and Cohen's Kappa, as well as AUC for the quarry class, are listed in Table 2, which suggest that the classifier archived a high accuracy using the parameters defined in Section 3.2.2.

Quarry Time Series
During the final step, the annual areal extent of the quarry was calculated.The quarry time series was produced after passing all images through the workflow that we described earlier spanning the period from 1984 to 2017, i.e., about 33 years of coverage.The image series had a gap in 2012 and, therefore, we actually attained a 34-year time series of quarry area.
The quarry area dynamics are presented as line charts (Figure 6).The left line chart indicates the area change of the quarry during the past 34 years.Generally, the area increased before 2013 and decreased during the following 3 years.In May 2017, the quarry experienced a great expansion again.We then portrayed the changes in the quarry boundaries in each period as follows.Each map delineates the quarry ranges before and after one period.In these maps, the blue area represents unchanged quarry regions, the red area represents extended quarry regions, and the cyan area represents recovered quarry regions.
Figure 7a,b shows the two stages of change in the first period, changes between 1984 and 1992 and changes between 1992 and 2002.In the early stage of quarry activities, since 1984, quarries had a scattered distribution.In the beginning, local people were motivated to exploit the rocks for the great profits from the building materials.They extracted stones and rocks without any planning and never considered the amount of prospecting for rocks in the quarry sites.Gradually they abandoned some smaller quarries and flocked to areas where plentiful rocks suitable for construction were extracted.By 1992, the scattered quarries in places that were rich in rocks were expanded and then merged into several large quarry sites.Others sites were abandoned and then re-covered with vegetation naturally in a short time for these quarries had only been minimally exploited.Then, from 1992 to 2002 in Figure 7b, the abandoned quarries continued their recovery with vegetation cover while the large quarry area kept expanding.By 2002, the scattered quarry sites can be distinguished as three main continuous quarry regions as shown in Figure 7f: Quarry A lies in the northwest, Quarry B lies in the south, and Quarry C lies in the northeast.During the 8 years between 2002 and 2010, each of the three quarry areas expanded rapidly, which is clearly illustrated by Figure 7c.The decrease mainly happened inside Quarry C; other small decreasing parts in Quarry B were primarily on the south edge of the quarry, which was exploited earliest.From 2010 to 2013 in Figure 7d, the area increased only slightly on the north side of Quarry B. Most quarry sites remained unchanged during this period.Figure 7e demonstrates that the quarry areas in all three regions decreased between 2013 and 2016, although a small expansion occurred in Quarry C and the western part of Quarry B. We should point out that the small decreasing parts inside Quarry B were newly built buildings.
From the line charts in Figure 6, two special years are identified for presenting very different tendencies: 2008 and 2017.Between May 2007 and August 2008, the quarry area had almost no change compared to its nearby years (Figure 8a).While the main trend was still toward a decrease in the three quarry regions from June 2016 to May 2017 (Figure 8b), a very large newly exploited area in the west of Quarry B resulted in an increment in the total area.Again, the small speckles inside Quarry B were a result of intense internal rebuilding.We then portrayed the changes in the quarry boundaries in each period as follows.Each map delineates the quarry ranges before and after one period.In these maps, the blue area represents unchanged quarry regions, the red area represents extended quarry regions, and the cyan area represents recovered quarry regions.
Figure 7a,b shows the two stages of change in the first period, changes between 1984 and 1992 and changes between 1992 and 2002.In the early stage of quarry activities, since 1984, quarries had a scattered distribution.In the beginning, local people were motivated to exploit the rocks for the great profits from the building materials.They extracted stones and rocks without any planning and never considered the amount of prospecting for rocks in the quarry sites.Gradually they abandoned some smaller quarries and flocked to areas where plentiful rocks suitable for construction were extracted.By 1992, the scattered quarries in places that were rich in rocks were expanded and then merged into several large quarry sites.Others sites were abandoned and then re-covered with vegetation naturally in a short time for these quarries had only been minimally exploited.Then, from 1992 to 2002 in Figure 7b, the abandoned quarries continued their recovery with vegetation cover while the large quarry area kept expanding.By 2002, the scattered quarry sites can be distinguished as three main continuous quarry regions as shown in Figure 7f: Quarry A lies in the northwest, Quarry B lies in the south, and Quarry C lies in the northeast.During the 8 years between 2002 and 2010, each of the three quarry areas expanded rapidly, which is clearly illustrated by Figure 7c.The decrease mainly happened inside Quarry C; other small decreasing parts in Quarry B were primarily on the south edge of the quarry, which was exploited earliest.From 2010 to 2013 in Figure 7d, the area increased only slightly on the north side of Quarry B. Most quarry sites remained unchanged during this period.Figure 7e demonstrates that the quarry areas in all three regions decreased between 2013 and 2016, although a small expansion occurred in Quarry C and the western part of Quarry B. We should point out that the small decreasing parts inside Quarry B were newly built buildings.
From the line charts in Figure 6, two special years are identified for presenting very different tendencies: 2008 and 2017.Between May 2007 and August 2008, the quarry area had almost no change compared to its nearby years (Figure 8a).While the main trend was still toward a decrease in the three quarry regions from June 2016 to May 2017 (Figure 8b), a very large newly exploited area in the west of Quarry B resulted in an increment in the total area.Again, the small speckles inside Quarry B were a result of intense internal rebuilding.

Uncertainties and Errors
For the classification step, we used a supervised learning method, which required a set of accurate and enriched samples.As we classified scenes of each year along the past 34 years, sampling was required for each scene and the accuracy of the results depends on the quality of the samples.Deriving operative samples from temporally adjacent samples and classifiers was demonstrated to work well in this study.However, further studies on other datasets or monitoring tasks are needed to confirm the conclusion.One more uncertainty is that when discussing the influence of social events and local policies on quarry changes, our analysis is based on the comparison of quarry areas in two jurisdictions considering some well-known events and policies.Maybe integrating some kind of quantization of the social events and policies into the analysis can present their close connection more clearly.

Social Discussion
Prior work has documented the feasibility of the Landsat archive for monitoring land cover changes and analyzing their dynamics in a high-resolution multi-temporal way [14].However, studies on quarry monitoring remain temporally limited.Whether using several decades of spanned space-borne imagery or higher resolution images focusing on a short time scale, prior studies have generally relied on a comparison of imagery from two to three time slices [1,5,6,9,10].Hence, with the focus on scarce observations so far, the socioeconomic reasons behind these quarry changes have not been thoroughly analyzed.In this study, we used a long time series data from the Landsat archive to fully investigate the annual change in the quarry area.In this section, as the annual monitored result of the quarry has been derived, we try to discuss the factors driving the quarry development from a socioeconomic viewpoint and demonstrate the potential benefits of annual monitoring in a socioeconomic study.

Comparison of Two Regions
The study area is located on the border of Beijing, Tianjin, and Hebei, making a comparison of quarry development among these areas possible.Since the development of quarry extraction is influenced deeply by human activities, it is also reasonable to discuss the quarry area extension and contraction with the different policies in the three regions, especially after the local government began to actively manage the mining.
Because the quarries in the Beijing region are actually exploited by people from Hebei Province, the study area was separated into two regions, one under the governance of Hebei Province and the

Uncertainties and Errors
For the classification step, we used a supervised learning method, which required a set of accurate and enriched samples.As we classified scenes of each year along the past 34 years, sampling was required for each scene and the accuracy of the results depends on the quality of the samples.Deriving operative samples from temporally adjacent samples and classifiers was demonstrated to work well in this study.However, further studies on other datasets or monitoring tasks are needed to confirm the conclusion.One more uncertainty is that when discussing the influence of social events and local policies on quarry changes, our analysis is based on the comparison of quarry areas in two jurisdictions considering some well-known events and policies.Maybe integrating some kind of quantization of the social events and policies into the analysis can present their close connection more clearly.

Social Discussion
Prior work has documented the feasibility of the Landsat archive for monitoring land cover changes and analyzing their dynamics in a high-resolution multi-temporal way [14].However, studies on quarry monitoring remain temporally limited.Whether using several decades of spanned space-borne imagery or higher resolution images focusing on a short time scale, prior studies have generally relied on a comparison of imagery from two to three time slices [1,5,6,9,10].Hence, with the focus on scarce observations so far, the socioeconomic reasons behind these quarry changes have not been thoroughly analyzed.In this study, we used a long time series data from the Landsat archive to fully investigate the annual change in the quarry area.In this section, as the annual monitored result of the quarry has been derived, we try to discuss the factors driving the quarry development from a socioeconomic viewpoint and demonstrate the potential benefits of annual monitoring in a socioeconomic study.

Comparison of Two Regions
The study area is located on the border of Beijing, Tianjin, and Hebei, making a comparison of quarry development among these areas possible.Since the development of quarry extraction is influenced deeply by human activities, it is also reasonable to discuss the quarry area extension and contraction with the different policies in the three regions, especially after the local government began to actively manage the mining.
Because the quarries in the Beijing region are actually exploited by people from Hebei Province, the study area was separated into two regions, one under the governance of Hebei Province and the other under the governance of Tianjin City, according to the boundary line between Hebei and Tianjin (Figure 9).Then, we separated the dynamics of the quarry area into two zones.
other under the governance of Tianjin City, according to the boundary line between Hebei and Tianjin (Figure 9).Then, we separated the dynamics of the quarry area into two zones.The amount of quarry area in Tianjin is quite small in the study area compared to that in Hebei Province, but is actually the biggest quarry site in Tianjin and most aggregates used in Tianjin buildings come from here.As it is so near to Hebei, its dynamics has a close relationship with that of the quarries in Hebei.Similarly, line charts are presented for a better observation of its dynamics (Figure 10).In order to compare them more easily, they are scaled in Figure 10b Behind these different tendencies, there are multiple socioeconomic reasons.In the early time, the government had no objection to the exploitation because of the great contribution of quarries to local revenue.Quarry activities were dominated by local people who were driven by the profits.Hence, the economy was the main factor affecting the quarries' expansion in the two regions.However, from 2008, the government started managing quarry activities.One big reason was The amount of quarry area in Tianjin is quite small in the study area compared to that in Hebei Province, but is actually the biggest quarry site in Tianjin and most aggregates used in Tianjin buildings come from here.As it is so near to Hebei, its dynamics has a close relationship with that of the quarries in Hebei.Similarly, line charts are presented for a better observation of its dynamics (Figure 10).In order to compare them more easily, they are scaled in Figure 10b, which indicates that their change trajectories before 2008 were very similar.From 2008 to 2016, the change trends of the quarries in Tianjin and Hebei developed in opposite directions.The quarry area in Tianjin decreased and then increased during this period, whereas that in Hebei grew first and then decreased.They both experienced a sharp increment in 2017.other under the governance of Tianjin City, according to the boundary line between Hebei and Tianjin (Figure 9).Then, we separated the dynamics of the quarry area into two zones.
. The amount of quarry area in Tianjin is quite small in the study area compared to that in Hebei Province, but is actually the biggest quarry site in Tianjin and most aggregates used in Tianjin buildings come from here.As it is so near to Hebei, its dynamics has a close relationship with that of the quarries in Hebei.Similarly, line charts are presented for a better observation of its dynamics (Figure 10).In order to compare them more easily, they are scaled in Figure 10b Behind these different tendencies, there are multiple socioeconomic reasons.In the early time, the government had no objection to the exploitation because of the great contribution of quarries to local revenue.Quarry activities were dominated by local people who were driven by the profits.Hence, the economy was the main factor affecting the quarries' expansion in the two regions.However, from 2008, the government started managing quarry activities.One big reason was Behind these different tendencies, there are multiple socioeconomic reasons.In the early time, the government had no objection to the exploitation because of the great contribution of quarries to local revenue.Quarry activities were dominated by local people who were driven by the profits.Hence, the economy was the main factor affecting the quarries' expansion in the two regions.However, from 2008, the government started managing quarry activities.One big reason was Beijing's Olympic Games in 2008, which caused a forceful restraint on quarry activities, arresting the extension of the quarry area from 2007 to 2008 in total.Specifically, the Tianjin government took stricter actions on regulating its quarry since January 2008 because a serious landslide inside a quarry site of our study area in Ji County, Tianjin, resulted in the death of 12 people (http://www.gov.cn/jrzg/2008-01/14/content_857398.htm).While these actions brought about an acute recovery of the quarry area, most of the truck drivers living on the quarry sites in Tianjin had to move to nearby sites in Hebei to continue their business, thereby causing an increase in the quarry area in Hebei despite the big inhibitory factor-Beijing Olympic Games in 2008.Overall, the quarry area in Tianjin decreased and the quarry area in Hebei increased while the total quarry area stayed nearly stable.After 2011, the economic transition in China required an environmentally friendly way in which to develop the economy.Inevitably, the total area of quarries decreased before an environmentally friendly developing pattern in the mining industry was found.Additionally, around 2013, CCTV-2 Half-Hour Economy reported that many local people contracted serious respiratory diseases because of the heavy dust in the air around the quarry sites in Hebei, and pointed out that this quarry region was the largest source of air pollutants on the east side of Beijing (http://finance.sina.com.cn/china/dfjj/20130621/220515877577.shtml).After that, all these quarry sites were forced by the government to close.Therefore, this time, the miners and truck drivers moved to sites in Tianjin, causing the quarry area in Tianjin to grow again.Therefore, we can see a different tendency in quarry changes in these two regions again and an overall decreasing trend in the study area.

Relationship between Quarry Area and GRP
It seems natural to relate quarry area increase and economic growth because the construction industry is one of the most profitable trades in China after reform and opening up.After we have acquired the annual monitoring results of the quarry, it also becomes possible to demonstrate this potential positive relationship.In this section, we discuss the relationship between the gross regional products (GRP) and quarry area in Langfang City using linear regression.
We first obtained the GRP in Langfang City from the Statistics Bureau of Hebei Province.The data and quarry area are depicted in Figure 11a where the line charts of GRP and quarry area express very similar developing trajectories that deviated after 2009.Moreover, from the analysis in Section 5.1, we know that in 2008 the government began to restrain quarry activities.We think this restraint shows its effects in the economy from 2009.Therefore, regression analyses between quarry area and GRP were conducted on data subsets of from 1990 to 2009 and from 2009 to 2016.The coefficient of determination between quarry area and GRP from 1990 to 2009 is 0.96531 (Figure 11b), evincing that the quarry area has good correlation with GRP in Langfang.The coefficient of determination for data after 2009 is 0.30004 (Figure 11c), which is much lower than that before, demonstrating that the correlation between quarry area and Langfang economy decreased after 2009.In summary, the above analysis indicates a very strong positive relationship between quarry area and GRP in Langfang City before 2009.
Beijing's Olympic Games in 2008, which caused a forceful restraint on quarry activities, arresting the extension of the quarry area from 2007 to 2008 in total.Specifically, the Tianjin government took stricter actions on regulating its quarry since January 2008 because a serious landslide inside a quarry site of our study area in Ji County, Tianjin, resulted in the death of 12 people (http://www.gov.cn/jrzg/2008-01/14/content_857398.htm).While these actions brought about an acute recovery of the quarry area, most of the truck drivers living on the quarry sites in Tianjin had to move to nearby sites in Hebei to continue their business, thereby causing an increase in the quarry area in Hebei despite the big inhibitory factor-Beijing Olympic Games in 2008.Overall, the quarry area in Tianjin decreased and the quarry area in Hebei increased while the total quarry area stayed nearly stable.After 2011, the economic transition in China required an environmentally friendly way in which to develop the economy.Inevitably, the total area of quarries decreased before an environmentally friendly developing pattern in the mining industry was found.Additionally, around 2013, CCTV-2 Half-Hour Economy reported that many local people contracted serious respiratory diseases because of the heavy dust in the air around the quarry sites in Hebei, and pointed out that this quarry region was the largest source of air pollutants on the east side of Beijing (http://finance.sina.com.cn/china/dfjj/20130621/220515877577.shtml).After that, all these quarry sites were forced by the government to close.Therefore, this time, the miners and truck drivers moved to sites in Tianjin, causing the quarry area in Tianjin to grow again.Therefore, we can see a different tendency in quarry changes in these two regions again and an overall decreasing trend in the study area.

Relationship between Quarry Area and GRP
It seems natural to relate quarry area increase and economic growth because the construction industry is one of the most profitable trades in China after reform and opening up.After we have acquired the annual monitoring results of the quarry, it also becomes possible to demonstrate this potential positive relationship.In this section, we discuss the relationship between the gross regional products (GRP) and quarry area in Langfang City using linear regression.
We first obtained the GRP in Langfang City from the Statistics Bureau of Hebei Province.The data and quarry area are depicted in Figure 11a where the line charts of GRP and quarry area express very similar developing trajectories that deviated after 2009.Moreover, from the analysis in Section 5.1, we know that in 2008 the government began to restrain quarry activities.We think this restraint shows its effects in the economy from 2009.Therefore, regression analyses between quarry area and GRP were conducted on data subsets of from 1990 to 2009 and from 2009 to 2016.The coefficient of determination between quarry area and GRP from 1990 to 2009 is 0.96531 (Figure 11b), evincing that the quarry area has good correlation with GRP in Langfang.The coefficient of determination for data after 2009 is 0.30004 (Figure 11c), which is much lower than that before, demonstrating that the correlation between quarry area and Langfang economy decreased after 2009.In summary, the above analysis indicates a very strong positive relationship between quarry area and GRP in Langfang City before 2009.

Conclusions
We used long time-series data from the Landsat archive to investigate the annual change in the quarry area.A 34-year annual monitoring result from 1984 to 2017 in a quarry region near Beijing was derived in this study using an efficient approach based on a modified random forest classification and quarry region analysis.A temporal smoothing strategy was proposed and applied to the multi-temporal result of classification for post-processing.The result exhibits a large growth in quarry area in the study region, namely from 4.2574 km 2 to 16.4969 km 2 , which is more than onetenth of the total area of the study region.The quarry area time series revealed a five-stage development with different growth patterns: (1) 1984-2002: the quarry area continually grew at a low speed, mainly below 0.5 km 2 per year (average 0.32 km 2 per year), with its area expanding from 4.2574 km 2 to 10.0584 km 2 ; (2) 2002-2010: the quarry area grew from 10.0584 km 2 to 15.9346 km 2 with an average annual growth rate higher than 0.7 km 2 per year; (3) 2010-2013: this period saw a modest increase in the quarry area with a total increment of 0.4451 km 2 during the 3 years; (4) 2013-2016: the quarry area declined for the first time from 16.3797 km 2 to 15.9658 km 2 ; (5) 2016-2017: a large growth in the quarry area recurred with the area in 2017 reaching 16.5002 km 2 .This study is the first to our knowledge to produce an annual quarry limits in one region for the past 34 years, which is more valuable than analyses on several observations as it can be used to reveal a potentially strong relationship between quarry changes and socioeconomic dynamics.For example, we did a comparison of quarry area dynamics between two nearby regions, from which their reactions to influential social events can be identified.With this time series of quarry area, we also proved a strong positive relationship between quarry area and local economic development.Moreover, as many ancient temples inside the study area were destroyed during the quarry activity, the time series of quarry area are of great benefit for the conservation of historic heritage.A further study may focus on a combined analysis of quarry area and other socioeconomic parameters, such population, highway mileage, and cement production.

Figure 1 .
Figure 1.Study area located at the boundary of Beijing, Tianjin, and Hebei; the upper left image is the Landsat imagery in 2016 depicting severe damage to the mountain body caused by mining activities; the lower left image is a digital surface model (DSM) map of the study area.Data comes from Advanced Land Observing Satellite (ALOS) DSM.(http://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm).

Figure 1 .
Figure 1.Study area located at the boundary of Beijing, Tianjin, and Hebei; the upper left image is the Landsat imagery in 2016 depicting severe damage to the mountain body caused by mining activities; the lower left image is a digital surface model (DSM) map of the study area.Data comes from Advanced Land Observing Satellite (ALOS) DSM.(http://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm).

FiltrationFigure 2 .
Figure 2. Representation of the methodological framework used in this study.

Figure 2 .
Figure 2. Representation of the methodological framework used in this study.

Figure 3 .
Figure 3. (a) DSM image of the study area; (b) Altitude mask of the study area.

Figure 4 .
Figure 4.An example of temporal smoothing on the confidence index sequence of one pixel.

Figure 3 .
Figure 3. (a) DSM image of the study area; (b) Altitude mask of the study area.

Figure 3 .
Figure 3. (a) DSM image of the study area; (b) Altitude mask of the study area.

Figure 4 .
Figure 4.An example of temporal smoothing on the confidence index sequence of one pixel.

Figure 4 .
Figure 4.An example of temporal smoothing on the confidence index sequence of one pixel.

Figure 5 .
Figure 5. Process illustration using 2002 as an example.(a) Edges of quarry; (b) Three defined zones; (c) Three types of zones in detail: zone (A) is the internal quarry region, which is fully accounted for in the area calculation; zone (B) is the margin of the quarry zone, which is partly accounted for in the

Figure 5 .
Figure 5. Process illustration using 2002 as an example.(a) Edges of quarry; (b) Three defined zones; (c) Three types of zones in detail: zone (A) is the internal quarry region, which is fully accounted for in the area calculation; zone (B) is the margin of the quarry zone, which is partly accounted for in the area calculation taking the confidence index as a weight coefficient; zone (C) is the region classified as non-quarry.
The right line chart represents its change speed at different time points, where each point means the change amount from the last year to this year.This chart insinuates different change rates at different time periods.According to its change patterns presented in the graph, five different change periods are distinguished: (1) The first period was from 1984 to 2002.In this period, the change speed kept at a low level, mainly under 0.5 km 2 per year, except for the speed in 1993.The speed changed slowly, taking roughly 4-5 years to rise and fall.(2) The next period with a high change speed continued from 2002 to 2010.Its speed was greater than 0.6 km 2 per year during this period, except in 2008.(3) The area then underwent a low-speed increase over the following 3 years from 2011 to 2013.(4) Its change rate dropped to negative values in the following 3 years (2013-2016), which means the quarry area went through a decline in these years.(5) Finally, from May 2016 to June 2017 sees a great areal growth in the quarry.

Figure 6 .
Figure 6.(a) Quarry area changes from 1984 to 2017; (b) the change rate of quarry area in each year.

Figure 6 .
Figure 6.(a) Quarry area changes from 1984 to 2017; (b) the change rate of quarry area in each year.

Figure 9 .
Figure 9.The study area is split into two regions according to the boundary between Hebei and Tianjin.

Figure 10 .
Figure 10.(a) Quarry area in two regions from 1984 to 2017; (b) Quarry area depicted with two y-axes for the two regions.

Figure 9 .
Figure 9.The study area is split into two regions according to the boundary between Hebei and Tianjin.

Figure 9 .
Figure 9.The study area is split into two regions according to the boundary between Hebei and Tianjin.

Figure 10 .
Figure 10.(a) Quarry area in two regions from 1984 to 2017; (b) Quarry area depicted with two y-axes for the two regions.

Figure 10 .
Figure 10.(a) Quarry area in two regions from 1984 to 2017; (b) Quarry area depicted with two y-axes for the two regions.

Figure 11 .
Figure 11.(a) Line charts of gross regional products (GRP) in Langfang and the quarry area in Langfang; (b) The scatter diagram of the regression variables before 2009; (c) The scatter diagram of the regression variables after 2009.

Table 1 .
Classification scheme with the number of training samples (number of sample pixels and sample areas in six selected years).

Table 2 .
Accuracy of the classification results on the ground truth data.