Greening rate in North Korea doubles South Korea

South and North Korea have experienced contrasting economic developments since the 1950s while sharing similar climates. Previous studies revealed an overall long-term greening trend across the Korean Peninsula based on greenness data from coarse-resolution satellite images. However, there has been no comprehensive comparison of the greenness patterns and the driving mechanisms between two countries due to the limitations of coarse-resolution satellite data. Here, we performed cross-calibration among Landsat sensors and adopted a phenology-based approach to generate Landsat annual maximum Normalized Difference Vegetation Index (NDVImax) time series for each pixel from 1986 to 2017. We found that over 1986–2017, the greening rate in North Korea was almost twice that of South Korea. Cropland in South Korea is the main source of the greening discrepancy. The expansion of agricultural facilities in the stable cropland area and urbanization in the cropland loss area of South Korea contributed 57% to the significant negative NDVImax trend, which was dominant over the forest NDVImax increase resulting from rising temperatures, CO2 fertilization effects and afforestation projects. However, in North Korea, CO2 fertilization effects in the stable cropland area and transition from grassland to cropland promoted an increase in NDVImax, despite decreasing NDVImax in forest areas due to deforestation. Our results highlight the need for delineating fine-scale land-use changes to advance our understanding of regional vegetation dynamics.


Introduction
The Northern Hemisphere has been warming rapidly since the last century (Lucht et al 2002, Slayback et al 2003. Warmer temperatures coupled with the fertilization effects of rising CO 2 (Wang et al 2021) and intensive land-use practices (Song et al 2018 are expected to alter vegetation composition and structure, contributing to increased vegetation productivity across the north-temperate zone. This enhancement of vegetation growth, reflected by an increase in vegetation greenness as seen from space (Ju and Masek 2016), is important in the global carbon cycle through a series of biogeochemical processes (Piao et al 2015). Vegetation greenness monitoring that relies on satellite remote sensing time series has reported diverse greening patterns and regimes among different regions (Chen et al 2014, Myers-Smith et al 2020, Choler et al 2021 across the northern hemisphere. With the increasing diversity of human practices (Song et al 2018), the study of greening in important ecological regions should be based on the unique trajectory of regional vegetation change. As an irreplaceable ecological barrier between China, Russia, and Japan, the Korean Peninsula, located in the center of Northeast Asia, has exhibited a pervasive greening trend by the coarse-resolution NDVI dataset Miller-Rushing 2014, Yun et al 2020). However, the currently widely used coarseresolution (e.g. 8 km) satellite observations limit our ability to further quantify and explain the vegetation greenness within the Korean Peninsula.
South and North Korea, which together occupy the Korean Peninsula, share similar climatic conditions (IPCC 2019) but exhibit diverse land use patterns due to different economic development models that have existed since the end of Korean War. After the 1990s when North Korea entered the 'arduous March' , natural forests continued to be cut down and reclaimed into farmland to alleviate severe fuel and food shortages, thus triggering extensive soil erosion and repeated flooding (Lee and Miller-Rushing 2014). Although North Korea has focused on land protection since the 2010s (Kim et al 2021), deforestation and forest degradation remain severe (Jin et al 2016). In contrast, in South Korea, 10 years afforestation programs have been launched by the national forestry department since the 1970s (Yun et al 2020), and this has increased forest cover and enhanced forest carbon sequestration. Still, the problems posed by rapid urbanization remain, where the expansion of impervious land surface degrades ecosystems and reduces the availability of arable land (Hwang et al 2022). These contrasting anthropogenic activities are intertwined with climate change and have resulted in widespread changes (Lee and Miller-Rushing 2014) in North and South Korea's ecosystems. However, the contribution of divergent land-use change (LUC) trajectories and climate disturbances to vegetation greenness in these two countries remains unclear.
Although dynamic vegetation models and satellite-derived models reveal that CO 2 (Zhu et al 2016), LUCs (Hansen et al 2013, Song et al 2018, temperature (Wu et al 2022) and snow cover changes (Piao et al 2003, Berner et al 2020 in the Northern Hemisphere act as the primary drivers of the greening trend, there are still many uncertainties in greening attribution. First, coarse-resolution vegetation imaging technologies that are widely used such as the Advanced Very High-Resolution Radiometer NDVI time series lack the spatial resolution (e.g. 8 km 2 ) to distinguish fine-scale changes in vegetation composition (McManus et al 2012, Ju and Masek 2016, Berner et al 2020, and this hinders the quantification and attribution of greening heterogeneity within and among different vegetation types (Sulla-Menashe et al 2018). Second, previous models' simulations lacking long-term high-resolution LUC information were unable to explicitly represent underlying complex land-use practices patterns , Liu et al 2019b, which are prone to underestimating the impacts of intensive and fragmented human practice (Zeng et al 2018) (e.g. small patches of land clearing and deforestation occurring at 10 s of m 2 ). Finally, rapid ecosystem changes need to be separated and accounted for in attribution approaches (Burrell et al 2020). Mixing stable vegetation communities with disturbed communities will tend to ignore the changes in ecological processes such as recovery from disturbance and increase attribution uncertainties (Sulla-Menashe et al 2018).
Here, we used the annual Landsat maximum growing season (May-September) NDVI (NDVI max ) to characterize greenness disparity and driving mechanisms between North and South Korea. We developed 30 m resolution annual land-cover classification maps to isolate land cover types and distinguish stable (no land cover change during the study period) areas from disturbed (land cover changes have occurred during the study period) areas. On this basis, we calculated the annual time series of NDVI max for each pixel from 1986 to 2017 (Berner et al 2020) and compared greenness change patterns among different land cover types in North and South Korea under stable and disturbed conditions. Moreover, we examined the extent to which greenness trends in North and South Korea were linked to climate variability and climate change (including temperature, precipitation, and radiation) and also to LUC using machine learning (see section 2).

Generating Landsat annual NDVI max time series
We developed a growing season maximum NDVI dataset annually (NDVI max ) from 1986 to 2017 over the entire Korean Peninsula using 30 m resolution NDVI data from Landsat satellites. All available Landsat 5, 7, and 8 (collection 1, level 1) surface reflectance images were processed on the GEE platform between 1 May and 30 September from 1986 to 2017. The cloud cover and cloud shadow pixels were excluded using the Quality Assessment band included in the Landsat collection 1 (Liu et al 2019a).
To correct the NDVI biases among the three Landsat sensors induced by the different spectral responses and radiometric resolution, we developed cross-calibration models using linear regression models that homogenized Landsat 5 and Landsat 8 to Landsat 7 (L7 NDVI = L5 NDVI × 0.9871 + 0.0403; L7 NDVI = L8 NDVI × 0.9317 + 0.0224) (see section 2.2 in section 2). The regression models possessed good fits (R 2 = 0.89 for harmonization between Landsat 5 and 7, R 2 = 0.87 for Landsat 7 and 8, with both P values < 0.0001, figure S1) and were therefore used for cross-calibrating annual NDVI max across the three Landsat satellites.
Additionally, unbalanced availability of satellite images before and after the 2000s may also introduce a spurious trend into the Landsat NDVI time series. To address the bias variability in the availability and timing of the clear Landsat observations per year, we estimated the annual Landsat NDVI max at the pixel level by modeling annual land surface phenology using a 17 years window (Berner et al 2020) of 1986-2017 considering the available observations were insufficient before 2000 at the pixel level, which forced us to pool enough observations April to October. Each phenological curve was developed by fitting a cubic spline using all available observations in a 17 years period. (B) Annual NDVImax estimation (blue point) derived from the median estimated NDVImax using all available growing season observations (red points) and phenological information (red curve). The vertical green lines represent the typical difference in NDVI between the specific day of each observation in the growing season and the day when the maximum NDVI was acquired. (C) Landsat image availability from 1986 to 2017 over the Korean Peninsula with cloud cover less than 80%. Before the year 1999 when only Landsat 5 operated, the number of acquired images was roughly less than 180. With the successive launches of Landsat 7 and Landsat 8 after 2000, the number of images with cloud cover less than 80% became increasingly available. to model phenological curves over the era of sparse observations.
We quantified land surface phenology for each grid cell in the extended growing season by fitting a cubic spline to all valid Landsat observations for each 17 years period from 1986 to 2017 (figure 1(A)). The phenological curves provided information regarding the typical differences in NDVI values between NDVI max and the day that each growing season observation was acquired ( figure 1(B)). By adjusting the individual growing season observations based on the timing of acquisition relative to the growing season maximum NDVI, the annual NDVI max for each grid cell was estimated (figure 2). This approach was first used to examine the greening trend in the Arctic tundra biome using random sampling sites across Arctic tundra (Berner et al 2020). More complete details regarding this method are provided in the supplementary material.

Developing land cover classification maps
We created high-resolution (30 m) annual landcover classification maps over 1986-2017 to isolate land cover types and to distinguish vegetation communities between stable and disturbed areas using imagery from the Landsat series of satellites (collection 1, level 1 surface reflectance images, including Landsats 5, 7, and 8).
To capture distinguishable spectral-temporal features for different land cover classes over the entire Korean Peninsula, we used the Continuous Change Detection and Classification (CCDC) algorithm (Zhu and Woodcock 2014) developed from all available cloud-screened Landsat images to construct the time series model for each 30 m pixel in the GEE platform. The spectral-temporal information from CCDC time-series models was acquired for distinguishing characteristic land cover types. We also prepared other ancillary data as potential inputs for classification. These ancillary data included the vegetation indices (i.e. NDVI, Enhanced Vegetation Index (EVI, Justice et al 1998), Soil-Adjusted Vegetation Index (Huete 1988)), Normalized Difference Water Index (McFeeters 1996), and topographic information such as elevation, slope, and aspect.
A total of 6729 samples that covers the entire Korean Peninsula were collected for training and validation (figure S2). We removed samples where  1986, 1987, 2001, 2016 and 2017. Prior to the 2000s, limited usable observations and extensive clouds and snow in the Korean Peninsula resulted in data deficiencies in some areas (black color). The phenology-based approach was adopted to estimate the NDVImax over the entire peninsula by developing the phenology curve. The whiskers of the box extended 1.5 times the interquartile range. the land-use type had changed from 1986 to 2017 identified by the CCDC algorithm to reduce classification errors caused by inconsistent spectral characteristics of the samples before and after the change in land use types. We then performed visual interpretation of each sample based on Landsat and high-resolution Google Earth imagery. We interpreted the samples as one of ten categories representing dominant land cover types in the study area, including cropland, deciduous forest, evergreen forest, grassland, mixed forest, shrubland, sparse vegetation, urban, water, and wetland (table S1). We then employed Random Forest (Breiman 2001), a non-parametric machine-learning classifier, to classify land cover types. To train the machine learning classifier, we randomly held out a selection (30%) of the land cover samples as the independent validation dataset. Random Forest classifier achieved a higher accuracy with an overall accuracy of 94.7% and a kappa coefficient of 0.929 when using the combination of spectral-temporal features derived from the CCDC algorithm, vegetation indexes, and topography information as the input features (table 1). We used the validation samples to create a confusion matrix to further assess the accuracy of the land cover classification results (table 2).

NDVI max trend attribution
To investigate the underlying drivers of the NDVI max change in the Korean Peninsula, we developed Random Forest (Breiman 2001), a machine learning algorithm, to predict the trend classes (i.e. greening, browning, and no significant trend) of Landsat NDVI max from 1986 to 2017. In this study, we used the time-series data of CO 2 concentration, climate factors (e.g. temperature, precipitation, and solar radiation) (see meteorological data in the supplementary material), and LUC drivers calculated above as the explanatory variables and NDVI max trend classes as the variables to be predicted. To capture both the contribution of climate trend and climate variability to vegetation greening (Qian et al 2019), we split the climatology factors (e.g. the per-pixel temperature, precipitation, and solar radiation) into trend and interannual variability components (Burrell et al 2020) and used them as independent climatic explanatory variables for vegetation greening prediction. We applied a 20 years moving average smoothing window (Burrell et al 2020) to remove the interannual climate variability, and the smoothed climate time series was preserved to determine the vegetation change driven by climate trend. We then used the smoothed data to detrend the originally observed climatology, and the variable of climate variability was thus determined by the difference between annually observed climatology and smoothed climatology. We then employed pair-wise correlations (Berner et al 2020) and excluded variables with high average absolute correlation to filter highly correlated variables (coefficients >0.75). In this manner, seven explanatory variables were preserved and used to establish the attribution model, including CO 2 concentration, temperature trend, temperature variability, precipitation trend, precipitation variability, radiation trend, and LUC.
To quantify the contribution of each explanatory variable, we calculated the Feature Importance based on the Mean Decrease Accuracy (Leroux et al 2017, Shi et al 2020. We then randomly partitioned the pixels of the study area into the training set (70%) and validation set (30%) to assess the performance of Random Forest model. To reduce Random Forest model uncertainties, the estimation steps were repeated 20 times, and the mean value of the determination coefficient (R 2 ) and standard deviation (SD) of model performance were calculated. The Random Forest model accounted for over 65% of the vegetation greenness in the Korean Peninsula. The entire construction and evaluation of the Random Forest model were implemented using the random-Forest (Liaw and Wiener 2002) and caret (Kuhn 2008) packages in R. Methods including calculation of Landsat NDVI max trends, partial correlation analysis and other datasets used in our study were described in the supplementary material.

Greening differences between North and South Korea
Our analysis of the Landsat NDVI max indicated significant greenness trends (Mann-Kendall test, P < 0.05) in greater than 44% of the Korean Peninsula over 1986-2017. Most (70%) of the significant trends were greening (significant positive), covering 31% of the Korean Peninsula, while 13% of the Korean Peninsula exhibited browning (significant negative) NDVI max trends (figures 3(A) and (D)). Although widespread greening was observed in the Korean Peninsula, there were differences in the magnitude and pattern of greenness changes between North and South Korea. The increase rate of annual mean NDVI max in North Korea was 0.008 per decade with the mean NDVI max of 0.792 (linear regression, p < 0.05, n = 32) over 1986-2017, which was nearly twice that of South Korea (0.005 per decade, mean NDVI max was 0.764) ( figure 3(C)). The averaged significant trend of North Korea (0.016 per decade) was also found to be twice the significant trend of South Korea (0.009 per decade). This greening trend difference between North and South Korea was confirmed by the temporal trend of advanced very high resolution radiometer (AVHRR) (0.031 ± 0.003 per decade in North Korea and 0.022 ± 0.004 per decade in South Korea Table 2. Confusion matrix of the land cover classification results was generated from the Random Forest classifier. The overall accuracy of the land cover maps with ten categories was 94.7%, and the largest confusion was between deciduous forest and mixed forest, and grassland and sparse vegetation.  In South Korea, both greening and browning areas were more extensive and accounted for 33% and 13%, respectively. In contrast, greening (29%) and browning (9%) have been observed in North Korea to a smaller extent ( figure 3(D)).

Greening discrepancies and their linkage between land-use types
To investigate which land cover type explained the greening discrepancy between North and South Korea, we compared the direction and magnitude of the NDVI max change in all stable and disturbed land cover types between North and South Korea based on the land cover classification maps (figures 4(A) and (B)). Our results revealed that the stable cropland in South Korea, which accounts for 17.5% of the whole nation ( figure 4(A)), was the most important factor for the lower NDVI max growth trend in South Korea (figure 5), with the browning trend achieving −0.012 ± 0.004 per decade ( figure 6(A)). In addition, the disturbed cropland in South Korea also exhibited significant negative NDVI max trends (−0.003 ± 0.002 per decade for cropland gain and −0.009 ± 0.004 per decade for cropland loss area) (figures 6(B) and (C)). Croplands in stable and disturbed areas together contributed 57% to the significant negative NDVI max trend in South Korea, and this was far higher than that contributed by the other land cover types (figures 5(B) and (C)). In contrast, croplands in North Korea exhibited a generally upward NDVI max trend in both stable and disturbed areas. In total, 71% (12.3 × 10 3 km 2 ) of the cropland in North Korea presented a greening trend, and of this, approximately three-quarters was identified as stable cropland (8.1 × 10 3 km 2 ). Contrary to the downward trend of NDVI max in croplands in South Korea, NDVI max in both stable (0.012 ± 0.005 per decade) and disturbed forest areas (0.014 ± 0.003 per decade for forest gain and 0.011 ± 0.004 per decade for forest loss area) in South Korea exhibited significant upward trends, with faster increase rates than those in North Korea (figures 6(D)-(F)). Especially in stable evergreen forest areas in South Korea, the NDVI max trend reached 0.02 per decade (figure S5). Forest area in South Korea was the largest contributor to greening in South Korea, which accounted for 63% of the positive trend over the whole country and exceeded the contribution of forests to the greening trend in North Korea (46%) (figures 5(B) and (C)). Compared to South Korea where only 14% of forest areas exhibited browning trends, the forests in North Korea experienced much stronger and more extensive browning (29%, 4.7 × 10 3 km 2 ).

Factors that drive greening discrepancies
We further explored the linkage between greening trend discrepancies and driving factors by constructing Random Forest model. The pixel-level results indicated that the driving mechanism differences between North and South Korea primarily manifested as discrepancies in the contributions of temperature variables and LUC to NDVI max change. Temperature variables, particularly temperature trends, explained approximately 33% of the greenness change in South Korea which surpassed that in North Korea (27%). However, the contribution of LUC to greenness change in North Korea was 8%, and this is larger than that in South Korea (3%), although the extent of land cover change in the two regions was similar, with proportions of 18% and 19%, respectively (figures 4(A) and (B)). The differences in the contribution of other variables to greenness change between North and South Korea were within 2%.
We then quantified the differences in the contribution of dominant factors to greening between North and South Korea. We found that temperature had dominant contributions over 77% of the greening pixels in South Korea, with high contribution values mainly located in the stable deciduous forest in the northern part and newly added evergreen forest in the southern part of South Korea (table 3, figures 7 and 4). Compared to South Korea, temperature as the dominant contributor occupied a smaller area (55%) in greening pixels of North Korea, followed by precipitation (21%) and radiation (12%) (figures 7(A)-(G)). The dependence of stable forests on temperature variability was stronger in South Korea than it was in North Korea (figures 8 and S6). However, the positive correlation between NDVI max and increasing atmospheric CO 2 concentration in the stable cropland of North Korea exceeded that in South Korea.

Effects of LUC on the greening discrepancy
We further compared the NDVI max trend in stable and disturbed areas between South and North Korea to investigate the differences in the specific driving mechanisms of LUC. The NDVI max trends in both stable and disturbed areas of North Korea were positive (0.009 ± 0.005 and 0.006 ± 0.003 per decade for stable and disturbed areas, respectively). In contrast, the disturbed areas in South Korea exhibited a slightly downward trend in NDVI max (figures 4(C) and (D)), regardless of the effects of climate variables and CO 2 .
The difference in NDVI max trends in disturbed areas between North and South Korea can be explained by quantifying the role of LUC in North and South Korea. LUC in South Korea provides dominant    contributions to over 4% of NDVI max greening, and this is half as much area as NDVI max browning (9%) (figure 7(H)). By contrast, in North Korea, LUC was the dominant driver of greening and browning trends over equal areas (each 9%). Land cover transformation paths in North and South Korea further elucidate the contribution disparities in greenness change.
Over the last three decades in North Korea, vegetation conversion from grassland to cropland (11%) and the forest (16%) have partly compensated for the impact of deforestation (40%) (table S3; figure 6), leading to an upward NDVI max trend in disturbed areas of North Korea. While in South Korea, cropland loss was dominant over the contribution of afforestation (grassland and shrubland transferred into forest) to vegetation greening (table S4), exhibiting significant negative NDVI max (−0.009 per decade) due to urbanization (figure 6(C)).

Discussion
Our study revealed that cropland was the primary contributor to the lower greening rate in South Korea, which exhibited significant browning trends in both disturbed and stable areas of South Korea. The expansion of urban impervious surface into croplands was the primary reason for cropland browning (Song et al 2018). This is consistent with a recent study which reported that urban expansion encroaching on croplands accounted for 72% of urban land growth in Asia (Liu et al 2019b). Over three decades, the urban area in South Korea has increased by 56% to reach 5.6 × 10 3 km 2 (table S4). Although urban areas only accounted for approximately 5% of South Korea, they contributed 16% of the vegetation browning (figure 5). In contrast, urban expansion areas in North Korea over the 30 years were small (table S3), which exerted a limited impact on greenness change. In addition to the disturbed cropland, the significant browning trend observed in stable croplands in South Korea is also worthy of attention. This could be explained by the expansion of agricultural facilities, especially the construction of greenhouses and agricultural infrastructure (figure 9). In the period of the 'White Revolution' (Boo et al 2013) that aimed at innovating and modernizing agricultural structure in South Korea, greenhouse achieved rapid growth in the cultivation area (increased 331%, from 1.3% of cultivated land in 1986 to 5.6% in 2000, figure 9(J)). This is consistent with the browning trends observed in South Korea's stable cropland ( figure 6). In addition, nationwide agricultural infrastructure, including agricultural housing facilities (Boo et al 2013), drainage and irrigation facilities, and livestock facilities has been extensively developed and expanded, with the proportion increasing from 2.9% of the cultivated land in 1993 to 5.1% in 2010 (figure 9(J)). It is noteworthy that although these greenhouses and agricultural infrastructure scattered in the cropland landscapes can be visually recognized in Google Earth high-resolution images (<5 m) (Wang et al 2022) (figure 9), they were undetectable by our 30 m Landsat-based land cover classification maps and were still classified as cropland, because many of these fragmented patches were smaller than the Landsat pixel scale. However, these scattered greenhouse and agricultural infrastructure could conspicuously reduce NDVI max levels (Hwang et al 2022). As shown in figure 9(D), after the construction of greenhouses, the maximum NDVI for pixels dropped from 0.8 to ∼0.3. On the other hand, traditional farming without greenhouses is still dominant in North Korea (Lee and Miller-Rushing 2014). Although North Korea has embarked on agricultural structural reforms in recent years (Kim et al 2021), their influences on vegetation greenness in croplands have not appeared clearly. It is worth noting that although the area of cropland exhibited a downward trend, the crop yield in South Korea was on a significant upward trend (figure S7), which was largely due to the fact that South Korea has benefited from advances in modern agricultural technology.
We found that the greater greening rate of forests in South Korea partially alleviated the discrepancy in the greening trend between North and South Korea. The NDVI max growth rate in the forest gain area of South Korea reached 0.014 per decade (figure 6(E)), which was higher than that in North Korea (0.01 per decade). We further revealed that forest age composition and management played important roles in NDVI max trends. In South Korea, the time period of forest gain was primarily concentrated prior to 2000, with the greening trend in most years before 2000 over 0.01 per decade. Forest age in disturbed forest areas in South Korea exhibited a positive correlation (r = 0.48, P < 0.05) with the NDVI max trend (figure 10). This implies that the greening rate was higher in earlier forest-gained areas which have grown for several decades. Though it is unclear how long it could increase NDVI max in the next decades as mature forests might not continuously exhibit increasing NDVI max . On the other hand, we found that North Korea did not show any significant correlation between forest age and NDVI max trends (figure 10). The discrepancy between South and North Korea may be related to forest management that has been applied in South Korea actively such as forest age, species, and composition managements (Engler et al 2014, Dong et al 2020. Although high-resolution NDVI max and land use classification maps were employed in our greening attribution research, it is noteworthy that other biological and physical factors, including nitrogen deposition (Zhu et al 2016), wildfire (McManus Kelly et al 2012), and non-climatic factors, particularly natural vegetation succession (Kennedy et al 2015), should be given further consideration to deepen our understanding of regional vegetation changes and climate feedbacks.

Conclusion
This study advanced our understanding of greening discrepancy and the driving mechanism between South and North Korea using the annual Landsat NDVI max time series. Our study revealed that the greening rate in North Korea almost doubles South Korea over 1986-2017. By developing high-resolution annual land use classification maps, we were able to show that the difference in vegetation greening rate between South and North Korea was linked to specific land use types. Cropland in South Korea, characterized by expansion of agricultural facilities in the stable cropland area and urbanization in the cropland loss area, was the main source of lower greening rate, although forest areas in South Korea showed a higher greening trend, benefiting from afforestation projects, forest management, and greater positive effects of temperature. By contrast, in North Korea, the expansion of cropland areas dominated over the pressure of deforestation, which led to an upward greening trend. Our research highlights that besides climate change, agricultural practices, forest management, and other LUCs, all played important roles in regional vegetation greenness change. Those complex components and their interactions should be considered in the regional carbon budget estimates.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// zenodo.org/search?page=1&size=20&q=7213038.

Data and code availability
The annual Landsat NDVI max dataset, land cover maps, and CO 2 datasets described in this study are available from the author through e-mail request. TerraClimate datasets are available on the Google Earth Engine platform (https://developers. google.com/earth-engine/datasets/catalog/IDAHO_ EPSCOR_TERRACLIMATE). Greenhouse proportion data, agricultural infrastructure area data, and annual cultivated land utilization ratio are publicly available at the Korean Statistical Information Service (KOSIS, https://kosis.kr/index/index.do) and the FAOSTAT database (FAO, Food and Agriculture Organization of the United Nations, www.fao.org/ faostat/en/#data/RL). All code used for preprocessing the data and creating the plot in this analysis is freely available on request from the author.