Delineation of Rice Productivity Projected via Integration of a Crop Model with Geostationary Satellite Imagery in North Korea

Abstract

Satellite images can be integrated into a crop model to strengthen the advantages of each technique for crop monitoring and to compensate for weaknesses of each other, which can be systematically applied for monitoring inaccessible croplands. The objective of this study was to outline the productivity of paddy rice based on simulation of the yield of all paddy fields in North Korea, using a grid crop model combined with optical satellite imagery. The grid GRAMI-rice model was used to simulate paddy rice yields for inaccessible North Korea based on the bidirectional reflectance distribution function-adjusted vegetation indices (VIs) and the solar insolation. VIs and solar insolation for the model simulation were obtained from the Geostationary Ocean Color Imager (GOCI) and the Meteorological Imager (MI) sensors of the Communication Ocean and Meteorological Satellite (COMS). Reanalysis data of air temperature were achieved from the Korea Local Analysis and Prediction System (KLAPS). Study results showed that the yields of paddy rice were reproduced with a statistically significant range of accuracy. The regional characteristics of crops for all of the sites in North Korea were successfully defined into four clusters through a spatial analysis using the K-means clustering approach. The current study has demonstrated the potential effectiveness of characterization of crop productivity based on incorporation of a crop model with satellite images, which is a proven consistent technique for monitoring of crop productivity in inaccessible regions.

1. Introduction

North Korea comprehensively depends on imports of staple crops primarily from China and the Russian Federation according to the report by The Food and Agriculture Organization, FAO and the World Food Program, WFP (FAO, 2011). The food shortage in the nation evolved into the famine in the mid-1990s and continued in the present day. Multiple interacting factors are attributable to the famine, including agricultural policy influenced by the political and economic situation, as well as weather conditions (Ireson, 2006). Global warming accompanied by abnormal and severe weather conditions has mainly affected the environments of the agricultural systems all over the world. Potential impacts of climate change on naturalresources have been specified that agriculture would predominantly sufferfrom a deficiency in water resources (Ko et al., 2012; Sowers et al., 2011). The degree of climate change in EastAsia is assumed more severe than the global trends (NIER, 2014; Schäfer, 2015). As North Korea is one of the most vulnerable countriesin thisregion, the agricultural environment in this country likely gets worse further due to the rapidly changing climate and unpredictable and severe weather conditions.

It was conveyed that the lack of advice and unsustainable national rural policies would cause the agricultural collapse in North Korea (Ireson, 2006). This deficiency in agricultural policy planning may adversely affect the future of the agricultural industry and cause a decline in gross agricultural production. A practical scheme to allow consistent monitoring and precise estimation ofstaple crop yieldsin North Korea is vital to help determine appropriate solutions to the food shortage issues for the grain aid organizations (Yeom et al., 2018). Precise and prompt assessment of the crop damage caused by natural disasters(e.g., crop defoliation, drought, and pest infestation)should benefit from well-made strategic design measures employed to meet the crop production demands (Doraiswamy et al., 2004). Nonetheless, it is challenging to secure immediate agricultural information on primary crop productivity due to the inaccessibility of the data sources for crop fields, which has stemmed from the past and current political issues in North Korea.

An assimilation technique of crop modeling with remote sensing information is a usefulscientific scheme for effective evaluation of crop growth and yield (Maas, 1992). While remote sensing has the advantage to observe crop growth conditions on time and with spatial variation, a crop model is a useful tool to describe the sequential crop growth and development during a growing season.Crop simulation providesfor pivotal tools for the continuous description of crop growth and their development by simulating the biophysical processes in the soil-plant-atmospheric system, offering significant positive foresight to support agricultural precision. However, crop modeling alone is an inadequate decision-support tool for generating general crop characteristics (Ko et al., 2015) and for estimating the regional crop yield when no spatial biophysical parameters are available.Most conventional crop models may be inadequate to simulate the specific spatial variations in seasonal crop growth.

Satellite sensor-based remote sensing can provide real-time crop growth information in considerable details to address the issue mentioned above. Satellite images can be useful information to observe crop growth to support on-site ground truth studies (Jones and Vaughan, 2010). In addition, a substantial benefit of remote sensing is that the spatial data can be achieved for any region on the Earth’s surface using satellite sources where the satellite track is reachable. However, a significant limitation ofthe opticalsatellite sensor-based images is the unfavorable atmospheric condition thatrestrictsthe timely retrieval of crop data, which is predominantly caused by cloud coverslinked to the periodic rainy monsoon season in many parts of Asia (Funk and Budde, 2009; Xiao et al., 2006). Utilizing the advantage of crop modeling and remote sensing techniques of each other could assist in negating their respective weaknesses. Crop models integrated with remote sensing imagesfromoperational satellites have been adopted to evaluate crop conditions and crop productivity at various geographical scales (Doraiswamy et al., 2005; Huang et al., 2016;Jeong et al., 2018; Kim et al., 2017; Zhao et al., 2013). A few different calibration procedures were applied to agree with model simulation and observation from remote sensing.

In this study, we used the GRAMI-rice model (Ko et al., 2015) that was further developed from a gramineous crop model, GRAMI (Maas, 1992) that employs remote sensing information for crop growth monitoring and yield mapping projects. GRAMI simulations of crop growth and yield using satellite data should carefully consider the spatiotemporal observation points of biophysical parameters to achieve reliable performance results for the region of interest. Temporal point observations alone are not adequate to reflect the actual crop phenology using indirect satellite indicators such as normalized difference vegetation index (NDVI). Among the biophysical and meteorological parameters of the crop model, vegetation indices(VIs)information isthe most critical factor to reflect actual crop growth and development, when compared to the meteorological elements such as solar insolation and air temperature. A continuous profile of the reflectance data of crop canopy is used to observe the current growth status of plants through a vegetation index indirectly (Hilker et al., 2008). Therefore, plant canopy is more sensitive to crop simulation performances in comparison with the other variables.

The GRAMI model is formulated with a ‘withinseason’ calibration procedure that agrees between simulation and observation based on a set of iterative mathematical processes(Maas, 1993).This calibration procedure systematically manipulates the initial conditions and parameters ofthe model that affect crop growth and development, in ordertominimize the error between simulation and observation from remote sensing. The ‘within-season’ calibration procedure allows infrequent observations to apply to calibrate the model. The optimization method formulated in GRAMIis a function to fill orsmooth noise and sparse remote sensing observations. Likewise, it isfeasible to use varioustemporalsmoothing techniquesto fill gaps and smooth noises in time-series remote sensing data, such as Fourier harmonics, threshold methods, and curve-fitting approaches (Atkinson et al., 2012; Chen et al., 2006;Jonsson and Eklundh, 2002;Roerink et al., 2000).

The rainy monsoon season during the summer influencesthe atmosphere in North Korea likewise the other East Asian nations (Qian et al., 2002). Thus, obtaining crop canopy data during the summer growing season is mostly hindered in the case of using polar orbit satellites such as the Moderate Resolution Imaging Spectroradiometer(MODIS)satellite (https:// modis.gsfc.nasa.gov, accessed on Dec. 10, 2018) despite the availability of Terra and Aqua sensors for improving temporal resolution. Yeom et al. (2015) introduced a potential solution to resolve this issue. They employed remotely sensed data from a high sequential resolution satellite to reproduce rice yield using the GRAMI-rice model. The satellite data were obtained explicitly from Geostationary Ocean Color Imager (GOCI) and Meteorological Imager (MI) of COMS (Choi et al., 2012).In the currentstudy, we also employed the GOCI-based NDVI profiles to monitor rice growth conditions and development in North Korea, in order to minimize the unfavorable weather conditions during the summer rainy season. It is also essential to secure spatial meteorological parameters of solar insolation and air temperature using satellite measurements and reanalysis data in order to reflect agricultural conditions to map crop productivity using the GRAMI model for inaccessible North Korea. As opposed toVIs,solarinsolation and airtemperature are physical variables that directly affect crop growth and development because they represent the energy source for photosynthesis and environmental control on morphogenesis of crops, respectively. Although there have been some efforts to use ground measurements to estimate the spatial distribution of biophysical parameters using a spatial interpolation method, this approach is restricted when considering the areas with sparse ground measurement networks. North Korea is such a region of the nation due to its complex terrain and inaccessibility. Therefore, crop growth and development must be simulated using biophysical parameters integrated with the remote sensing information and reanalysis data.

The objective of this study was to propose a delineation approach to monitor and determine the agronomic environment adequately in the unapproachable crop fields of North Korea based the GRAMI model integrated with biophysical and meteorological parametersfrom theCOMS GOCI and MI images. The implications of the study far reach potential to produce viable solutionsin the area of crop modeling for enhancing the agricultural quality and yield with modern day techniques, while addressing data sparsity constraints using satellite-based remotely sensed products integrated into the GRAMI model.

2. Study Region and Data

1) Characteristics of the study area

The monitoring of the agricultural system environmentsin North Korea was accomplished using a novel approach ofintegrating the GRAMImodel with satellite-based biophysical parameters. North Korea, or the Democratic People’s Republic of Korea (Fig. 1a), resides in East Asia and, bordering China, Russia, and South Korea. Agriculture in North Korea comprises ~ 25.2 % of the total gross domestic product (GDP), being the third most important economic factor (CIA, 2017).Arable lands embrace approximately 1.96 M ha (19.5 %) of the total area of the country (12.0 M ha), with comprising 0.32 M ha of perennial crops such as mulberry and fruit (CIA, 2017; Zhang et al., 2017). Rice is the most important staple crop in the country, followed by corn,soybean, and potato (FAO, 2015).In this study, we focused on rice because of the limited ground truth data of the other staple crops required to calibrate the model. Residing between 37°N and 43°N, North Korea has an extended winter period and a few frost-free days. Due to the severe and stretched winter environment, the single cropping system dominates the cultivation.

Fig. 1. A map of North Korea (a) with the provincial boundaries and geographic locations (presented in red) of Paju (b) and Cheorwon (c) considered for validation of the GRAMI model in this study in South Korea.

The prevailing climate in North Korea is humid continental while it experiences the impact of the combined climate atmosphere from the dry land and the ocean according to the Koppen-Geiger climate classification scheme (Peel et al., 2007). The summer istypically the hottest,most humid, and the wettest time of year as the southern and southeastern monsoon winds convey moist airfrom the Pacific Ocean.Almost 60 % of all rainfall is concentrated between June and September.Therefore, it obstructsto obtain the surface spectralreflectance data ofthe crop canopy due to longlasting clouds during the summer growing season.

2) Satellite data

Biophysical and meteorological parameters of VI and solar insolation were retrieved from the COMS data, with an aim to acquire spatiotemporal information on paddy productivity. COMS is a geostationary satellite developed by the Korea Aerospace Research Institute (KARI) and launched on June 27, 2010, stationed at a longitude of 128.2°E.COMS is mounted with two individual payloads of the GOCI and the MI sensors (Table 1). The core objectives of GOCI include detecting,monitoring, and predicting short-term biophysical phenomena; analyzing bio-geochemical variables and cycles; and detecting Yellow dust and land classified information. GOCI with the eight solar spectral wavebands observesthe Korean peninsula and its surrounding regions eight times a day. GOCI was mainly planned to monitor the ocean environments, using the spectral waveband of a Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) (Wang, 1999).

Table 1. The detailed characteristics of COMS GOCI, COMS MI and MODIS sensors for estimating crop productivity of North Korea

Meanwhile, the high temporal resolution of the observations and vegetation-sensitive spectral bands of GOCI allow adapting to land surface-based applications such as monitoring crop information (Yeom et al., 2012;Yeom and Kim, 2015).The MIsensor ofCOMS observes the near-real-time weather conditions using five dominant multispectral wavebandsfrom visible to infrared (IR) wavelengths. In the current study, the GOCI was used to obtain the cumulative VIssince the set-up is equipped with useful vegetation bands of Red and near infrared (NIR).As the IR channels of MI can be beneficialforinterpreting the intricate cloud effects, the solar insolation data were obtained from the MI.

3. Study Methods

1) GRAMI-rice crop model

The GRAMI-rice model, employed in this study, is formulated to assimilate remotely sensed information, to evaluate the potential crop productivity based on combining advantages of each technique of crop modeling and remote sensing (Ko et al., 2015). This model was initially designed to be able to receive remote sensing data as an input to execute the ‘withinseason’ calibration procedure by Maas (1993). In this procedure, the simulated vegetation indices (VIs) of crop canopies or LAI are compared with the corresponding observed values to allow agreement within a presetrange of errors. Four different parameters (L0, a, b, and c) are employed in the GRAMI-rice model to frame the mathematical procedures of crop growth. In the current study, the Bayesian method is adopted to obtain these parameters based on a prior distribution chosen according to the estimatesfrom the previousstudies.We assumed that a log-log regression model with a slope of approximately two-thirds (2/3) could describe the relationship between reflectance and leaf area index (LAI). Based on thistheory, the log-log linear regression models were employed to formulate the relationships between five VIs and the LAIs. The detailed description of this process is as follows. The empirical model below is employed for each VI, labeled l = 1, 2, 3, 4, and 5, respectively.

\(\log \left(\mathrm{V} \mathrm{I}_{\mathrm{t}}\right)=\alpha_{V I}+\beta_{V I} \log \left(L A I_{t}\right)+\epsilon_{t}\)       (1)

where α_VI, β_VI, and є_t ~ N(0, σ²_VI) represent the intercept, slope, and error of the linear regression model, respectively. The evolution of the LAI for each pixel is explained by the GRAMI-rice model using four parameters as θ = (L₀, a, b, and c) that are assumed to be generated from the prior distribution ψ ~ N(μ, D). In these parameters, the transformations are used to guarantee that allfour parameters(L₀, a, b, and c)range from 0 to 1 as follows:

\(\begin{array}{c} \Psi=\left(\psi_{1}, \psi_{2}, \psi_{3}, \psi_{4}\right)= \\ \left(\log \frac{a}{1-a}, \log \frac{b}{1-b}, \log \frac{c}{1-c}, \log \frac{L_{0}}{1-L_{0}}\right) \\ \theta=\theta(\psi)= \\ \left(\frac{e^{\psi_{1}}}{1+e^{\Psi_{1}}}, \frac{e^{\psi_{2}}}{1+e^{\Psi_{2}}}, \frac{e^{\Psi_{3}}}{1+e^{\psi_{3}}}, \frac{e^{v_{4}}}{1+e^{v_{4}}}\right) \end{array}\)       (2)

Both the regression coefficients (α_l , β_l , σ_l² ) and the hyper-parameters(μ, D) are then obtained fromthe data collected in previous studies (Kim et al., 2017; Ko et al., 2015).These contain both theVIs and themeasured LAI. The parameter μ is quantified using the default values(L0 = 0.2, a = 3.25 × 10^-1 , b = 1.25 × 10^-3 , and c = 1.25 × 10^-3). Parameter D is chosen once a diagonal matrix with all diagonal elements corresponds to 0.5.

The following four processes are employed to obtain θ for each pixel.

Step 1 is to set μ as the initial guess of ψ for each pixel. In Step 2, define LAIt = G˜(t; ψ) = G(t; θ(ψ)) and consider the objective function as follows:

\(\begin{array}{c} \sum_{l=1}^{5}\left\{\frac{1}{\sigma_{1}^{2}} \sum_{t=1}^{n}\left(\log V I_{l t}-\alpha_{1}-\beta_{1} \log \tilde{G}(t ; \psi)\right)^{2}\right\}+ \\ (\psi-\mu)^{\prime} D^{-1}(\psi-\mu) \end{array}\)       (3)

where the parameter ψ is estimated by minimizing the above function, and the optimization is performed using the POWELL optimization routine (Press et al., 1992). Step 3 isto generate the simulated curve for each pixel from the estimated ψ in Step 2. Lastly in Step 4, update μ, D asthe sample means and sample variances of the estimates in Step 2.

A Crop Information Delivery System (CIDS) was formulated by Ko et al. (2015), as an extended version of the GRAMI-rice model that uses remote sensing images to project pixel-based crop growth and yield maps(Fig. 2a).CIDS employs both pixel-based remote sensing data and climate data as the system’s inputs where climate data are used either a single weather station or from multiple weather stations (pixels) depending on the situation. The GRAMI-rice model then applied to simulate crop growth in each pixel using both types of input data.

 

Fig. 2. Schematic diagram of modeling of spatiotemporal crop productions using the GRAMI model. (LAI = leaf area index; PAR = photosynthetic active radiation; VI = vegetation index).

The four processes (Fig. 2b) involved in simulating daily rice growth are:(1)interception and absorption of the incidentsolarradiation by the leaves;(2) calculation of growing degree-days (GDD); (3) production of the new dry mass by the leaf canopy; and (4) determination ofthe LAI partitioning ofthe new dry mass.The details of these procedures and related equations have been described in earlier studies by Ko et al. (2015) and Jeong et al. (2018). In this study, we have applied the same initial conditions and parameter values used in these studies to calibrate the GRAMI-rice model.

2) Cumulative crop NDVI based on the semi-empirical BRDF model

The red and the near infrared (NIR) channels of GOCI were used to calculate the cumulative NDVI profiles that were used for the input variable of the GRAMI-rice model. It is essential to secure stable NDVI profiles for consistent simulation of crop yield in monsoon areas asthey act as a driving input variable for the GRAMI-rice model. While MODIS-based NDVI is a satellite product most frequently adopted, it can be challenging to obtain continuous NDVI profiles from MODIS during the summer monsoon season in the presentstudy area. Missing values due to clouds or aerosol effects were interpolated and simulated using novel state-of-the-art methodologies of polynomial, function-fitting,filtering, and interpolating approaches. These techniques have generated reasonable estimates of vegetation phenology using optical satellites such as MODIS (Gao et al., 2008; Poggio et al., 2012). However,many uninterruptedmeasurements ofthe crop area should be made to reflect the actual vegetation conditions more accurately. In the current study, we used VIs data from GOCI on the geostationary COMS instead of the MODIS-based product to address this issue.

The tabulated valuesfrom the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) were employed to performatmospheric correction procedure before estimating the NDVI profiles, according to the methodology proposed by Zhao et al. (2000). Surface reflectance data of GOCI were atmospherically corrected using input values ofthe aerosol optical depth (AOD), water vapor, and the total ozone from the MODIS products of MOD04, MOD05, and MOD07. When these MODIS products were unavailable mainly because ofthe cloud contaminations, we substituted the atmospheric constituents such as AOD or water vapor from the Korea MeteorologicalAdministration (KMA) for the respective MODIS derived values (Kim et al., 2008).

The surface anisotropy effects of the geostationary satellites arise predominantly due to the changing sun position. We used the semi-empirical bidirectional reflectance distribution function (BRDF) to obtain angular independent crop NDVI profiles from GOCI. The Ross-Thick Li-SparseReciprocal(RTLSR) model was used to correspond to the Roujean BRDF model that is a linear combination of three primary scattering components ofisotropic scattering, volumetric scattering, and geometric scattering (Ross, 1981; Roujean et al., 1992):

\(\begin{aligned} \rho\left(\theta_{s}, \theta_{v}, \emptyset\right)=& f_{\text {iso}}+f_{\text {geo}} k_{\text {geo}}\left(\theta_{s}, \theta_{v}, \otimes\right)+\\ & f_{\text {vol}} k_{\text {vol}}\left(\theta_{s}, \theta_{v}, \varnothing\right) \end{aligned}\)       (4)

where f_iso, f_geo, and f_vol are the spectrally dependent model parameters, fiso is the Lambertian reflectance in the nadir direction, f_geo is the coefficient of the LiSpare-Reciprocal geometric kernel k_geo, and f_vol is the coefficient of the Ross-Thick volumetric kernel k_vol. The parameters of k_geo and k_vol are relevant kernel functions of the viewing zenith angle θ_v, the solar zenith angle θ_s, and the relative azimuth angle ø.These parameters provide shapes for the volumetric and geometric-optical scattering of BRDFs.

The kernel coefficients of the BRDF model were estimated independently for each pixel location by the inversion of Eq. (4) using satellite observations with sensor-target-solar geometry sensed during the 16-day composite period (Lucht et al., 2000; Schaaf et al., 2002; Schaaf et al., 2011). The Nadir BRDF-Adjusted Reflectance (NBAR) for the MODIS vegetation sensitive bands was estimated fromEq.(4)to normalize the nadir view angles and the mean solar zenith angle. BRDFAdjustedReflectance (BAR)for GOCI was also estimated from Eq.(4) based on itsfixed viewing angle and mean solar zenith angle without adjusting the viewing zenith angle to the nadir direction. This estimation was carried out because the GOCIsensor is unable to gain the nadir direction angularsampling over the study area (Yeom and Kim, 2013).

MODIS NBAR products (MCD43) in use were retrieved every eight days due to the archival limitations (Schaaf et al., 2002; Schaaf et al., 2011). In the present study, theBARvaluesfromthe GOCIsensor ofCOMS were also estimated by a daily rolling strategy over a 16-day retrieval period to obtain a more intuitive interpretation of phenology characteristics and to capture more subtle details (Ju et al., 2010; Shuai et al., 2013).

3) Insolation from COMS MI based on the physical model

The incident solar radiance on the surface (insolation) that is a primary meteorological input variable for the GRAMI model was estimated using the COMS MI to reflect the energy source of photosynthesis on crop canopies. A pixel-based physical model dedicated by Kawamura et al. (1998) was adjusted using instantaneoussatellite observations and atmospheric information, out of various satellitebased insolation models (Nunez, 1993; Otkin et al., 2005). The adjustment was practiced because of difficulty to interpret the radiation influences of atmospheric components and clouds under the hemispherical sphere condition due to complex physical features and time-consuming scheming (Kawamura et al., 1998; Yeom et al., 2016).

The design ofthe satellite-based solarinsolation with COMS MI was adopted from the Kawamura physical model(Kawamura et al., 1998). However, themodified model had an improved cloud factor as it considered visible reflectance of the satellite and the solar zenith angle instead of the brightness temperature because the pass depth of cloudsis more sensitive to the amount of irradiance attenuation (Yeom et al., 2012). The following formulas are details of the physical model (Kawai and Kawamura, 2005; Tanahashi et al., 2001; Yeom et al., 2012).

\(S_{T}=S_{I}+S_{R}+S_{A}\)       (5)

\(S_{I}=S\left(\tau o \tau_{R}-\alpha_{w}\right) \tau_{A}\)        (6)

\(S_{R}=S \tau_{0}\left(0.5\left(1-\tau_{R}\right)\right) \tau_{A}\)       (7)

\(S_{A}=S\left(\tau o \tau_{R}-\alpha_{w}\right) F_{C} \omega_{O}\left(1-\tau_{A}\right)\)       (8)

\(S=I\left(d_{M} / d\right)^{2} \cos \theta\)       (9)

where S_T, S_I, S_R, and S_A are the total solar insolation, direct irradiance, diffuse irradiance due to Rayleigh scattering, and the diffuse irradiance due to scattering by aerosols, respectively. Further details of these parameters are described in Table 2 (Kawai and Kawamura, 2005; Yeom et al., 2016).

Table 2. Description of parameters used for estimation of satellite-based insolation

4) Air temperature from KLAPS based on the numerical model

Air temperature data to simulate the rice production were obtained from the Korea local analysis and prediction system, KLAPS (https://data.kma.go.kr/ cmmn/main.do, accessed on Dec. 12, 2018). The KLAPS was developed to predict weather conditions of the Korean peninsula with a grid resolution of 5 km, up to 24 times a day for 12 hours. The KLAPS simulates reanalysis data with a high resolution of 1.5 km based on its analysis system using all potential measured weather data from the region ofinterest(Kim et al., 2002). The KLAPS also acclimates the data assimilation part of the local analysis and prediction system (LAPS) developed by the US National Oceanic and Atmospheric Administration/Forecast Systems Laboratory (NOAA/FSL). The KLAPS is classified into both data collection and analysis components.The analysis regime is composed of 3-dimensional wind, temperature, humidity, cloud, precipitation, and soil analysis procedures as well as the surface analysis procedures. More details of the KLAPS can be referenced from the earlier studies (Albers, 1999; Albers et al., 1996; McGinley et al., 1991).

5) Classification of paddy fields

Paddy fields in North Korea were classified using several geographical variables and the GOCI NDVI phenological indices(Table 3).The paddy classification index (PCI) was calculated according to the following formula using the median values of the 4-year daily GOCI NDVI data (2011 to 2014).

\(\mathrm{PCI}=\frac{N D V I_{\text {hanest}}-N D V I_{\text {tranglant}}}{N D V I_{\text {hanest}}+N D V I_{\text {tranglant}}}\)        (10)

Table 3. Phenological parameters and topographical variables for classification of paddy fields

where NDVI_harvest andNDVI_transplant represent the NDVI_values at harvest (~ day of year 270) and the NDVI _values at transplanting (~ day of year 140).Geographical variables were calculated using the topography analysis modules as per the SAGA GIS program (Conrad et al., 2015).

The Random Forest (RF) algorithm was used to classify the paddy field (Breiman, 2001). RF is an ensemble classifier that combines a series of classification trees and an operative technique for remote sensing applications owing to the simplicity of model construction and the reliable accuracy of its classification (Belgiu and Drăguţ, 2016). The R package “caret” was used to tune the size ofthe variable subset of RF, described by Kuhn and Johnson (2013). The number of decision trees was set to 1,000. The recursive feature elimination (RFE) was used for variable selection purposes, which calculates the variable importance and removes the least significant variables one by one over again until the best one remains (Guyon and Elisseeff, 2003). The RFE algorithm was implemented by the functions of the R package “caret” (Kuhn and Johnson, 2013). A resampling method for the RFE procedure was performed with repetition five times with crossvalidation at ten times. The variables selected through RFE were in the following order: wetness index, elevation, paddy classification index, maximum peak value, maximum fading date, maximum peak date, maximum fading rate, slope degree, and maximum growth rate. A reference dataset for classification of paddy fields was assembled for calibration using the land use map of two Northern provinces of Kangwon and Gyeonggi in South Korea (data not shown). Since North Korea is unreachable for ground truth data, validation was performed using a digital map of paddy field areas in South Hwanghae, which was independently obtained (data not shown). Three satellite images with high spatial resolution from each of KOMPSAT-3, KOPMSAT-3A, and RapidEye were employed for this purpose, being classified using an on-screen digitization method.

6) Cluster Analysis and statistical evaluation

The study area of North Korea wassplit into 1,300 × 1,500 pixels.The determined paddy areasin pixel were available in only approximately 1% of pixels: 19,531 in 2011, 19,521 in 2012, 19,540 in 2013, and 19,472 in 2014. Approximately 19,409 pixels have observations in all four years. These exhibited the variation in maturity dates, which ranged from115 to 184 days after planting (DAP). The distribution of the maturity dates is shown in Fig. 3 for each year.

Fig. 3. Histograms of the maturity date in 2011 (a), 2012 (b), 2013 (c), and 2014 (d).

Cluster analysis was used to classify the riceproducing regions where the classification was based on the simulated LAI values from DAP = 1 to DAP = 184. Simulated values between the maturity date and the maximum DAPwere filled with zeroes. Both yearby-year analysis and combined four-year analysis were performed. In the year-by-year analysis, the simulated LAI of each DAP was considered as a variable, and K-means clustering was applied to these 184 variables. The elbow ofthe scree plot is 4 for all years,suggesting that 4 clusters are sufficientforthe classification. Mean simulated-LAI curve for each year are presented in Fig. 4. In the combined four-year analysis, the fouryear-averaged LAI curves of all 19,409 pixels (with observations) were first obtained. K-means clustering was then applied to the 184 simulated LAI values on the averaged LAI curves (Fig. 5). The scree plot, density plots, and mean simulated-LAI curve for the four-year averages are also reported. The procedure of K-mean clustering was performed as the following three steps. In step 1, K initial cluster centroids were specified. In step 2, the cluster was assigned, whose centroid is nearest in terms of Euclidean distance for each observation, followed by updating the centroid. In the last step, step 2 was repeated until no more reassignments take place.

Fig. 4. Simulated curves of seasonal changes of mean leaf area index versus days after planting for each of four clusters (1-4) from 2011 to 2014.

Fig. 5. Histograms of the maturity date (a), scree plot against number of clusters (b), proportions of the four clusters in the whole country (c), and four year (2011-2014) averaged simulated curves of seasonal changes of mean leaf area index (LAI) versus days after planting for each of four clusters (d).

Several statistical analysis methods were adopted to evaluate the reliability of the study results. Yield data were analyzed by comparing simulations and observations according to a two-sample t-test using R software version 3.5 (https://www.r-project.org/, accessed on Dec. 20, 2018). Two statistical indices, root mean square difference (RMSD) and model efficiency (ME) (Nash and Sutcliffe, 1970), were also used for model evaluation. The physical meaning of the ME is that it is a normalized statistic that determines the relative magnitude of the residual variance compared to the observed data variance. ME indicates how well the plot of observed versussimulated data fitsthe 1:1 line. ME values range from –∞ to 1; the closer to 1, the more accurate the model. When values are close to zero, the model predictions are less than or as accurate as the observed mean.

4. Results and Discussion

1) Classification of paddy fields and simulation of rice yield

The totalrice field area of North Korea was estimated to be 491,325 ha as an average (Fig. 6 and Table 4). The average value was determined and applied in the currentstudy, considering the limited availability ofthe input information for further detailed year-to-year classification. The overall accuracy of the South Hwanghae validation map was 0.89, with a Kappa coefficient of 50.1. We found that areas classified as paddy fields (491,325 ha) in North Korea generally agreed with reportsfrom FAO (2015) and Zhang et al. (2017) for similar years.According to the FAO report, paddy field areas decreased from 571,360 to 525,000 ha over the study period (2011 to 2014). However, we could not determine the annual variation in the paddy field area due to the limited satellite data used in this study. While this may be a limitation of our results, we minimized potential errors in productivity projections by removing pixels assigned to non-paddy fields during the simulation process. For accurate classification, we found that elevation, PCI, and wetness index are the three most influential parameters (see Table 3) used in this study.

Table 4. The error matrix associated with the paddy area classification map (Fig. 6)

Fig. 6. Classified paddy field areas in North Korea.

Before applying remotely sensed satellite images to simulations of rice yield in North Korea, we calibrated the GRAMI-rice model using the experimental field site inCheorwon, South Korea,followed by validation using geospatial data at two administrative districts bordering North Korea (Yeom et al., 2018). Fig. 7 shows the projection of the geographical distribution of paddy yield in North Korea during the crop-growing season from 2011 to 2014. The annual variation in the paddy productivity reflects an excellent spatial representation.Rice paddy areas with high productivity are mainly located in the Hwanghae province, within the areas classified asthe paddy fieldsin North Korea. According to the report by FAO (2015), annual paddy yields have exhibited an increase until the year 2013, before declining in the year 2014. Looking into the climate conditions for the study period, the average temperatures during the reproductive growth stage (June toAugust) of rice become elevated from 2011 to 2014 in the classified paddy areas of North Korea. The average temperatures were 14.8°C in 2011, 16.3°C in 2012, 16.8°C in 2013, and 17.8°C in 2014.We assume that the elevated temperatures are attributable to the increase in yield from 2011 to 2013. The temperature effect on yield agreesto the earlier report in the current latitudinal zone by Kim et al. (2013). The reason for the yield reduction in 2014 contrariwise appears to be related to the potential shortage of the reservoirs for irrigation to paddy fields. FAO (2015) reported comparatively less volume of water in irrigation reservoirs in 2014 (1,000,000 m³) than the quantities in 2012 (3,746,000 m³) and 2013 (3,644,000 m³). Meanwhile, Fig. 7 indicatesthat the paddy productivity in Hwanghae and Pyongan provinces has decreased in the year 2014 compared to the previous years, indicating that the proposed method in this paper has simulated the spatiotemporal productivity of rice paddies efficiently. Statistically, it was evident that the spatiotemporal variability of paddy yield in North Korea was well reproduced, with an RMSE value of 0.60 t ha^-1 and statistically significant (at p = 0.101) between the simulation and observation dataset (using the values based on FAO reports) with a two-sample t-test (Table 5). While the simulated paddy yield in the administrative provinces of N. Hamgyong and S. Pyongan showed significantly different valuesfrom the observations (at p < 0.05), the other seven provinces under investigation registered statistically significant agreement between the simulations and observations (at p > 0.05).The observed yieldsin North Korea based on the FAO reports generally ranged within the standard deviations (SD) of the corresponding simulated mean yields (Fig. 8). The SD values of several data points farther than the 1:1 ratio are attributed to the simulated yields in N. Hamgyong and S. Pyongan.

Fig. 7. Projections of paddy yields in North Korea for the five years in 2011 (a), 2012 (b), 2013 (c), and 2014 (d). Values at the bottom of each map represent the yearly mean yield ± standard deviation.

Fig. 8. Simulated and observed yields of paddy rice for the nine administrative provinces from 2011 to 2013 and the entire North Korea (NK) from 2011 to 2014 (refer to Table 5). Vertical error bars represent the standard deviations of the simulated mean yields in NK.

Table 5. Comparisons between simulated and observed paddy yields in the nine administrative provinces for three years (2011, 2012, and 2013) and the entire North Korea for four years (2011-2014). Comparison statistics applied were root mean square error (RMSE) and two-sample t-test

Simulated values of crop growth obtained with our model showed statistically significant agreement with the observations as well (data not shown). It was previously demonstrated that incorporation of the GRAMI model with remote sensing data could be applied tomonitor growth and yield of paddy rice using operational satellites (Jeong et al., 2018; Kim et al., 2017; Padilla et al., 2012).In thisstudy, we showed that the GRAMI-rice version is also applicable to regional crop monitoring projects by successfully simulating paddy productivity in North Korea (Fig. 7).While there are some weaknesses in the model, such as strong dependence on remote sensing data in order to carry out simulation tasks with limited input parameters and/or variables, it has significant applications for inaccessible or data-sparse regions of interest. In such a region, it can be almost impossible to monitor or simulate crop productivity without using operational satellite-based remote sensing data. For example, Lee et al. (2016)showed thatsatellite data could determine the impacts of drought on croplands over North Korea.

In addition, the strong dependence of the model on remote sensing data can also be an advantage in more or fewer means. First, the input requirements are simple, using actual observations that represent the environmental conditions. Second, the optimization method improves simulation performance. Third, satellite-based remote sensing is available for any region of interest on the Earth. Meanwhile, the shortcomings include restricted observations during the crop-growing season and inadequate illustration of remote sensing information. These can eventually cause disagreement between simulations and observations.

2) Cluster analysis of rice productivity

The regional characteristics of the climate and environment in North Korea were grouped into four categories using the K-means clustering (Fig. 9).It was evident from this study that the relative proportions of the four groups changed between the years 2011 (Fig. 10) and 2014 (Fig. 11). The peak values changed each year as well (Figs A1 and 2). In fact, in the year 2011, low production categories 1 and 2 exceeded groups 3 and 4, but after 2012, the situation appears to be improved. It istherefore interpreted that the present remote sensing and modeling approach can allow a broad representation of the dynamics of climatic conditions in any study region over time.

Fig. 9. Delineation (a) and distributions (b-e) of the four clusters of NorthKorea paddy grain production regions, classified according to the four year (2011-2014)-averaged simulated curves of seasonal changes of leaf area index (LAI) using K-means cluster analysis with K = 4 and clusters: C = 1 (b), C = 2 (c), C = 3 (d), and C = 4 (e).

Fig. 10. Delineation (a) and distributions (b-e) of the four clusters of North Korea paddy grain production regions, classified according to the year 2011 based on simulated curves of seasonal changes of leaf area index (LAI) using K-means cluster analysis with K = 4 and clusters: C = 1 (b), C = 2 (c), C = 3 (d), and C = 4 (e).

Fig. 11. Delineation (a) and distributions (b-e) of the four clusters of North Korea paddy grain production regions, classified according to the year 2014 based on simulated curves of seasonal changes of leaf area index (LAI) using K-means cluster analysis with K = 4 and clusters: C = 1 (b), C = 2 (c), C = 3 (d), and C = 4 (e).​​​​​​​

Crop productivity differs significantly with the prevailing climatic conditions and the soil types. According to the Koeppen-Geiger classification system (Peel et al., 2007), North Korea encompasses five different climate types: Dfa, Dfb, Dfc, Dwa, and Dwb (Lee et al., 2016). However, using the K-means grouping, the regional characteristics of the North Korean climate and environment can be summarized as four categories instead of five climate types. As soil type data is not available in North Korea, the classification of rice-producing regions using environmental data alone becomes a very challenging issue.Cluster analysisis a descriptive technique aiming to group observations into homogenous classes or clusters. While the solution is not technically unique and dependent upon the choices of the analyst, the technique can allow us to determine groups of observations internally categorized by a high level of consistency. This methodology should provide an option to find out a way of characterizing crop productivity into different classes. An essential contribution ofthisstudy was, therefore, to demonstrate how remote sensing data and the interpolated data from the GRAMI-rice model could be utilized to outline the agricultural regions of a developing country.

5. Conclusions

productivity in North Korea based on simulation of rice yield using the GRAMI-rice model integrated with satellite-based biophysical and meteorological parameters. North Korea requires international aid for the shortage of food although the government is internationally isolated. The results of the study can help find out a new approach for monitoring of agricultural production conditions and design of a decision supportsystem mainly unreachable regions as well as make policy decisionsfor the support of North Korea more efficiently. Meanwhile, there have been few studies performed to assess the agricultural environment in North Korea due to difficultiesto obtain scientific information on agricultural systems such as the productivity of crops of the country. Remote sensing can be a useful tool to monitor spatial areas of inaccessible regions (e.g., North Korea). Remotely sensed information could be applied to improve the possible weakness of assembling data to advance novel mechanisms for establishing or strengthening agricultural precision technologies and addressing production issues(e.g., drought impacts on crops).The accurate inferences of the methodology accomplished in the study should provide constructive significances. The technique adopted could progress the science of remote sensing of the agricultural systems to address agricultural uncertainty and precision issues.Itmay also provide useful methods to deliver several confronting concernsin the face of climate change,such as drought risks, early warning systems, and capital costreduction against measurement-based approaches in developing countries.