Spatial Patterns and Driving Forces of Greenhouse Land Change in Shouguang City , China

In the feature extraction step, preprocessed images were subject to a Kauth-Thomas transformation (also called a ‘tasselled cap transformation’), computation of normalized difference vegetation index (NDVI), a Karhunen-Loeve transformation, and a minimum noise fraction (MNF) rotation. Kauth-Thomas transformation is a linear transformation that rotates coordinate space. After rotation, coordinate axes point in the direction most closely related to ground features. Thus, Kauth-Thomas transformation is able to identify the features of vegetation and soil in multispectral space. Three sub-images representing brightness, greenness, and moisture were separately generated after transformation. The NDVI denotes reflectivity of vegetation in near-infrared and red bands. This index is the optimal indicator for vegetation growth status and coverage, as it can eliminate the impact of altitude angle, satellite observation angle, terrain, and radiation change induced by cloud shadows or other atmospheric conditions. The formula used to calculate NDVI is as follows:


Introduction
Agriculture is essential to sustain population levels as well as for the development of material (e.g., industry) and non-material production sectors [1,2].Although countries around the world vary in their speed of agricultural development, all pursue modern approaches because they enable both high productivity and economic benefits [3,4].Greenhouse use has emerged in response to these developments, as an important tool in modern agriculture [5] that can enhance both the quality and output of products by using artificial facilities to alter horticultural environments [6].Because greenhouses are effective ways to utilize water and land resources, as well as enhance output rates [7], it is not surprising that their usage has comprised an increasing proportion of global agriculture [8].
Greenhouse area has been expanding around the world [9].Over the last two decades, greenhouse area worldwide has increased five-fold and this expansion continues to accelerate [10].Usage globally is mainly concentrated in Asia, with China accounting for almost half of the overall worldwide total [11].To meet accelerating demands for vegetables, the greenhouse area in China doubled between 2006 and 2010 [12].Europe boasts the second largest greenhouse area globally, encompassing over 20% of the total [13]; in this region, greenhouses are mostly concentrated in Italy, Spain, and the Netherlands.In contrast, the greenhouse area in North and South America accounts for just a small proportion of the global total; and though greenhouse production is an important component of the U.S. horticultural sector, total greenhouse area is just 9100 hectares, of which 8000 hectares is used for ornamental production [14].
Most countries around the world do not collect regular statistics on greenhouse area.Indeed, of the few that do, just a handful collates yearly figures, including the Netherlands, Germany, and Canada.Thus, detecting greenhouse land cover using RS data has become a research target; very high resolution satellite images have been utilized in previous studies to identify greenhouse land cover.For example, Sönmez and Sari [15] manually identified the boundaries of greenhouses in the Antalya region of Turkey using IKONOS satellite images, achieving a classification accuracy of 96%, while Aguera et al. [16] used QuickBird satellite images to detect greenhouses in southeastern Spain.In earlier work, Aguera et al. [17] and Aguera and Liu [18] developed an automatic greenhouse extraction method based on IKONOS and Quickbird satellite images, achieving classification accuracies of 90% and 84.20%, respectively.Using the same satellite images as the source of data, Carjaval et al. [19] then applied an artificial neural network classifier to identify greenhouse land cover in southeastern Spain.
In addition to IKONOS and Quickbird images, researchers have also investigated the use of other satellite data sources to map greenhouse land cover.Tarantino and Figorito [20], for example, showed that traditional classification techniques face challenges when used to identify greenhouse land cover in complex areas.These workers therefore applied an object classifier to detect greenhouse land cover using aerial images, achieving an accuracy of 90.25%.Building on this, Koc-San and Sonmez [21] used WorldView-2 satellite images to map greenhouse land cover in Antalya, Turkey.These images are quite promising for the detection of greenhouse land cover.
However, while it is feasible to detect greenhouse land cover using very high resolution satellite images, mapping issues remain.In the first place, high resolution images are generally not free to acquire, and often incur high cost.Thus, for a country such as China that has a huge area, it has been very difficult to monitor large-scale greenhouse land cover using high-resolution satellite images.Secondly, time series encapsulated by high-resolution satellite images are usually short, limiting long-time detection of greenhouse land cover changes.Finally, most of the previous research has focused on determining the current status of greenhouse land cover, and has ignored analysis of underlying drivers.Such an analysis will contribute to a deeper understanding of the expansion of greenhouse land cover and provide a scientific basis for future forecasting and land use policy decisions.
The aim of this paper is to map and analyze changes in greenhouse land cover within the city of Shouguang, China, using free-to-download Landsat images that encapsulate a continuous time series of over 40 years.We also attempt to discern the factors driving greenhouse land cover change in this region by applying a logistic regression model.

Study Area
The city of Shouguang (118 • 32 E-119 • 10 E, 36 • 41 N-37 • 19 N) is situated on the coastal plain of north-central Shandong Province and covers an area of 2200 km 2 (Figure 1).This city has a warm-to-temperate continental monsoonal climate, characterized by a dry spring, a hot and rainy summer, a cool autumn, and a dry and cold winter with little snow.The altitude of the city varies between 1 m and 110 m, decreasing from south to north.Shouguang is the cradle of Chinese winter greenhouses and enjoys the reputation as the "home of vegetables".This city is the largest production base and distribution center of commercial vegetables in China.

Figure 1.
Location and classification verification for Shouguang, using ground reference data from Google Earth.We selected 50 samples classified as greenhouse land in 2015 and verified them using Google Earth images.Of these samples, the six highlighted in red were incorrectly classified as greenhouse land, while of the ten enlarged boxes, 1 and 8 actually comprise arable land.The other six comprise greenhouse land.

Data Sources
We used RS images to extract greenhouse land cover in Shouguang during 2000, 2010, and 2015.All these images consist of ETM data downloaded from the Landsat 8 data sharing system [22].The properties of these images are listed in Table S1.As one ETM image scene covers an area of 185 km × 185 km, a total of three images were downloaded in order to detect greenhouse land cover in Shouguang.The meteorological data used in this paper were sourced from the website of the China Meteorological Administration [23], while digital elevation model (DEM) data at a spatial resolution of 30 m were extracted from the Geospatial Data Cloud [24].Road and settlement data were collected from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [25].

Greenhouse Land Cover Mapping
Five steps were followed in order to map greenhouse land cover in Shouguang: image pre-processing, feature extraction, supervised classification, post-processing, and thematic mapping (Figure 2).Of these, the aim of the first step, image pre-processing, is to improve image data while suppressing distortion [26]; however, as ETM images have already been geometrically corrected and radiometrically calibrated before download, we only performed image mosaic and masking processes, as well as de-striping, cloud removal, and atmospheric correction using the software packages ENVI 5.1 (Exelis Visual Information Solutions: Broomfield, CO, USA, 2013) and Arcgis Location and classification verification for Shouguang, using ground reference data from Google Earth.We selected 50 samples classified as greenhouse land in 2015 and verified them using Google Earth images.Of these samples, the six highlighted in red were incorrectly classified as greenhouse land, while of the ten enlarged boxes, 1 and 8 actually comprise arable land.The other six comprise greenhouse land.

Data Sources
We used RS images to extract greenhouse land cover in Shouguang during 2000, 2010, and 2015.All these images consist of ETM data downloaded from the Landsat 8 data sharing system [22].The properties of these images are listed in Table S1.As one ETM image scene covers an area of 185 km × 185 km, a total of three images were downloaded in order to detect greenhouse land cover in Shouguang.The meteorological data used in this paper were sourced from the website of the China Meteorological Administration [23], while digital elevation model (DEM) data at a spatial resolution of 30 m were extracted from the Geospatial Data Cloud [24].Road and settlement data were collected from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [25].

Greenhouse Land Cover Mapping
Five steps were followed in order to map greenhouse land cover in Shouguang: image pre-processing, feature extraction, supervised classification, post-processing, and thematic mapping (Figure 2).Of these, the aim of the first step, image pre-processing, is to improve image data while suppressing distortion [26]; however, as ETM images have already been geometrically corrected and radiometrically calibrated before download, we only performed image mosaic and masking processes, as well as de-striping, cloud removal, and atmospheric correction using the software packages ENVI 5.1 (Exelis Visual Information Solutions: Broomfield, CO, USA, 2013) and Arcgis 10.3 (Esri: New York, NY, USA, 2014).As a result of this image pre-processing, clearer and more accurate images of the Shouguang administrative region were obtained (see Figure S1).
10.3 (Esri: New York, NY, USA, 2014).As a result of this image pre-processing, clearer and more accurate images of the Shouguang administrative region were obtained (see Figure S1).The second step, feature extraction, aims to generate new fused images that can be used to increase classification accuracy at later stages.This step includes three sub-steps: spectral feature extraction, textural feature extraction, and image fusion (Figure 2).Of these, spectral features are the unique spectral reflection and radiation properties of different ground features, the result of differences in material compositions and structure.Thus, using spectral features, it is possible to identify greenhouse land cover in images, although noise that can influence image interpretation is usually present in original images.The aim of spectral feature extraction is therefore to extract useful spectral information via analysis or calculation.To do this, four methods (i.e., calculation of normalized difference vegetation index (NDVI); extraction of first and second spectral components; brightness, greenness, and wetness of images; and principal component analysis (PCA)) were applied using the software ENVI5.1 (see Supplementary Materials for full methodological details).The next sub-step involves the extraction of textural features, a set of metrics computed in image processing designed to quantitatively describe the perceived texture of an image.This provides information on the spatial arrangement of color or intensities in an image [27][28][29]; extraction was performed using a PCA analysis, performed during the spectral feature extraction sub-step.Thus, textural features of images were analyzed in this study using eight common indexes [30,31] encompassing mean, variance, homogeneity, contrast, dissimilarity, entropy, angular second moment, and correlation.Following these two sub-steps, spectral and textural features were then fused to create new satellite images which were used as the basis for greenhouse land cover classification.
In the third processing step, we used a SVM for image classification.In this context, a SVM is a machine-learning algorithm based on statistical learning theory, Vapnik-Chervonenkis dimensions, and structural risk minimization [32].This approach enables the management of high-resolution, multiple-band images, as well as large-sized segmented satellite data [33].Indeed, because SVMs have favorable generalization capacities, they are often used to solve small-sample nonlinear and high-dimensional pattern identification problems, as well as for prediction and comprehensive evaluation [34].SVMs work by generating a hyperplane that represents an optimal separation of The second step, feature extraction, aims to generate new fused images that can be used to increase classification accuracy at later stages.This step includes three sub-steps: spectral feature extraction, textural feature extraction, and image fusion (Figure 2).Of these, spectral features are the unique spectral reflection and radiation properties of different ground features, the result of differences in material compositions and structure.Thus, using spectral features, it is possible to identify greenhouse land cover in images, although noise that can influence image interpretation is usually present in original images.The aim of spectral feature extraction is therefore to extract useful spectral information via analysis or calculation.To do this, four methods (i.e., calculation of normalized difference vegetation index (NDVI); extraction of first and second spectral components; brightness, greenness, and wetness of images; and principal component analysis (PCA)) were applied using the software ENVI5.1 (see Supplementary Materials for full methodological details).The next sub-step involves the extraction of textural features, a set of metrics computed in image processing designed to quantitatively describe the perceived texture of an image.This provides information on the spatial arrangement of color or intensities in an image [27][28][29]; extraction was performed using a PCA analysis, performed during the spectral feature extraction sub-step.Thus, textural features of images were analyzed in this study using eight common indexes [30,31] encompassing mean, variance, homogeneity, contrast, dissimilarity, entropy, angular second moment, and correlation.Following these two sub-steps, spectral and textural features were then fused to create new satellite images which were used as the basis for greenhouse land cover classification.
In the third processing step, we used a SVM for image classification.In this context, a SVM is a machine-learning algorithm based on statistical learning theory, Vapnik-Chervonenkis dimensions, and structural risk minimization [32].This approach enables the management of high-resolution, multiple-band images, as well as large-sized segmented satellite data [33].Indeed, because SVMs have favorable generalization capacities, they are often used to solve small-sample nonlinear and high-dimensional pattern identification problems, as well as for prediction and comprehensive evaluation [34].SVMs work by generating a hyperplane that represents an optimal separation of linearly separable classes in a decision boundary space [35] (see Supplementary Material for a full mathematical description).In brief, a SVM forecasts the destination values of test data by developing a model using training data.The basic process of running a SVM includes steps to select a training data set, compute separability, training using a SVM algorithm, and the acquisition of classification results.Thus, for this study, land use in Shouguang was divided into seven types encompassing greenhouses, other arable, forested, grassland, water area, construction, and unused.Training samples of these different land use types were then selected using pre-processed ETM images (see Supplementary Material for methods); separability was calculated to be 1.8, which indicates that ground features are easily distinguishable from one another.We therefore used a SVM to classify land use in Shouguang, generating three land use maps (Figure 3).Training samples of these different land use types were then selected using pre-processed ETM images (see Supplementary Material for methods); separability was calculated to be 1.8, which indicates that ground features are easily distinguishable from one another.We therefore used a SVM to classify land use in Shouguang, generating three land use maps (Figure 3).As part of post-processing, it is important to assess accuracy.To do this, the most popular technique applies an error, or confusion, matrix [33,[36][37][38].Based on this matrix, we calculated both overall accuracies and Kappa Statistics (see Supplementary Material for full method) using validation data [36,[38][39].A total of 100 samples that differed from training classification samples were selected to assess classification accuracy.The results of this step show that classifications of Shouguang land use types for 2000, 2010, and 2015 all attained a satisfactory level of accuracy: 98.17%, 98.54%, and 96.60%, with Kappa coefficients of 0.9738, 0.9813, and 0.9548, respectively.In addition, as ground referencing is an essential component of any analysis that relies on RS data, we also visually interpreted Google Earth data at a resolution of 0.5 m to further verify our classifications.However, as Google Earth data for 2000 and 2010 were unavailable, we were only able to apply this method to verify classification accuracy for 2015 (see Figure 1).A total of 50 samples for 2015 were selected, of which 44 were found to be correctly classified as greenhouse land.This translates to a ground referenced classification accuracy of 88%.

Driving Force Analysis of Greenhouse Land Cover Change
To explore spatial patterns in greenhouse land cover change in the city of Shouguang, as well as the possible reasons for these changes, we investigated 13 potential "driving forces", including the distance to the nearest coastline (coastline), the distance to the nearest highway (highway), the distance to the nearest railway (railway), the distance to the nearest provincial road (provroad), the distance to the nearest ordinary road (road), the distance to the nearest urban settlement (urban), the distance to the nearest rural settlement (rural), and the distance to the nearest river (river), as well as elevation, slope, aspect, soil water content (swc), and precipitation (∆precip, encapsulating As part of post-processing, it is important to assess accuracy.To do this, the most popular technique applies an error, or confusion, matrix [33,[36][37][38].Based on this matrix, we calculated both overall accuracies and Kappa Statistics (see Supplementary Material for full method) using validation data [36,38,39].A total of 100 samples that differed from training classification samples were selected to assess classification accuracy.The results of this step show that classifications of Shouguang land use types for 2000, 2010, and 2015 all attained a satisfactory level of accuracy: 98.17%, 98.54%, and 96.60%, with Kappa coefficients of 0.9738, 0.9813, and 0.9548, respectively.In addition, as ground referencing is an essential component of any analysis that relies on RS data, we also visually interpreted Google Earth data at a resolution of 0.5 m to further verify our classifications.However, as Google Earth data for 2000 and 2010 were unavailable, we were only able to apply this method to verify classification accuracy for 2015 (see Figure 1).A total of 50 samples for 2015 were selected, of which 44 were found to be correctly classified as greenhouse land.This translates to a ground referenced classification accuracy of 88%.

Driving Force Analysis of Greenhouse Land Cover Change
To explore spatial patterns in greenhouse land cover change in the city of Shouguang, as well as the possible reasons for these changes, we investigated 13 potential "driving forces", including the distance to the nearest coastline (coastline), the distance to the nearest highway (highway), the distance to the nearest railway (railway), the distance to the nearest provincial road (provroad), the distance to the nearest ordinary road (road), the distance to the nearest urban settlement (urban), the distance to the nearest rural settlement (rural), and the distance to the nearest river (river), as well as elevation, slope, aspect, soil water content (swc), and precipitation (∆precip, encapsulating the change between 2000 and 2015).As precipitation data for 2015 were unavailable, data for 2012 were used in this study.
We selected driving forces for analysis by relying on expert knowledge and previous research.In the first place, we assumed that changes in greenhouse land cover are closely related to topography [40][41][42], as the key factors of elevation, slope, and aspect influence the complexity and cost of construction.In addition, climatic factors, such as precipitation, influence the growth of vegetables and are thus important in decisions about greenhouse construction [9,[43][44][45][46], while distance factors, such as the distance to the nearest road or railway, will influence transportation convenience and thus the final cost of vegetables wholesale [47][48][49].Additional distance factors, such as the distance to the nearest settlement or river, can also influence the management of vegetables planted in greenhouse facilities [50][51][52][53].
We analyzed the relationship between spatial patterns of greenhouse land cover and our driving forces using binary logistic regression [54].Although all driving force data were transformed to a grid at 30 m resolution for this analysis, it should be noted that if all pixels are selected as regression samples, serious spatial autocorrelation will occur and no meaningful results will be acquired [40,41,55,56].Although the effects of spatial autocorrelation on the results of regression analyses can be minimized by using just a sub-sample of pixels [56,57], the appropriate proportion has varied in previous studies between 15% and 40% [40,55,58].In this study, we used an approximate median sample proportion (30%) based on previous research, and analyzed variation in driving forces from the contrasting perspectives of land expansion and contraction.
In a binary logistic regression, independent variables are used as forecasted values to compute the occurrence probability of an event.Thus, this probability can help to quantitatively discern the relationship between land use distribution and various influencing factors [59].In this study, the dependent variable (Y) denotes greenhouse land expansion or contraction; this binary value, 1 or 0, denotes the presence or absence of expansion and contraction [60,61].The probability (P i ) describes the extent to which Y will change to 1 (see Equation (1)); thus, P i is very close to either 0 or 1, as follows: In this expression, P i (Y = 1|X 1 , X 2 , . . ., X n ) denotes the probability of Y given by X i (where i is the driving force variable) and thus the change from non-greenhouse to greenhouse land cover (expansion) and vice versa (contraction).At the same time, β i is the regression coefficient for each driving force variable, while 1-P i is the probability of no greenhouse land expansion or contraction.Thus, subsequent to Logit transformation, the equation for the binary logistic regression model is as follows: Independent coefficients can be assessed in this expression using maximum likelihood estimation; thus, analyzing the significance of this regression as well as the value for each coefficient, it is possible to interpret how each driving force affects greenhouse land expansion or contraction.
Relative operating characteristics (ROC) analysis can then be used to further assess the explanatory power of logistic regression models [41,62].In this context, ROC values indicate that spatial patterns of most land use types can be reasonably explained by the independent variables [42].ROC analysis is based on a rate curve of true versus false positives for a range of threshold values that are then used to classify probabilities into two classes [55].In general, the accuracy of a model can be considered to be either low, credible, or high when ROC values are lower than 0.7, between 0.7 and 0.9, and higher than 0.9, respectively [41,55,63,64].ROC values for the expansion and contraction models reported in this study are 0.786 and 0.844 (Figure 4), respectively, suggesting that our analysis of driving forces is credible.
values that are then used to classify probabilities into two classes [55].In general, the accuracy of a model can be considered to be either low, credible, or high when ROC values are lower than 0.7, between 0.7 and 0.9, and higher than 0.9, respectively [41,55,63,64].ROC values for the expansion and contraction models reported in this study are 0.786 and 0.844 (Figure 4), respectively, suggesting that our analysis of driving forces is credible.This ROC curve was used to test the degree of fit of each binary logistic regression; the area ratio under the curve corresponds to a ROC value between 0 and 1.Thus, larger is the value, the higher is the explanatory power of the logistic regression model.

Analysis of the Distribution of Greenhouse Land Use around Urban Areas
To determine the spatial relationship between greenhouse land and urban areas, we created 13 buffers around each urban settlement, at 2 km intervals.Thus, by analyzing the distribution of greenhouse land in each buffer, we investigated associations between this land use type and urban settlements.

Land Use Changes
The results of this study show that in 2000 the main land use types within the city of Shouguang were arable (32.38%), construction (17.26%), water areas (21.60%), and greenhouse land (22.60%) (Figure 3; Table 1).The remaining three land use types accounted for just 6.16% of the total area.Between 2000 and 2015, changes in land use within Shouguang included the rapid expansion of greenhouse land (50.51%), forests (46.46%), and construction land (46.78%), as well as drastic contractions in areas of other arable land (-49.27%) and grassland (-76.98%).

Changes in Greenhouse Land Use
Results show that, although greenhouse land use expanded continuously between 2000 and 2015, the rate of expansion was markedly higher in the later period (i.e., 8.71%/year between 2010 and 2015) compared to earlier (i.e., 0.49%/year between 2000 and 2010) (Table 2).During the earlier period, a total of 223.25 km 2 of land was converted to greenhouse use, while 198.28 km 2 of this This ROC curve was used to test the degree of fit of each binary logistic regression; the area ratio under the curve corresponds to a ROC value between 0 and 1.Thus, larger is the value, the higher is the explanatory power of the logistic regression model.

Analysis of the Distribution of Greenhouse Land Use around Urban Areas
To determine the spatial relationship between greenhouse land and urban areas, we created 13 buffers around each urban settlement, at 2 km intervals.Thus, by analyzing the distribution of greenhouse land in each buffer, we investigated associations between this land use type and urban settlements.

Land Use Changes
The results of this study show that in 2000 the main land use types within the city of Shouguang were arable (32.38%), construction (17.26%), water areas (21.60%), and greenhouse land (22.60%) (Figure 3; Table 1).The remaining three land use types accounted for just 6.16% of the total area.Between 2000 and 2015, changes in land use within Shouguang included the rapid expansion of greenhouse land (50.51%), forests (46.46%), and construction land (46.78%), as well as drastic contractions in areas of other arable land (-49.27%) and grassland (-76.98%).

Changes in Greenhouse Land Use
Results show that, although greenhouse land use expanded continuously between 2000 and 2015, the rate of expansion was markedly higher in the later period (i.e., 8.71%/year between 2010 and 2015) compared to earlier (i.e., 0.49%/year between 2000 and 2010) (Table 2).During the earlier period, a total of 223.25 km 2 of land was converted to greenhouse use, while 198.28 km 2 of this existing land use type was converted to other uses.This inter-conversion resulted in a net increase in greenhouse land of just 4.85 km 2 , while over the second five-year period, this reached 287.79 km 2 , not much higher than during the earlier period (223.25 km 2 ).However, the net decrease in land used for greenhouses fell as low as 52.75 km 2 during this period, resulting in a more rapid overall increase in this land use type.The inter-conversion of land use types during the earlier period mainly comprised other arable land and greenhouse land, while this transition was mainly between greenhouse and construction land during the later period.Results show that between 2000 and 2010, expansion in greenhouse land was mainly concentrated in the south of the city, with a modest increase in the northwest, while contraction mainly took place in central region (Figure 5a).However, between 2010 and 2015, expansion in greenhouse land mainly took place in the central and southern regions (Figure 5b), while the southern and northwestern sectors, as well as some sites in the east, experienced slight contraction.
existing land use type was converted to other uses.This inter-conversion resulted in a net increase in greenhouse land of just 4.85 km 2 , while over the second five-year period, this reached 287.79 km 2 , not much higher than during the earlier period (223.25 km 2 ).However, the net decrease in land used for greenhouses fell as low as 52.75 km 2 during this period, resulting in a more rapid overall increase in this land use type.The inter-conversion of land use types during the earlier period mainly comprised other arable land and greenhouse land, while this transition was mainly between greenhouse and construction land during the later period.Results show that between 2000 and 2010, expansion in greenhouse land was mainly concentrated in the south of the city, with a modest increase in the northwest, while contraction mainly took place in the central region (Figure 5a).However, between 2010 and 2015, expansion in greenhouse land mainly took place in the central and southern regions (Figure 5b), while the southern and northwestern sectors, as well as some sites in the east, experienced slight contraction.

The Relationship between Greenhouse Land Use and Urban Settlements
The results show that distance to the nearest urban settlement is generally negatively correlated with area of greenhouse land (Figure 6).The total land area used for greenhouses was A decrease denotes a change from greenhouse land to another land use type, while an increase denotes the opposite.

The Relationship between Greenhouse Land Use and Urban Settlements
The results show that distance to the nearest urban settlement is generally negatively correlated with area of greenhouse land (Figure 6).The total land area used for greenhouses was more than 100 km 2 within the three buffers closest to urban settlements (i.e., between 0 km and 2 km, 2 km and 4 km, and 4 km and 6 km).The area of this land use type also gradually declines with increasing distance from urban settlement (Figure 6); compared to 2000 and 2010, the area of greenhouse land in 2015 was increasingly concentrated within the 10 km buffer zone.
Sustainability 2017, 3, 359 9 of 15 more than 100 km 2 within the three buffers closest to urban settlements (i.e., between 0 km and 2 km, 2 km and 4 km, and 4 km and 6 km).The area of this land use type also gradually declines with increasing distance from urban settlement (Figure 6); compared to 2000 and 2010, the area of greenhouse land in 2015 was increasingly concentrated within the 10 km buffer zone.

Expansion-Inducing Forces
The results of our binary logistic regression (Table 3) show that five variables significantly influenced the expansion of greenhouse land use between 2000 and 2015; these were distances to the nearest coastline and rural settlement, as well as elevation, slope, and precipitation.Specifically, positive and negative values for regression coefficients show that distance to the nearest coastline and slope both had positive effects on expansion, while the other three driving forces exerted negative effects.The positive effects on expansion of distance to the nearest coastline and slope mean that greenhouse land use has preferentially expanded in areas far from the coastline zone, as well as in areas with relatively high slopes.In contrast, the negative effects of distance to the nearest rural settlement, elevation, and changes in precipitation between 2000 and 2012 have also caused greenhouse land use to expand preferentially in areas close to rural settlements and in regions characterized by low elevation and small changes in precipitation.The results of our binary logistic regression (Table 3) show that five variables significantly influenced the expansion of greenhouse land use between 2000 and 2015; these were distances to the nearest coastline and rural settlement, as well as elevation, slope, and precipitation.Specifically, positive and negative values for regression coefficients show that distance to the nearest coastline and slope both had positive effects on expansion, while the other three driving forces exerted negative effects.The positive effects on expansion of distance to the nearest coastline and slope mean that greenhouse land use has preferentially expanded in areas far from the coastline zone, as well as in areas with relatively high slopes.In contrast, the negative effects of distance to the nearest rural settlement, elevation, and changes in precipitation between 2000 and 2012 have also caused greenhouse land use to expand preferentially in areas close to rural settlements and in regions characterized by low elevation and small changes in precipitation.
Regions adjacent to the coastline zone are usually characterized by serious soil salinization and are thus unsuitable for greenhouse land use expansion.Greenhouses also require a large amount of labor; thus, areas close to rural settlements can be convenient for vegetable cultivation and for attracting construction workers.The low-elevation areas of Shouguang are relatively warm and thus are suitable for vegetable cultivation.These areas have witnessed expansions in greenhouse land use.At the same time, however, it is interesting to note that this land use type in Shouguang has been preferentially expanded in hilly areas rather than flat ones because of the scarcity of flat arable land.Shouguang is a densely-populated region with a long history of greenhouse land use; however, expansion of this land use type in flat areas is problematic, because many are already covered with greenhouses.In addition, because the Shouguang region is relatively flat overall, even sloping areas do not necessarily limit greenhouse construction.This land use type has also tended to expand in areas characterized by slight changes in precipitation, because stable climatic condition reduce the risk of vegetable cultivation disasters.

Contraction-Inducing Forces
Regression results (Table 4) show that distances to the nearest coastline, railway, urban settlement, river, and road all exert significant effects on the contraction of land used for greenhouses.Indeed, with the exception of distance to the nearest railway, regression coefficient values show that the other four factors all exert negative effects on greenhouse land contraction.The area of land used for greenhouses has preferentially contracted close to coastlines, urban areas, rural settlements, and roads, as well as in regions far from railways.Abbreviations as in Table 3.
Greenhouse land use has also preferentially contracted in regions adjacent to the coastline, as these areas are not suitable for vegetable cultivation because of their high soil salt content.In addition, land area used for greenhouses has also tended to contract in regions adjacent to urban areas and rural settlements, as land has a high probability of being occupied by these usage types.Because greenhouses also tend to generate soil and water pollution, this land use type has tended to contract in regions close to rivers, as well as far from railways because of transportation issues.

Discussion
In addition to the spatial driving factors discussed above, increases in greenhouse land use area in Shouguang are closely related to policy.Because this city is a major Chinese vegetable center, it is necessary for Shouguang to continuously enlarge its area of cultivation to satisfy domestic and foreign demands for produce.Thus, to achieve this goal, the local government has enacted a series of policies to encourage the construction of greenhouses.For example, because construction costs are relatively high compared to the low income of farmers in this region, the government provides subsidies for greenhouse-building and has reduced the fees related to construction.Banks in the city are also encouraged to make loans to farmers to build greenhouses, while the technical support provided by the government and other social organizations is a popular way to help farmers improve greenhouse management.
Changes in dietary patterns present another potential factor that has induced expansion in the land used for greenhouses.It has been reported that dietary patterns globally are changing from low to high calorie intakes [65], and this leads to increased consumption of vegetables compared to cereals.Between 2000 and 2013, the annual cereal consumption of Chinese residents increased by 19.75%, while that of vegetables increased by 48.43% [66], a trend that is projected to continue.According to the Guideline of Food and Nutrition Development in China [67], the vegetable consumption of Chinese residents is projected to increase by 44.48% by 2020 compared to consumption in 2014, while the consumption of grain is expected to decreased by 4.26%.These changes in Chinese dietary patterns are likely to induce further conversions from arable land used for the cultivation of grain to greenhouse land for vegetable production.
The global food surplus may also affect greenhouse land expansion in Shouguang.For example, food availability in 2010 was estimated to be 20% higher than required globally [68], leading to a waste problem.Thus, as a country that suffers from a severe food deficit, China has long emphasized grain production; following 12 years of continuously increasing output, grain production in China reached 620 million tons in 2015 [69] and led to large volumes of food waste.It has been reported that household food waste in China equates to roughly 2.5% of grain output each year [70].Thus, this Chinese food surplus provided another potential driving force for the reduction of grain production on arable land and the expansion of vegetable production in greenhouses.
In order to reduce the effects of spatial autocorrelation, we randomly selected 30% of our samples and performed a binary logistic regression based on previous results [40,41,55,56,58].However, the use of this sampling method to assess the robustness of regression results is problematic; regression of a different set of samples might lead to slight changes in results and influence the analysis of driving forces.More robust assessment of regression results using a wider range of samples should be emphasized in future work.

Conclusions
This research focused on detecting temporal and spatial patterns in greenhouse land use as well as the forces driving this change in the rapidly developing Chinese city of Shouguang.The results of this study show that greenhouse land use in Shouguang increased by 50.51% between 2000 and 2015.However, 90.39% of this expansion occurred in the period between 2010 and 2015.The area of land used for greenhouses in Shouguang overtook other arable land in 2010 and had reached 2.07 times this level by 2015.Grain production in this city is gradually being replaced by greenhouse vegetables.
Our comparisons show that greenhouse area in Shouguang has preferentially expanded in areas far from the coastal zone, in regions close to rural settlements, and in those that are characterized by low elevations and slight precipitation.In addition, since the availability of flat areas for greenhouse construction is limited, new increases in this land use have tended to be in the hilly regions of Shouguang.At the macroscale, encouragements due to government policy, changes in dietary patterns, and the global food surplus may also have forced expansions in greenhouse land.

Figure 1 .
Figure 1.Location and classification verification for Shouguang, using ground reference data from Google Earth.We selected 50 samples classified as greenhouse land in 2015 and verified them using Google Earth images.Of these samples, the six highlighted in red were incorrectly classified as greenhouse land, while of the ten enlarged boxes, 1 and 8 actually comprise arable land.The other six comprise greenhouse land.

Figure 2 .
Figure 2. Flow chart illustrating greenhouse land mapping.The process of greenhouse land cover mapping was divided into five steps, as shown on the right of this figure.Steps one and two focus on the processing of satellite images, while steps three and four are classification.

Figure 2 .
Figure 2. Flow chart illustrating greenhouse land mapping.The process of greenhouse land cover mapping was divided into five steps, as shown on the right of this figure.Steps one and two focus on the processing of satellite images, while steps three and four are classification.
Sustainability 2017, 3, 359 5 of 15 linearly separable classes in a decision boundary space[35] (see Supplementary Material for a full mathematical description).In brief, a SVM forecasts the destination values of test data by developing a model using training data.The basic process of running a SVM includes steps to select a training data set, compute separability, training using a SVM algorithm, and the acquisition of classification results.Thus, for this study, land use in Shouguang was divided into seven types encompassing greenhouses, other arable, forested, grassland, water area, construction, and unused.

Figure 3 .
Figure 3. Maps to show the distribution of land use in: 2000 (a); 2010 (b); and 2015 (c).The arable land use category includes paddies and dry land, while the forest category includes forests, shrubs, and woods.The grassland category includes grasses of different densities, while the water area category includes rivers, lakes, reservoirs, and ponds.The construction land use category includes urban areas and rural settlements, while the unused land category includes sandy land, as well as gobis and salinas.

Figure 3 .
Figure 3. Maps to show the distribution of land use in: 2000 (a); 2010 (b); and 2015 (c).The arable land use category includes paddies and dry land, while the forest category includes forests, shrubs, and woods.The grassland category includes grasses of different densities, while the water area category includes rivers, lakes, reservoirs, and ponds.The construction land use category includes urban areas and rural settlements, while the unused land category includes sandy land, as well as gobis and salinas.

Figure 4 .
Figure 4. ROC curve for: expansion (a); and contraction (b) regression models.This ROC curve was used to test the degree of fit of each binary logistic regression; the area ratio under the curve corresponds to a ROC value between 0 and 1.Thus, larger is the value, the higher is the explanatory power of the logistic regression model.

Figure 4 .
Figure 4. ROC curve for: expansion (a); and contraction (b) regression models.This ROC curve was used to test the degree of fit of each binary logistic regression; the area ratio under the curve corresponds to a ROC value between 0 and 1.Thus, larger is the value, the higher is the explanatory power of the logistic regression model.

Figure 5 .
Figure 5. Changes in greenhouse land use: between 2000 and 2010 (a); and between 2010 and 2015 (b).A decrease denotes a change from greenhouse land to another land use type, while an increase denotes the opposite.

Figure 5 .
Figure 5. Changes in greenhouse land use: between 2000 and 2010 (a); and between 2010 and 2015 (b).A decrease denotes a change from greenhouse land to another land use type, while an increase denotes the opposite.

Figure 6 .
Figure 6.Areas of land used for greenhouses within various buffers at 2 km intervals around urban areas.The area of greenhouse land within these buffers gradually decreases closer to urban areas.

Figure 6 .
Figure 6.Areas of land used for greenhouses within various buffers at 2 km intervals around urban areas.The area of greenhouse land within these buffers gradually decreases closer to urban areas.

Table 2 .
Changes in area (km 2 ) of land used for greenhouses between 2000 and 2015.
Abbreviations: GHL, greenhouse land; Unchanged GHL, GHL that did not undergo land use conversion; Increase, newly added GHL converted from other land use types; Decrease, GHL converted into other land use types.

Table 2 .
Changes in area (km 2 ) of land used for greenhouses between 2000 and 2015.
Abbreviations: GHL, greenhouse land; Unchanged GHL, GHL that did not undergo land use conversion; Increase, newly added GHL converted from other land use types; Decrease, GHL converted into other land use types.

Table 3 .
Analysis of forces driving expansion in greenhouse land use between 2000 and 2015.

Table 3 .
Analysis of forces driving expansion in greenhouse land use between 2000 and 2015.
Abbreviations: Swc, soil water content in 2008; ∆precip, changes in precipitation between 2000 and 2012; B, regression coefficient values; S.E., standard error; Wald, statistical value that tests the influence of an independent variable on a dependent one; df, degree of freedom; Sig., P-value for which "**" and "***" denote a significant coefficient at an alpha level of 0.01 and 0.001, respectively; Exponent (B), odds-ratio for each variable.

Table 4 .
Analysis of forces driving contraction in greenhouse land use between 2000 and 2015.