The Intensity Analysis of Production Living Ecological Land in Shandong Province, China

Due to the limitedness of land, the coordinated development of production, living, and ecological (PLE) land is essential for sustainable development. A clear understanding of PLE land change is necessary given the increased human activities, especially in developing regions. This study first reclassified remotely sensed landuse maps in Shandong province into PLE land. Then the spatiotemporal change of PLE land between 2000 and 2015 was analyzed using spatial trajectory and intensity analysis methods. The results show that the rate of PLE land change in the interval of 2000–2005 was the highest, and it kept decreasing during 2005–2010 and 2010–2015. The overall quantity component accounts for more difference than the exchange and shift components for all intervals. At the category level, the largest quantity component of change was PE land loss, followed by LP land gain. LP land gain targeted PE land intensively in each interval. The loss of E land was mainly to PE land in terms of size, and to PL land in terms of intensity. The findings of this study contribute to a deeper understanding of the spatiotemporal transitions of PLE land in Shandong province, which could help policy making for PLE land regulation.


Introduction
Given current economic development and rapid urbanization and their effects on the ecological environment, the balance of land multi-function has attracted increasing attention from scholars, planners, and policy makers [1]. How Production, Living, and Ecological land (PLE) is being spatiotemporally adjusted in the process of urbanization deserves to be studied. Optimizing PLE land has become an important principle in spatial planning and sustainable development strategy drafting [2]. However, many previous studies focused on the spatiotemporal change of remotely-sensed land-use types [3], with inadequate analysis on the change of multifunction-based land patterns. A comprehensive investigation on PLE land transition is more necessary than ever.
PLE land was emphasized in the report of the 18th National Congress of the Communist Party of China in 2012 for promoting the coordinated and sustainable development of economy and ecology [4]. Land can provide a variety of products and services for humans. Collectively, the production function refers to ability of providing agricultural, forestry and industrial products. The living function refers to the role of human habitation and survival. The ecological function is the capacity to improve the ecological environment and guarantee ecology security. In the past, people paid more attention to the production and living function of land. However, the ecological function is also very important for rational and sustainable utilization of land. Recently, studies on PLE land have increased. There are

Data
The land use data for 2000,2005,2010, and 2015, with 30 m × 30 m spatial resolution, were obtained from the Resource and Environment Data Cloud Platform [41]. The dataset was produced by the Chinese Academy of Sciences (CAS) based mainly on the Landsat TM/ETM+/OLI and GF-2 data using the human-machine interactive method. More details of the interpretation technology were presented in a previous study [42]. The interpretation results were validated against extensive field surveys, and the average classification accuracy was higher than 90% [43]. The land use types include 6 primary classes and 25 secondary classes. For example, the cropland includes paddy and non-irrigated uplands. The forest type includes high density forest, low density forest, shrub land, and orchards. The grassland includes high-, medium-, and low-density grassland. The water area includes rivers, canals, lakes, permanent glaciers, etc. The built-up land includes urban, rural settlements, and industrial land. The unused land includes desert, Gobi, wetland, bare land, rock, etc. This land use dataset has been widely used in land use change related research at the national and local levels [3,43].

Data
The land use data for 2000,2005,2010, and 2015, with 30 m × 30 m spatial resolution, were obtained from the Resource and Environment Data Cloud Platform [41]. The dataset was produced by the Chinese Academy of Sciences (CAS) based mainly on the Landsat TM/ETM+/OLI and GF-2 data using the human-machine interactive method. More details of the interpretation technology were presented in a previous study [42]. The interpretation results were validated against extensive field surveys, and the average classification accuracy was higher than 90% [43]. The land use types include 6 primary classes and 25 secondary classes. For example, the cropland includes paddy and non-irrigated uplands. The forest type includes high density forest, low density forest, shrub land, and orchards. The grassland includes high-, medium-, and low-density grassland. The water area includes rivers, canals, lakes, permanent glaciers, etc. The built-up land includes urban, rural settlements, and industrial land. The unused land includes desert, Gobi, wetland, bare land, rock, etc. This land use dataset has been widely used in land use change related research at the national and local levels [3,43].

Classification of Production-Living-Ecological Land
Some land use types often serve multiple functions [6]. For example, beyond the primary function of food production, cropland can also provide ecological benefits [44]. Urban area serves both the function of human living and production. Each land use type has a primary and secondary function. Zhang et al. (2017) put forward a Production Living Ecological land classification system for China, which consists of four categories: productive ecological land (PE), ecological productive land (EP), Living productive land (LP), and ecological land (E) [8]. The PE land includes cropland and orchard, which have important production functions, but also have some ecological functions. The EP land includes pasture, timberland, and aquaculture area. The LP land includes urban built-up and rural settlement. The E land includes ecological regulate regions, such as forest and grassland where are non-timber and non-pasture, and ecological conservation regions that are conserved as unused land.
In this study, we adopted the classification system of Zhang et al. [8]. We also made some modifications according to the local land use characteristics of the study area. We added the productive living land (PL) in this study. The reason is that the proportion of industrial built-up land in Shandong province is much larger than that within the national scale. The primary function of the industrial land is production, and they also serve the function of human living. Therefore, we classified the industrial built-up land as PL type, which are included in LP land in the classification system of Zhang et al. [8]. As a result, we reclassified the land use data of Shandong province into five function types: PE, EP, LP, PL, and E land based on the classification principle mentioned above. Then the spatiotemporal analysis on the PLE land between 2000 and 2015 was performed to describe PLE land succession on a pixel-by-pixel basis.

Intensity Analysis
Intensity analysis is a method of comparing maps with the same categories, spatial extent, and resolution [18][19][20]. For example, it could be used to analyze the land use maps at different time points and shows the size and intensity of land use change at the interval, category, and transition levels. Similarly, the spatiotemporal changes of PLE land between 2000 and 2015 were measured using the intensity analysis method in this study.
Intensity analysis quantifies the changes of size and intensity at increasingly detailed levels. Firstly, the interval level measures the total scale and speed of change in different time intervals. Equations (1) and (2) are as follows: where S t is the annual change rate during the period [Y t ,Y t+1 ]; Y t is the year at time point t. Subscripts i and j denote the categories, and the total number of categories is J, which is five in this study (PE, EP, PL, LP, and E land). C tij is the area that changes from category i in year Y t to j in year Y t+1 .
where U is the uniform value, which is the mean change rate. T is the number of time points, which is four in this study (2000, 2005, 2010, and 2015). If S t is less than U, then the change rate of that period is relatively slow. Likewise, if S t is larger than U, then the change rate of that period is relatively fast. Furthermore, the overall change could be separated into three components, which are quantity, interval, and shift [19,20]. Quantity difference is the total amount of change between maps of different time points (Figure 2a); exchange is the spatial allocation difference that pairwise transitions cause Sustainability 2020, 12, 8326 5 of 17 ( Figure 2b); shift is the spatial allocation difference that non-pairwise transitions cause (Figure 2c). The equations to compute these three components are presented in Pontius et al. [19,20].
Sustainability 2020, 12, x FOR PEER REVIEW 5 of 17 The second level is category level, which produces the change trends for each individual category. For each period, the annual gain and loss of each category are calculated by equations (3) and (4): where Gaintj is annual gross gain intensity of category j and Lossti is annual loss intensity of category i for the period [Yt, Yt+1]. St computed by Equation (1) where Changetj is gross change of category j for the time interval t. Changetj includes three components, which are the quantity component Quantitytj and the allocation components Exchangetj and Shifttj. The third level is transition level, which measures the size and intensity that a category gains from other categories during a certain time interval. The intensities of transition to category n from each category i (i ≠ n) are calculated using Equation (9): where Ctin is the area that changes from category i in the year Yt to category n in the year Yt+1. Equation (9) thus gives J-1 transition intensities each period. The uniform transition intensity from other categories to category n during each period (Wtn) is computed using Equation (10).
Then if Rtin is less than Wtn, the gain of category n avoids j during the certain interval. Alternatively, if Rtin is larger than Wtn, the gain of category n targets j during that interval. The second level is category level, which produces the change trends for each individual category. For each period, the annual gain and loss of each category are calculated by Equations (3) and (4):

Spatiotemporal Distributions of Prodution-Living-Ecological Land Changes
where Gain tj is annual gross gain intensity of category j and Loss ti is annual loss intensity of category i for the period [Y t , Y t+1 ]. S t computed by Equation (1) is the uniform intensity for each period. If the Gain tj or Loss ti is less than S t , then the change of category j or i during that period is dormant. On the contrary, if the Gain tj or Loss ti is larger than S t , then the change of category j or i is active. For the time interval t, the gross change of category j and its three components are computed by Equations (5)-(8) [45]. Change tj = Gain tj + Loss tj (5) Quantity tj = Gain tj − Loss tj (6) Shi f t tj = Change tj − Quantity tj − Exchange tj (8) where Change tj is gross change of category j for the time interval t. Change tj includes three components, which are the quantity component Quantity tj and the allocation components Exchange tj and Shift tj . The third level is transition level, which measures the size and intensity that a category gains from other categories during a certain time interval. The intensities of transition to category n from each category i (i n) are calculated using Equation (9): where C tin is the area that changes from category i in the year Y t to category n in the year Y t+1 . Equation (9) thus gives J − 1 transition intensities each period. The uniform transition intensity from other categories to category n during each period (W tn ) is computed using Equation (10). Then if R tin is less than W tn , the gain of category n avoids j during the certain interval. Alternatively, if R tin is larger than W tn , the gain of category n targets j during that interval.

Spatiotemporal Distributions of Prodution-Living-Ecological Land Changes
The spatial patterns of each category are shown Figure 3. The main type in the study area is PE land, which was 67.7% of the total area in 2000. The production function of PE land here is agricultural provision. The EP and E lands are mainly located in the upland area. Over the period from 2000 to 2015, the expansion of LP land was remarkable. The area of LP land increased rapidly from 17,532 km 2 to 20,102 km 2 . Most of the newly expanded LP land was in the neighborhood of the original LP land. Due to the expansion of LP land, other categories were encroached by it to varying extents. A large proportion of the expansion of LP land was taking place in the PE land ( Figure 4). Consequently, PE land suffered the largest decrease size. The reason for this may be because the PL lands and the PE lands are frequently located in proximity to each other. The spatial patterns of each category are shown Figure 3. The main type in the study area is PE land, which was 67.7% of the total area in 2000. The production function of PE land here is agricultural provision. The EP and E lands are mainly located in the upland area. Over the period from 2000 to 2015, the expansion of LP land was remarkable. The area of LP land increased rapidly from 17,532 km 2 to 20,102 km 2 . Most of the newly expanded LP land was in the neighborhood of the original LP land. Due to the expansion of LP land, other categories were encroached by it to varying extents. A large proportion of the expansion of LP land was taking place in the PE land ( Figure 4). Consequently, PE land suffered the largest decrease size. The reason for this may be because the PL lands and the PE lands are frequently located in proximity to each other.    Figure 5a represents the uniform intensity value, which was calculated using Equation (2), and it was the mean change rate of the three time intervals. Furthermore, the overall annual change of each time interval was divided into three components: quantity, exchange, and shift ( Figure 5b). The exchange and shift are two components of allocation difference, where exchange is the allocation difference that caused by Sustainability 2020, 12, 8326 7 of 17 pairwise transitions, and shift is the allocation difference caused by non-pairwise transition ( Figure 2). The quantity component accounted for more difference than allocation components (exchange and shift) for all the three periods. Exchange was larger than shift during all the three periods.

Interval Level
The results of interval level intensity analysis are shown in Figure 5 Figure 5(a) represents the uniform intensity value, which was calculated using Equation (2), and it was the mean change rate of the three time intervals. Furthermore, the overall annual change of each time interval was divided into three components: quantity, exchange, and shift ( Figure 5 (b)). The exchange and shift are two components of allocation difference, where exchange is the allocation difference that caused by pairwise transitions, and shift is the allocation difference caused by nonpairwise transition ( Figure 2). The quantity component accounted for more difference than allocation components (exchange and shift) for all the three periods. Exchange was larger than shift during all the three periods.

Interval Level
The results of interval level intensity analysis are shown in Figure 5 Figure 5(a) represents the uniform intensity value, which was calculated using Equation (2), and it was the mean change rate of the three time intervals. Furthermore, the overall annual change of each time interval was divided into three components: quantity, exchange, and shift ( Figure 5 (b)). The exchange and shift are two components of allocation difference, where exchange is the allocation difference that caused by pairwise transitions, and shift is the allocation difference caused by nonpairwise transition (Figure 2). The quantity component accounted for more difference than allocation components (exchange and shift) for all the three periods. Exchange was larger than shift during all the three periods.

Category Level
The areas of the annual gain or loss of each category during the three periods are revealed in Figure 6. The net change of each category was also derived. The PE land had the largest loss, and LP land had the largest gain during all of the three time intervals. The gain size of PL land was also relatively large. The loss in size of E land was big during the first time interval, then decreased obviously during the middle and the latter interval. The LP land loss sizes were small across the three periods. The sizes of gain and loss of all the categories were largest during the period 2000-2005, and were the smallest during the period 2010-2015. Figure 6. The net change of each category was also derived. The PE land had the largest loss, and LP land had the largest gain during all of the three time intervals. The gain size of PL land was also relatively large. The loss in size of E land was big during the first time interval, then decreased obviously during the middle and the latter interval. The LP land loss sizes were small across the three periods. The sizes of gain and loss of all the categories were largest during the period 2000-2005, and were the smallest during the period 2010-2015.   land had the largest gain during all of the three time intervals. The gain size of PL land was also relatively large. The loss in size of E land was big during the first time interval, then decreased obviously during the middle and the latter interval. The LP land loss sizes were small across the three periods. The sizes of gain and loss of all the categories were largest during the period 2000-2005, and were the smallest during the period 2010-2015.    (Figure 8a). A closer look at the processes revealed that the PE land includes cropland and orchard, and PE land includes pasture and timberland; therefore, a large part of allocation changed during the former two periods, involving the transition from cropland and orchard to pasture and timberland, and from pasture and timberland to cropland and orchard. The sign of net loss and net gain associated with each category denoted the loss and gain of the quantity component during each interval, respectively. LP and PE land were mainly responsible for the quantity change component during all the three time intervals, while LP and PL land were net gaining categories, and PE, EP, and E land were net losing categories. The process was an expansion of LP land (such as urban area) in the study area, where usually the available land was PE (such as cropland) and E land.

Transition Level
The quantity changes of PE land during the three intervals were net loss. The annual loss amount of PE land is revealed in Figure 9a, while the intensities for the categories to which PE land loses are revealed by Figure 9b. LP land had the largest transition size, followed by PL land. The transition intensity to PL was larger than to LP land, because PL accounted for a larger percentage of the study area than LP land. The red vertical dashed line in Figure 9b indicates the uniform change intensity for PE land loss, which was calculated using Equation (6). The bars of LP land and PL land ended beyond the uniform line. This reveals that when PE land was lost, it tended to lose intensively to LP and PL land more than other categories. That is to say, the loss of PE land targeted LP and PL land but avoided EP and E land during all the time intervals.
The quantity changes of PE land during the three intervals were net loss. The annual loss amount of PE land is revealed in Figure 9(a), while the intensities for the categories to which PE land loses are revealed by Figure 9(b). LP land had the largest transition size, followed by PL land. The transition intensity to PL was larger than to LP land, because PL accounted for a larger percentage of the study area than LP land. The red vertical dashed line in Figure 9(b) indicates the uniform change intensity for PE land loss, which was calculated using Equation (6). The bars of LP land and PL land ended beyond the uniform line. This reveals that when PE land was lost, it tended to lose intensively to LP and PL land more than other categories. That is to say, the loss of PE land targeted LP and PL land but avoided EP and E land during all the time intervals.   of PE land is revealed in Figure 9(a), while the intensities for the categories to which PE land loses are revealed by Figure 9(b). LP land had the largest transition size, followed by PL land. The transition intensity to PL was larger than to LP land, because PL accounted for a larger percentage of the study area than LP land. The red vertical dashed line in Figure 9(b) indicates the uniform change intensity for PE land loss, which was calculated using Equation (6). The bars of LP land and PL land ended beyond the uniform line. This reveals that when PE land was lost, it tended to lose intensively to LP and PL land more than other categories. That is to say, the loss of PE land targeted LP and PL land but avoided EP and E land during all the time intervals.   The quantity changes of PE land were net loss during the first interval, and net gain during the latter two intervals. The annual transition sizes and intensities for the categories from which EP land gains are shown in Figure 11. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. If one category's intensity bar ends on the left of the uniform line, then the transition of this category avoids EP land, and vice versa. The transition size from EP to PE land was largest during the first time interval. The transition intensity reveals that the loss of EP land was targeting PL and E land during all the intervals. The losses intensity of EP land to LP land was also larger than the uniform line during the latter two periods.
The quantity changes of PE land were net loss during the first interval, and net gain during the latter two intervals. The annual transition sizes and intensities for the categories from which EP land gains are shown in Figure 11. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. If one category's intensity bar ends on the left of the uniform line, then the transition of this category avoids EP land, and vice versa. The transition size from EP to PE land was largest during the first time interval. The transition intensity reveals that the loss of EP land was targeting PL and E land during all the intervals. The losses intensity of EP land to LP land was also larger than the uniform line during the latter two periods. For transition to EP land, the annual transition sizes and intensities for the categories from which EP land gains are presented in Figure 12. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. The transition sizes from PE to EP land are the largest during all of the three time intervals. However, the transition intensity reveals that the gain of EP land was avoiding PE land during all the intervals. Therefore, the large transition area from PE land to EP land could be attributed to the large percentage of PE land with the study area. The gain of EP land was targeting E land more than PE land in terms of transition intensity. The interval level and categorical level analysis above show that the gain of LP land was largest both in terms of size and of intensity, followed by PL land. At the transition level, the annual transition size and the intensity for the categories from which LP land gained are shown in Figure 13. LP land gained first from PE land, and then from PL and EP land in terms of the annual transition For transition to EP land, the annual transition sizes and intensities for the categories from which EP land gains are presented in Figure 12. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. The transition sizes from PE to EP land are the largest during all of the three time intervals. However, the transition intensity reveals that the gain of EP land was avoiding PE land during all the intervals. Therefore, the large transition area from PE land to EP land could be attributed to the large percentage of PE land with the study area. The gain of EP land was targeting E land more than PE land in terms of transition intensity.
The quantity changes of PE land were net loss during the first interval, and net gain during the latter two intervals. The annual transition sizes and intensities for the categories from which EP land gains are shown in Figure 11. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. If one category's intensity bar ends on the left of the uniform line, then the transition of this category avoids EP land, and vice versa. The transition size from EP to PE land was largest during the first time interval. The transition intensity reveals that the loss of EP land was targeting PL and E land during all the intervals. The losses intensity of EP land to LP land was also larger than the uniform line during the latter two periods. For transition to EP land, the annual transition sizes and intensities for the categories from which EP land gains are presented in Figure 12. The red vertical dashed line indicates the uniform change intensity for the gain of EP land during each interval. The transition sizes from PE to EP land are the largest during all of the three time intervals. However, the transition intensity reveals that the gain of EP land was avoiding PE land during all the intervals. Therefore, the large transition area from PE land to EP land could be attributed to the large percentage of PE land with the study area. The gain of EP land was targeting E land more than PE land in terms of transition intensity. The interval level and categorical level analysis above show that the gain of LP land was largest both in terms of size and of intensity, followed by PL land. At the transition level, the annual transition size and the intensity for the categories from which LP land gained are shown in Figure 13. LP land gained first from PE land, and then from PL and EP land in terms of the annual transition The interval level and categorical level analysis above show that the gain of LP land was largest both in terms of size and of intensity, followed by PL land. At the transition level, the annual transition size and the intensity for the categories from which LP land gained are shown in Figure 13. LP land gained first from PE land, and then from PL and EP land in terms of the annual transition area. However, the fact that the LP land gained more intensively from PL land than from PE land during the first two time intervals is revealed in Figure 13b. The reason is that the initial size of PE land was larger than PL land. Then the situation changed for the interval of 2010-2015, during which the bar for PL land was to the left of the uniform line. PE land was the only category that ended beyond the uniform line during all of the three time intervals.
area. However, the fact that the LP land gained more intensively from PL land than from PE land during the first two time intervals is revealed in Figure 13(b). The reason is that the initial size of PE land was larger than PL land. Then the situation changed for the interval of 2010-2015, during which the bar for PL land was to the left of the uniform line. PE land was the only category that ended beyond the uniform line during all of the three time intervals. The results for the transitions to PL land are represented in Figure 14. PL land gained mainly from PE and E land in terms of annual transition size. However, the bar for PE land transition intensity extended before the uniform line, which revealed that PE's start size was more crucial than PE's transition intensity to account for the area of the transition. Both the transition size and the transition intensity of E land to PL land were large, except that the gain of PL land avoided E land during 2010 and 2015. The gain of PL targeted PE land during 2010-2015. The quantity changes of E land were net loss during all the intervals. The results for the loss of E land are shown in Figure 15. The transition to PE land was much larger than to other categories during the first time interval. In terms of intensity, we see that the bar for EP and PL land extended to the right of the uniform line during all the intervals. Concerning the status of targeting and avoiding with regard to E land loss, all categories were stationary across the three periods. The results for the transitions to PL land are represented in Figure 14. PL land gained mainly from PE and E land in terms of annual transition size. However, the bar for PE land transition intensity extended before the uniform line, which revealed that PE's start size was more crucial than PE's transition intensity to account for the area of the transition. Both the transition size and the transition intensity of E land to PL land were large, except that the gain of PL land avoided E land during 2010 and 2015. The gain of PL targeted PE land during 2010-2015.
area. However, the fact that the LP land gained more intensively from PL land than from PE land during the first two time intervals is revealed in Figure 13(b). The reason is that the initial size of PE land was larger than PL land. Then the situation changed for the interval of 2010-2015, during which the bar for PL land was to the left of the uniform line. PE land was the only category that ended beyond the uniform line during all of the three time intervals. The results for the transitions to PL land are represented in Figure 14. PL land gained mainly from PE and E land in terms of annual transition size. However, the bar for PE land transition intensity extended before the uniform line, which revealed that PE's start size was more crucial than PE's transition intensity to account for the area of the transition. Both the transition size and the transition intensity of E land to PL land were large, except that the gain of PL land avoided E land during 2010 and 2015. The gain of PL targeted PE land during 2010-2015. The quantity changes of E land were net loss during all the intervals. The results for the loss of E land are shown in Figure 15. The transition to PE land was much larger than to other categories during the first time interval. In terms of intensity, we see that the bar for EP and PL land extended to the right of the uniform line during all the intervals. Concerning the status of targeting and avoiding with regard to E land loss, all categories were stationary across the three periods. The quantity changes of E land were net loss during all the intervals. The results for the loss of E land are shown in Figure 15. The transition to PE land was much larger than to other categories during the first time interval. In terms of intensity, we see that the bar for EP and PL land extended to the right of the uniform line during all the intervals. Concerning the status of targeting and avoiding with regard to E land loss, all categories were stationary across the three periods. Sustainability 2020, 12, x FOR PEER REVIEW 13 of 17

Discussion
For all of the three time intervals, the largest quantity component of change among the five categories was for LP, which presented net gain. The second was for PE, which presented net loss. The further categorical results show that the largest gain of LP land was from PE land, and the largest loss of PE land was to LP land in terms of size, and to PL land in terms of intensity. The findings provide evidences to justify the concern of cropland loss because of explosive urban sprawl. The conversion of PE to LP and PL land was likely due to the urbanization, i.e., the transitions from cropland to urban area and industrial land, which were common in the study area.
Since 2000, the government had realized the issue of food insecurity due to the loss of cropland caused by rapid urbanization. The cropland protection policies were initiated due to concern for food security [46]. The cropland requisition-compensation balance policy requires a quantity balance between the losses in cropland that transfer to built-up land and the gains of cropland elsewhere [47]. In this study, we found that the gain of PE land was mainly from E land, and large size of E land transferred to PE land especially during the first interval. The potential negative effects of the cropland requisition-compensation balance policy on ecological environment deserve attention from the policymakers and public, especially in ecologically fragile regions. Meanwhile, the loss intensity from E land to PL land was the largest. This finding implies that the tradeoff and synergy between industrial production and ecosystem regulating and supporting services should be evaluated and considered in the land management.
The PLE land change was the fastest during 2000-2005; then, the change rate decreased in the latter two periods. The results suggest that the PLE land change rate in Shandong province decelerated steadily from 2000 to 2015. A few reasons for this are that the accumulation of negative effects on the environment, the measures the government has put forward to balance the positive effect on the economy, and the negative effect on food production have caused public concern. Land use policies and management including land consolidation, urban growth boundaries, and an ecological civilization program had been implemented to improve sustainability [48]. Another underlying reason is that the growth rate of GDP declined after 2010.
This study also has some limitations. First, there are some uncertainties in the PLE land classification system. The classification systems and the concepts for PLE land have not been unified until now. There are still some controversies in linking the PLE classification systems with the remotely sensed land cover classification systems. For example, the unused lands include the desert with vegetation cover lower than 5%, the saline and alkaline land, and bare rock, which are not suitable for production and living. Whether the unused lands should be included in the ecological land is indefinable. In this study, we classified unused lands as E land considering their ecological functions. Further research into the classification system of PLE land using scientific methods and techniques is needed. Second, the spatial resolution of the land use data used in this study is 30 × 30

Discussion
For all of the three time intervals, the largest quantity component of change among the five categories was for LP, which presented net gain. The second was for PE, which presented net loss. The further categorical results show that the largest gain of LP land was from PE land, and the largest loss of PE land was to LP land in terms of size, and to PL land in terms of intensity. The findings provide evidences to justify the concern of cropland loss because of explosive urban sprawl. The conversion of PE to LP and PL land was likely due to the urbanization, i.e., the transitions from cropland to urban area and industrial land, which were common in the study area.
Since 2000, the government had realized the issue of food insecurity due to the loss of cropland caused by rapid urbanization. The cropland protection policies were initiated due to concern for food security [46]. The cropland requisition-compensation balance policy requires a quantity balance between the losses in cropland that transfer to built-up land and the gains of cropland elsewhere [47]. In this study, we found that the gain of PE land was mainly from E land, and large size of E land transferred to PE land especially during the first interval. The potential negative effects of the cropland requisition-compensation balance policy on ecological environment deserve attention from the policymakers and public, especially in ecologically fragile regions. Meanwhile, the loss intensity from E land to PL land was the largest. This finding implies that the tradeoff and synergy between industrial production and ecosystem regulating and supporting services should be evaluated and considered in the land management.
The PLE land change was the fastest during 2000-2005; then, the change rate decreased in the latter two periods. The results suggest that the PLE land change rate in Shandong province decelerated steadily from 2000 to 2015. A few reasons for this are that the accumulation of negative effects on the environment, the measures the government has put forward to balance the positive effect on the economy, and the negative effect on food production have caused public concern. Land use policies and management including land consolidation, urban growth boundaries, and an ecological civilization program had been implemented to improve sustainability [48]. Another underlying reason is that the growth rate of GDP declined after 2010.
This study also has some limitations. First, there are some uncertainties in the PLE land classification system. The classification systems and the concepts for PLE land have not been unified until now. There are still some controversies in linking the PLE classification systems with the remotely sensed land cover classification systems. For example, the unused lands include the desert with vegetation cover lower than 5%, the saline and alkaline land, and bare rock, which are not suitable for production and living. Whether the unused lands should be included in the ecological land is indefinable. In this study, we classified unused lands as E land considering their ecological functions. Further research into the classification system of PLE land using scientific methods and techniques is needed. Second, the spatial resolution of the land use data used in this study is 30 × 30 m. The resolution is higher than many other studies of similar space scale. However, some typical ecological land was still ignored because their areas are relatively small. For example, some urban ecological land may be included in LP land, and some farm protection ecological land may be included in PE land, etc. To reduce error caused by mixed pixel, higher resolution land use data should be used. In this study, the study period is from 2000 to 2015. The reason is that the accuracy of the land use data before 2000 was too low to reveal the land use processes. In the future, more accurate land use data in 2020 and before 2000 would be derived using more advanced techniques. The time scale could then be refined.
This study investigated the historical PLE land changes; however, the projection of PLE land dynamics in the future is also a key issue for optimal allocation. A multitude of land use models have been developed for scenario projections, like CLUMondo [49], LandCA [50], LANDSCAPE [51], FLUS [52], and so forth. The next step would be to predict PLE land distribution under different scenarios based on land use dynamic models. The driving factors of PLE land change should also be analyzed quantitatively using statistical or machining learning methods in the future.

Conclusions
To sum up, this study presents a full picture of PLE land changes in Shandong province between 2000 and 2015. We first defined the PLE land classification system based on previous studies and the local characteristics of the study area. Then we assessed spatiotemporal change in PLE land by spatial trajectory analysis and intensity analysis. Overall, the results show that the PLE land, especially the PE and LP land, experienced notable changes during the studied period. In terms of the overall change level, the PLE land change was fastest during the first time interval 2000-2005, and was the slowest during 2010-2015. The quantity component is the largest among the three components during all intervals. At the categorical level, the size of PE land loss was the largest, followed by the gain of LP land. In terms of intensity, active transitions include the gains of LP and PL land during all intervals and the losses of E land during the first two intervals. The transition level results reveal that LP land's gain targeted PE land during all intervals. The largest transition to PL was also from PE land. Large size of E land transferred to PE land during the first interval, and E land's loss targeted PL land in all intervals. As noted above, policy making for PLE land coordination should be put forward. For example, the land use policies of the study area should take into account the protection of PE land. The appropriate restriction of LP land expansion would be helpful in preventing the conflicts with PLE land. Meanwhile, there should be concerted effort by government to supervise for effective enforcement of policies that protect E land from being occupied by PE land and PL land. It is essential to coordinate the PLE land for socioeconomic development, food provision, and ecological protection and abide by the rules of food security and ecology. The methods used in this research could be applied in the studies of other provinces or regions easily, however, it would be necessary to modify the PLE land classification principles according to the local conditions.