The spatial information entropy has been widely used in studies such as geography, ecology and life sciences to describe or explain the heterogeneity of data or the distribution of “things” in space due to its ability to characterize the amount of information, uncertainty, heterogeneity and other concepts, in this study we applied it to quantify the global change of GCL clustering in Shandong over the past 30 years, which can be defined as (Leibovici et al., 2014): H ( X ) = ∑ i = 1 I p ( x i ) log 1 p ( x i ) where X is the discrete variable containing the corresponding value x i (in this study is divided into x 0 and x 1 , representing other land use types and GCL, respectively); p ( x i ) is the probability mass function (PMF) of X , and its result is based on PMF quantifying the average information of X . The lower Shannon information entropy value indicates that the spatial distribution of GCL is more clustered, and the higher Shannon information entropy value indicates that the spatial distribution of GCL is more discrete.

In order to further reflect the spatial heterogeneity of GCL clustering, the density of GCL in each year was classified based on the kernel density estimation (KDE). KDE is a non-parametric test method used to estimate unknown density functions, which has been widely used in natural disasters, public health, industrial spatial layout hotspot analysis and detection, etc. Its calculation formula is as follows (Cai et al., 2013): f ^ ( x , y ) = 3 n h 2 π ∑ i = 1 n ( 1 − ( x − x i ) 2 + ( y − y i ) 2 h 2 ) 2 where f ^ ( x , y ) is the kernel density of GCL in the grid to be estimated, h is the kernel density estimation bandwidth (set to a 30-year average bandwidth of 42,926 m that guaranted the comparability of each KDE result over the past 30 years), x i , y i are the coordinates of the grid centroids to be estimated, n is the number of all sample points within the bandwidth, and x , y are the coordinates of the central sample points within the bandwidth.

One of the significant advantages of remote sensing-based data is the ability to continuously observe the same area and thus obtain continuous information (Maus et al., 2016). Therefore, we proposed a time-series segmentation and sliding algorithm (TSS) based on the annual remote sensing mapping products of GCL, aiming to analyze the spatiotemporal dynamics of GCL in cumulative duration, demolition frequency and expansion trajectory.

As shown in Figure 2, the annual GCL maps were stacked to obtain the classification labels of each pixel from 1989 to 2018, and the temporal label sequences that need to be analyzed were extracted based on the temporal segmentation window N. For the cumulative duration analysis, the cumulative duration of GCL for each pixel at that period can be obtained by accumulating all labels of the segmented sequence. For the demolition frequency analysis, the time sliding window n was further set to 2, and the frequency of “01” or “10” combinations at that period can be obtained by sliding backward from the starting year in steps of 1. For the expansion trajectory analysis, we set the time sliding window n to 1, slide backward from the starting year in steps of 1, and assign the year information to each pixel when “1” appears for the first time.

Regarding the expansion intensity, previous studies have focused on urban expansion as represented by built-up or impervious areas. However, the major difference between GCL and urban is that urban expansion is generally considered to be irreversible (Li et al., 2015). In terms of GCL, due to the fluctuation of demolition and construction, the information between each year would be lost if the urban expansion intensity calculation rule is followed. Therefore, we proposed an annual expansion index (AEI) based on annual growth to reflect the true expansion intensity over a certain period, which was calculated as follows: AEI = ∑ i = 1 n GC L i n where G C L i represents the annual GCL growth from year i to year i + 1 in the specific grid, and n is the study period.

Spatially stratified heterogeneity is a universal characteristic of geographic data, and for geographic data in different layers, if each layer of independent variables has an influence on the dependent variable, its spatial distribution should be similar (Wang et al., 2010). Geodetector is such a statistical method for detecting the spatial heterogeneity and revealing the driving factors behind it based on this assumption (Wang et al., 2016). In this study, we applied the factor detection in Geodetector to explain the local driven mechanism between GCL expansion and the selected variables, which is calculated as follows: q = 1 − ∑ h = 1 L N h σ h 2 N σ 2 where L is the total number of stratifications of the dependent variable Y or the independent variables X, N h and N are the number of units in stratum h and the whole region, respectively, σ h 2 and σ 2 are the variance of the dependent variable Y in stratum h and the whole region, respectively. q ∈ [ 0,1 ] , and the closer the value of q is to 1, the stronger explanatory power of this independent variable for the spatial heterogeneity of the dependent variable (Zhang et al., 2019).

Granger causality test is to determine whether a time-series change in one variable is caused by a time-series change in another variable (Freeman, 1983). In this study, we first used the Augmented Dickey-Fuller test (ADF) to test whether the time-series variables were stationary, which is based on the following unit root tests (Krämer, 1998): Δ X t = α + β t + δ X t − 1 + ∑ i = 1 m β i Δ X t − i + ϵ t where the ADF tests whether the series xt rejects the original hypothesis by judging the estimated value of δ. When δ is less than the critical value of the correlation test, the series Xt is a stationary series, otherwise, the series Xt is a non-stationary series, and a further differential process is required to judge the stationarity of the differenced series.

For variables that are still non-stationary and belong to the same-order of single integer after differencing, further cointegration test based on the Engle-Granger method is required to determine whether there is a long-term stable relationship between the variables, which is conduct by establishing the least squares regression equation as following (Lee and Lee, 2015): Y t = b 0 + b 1 X t + ε t when ε is less than the critical value of the correlation test, the residual series ε t remains a cointegration relationship between variables, otherwise, there is no cointegration relationship between variables in the residual series ε t .

After testing the stationarity of the variables and the cointegration relationship between the variables, Granger causality tests were carried out on the time-series variables that were stationary and had cointegration relationships, and the following two regression models were estimated: Y t = β 0 + ∑ i = 1 m β i Δ Y t − i + ∑ i = 1 l α i X t − i + μ t X t = δ 0 + ∑ i = 1 m δ i Δ X t − i + ∑ i = 1 l λ i Y t − i + ν t where the original hypothesis is that Y t is not the Granger cause of X t and X t is not the Granger cause of Y t , and the F-test statistic is performed with the sum of squares of residual μ t and the sum of squares of residual ν t , and if the F statistic is higher than the critical value corresponding to the significance level, the original hypothesis is rejected, which means that there is Granger causality.

The spatial information entropy was calculated for all years using the annual maps of GCL in Shandong (Figure 3A). The unit area entropy value of GCL exhibited a decreasing trend in general, with the maximum value appearing in 1990 and the minimum value appearing in 2017, which can be divided into three stages: 1) from 1989 to 1996, the average unit area entropy value was 0.044; 2) from 1996 to 2005, the average unit area entropy value was 0.036; 3) from 2005 to 2018, the average unit area entropy value was 0.031. It indicated that the spatial clustering of GCL in Shandong over past 30 years was gradually increasing in general, where the highest level of spatial clustering was in 2017 and the lowest was in 1990. Combined with the changes in its total area, it can be seen that entropy value per unit area decreased continuously with the expansion of greenhouses in general, which means that the spatial clustering has a positive synergistic relationship with the change of total area of GCL.

We further explored the spatial heterogeneity of GCL clustering over 30 years by using the KDE method to partition the spatial distribution density (Figure 3B). The percentage of low density zones decreased from 96.82% in 1989 to 66.38% in 2018, with a decrease of 30.44% and an average annual decrease of 1.05%; the medium-low density zones increased from 0.54% in 1989 to 23.51% in 2018, with an increase of 22.97% and an average annual increase of 0.79%; the medium density zones increased from 0.13% in 1995 to 7.80% in 2018, with an increase of 7.67% and an average annual increase of 0.33%; the medium-high density zones increased from 0.30% in 1998 to 1.36% in 2018, with an increase of 1.06% and an average annual increase of 0.05%; and the high density zones increased from 0.43% in 2001 to 0.95% in 2018, with an increase of 0.47% and an average annual increase of 0.03%. It indicated that despite the GCL expansion in Shandong was toward continuous clustering, but mainly dominated by the medium-low and medium densities.

Based on the TSS algorithm, we calculated the cumulative duration of GCL over the past 30 years. As shown in Figure 4, GCL with the cumulative duration of 1 year accounted for 55.32% and showed a scattered distribution pattern mainly around the central mountainous area, and when the cumulative duration reached more than 2 years, its percentage rapidly decreased to less than 18.03%. And it can be found that the areas where the cumulative duration of GCL was more than 15 years were mainly distributed in northwestern Weifang, southwestern Linyi and western of Liaocheng, showing a triangular distribution pattern, which is similar to the pattern of the spatial clustering, indicating that the higher level of spatial clustering can extend the cumulative duration of GCL. In addition, it can be seen from the detail enlargements in Figures 4A–C that there was a certain spatial continuity in the extension of the cumulative duration.

Then the demolition frequency of GCL over the past 30 years was also calculated based on the TSS algorithm. As shown in Figure 5, the demolition frequencies of GCL in the range of 0.8∼1.0 accounted for 60.65% during the 30-year period, and showed a scattered distribution pattern mainly around the central mountainous area. The demolition frequencies in the range of 0.6∼0.8–0.4∼0.6 accounted for about 13%, and when the demolition frequency dropped to 0.4, its percentage decreased rapidly to below 12.33%. At the same time, it can be found that the areas with demolition frequency below 0.2 were mainly distributed in northwestern Weifang, southwestern Linyi and western Liaocheng, which is similar to the patterns of the spatial clustering as well as the cumulative duration, indicating that the higher level of spatial clustering can not only ensure the duration but also the continuity of GCL. In addition, from the detailed enlarged Figures 5A–C, it can be seen that there is a certain continuity in the change of demolition frequency, in which the area with low demolition frequency tends to form the area with relatively low demolition frequency around it.

Based on the TSS algorithm, we also identified the year in which GCL first appeared in Shandong province over the past 30 years. As shown in Figure 6, the expansion trajectory of GCL roughly formed a circular expansion pattern around the central mountainous area, among which the expansion trajectories were more disorderly in Dezhou, Jinan, Binzhou and Dongying, and more orderly in Weifang, Liaocheng, Linyi, Zaozhuang and Jining. Combined with its spatial clustering pattern, it can be seen that the expansion trajectories were more orderly in the GCL concentration areas. In terms of each period, in period 2004∼2008 and period 2014∼2018 accounted for the most percentage of GCL expansion with 25.18% and 20.62%, respectively, followed by period 1999∼2003 (18.92%) and period 2009∼2013 (18.40%), and the least was period 1989∼1993 (7.22%) and period 1994∼1998 (9.67%). In addition, from the detail enlarged Figures 6A–C, in northwestern Weifang, southwestern Linyi and western Liaocheng showed an obvious pattern of circle outward expansion trajectories.

We further calculated the GCL expansion intensity in the above-mentioned periods at the 1 km grid scale, and classified it into five types: slow expansion, low-speed expansion, medium-speed expansion, high-speed expansion, and rapid expansion based on the natural break method (De Smith et al., 2007). As shown in Figure 7, the areas with rapid expansion intensity were mainly located in the northwestern plains represented by Weifang in period 1989∼1999 and continued to expand in period 1999∼2004, the southern plains represented by Linyi began to appear in the rapid expansion intensity in period 1999∼2004 as well. In period 2004∼2014, the areas with rapid expansion intensity shrank to the plain areas of northwestern and southern Shandong, and appeared in the plain areas represented by Liaocheng in the west. After that, the new areas with rapid expansion intensity appeared in the southeastern of Weifang and the junction of Jining and Heze from 2014 to 2018.

In order to reveal the differences in the local driven mechanism of GCL at different periods, we first divided the past 30 years into three different periods based on the annual change rate (ACR) of the total area of GCL, which were the stable period (1989∼1996, with an ACR of −2.28 kha/year), the growth period (1997∼2005, with an ACR of 23.09 kha/year) and the fluctuating period (2006∼2018, with an ACR of 14.03 kha/year), respectively. Meanwhile, the ACR of the selected local dynamic variables were calculated and identified as “ACR_XX” within each period. For those local stable variables that vary little over time, were also quantified based on factor detection at different periods.

Regarding the results of factor detection for the local dynamic variables in each period (Figure 8A), the top three q-values were ACR_BER (0.1164), ACR_EI (0.1153) and ACR_AE (0.1055) in the stable period, ACR_LRS (0.0389), ACR_RO (0.0386) and ACR_AM (0.0354) in the growth period, and ACR_AE (0.0365), ACR_AM (0.0328) as well as ACR_BER (0.0316) in the fluctuating period. It indicated that 1) budget expenditure for rural development reflecting government support and average earning of local famers reflecting farmers’ willingness were both the major driving factors for farmers to promote GCL in the stable as well as the fluctuating period, 2) effective irrigated area reflecting agricultural production conditions only played a major driving role in the stable period, 3) local retail sales reflecting local consumer demand and road mileage reflecting transportation conditions played a key driving role in the growth period.

Regarding the results of factor detection for the local stable variables in each period (Figure 8B), the top three q-values were ST (0.0410), AVE (0.0258) and RIV (0.0073) in the stable period, AVE (0.0147), ST (0.0093) and DU (0.0040) in the growth period, and AVE (0.0108), RIV (0.0069) as well as ST (0.0410) in the fluctuating period. It indicated that 1) elevation and soil type were always the main driving force for GCL expansion in all local stable variables over time, 2) the driving force of the distance to rivers exhibited a significant performance in both the stable and fluctuating periods among the location aspect, 3) compared with elevation, the driving force of slope showed an insignificant performance in all periods.

In order to reveal the driven mechanism of the supply and demand of vegetables in external markets on GCL expansion in Shandong, the following target and explanatory time-series were selected: total area of GCL in Shandong province ( Y ); vegetable production in Guangxi province ( X 1 ), Hainan province ( X 2 ), Sichuan province ( X 3 ), Hunan province ( X 4 ), Hubei province ( X 5 ), Guangdong province ( X 6 ), Jiangsu province ( X 7 ); Total value of agricultural exports in Shandong province ( X 8 ); consumption demand of vegetable in Beijing ( X 9 ) and Shanghai ( X 10 ). To eliminate the effects caused by heteroskedasticity of each time-series, which have been treated with logarithm and denoted as “ ln ⁡ X n ”.

As shown in Supplementary Material S2, the ADF values of Y , X 3 , X 4 , and X 7 in the original series were greater than all critical values, indicating that these four original series were non-stationary. The ADF values of X 1 , X 2 , X 5 , and X 6 in the original series were less than the critical values at 1% or 5% significance level, indicating that these four original series were stationary. The results of further differencing of each series showed that series X 2 was stationary after the second-order differencing and the rest series were stationary after the first-order differencing, which means that the target series Y and the explanatory series X 1 , X 3 , X 4 , X 5 , X 6 , and X 7 belong to the first-order of single integer after differencing, which satisfied the premise of the cointegration test, while the target series Y and the explanatory series X 2 belong to the second-order of single integer after differencing, which also satisfied the premise of the cointegration test. Further based on the Engle-Granger method, the target series was modeled with each explanatory series to determine whether there was a cointegration relationship between them. As shown in Supplementary Material S3, the ADF test statistics between the target series and each explanatory series were all less than the critical values at all significance levels, indicating that there was a cointegration relationship between the target series and each explanatory series, and the Granger causality test between them can be further analyzed.

Therefore, the target series after differencing was selected for Granger causality tests with each explanatory series after the same-order differencing at lag 1, lag 3, and lag 5 (Supplementary Material S4). From the perspective of the supply of vegetables in the external markets, vegetable production change in Guangxi, Hainan, Hunan, and Hubei province were a one-way Granger causality of the development of GCL in Shandong province, vegetable production change in Sichuan province was a two-way Granger causality of the development of GCL in Shandong province, and no Granger causality between vegetable production change in Guangdong province and the development of GCL in Shandong. From the perspective of the consumption demand of vegetable in the external markets, vegetable consumption in Beijing was a one-way Granger causality of the development of GCL in Shandong, and no Granger causality between both vegetable consumption in Shanghai as well as abroad and the development of GCL in Shandong.