Regional Crop Gross Primary Productivity and Yield Estimation Using Fused Landsat-MODIS Data

Abstract

Accurate crop yield assessments using satellite remote sensing-based methods are of interest for regional monitoring and the design of policies that promote agricultural resiliency and food security. However, the application of current vegetation productivity algorithms derived from global satellite observations is generally too coarse to capture cropland heterogeneity. The fusion of data from different sensors can provide enhanced information and overcome many of the limitations of individual sensors. In thitables study, we estimate annual crop yields for seven important crop types across Montana in the continental USA from 2008–2015, including alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat. We used a satellite data-driven light use efficiency (LUE) model to estimate gross primary productivity (GPP) over croplands at 30-m spatial resolution and eight-day time steps using a fused NDVI dataset constructed by blending Landsat (5 or 7) and Terra MODIS reflectance data. The fused 30-m NDVI record showed good consistency with the original Landsat and MODIS data, but provides better spatiotemporal delineations of cropland vegetation growth. Crop yields were estimated at 30-m resolution as the product of estimated GPP accumulated over the growing season and a crop-specific harvest index (HI_GPP). The resulting GPP estimates capture characteristic cropland productivity patterns and seasonal variations, while the estimated annual crop production results correspond favorably with reported county-level crop production data (r = 0.96, relative RMSE = 37.0%, p < 0.05) from the U.S. Department of Agriculture (USDA). The performance of estimated crop yields at a finer (field) scale was generally lower, but still meaningful (r = 0.42, relative RMSE = 50.8%, p < 0.05). Our methods and results are suitable for operational applications of crop yield monitoring at regional scales, suggesting the potential of using global satellite observations to improve agricultural management, policy decisions and regional/global food security.

1. Introduction

Accurate quantification of crop yield at regional to global scales is important in supporting policy- and decision-making in agriculture [1,2,3,4]. Numerous approaches have been developed to estimate crop yield for various cropping systems [1,3,4,5,6,7]. Agricultural surveys provide a reliable way to estimate crop yield, but are less effective over larger regions due to excessive time and budget constraints [8]. In recent decades, satellite remote sensing has been employed for agricultural applications, including crop yield monitoring [9,10,11]. Current operational satellite records, including Landsat and MODIS (MODerate resolution Imaging Spectroradiometer), are sensitive to photosynthetic vegetation cover and provide frequent observations with global coverage and consistent sampling, as well as relatively long-term overlapping records.

Traditional crop yield estimation methods have used empirical relationships between vegetation biomass and remote sensing spectral vegetation indices to estimate yields [3,11,12]. For example, crop yield derived from MODIS NDVI data from 2000 to 2006 in the Canadian Prairies for barley, canola, field peas and spring wheat accounted for 48 to 90%, 32 to 82%, 53% to 89% and 47 to 80% of the variability in reported crop yield from Statistics Canada, respectively [13]. However, these empirical models are fundamentally simple and specific to the limited areas and conditions from which they were developed and cannot easily be extended to other areas.

Another approach involves estimating crop yield as the product of vegetation gross primary productivity (GPP) and an empirical harvest index (HI) specific to different crop types. GPP, representing the total carbon uptake by plant photosynthesis, can be estimated at spatial and temporal scales suitable for cropland applications using a light use efficiency (LUE) model driven by remote sensing inputs [14,15,16]. Two global operational GPP products are currently produced using the satellite data-driven LUE model logic, including the NASA MODIS MOD17 and SMAP (Soil Moisture Active Passive) Level 4 Carbon (L4C) products [15,17]. The MOD17 product provides continuous GPP estimates with eight-day temporal fidelity and 500-m spatial resolution (Version 6) spanning all global vegetated ecosystems and extending from 2000 to the present. The L4C product uses a similar LUE model framework driven by combined satellite information from MODIS and SMAP sensors to estimate GPP and underlying environmental constraints to vegetation growth, including soil moisture-related water supply controls; the L4C product is derived globally from 2015 to the present and provides daily temporal fidelity and 1 to 9-km resolution. However, while operational GPP products derived from global satellite observations provide consistent and frequent temporal sampling, the spatial scale of these products may be too coarse for many agricultural applications; the global LUE algorithm parameterizations for croplands used in the MOD17 and L4C products also only distinguish general crop functional types (e.g., cereal vs. broadleaf), which can degrade GPP accuracy for agricultural ecosystems [17,18,19,20,21]. Alternatively, GPP products derived using LUE model parameterizations that distinguish a greater number of crop types, and with finer spatial resolution (e.g., 30-m) and suitable temporal fidelity (e.g., eight-day), may overcome many of the above limitations while enhancing the utility of these data for agricultural applications. The Landsat TM and ETM+ sensors on the Landsat 5 and 7 platforms provide 30-m resolution imagery that is well suited for capturing surface spectral reflectance heterogeneity at the level of individual agricultural fields [22,23]. However, Landsat has limited temporal coverage due to a long revisit cycle (16-day), data loss from atmosphere aerosol and cloud contamination and failure of the Landsat 7 sensor Scan Line Corrector (SLC) in May 2003 [23,24]; these factors contribute to degraded Landsat utility for cropland monitoring [23,25]. MODIS provides similar spectral information as Landsat, but with more frequent eight-day composite global observations of surface conditions, albeit at a coarser (250 to 1000 m) spatial resolution that is less suitable for heterogeneous agricultural landscapes [22]. Data fusion methods have been used to reduce the constraints of single sensor remote sensing by blending similar spectral information from Landsat and MODIS to generate harmonized multi-sensor observations, providing both relatively fine-scale spatial resolution and frequent temporal sampling [22,25,26,27,28,29,30]. The spatial and temporal adaptive reflectance fusion model (STARFM), developed by [25], has been widely used for blending surface reflectance data from Landsat and MODIS. The STARFM approach was modified and improved by [23] for more complex situations. However, the STARFM model is computationally expensive, which can impose a constraint on regional applications. Alternatively, a rule-based piecewise regression model based on MODIS and Landsat data was developed by [31] to derive 30-m NDVI (Normalized Difference Vegetation Index) maps over northeastern Colorado. A simple linear relationship between NDVI and fPAR (fraction of photosynthetically-active radiation) was also used to downscale 250-m MODIS fPAR data to 30-m resolution using Landsat TM NDVI data [32]. Similar vegetation information may be used as primary inputs to a satellite-based LUE model to derive GPP with enhanced spatial and temporal resolution suitable for agricultural applications.

Although GPP represents the total amount of carbon accumulation in vegetation biomass from photosynthesis, most agricultural applications are concerned with estimating crop grain yield, which generally represents a much smaller portion of GPP. Crop yield can be estimated as the product of GPP and an empirical harvest index (HI) that defines the conversion ratio between crop GPP and grain yield, with HI varying for different crop types and environmental conditions [33,34,35,36,37]. MODIS GPP was converted to wheat yield for the 2001 and 2002 growing seasons in Montana and North Dakota by [38]. The resulting county-level estimated wheat yields from this study showed relatively low correlations with observed wheat yield in Montana (R² value of 0.46 and 0.33) and North Dakota (R² of 0.06 and 0.16), respectively; however, the model results explained 67% and 33% of the total variance in observed wheat yield over both states [38]. Inaccurate HI settings may propagate to significant errors in the final crop yield calculations; however, the integration of regional survey data can improve HI accuracy and reduce uncertainties in crop yield predictions. The U.S. Department of Agriculture (USDA) National Agricultural Statistics Service (NASS) provides annual crop yield and harvested area information at the county level for specific crop types within the continental USA, which can be used as a benchmark for crop model calibration (e.g., HI) and validation.

Google Earth Engine (GEE) is a cloud-based platform designed for efficient planetary-scale geospatial analysis. GEE has a large data catalog co-located with massive CPU that allows for interactive data exploration, providing an easy and quick way to process and analyze data [39]. GEE has been used for mapping crop yields with satellite data spanning multiple states and years (2000 to 2013) in the Midwestern United States [40]. GEE has also been used for other hydrological and ecological applications, including mapping global forest change [41], detecting global surface water change [42] and generating a dynamic Landsat-based NDVI product for the Conterminous United States [43].

In this paper, we show that a global satellite data-driven productivity model (i.e., the MODIS MOD17 algorithm) can be adapted for regional crop assessment. Our approach is significant because it enables crop monitoring over large regions through the modification and application of existing operational satellite remote sensing data and models. We demonstrate the model approach by characterizing the spatial extent and annual variability in cropland productivity (GPP) and yield over the state of Montana in the continental USA. Montana’s agriculture is mostly rainfed, and production is very sensitive to climate variability. Montana is also a major producer of wheat, barley and pulse crops and is representative of intermountain western agriculture.

We used a satellite data-driven LUE modeling approach based on the MODIS MOD17 algorithm to derive relatively high spatial (30-m) resolution GPP and crop yield predictions across the region and spanning the recent (2008 to 2015) satellite record. A data fusion approach was used for blending 30-m Landsat 5 and Landsat 7 NDVI with 250-m eight-day MODIS NDVI observations to produce harmonized 30-m eight-day fused NDVI data record for Montana extending from 2008 to 2015. The 30-m NDVI record was used as a primary LUE model input to estimate fPAR and GPP at a similar 30-m resolution and eight-day time step for the major Montana crop types. Both the LUE modeling framework and satellite NDVI data fusion were implemented within the GEE framework for this investigation. We derived 30-m annual crop yield maps encompassing Montana cropland areas using the resulting 30-m GPP record and crop-specific HI coefficients calibrated and validated using NASS county-level crop production data and field-level yield data. These results were used to assess regional patterns and anomalies in cropland productivity and yield over the multi-year satellite record. A detailed summary of the methods and results from this study is described below.

2. Materials and Methods

2.1. Study Area

Montana (MT, 44°~49°N, 104°W~116°W) encompasses 380,832 km² of the northwest continental United States and is the largest landlocked state in the country. Agriculture is the largest economic sector in MT [44]. The Cropland Data Layer (CDL) from National Agricultural Statistics Service (NASS) provides cropland classification maps over the continental U.S. domain from 2008 to the present; the CDL is a geo-referenced and crop-specific land cover data layer derived using relatively fine-scale satellite imagery (e.g., Landsat 4/5/7/8, Terra MODIS NDVI, etc.) with a decision tree classifier. Croplands encompassed 29.60 ± 1.98 percent of the entire MT land area over the eight-year (2008 to 2015) study period based on the 30-m CDL product. The MT croplands are dominated by seven major crop types, including maize, barley, durum wheat, spring wheat, winter wheat, alfalfa and peas, which were selected in this study for subsequent GPP and crop yield predictions. The MT croplands are mainly distributed in the northcentral and northeastern portions of the state and are more dispersed in central and southern areas for the 2008 to 2015 study period (Figure 1a). Alfalfa is dispersed across MT, but is concentrated in the northwest and central areas of the state. Barley is also grown in most counties, but is more prevalent in the north central portion of MT, together with spring wheat and winter wheat. Some spring wheat is also located in the northeast, where durum wheat is also planted. Peas are planted near barley, spring wheat, winter wheat and durum wheat. Maize is mainly distributed in the southeastern, northeastern and south central portions of the state over the 2008 to 2015 record. Crop types and planted areas show strong annual variations in MT from 2008 to 2015 (Figure 1b). For example, the total cropland area of the seven major MT crops ranged from an annual minimum of 4.13 × 10⁶ hectares in 2009 to a maximum of 6.34 × 10⁶ hectares in 2014; whereas individual crop types ranged from 0.7% (maize) to 30.7% (spring wheat) of the total planted area in MT over the study period.

MT is located in the northern temperate latitudes and has a semi-arid continental climate with an extended winter frozen season, where the potential growing season extends from the end of snow cover depletion in spring to the onset of persistent frozen temperatures in autumn. Based on the reported planting and harvest dates for field crops [45], the planting dates for the major MT crop types (except winter wheat) range from early April to mid-June, whereas harvest dates extend from late July to early October. Winter wheat is usually planted in early September and harvested during the following July to August period. Thus, for this study, we assume an annual period for seasonal crop development, harvest and yield from April to September over the 2008 to 2015 record.

2.2. Satellite Based Crop GPP and Yield Modeling

2.2.1. Developing 30-m MODIS-Landsat Fused NDVI Maps for MT from 2008 to 2015

All data processing and modeling work in this study was conducted in the Google Earth Engine (GEE) framework, which has an extensive library of useful datasets and functions, while providing an efficient and quick way to process and analyze large geospatial datasets. The geospatial data used in this study are summarized in

Table 1. Descriptions of data used in this study.

Data	Time Period	Spatial Resolution	Temporal Resolution
Landsat 5 surface reflectance	2008 to 2011	30-m	16 day
Landsat 7 surface reflectance	2012 to 2015	30-m	16 day
MOD09Q1	2008 to 2015	250-m	8 day
MOD17A2H GPP	2008 to 2015	500-m	8 day
Flux tower based GPP	2000 to 2006	1-km	daily
Cropland Date Layer	2008 to 2015	30-m	annual
Gridded Surface Meteorological Dataset	2008 to 2015	4-km	daily
USDA NASS crop yield/production data	2008 to 2015	County	annual
Crop yield field measurements	2008 to 2015	10-m	annual

. The 30-m resolution surface spectral reflectance data records generated from the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS, [46]) from Landsat 5 (2008 to 2011) and Landsat 7 (2012 to 2015) were extracted and processed over the MT domain. Through LEDAPS, lower order spectral imagery from Landsat 5 (or 7) were calibrated, converted to top-of-atmosphere reflectance and then atmospherically corrected using the MODIS/6S methodology to produce the surface reflectance product [46]. Cloud screening was applied to the reflectance data to remove cloud-contaminated pixels according to cloud mask information provided by the surface reflectance data products [47]. The surface reflectance data for Landsat Bands 3 (red: 0.63 to 0.69 μm) and 4 (near-infrared: 0.76 to 0.90 μm) were then selected to calculate the 30-m resolution NDVI for each image sequence. Simultaneously, the MODIS 250-m 8-day global surface reflectance product (MOD09Q1), after masking out cloud-contaminated pixels, was used to derive the MODIS NDVI record over the MT domain and 2008 to 2015 study period. Both Landsat and MODIS NDVI data were reprojected to a consistent geographic projection using the WGS84 datum. Since the MODIS (250-m 8-day) and Landsat (30-m) records spanned the same MT domain and period, we assumed that the NDVI retrievals from the different sensor records are consistent and comparable. Therefore, a pixel-wise linear regression model was applied to blend the Landsat and MODIS NDVI records to construct a fused NDVI record with consistent 30-m spatial resolution and 8-day temporal fidelity over the MT domain and for each year of record from 2008 to 2015 (https://code.earthengine.google.com/1fb49cd39a59bf1d6d653bc04499c690). The linear regression model describing the relationship between Landsat and MODIS NDVI was developed and applied on a pixel basis. The fused 30-m NDVI record was derived using the resulting slopes and intercepts from the linear regression model and the MODIS NDVI data. The underlying assumption of Landsat and MODIS spectral consistency used to derive the fused NDVI record was determined to have a suitable signal-to-noise ratio for our agricultural application, but may be insufficient for distinguishing subtler environmental trends.

2.2.2. Cropland GPP Estimation

GPP for the seven major crop types in MT during 2008 to 2015 was calculated using an LUE framework similar to the MODIS MOD17 model [14]:

$$\mathrm{GPP}=LUEmax\times f\left(VPD\right)\times f\left(T\right)\times fPAR\times SWrad\times 0.45$$

(1)

$$f\left(T\right)=\{\begin{array}{cc}0& T\le T_{min}\\ \frac{T-T_{min}}{T_{max}-T_{min}}& T_{min}<T<T_{max}\\ 1& T\ge T_{max}\end{array}$$

(2)

$$f\left(VPD\right)=\{\begin{array}{cc}0& VPD\ge VP{D}_{max}\\ \frac{VP{D}_{max}-VPD}{VP{D}_{max}-VP{D}_{min}}& VP{D}_{min}<VPD<VP{D}_{max}\\ 1& VPD\le VP{D}_{min}\end{array}$$

(3)

where LUEmax is a prescribed maximum light use efficiency (g C MJ⁻¹), which is specific for different biome types (

Table 2. LUE model Biome Properties Look-Up Table (BPLUT) used for the MT cropland GPP calculations.

Crop Type	LUEmax (g C MJ−1)	VPDmax (Pa)	VPDmin (Pa)	Tmax (°C)	Tmin (°C)
Cereal crop	2.55	6940	1	45.85	−23.15
Broadleaf crop	2.5	7000	1500	27.85	−2.15

); photosynthetically-active radiation (PAR) is estimated as 45% of shortwave radiation (SWrad × 0.45; [48]); VPD is the mean daily vapor pressure deficit (Pa); T is the minimum daily surface air temperature (°C); VPD_max, VPD_min, T_max and T_min represent respective maximum and minimum VPD and temperature (T) levels for plant photosynthesis, which vary for different plant types (Table 2); f(VPD) and f(T) are dimensionless (0 to 1) scalar linear ramp functions ranging between optimal (1) and fully constrained (0) levels under unfavorable atmospheric moisture deficit (VPD) and minimum daily temperature conditions; the product of f(VPD) and f(T) describes the daily LUE and GPP reduction from potential (LUEmax) conditions due to environmental stress. The LUE model environmental response functions are prescribed for different biome types within a global Biome Properties Look-Up Table (BPLUT) that distinguishes two major cropland functional types for cereal and broadleaf growth forms (Table 2). For estimating GPP, the seven major MT crop types were grouped into each of these two general plant functional types as either cereal (maize, barley, durum wheat, spring wheat, winter wheat) or broadleaf (peas, alfalfa) BPLUT categories.

The daily minimum air temperature (T) and shortwave radiation (SWrad) inputs to the LUE model were obtained from a Gridded Surface Meteorological dataset (Gridmet, [49]), which provides relatively high spatial resolution (~4 km) daily surface meteorological data from 1979 to 2017 over the continental United States. The daily VPD inputs were derived using relative humidity and air temperature data from Gridmet. The daily surface meteorology for each 30-m pixel used for GPP estimation was subsampled from the overlying 4-km resolution Gridmet cell. The BPLUT parameters (VPD_max, VPD_min, T_min and T_max) for the cereal and broadleaf crop types were obtained from the Soil Moisture Active Passive (SMAP) L4C model GPP product [16]. The L4C GPP product uses a similar LUE model framework and MODIS vegetation observations as the MOD17 GPP product, but the L4C BPLUT parameters are calibrated using historical tower eddy covariance CO₂ flux observations from a global tower carbon-flux network (FLUXNET) that includes multiple cropland tower sites.

The fraction of photosynthetically-active radiation (fPAR) absorbed by crops was calculated from the 30-m 8-day fused NDVI as [1,50]:

$$fPAR=\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\ast \left(fPA{R}_{max}-fPA{R}_{\min }\right)$$

(4)

where fPAR_max and fPAR_min, representing the maximum and minimum fPAR values of vegetation in the domain during the growing season, were assigned respective values of 0.95 and 0.01, following [1] and [50]. NDVI_max and NDVI_min represent the NDVI values corresponding to 98% and 2% of the NDVI frequency distributions for all seven MT crop types each year (April to September) over the entire study period.

Tower eddy covariance CO₂ flux measurement-based daily GPP observations were used to evaluate the satellite-based GPP estimates from this study. However, since no available flux tower sites were located in MT cropland areas, an alternative tower site representing a natural grassland from Fort Peck, MT (48.3077°N, 105.1019°W; US-FPe), was selected for the comparison. The US-FPe tower measurement-based daily GPP observations were compared against overlying satellite-based GPP estimates representing a cereal crop (Table 2) for this location. The US-FPe tower and model data records represented different time periods (2000 to 2006 vs. 2008 to 2015), so we compared GPP 8-day climatology records derived from the model (GPP_M) and flux tower observations (GPP_F). We also compared GPP_M at the US-FPe tower location with the MODIS GPP Version 6 (MOD17H2) product, which has 500-m resolution and 8-day temporal fidelity. The GPP_M and MOD17H2 8-day climatology records were compared for randomly selected locations representing different MT crop types over the 2008 to 2015 period. Both the GPP_M and MOD17H2 products were extracted as the mean values within 1 × 1 km² windows centered over selected locations representing uniform crop type conditions to minimize potential geolocation errors.

2.2.3. Crop Production and Yield Estimation

Crop yield was calculated as the product of the estimated cumulative 8-day crop GPP for the growing season between April and September for each year of the 2008 to 2015 study period and an empirical harvest index (HI) specific to each of the seven MT crop types. Crop production was obtained by integrating the estimated yield over the area planted for each crop type at both county and MT state levels. Modeled crop production (P_M) was compared with respective MT county-level annual production values for each crop type reported by the USDA. The range of HI values for each crop type was assembled from the literature and summarized in

Table 3. Estimated annual harvest index (HI) values for the seven major MT crop types.

Crop Type	HI from Literature *	Calibrated HIGPP
Alfalfa	0.07 to 0.18	0.55
Barley	0.30 to 0.62	0.42
Maize	0.25 to 0.58	0.44
Durum Wheat	0.31 to 0.43	0.22
Peas	0.33 to 0.59	0.28
Spring Wheat	0.31 to 0.53	0.24
Winter Wheat	0.33 to 0.53	0.35

* HI from the literature includes HI values derived for converting biomass, gross and net primary production to crop yield.

; this range includes HI values used to convert aboveground biomass, GPP or net primary production to crop yield. However, in this study, we only consider HI as the conversion ratio from GPP to crop yield (HI_GPP). The Monte Carlo Markov chain (MCMC) method was used to calibrate HI_GPP (Table 3) by minimizing the root mean square error between modeled crop production (P_M) and two-thirds of the NASS reported county scale crop production (P_N) records for each crop type in MT.

We also compared the 30-m model crop yield estimates (Yield_M) against independent field-scale crop yield observations (Yield_F) obtained from 10 m by 10 m plots at four different MT farms. The Yield_F observations were obtained from GPS linked yield monitors on combine harvesters and represent four different crop types, including barley, peas, spring wheat and winter wheat. Here, the collocated spatial mean Yield_M results and Yield_F observations for individual crop types were compared over coarser 90 m by 90 m windows to reduce potential geolocation errors for each year of record from 2008 to 2015.

2.3. Statistical Metrics

The correlation coefficient (r) was used to characterize the correspondence between the model results (Mod) and validation (Val) datasets at a 0.05 p-value significance threshold. Bias, calculated as the difference between Mod and Val (bias = Mod-Val), and the root mean square error (RMSE) were used to evaluate model performance in relation to the GPP and yield observations. The model bias and RMSE metrics were also expressed as a relative percentage of the validation datasets.

3. Results

3.1. NDVI Fusion

An example comparison showing the baseline Landsat and MODIS NDVI records relative to the fused 30-m NDVI record developed from this study is shown for the Valley County sub-region of MT (Figure 2). The Landsat 5 (or 7) NDVI record has extensive gaps in spatial and temporal coverage relative to the MODIS record due to less frequent Landsat temporal sampling (16-day revisit cycle) and data loss from atmosphere aerosol and cloud contamination effects, as well as the SLC-off issues in Landsat 7 (Figure 2a). Nevertheless, the Landsat data provide approximately eight-fold improved spatial resolution NDVI observations relative to MODIS, which enhances the delineation of agricultural fields and heterogeneous crop types. However, the MODIS (MOD09Q1) NDVI record (Figure 2b) provides complete spatial coverage and continuous eight-day sampling. The resulting fused NDVI record (Figure 2c) benefits from the combined qualities of both Landsat and MODIS, while minimizing the limitations of the individual sensor records. For Valley County and other MT areas, the fused NDVI 30-m record shows similar spatial patterns and magnitudes as the component Landsat and MODIS NDVI records. The fused NDVI record also provides complete spatial coverage with continuous eight-day sampling at 30-m resolution over the entire MT study domain and 2008 to 2015 record.

For the entire MT domain, the Landsat 5 and MODIS NDVI records were favorably correlated (r = 0.78 ± 0.22, p < 0.05) for the 2008 to 2011 portion of the study period overlapping with Landsat 5. Similar favorable correlations were also found between Landsat 7 and MODIS NDVI records for the 2012 to 2015 portion of the study period (r = 0.77 ± 0.21, p < 0.05). These results indicate suitable conditions for deriving a fused NDVI record using the linear regression model relationship between the Landsat and MODIS NDVI records.

One pixel for each crop type, including alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat, was randomly selected within the MT domain and over the 2008 to 2015 study period to illustrate the NDVI seasonal variations from Landsat 5 (LS5), Landsat 7 (LS7), MODIS and the fused NDVI records derived from MODIS and Landsat 5 (LS5_fused) and Landsat 7 (LS7_fused) inputs (Figure 3). The NDVI time series from the different sensor products in each plot in Figure 3 have consistent geographic coordinates and overlapping time periods, but have varying footprints due to the different pixel sizes represented from each data record. All of the NDVI records show similar seasonality for each of the seven crop types, with seasonal maximum NDVI values generally occurring in the summer (June to August) and early spring (February to April) and seasonal minimum NDVI values occurring in late autumn (September to October) and winter (December to January) periods. The LS5 and LS7 records have relatively sparse growing season temporal coverage from 2008 to 2015 due to cloud contamination and long revisit cycle (16-day), while MODIS NDVI and the two fused NDVI records provide continuous eight-day time series spanning the entire study period. The LS5 and MODIS NDVI records were well correlated for each MT crop type during the overlapping period from 2008 to 2011 (0.91 ≤ r ≤ 0.98, p < 0.05). The LS7 and MODIS NDVI records for the seven MT crop types also showed strong correspondence for the overlapping period from 2012 to 2015 (0.93 ≤ r ≤ 0.99, p < 0.05). Strong correlations were found between the LS5 and LS5_fused (r ≥ 0.91, p < 0.05) and LS7 and LS7_fused (r ≥ 0.95, p < 0.05) records. The MODIS NDVI record was also strongly correlated with the LS5_fused (r ≥ 0.99, p < 0.05) and LS7_fused (r ≥ 0.97, p < 0.05) results. These results indicated that the LS5_fused and LS7_fused records preserve similar variability as the native LS5, LS7 and MODIS NDVI records. There was no significant difference between LS5_fused and LS7_fused NDVI records from the different crop types for the 2008 to 2011 period of overlapping LS5 and LS7 observations. These results indicated the reasonable accuracy and performance for our agricultural application, allowing production of a continuous (2008 to 2015) 30-m, eight-day NDVI record for MT derived using the linear regression model developed between overlapping 250-m MODIS and 30-m Landsat 5 (2008 to 2011) and Landsat 7 (2012 to 2015) NDVI records.

3.2. Regional GPP Estimation over Montana

The fused NDVI record was used with Gridmet daily meteorological data and the yearly NASS CDL as primary inputs to the LUE model (Equations (1) to (4)) to derive 30-m eight-day GPP (GPP_M) across MT from April to September for the seven major crop types and eight-year (2008 to 2015) study period. The mean GPP_M (2008 to 2015) results across the MT domain showed distinct spatial patterns corresponding to the different crop types. The GPP_M results showed higher productivity (>4 g C m⁻² day⁻¹) in the northwest and central areas, where alfalfa and barley were located. Higher GPP_M values were also found in the southeast portion of MT associated with maize production. The central and northeast areas of MT showed slightly lower GPP_M (<2.8 g C m⁻² day⁻¹) associated with durum wheat, spring wheat and winter wheat production.

The GPP_M mean eight-day climatology derived from April to September over the 2008 to 2015 study period is compared with the MODIS GPP (MOD17A2H) mean eight-day climatology for the same period from randomly-selected pixels representing the different MT crop types in Figure 4. Here, both GPP_M and MODIS GPP were extracted as the spatial mean values within 1 × 1 km² windows representing each major crop type. The GPP_M results showed consistent seasonal variations, but different magnitudes among the different MT crop types relative to the MOD17A2H global GPP product (Figure 4). The GPP_M for alfalfa showed stronger annual variations than the other crop types during the 2008 to 2015 study period, indicated by a larger temporal standard deviation. The GPP_M results also showed less temporal variability than MODIS GPP for barley, peas, spring wheat and winter wheat, relative to the other crop types. At US-FPe, the GPP_M mean eight-day climatology derived from the 2008 to 2015 record was compared with the GPP_F during the growing season (r = 0.77), but with a relatively large bias (1.73 g C m⁻² day⁻¹) and RMSE (2.16 g C m⁻² day⁻¹, Figure 4); these results are consistent with the GPP_M depiction of a cereal cropland for this site relative to the GPP_F observations representing a less productive natural grassland [51]. The GPP_M results also showed stronger correspondence with the MODIS GPP record (r = 0.97, p < 0.05) at the US-FPe site than with the GPP_F, including smaller positive GPP_M bias.

The eight-day 30-m GPP_M results were compared with the eight-day MODIS (MOD17A2H) 500-m GPP record over the entire 2008 to 2015 study period for each selected MT crop type location, and summarized in

Table 4. Summary of statistical metrics describing relations between the modeled GPP (GPP_M) and Version 6 MODIS GPP product (MOD17A2H) for the seven crop types in Montana from April to September during 2008 to 2015 (all of the correlations are significant with p < 0.05).

Crop Type	r	Bias	RMSE
		(g C m−2 day−1)	(g C m−2 day−1)
Alfalfa	0.81	0.45	1.48
Barley	0.87	−0.14	1.51
Maize	0.38	2.69	3.92
Durum wheat	0.89	0.80	1.47
Peas	0.91	0.74	1.44
Spring wheat	0.92	0.82	1.39
Winter wheat	0.90	1.20	1.66

. The eight-day GPP_M results showed generally favorable correspondence (0.38 ≤ r ≤ 0.92) with the MODIS GPP record for all seven major MT crop types.). For alfalfa, the GPP_M results were relatively consistent with the MODIS GPP record from April to mid-June, but with a more productive bias between mid-July and early September (e.g., Figure 4). The seasonal variations and magnitudes of GPP_M for barley were relatively close to the MODIS GPP record (r = 0.87; bias = −0.14 g C m⁻² day⁻¹). Significant differences in magnitude between GPP_M and MODIS GPP occurred during the peak growing season (mid-June to August) for maize, peas and the three wheat species (e.g., Figure 4). However, these crop types all showed strong correlations between GPP_M and MODIS GPP during the 2008 to 2015 growing seasons (r ≥ 0.89, Table 4). Moreover, the GPP_M results showed the lowest correlations (r = 0.38, p < 0.05) and largest positive bias (2.69 g C m⁻² day⁻¹) with MODIS GPP for maize compared to the other six MT crop types. These results are consistent with a generally low LUE and productivity bias in the MODIS GPP record for croplands [17,51,52].

3.3. Crop Yield Monitoring in Montana

Using the calibrated HI_GPP values for each crop type (Table 3), crop yield (Yield_M) was estimated for each MT crop type over the 2008 to 2015 study period. Crop production was then calculated as crop yield multiplied by the total planted area each year for each crop type across MT. The comparisons of the estimated crop production (P_M) with reported annual county-level NASS crop production (P_N) values for the 2008 to 2015 record and MT domain are summarized in

Table 5. Comparisons between mean county-level crop production from the USDA NASS survey (P_N) and mean modeled crop production (P_M) for each county in Montana from 2008 to 2015 (all of the correlations are significant with p < 0.05).

Crop Type	Mean of PN (103 Ton)	Mean of PM (103 Ton)	r	Bias (103 Ton)	Relative Bias	RMSE (103 Ton)	Relative RMSE
Alfalfa	67.22 ± 48.97	59.81 ± 50.25	0.85	−7.41	−11.0%	28.44	42.3%
Barley	33.19 ± 44.81	31.50 ± 42.77	0.96	−1.69	−5.1%	12.22	36.8%
Maize	10.27 ± 8.19	9.91 ± 7.84	0.91	−0.36	−3.5%	3.46	33.6%
Peas	10.70 ± 13.87	10.94 ± 14.39	0.98	0.24	2.2%	3.10	29.0%
Durum Wheat	39.55 ± 68.91	35.50 ± 72.21	1.00	−4.05	−10.2%	8.67	21.9%
Spring Wheat	59.41 ± 70.94	52.49 ± 70.45	0.97	−6.91	−11.6%	19.05	32.1%
Winter Wheat	70.84 ± 110.11	69.87 ± 107.93	0.99	−0.97	−1.4%	14.87	21.0%
All Crops	51.08 ± 69.22	47.10 ± 68.20	0.96	−3.98	−7.8%	18.89	37.0%

and Figure 5; these results indicated strong correspondence between P_M and P_N at the MT county level (r ≥ 0.85, p < 0.05). The P_M results also showed large spatial and annual variability in MT crop yields, which were generally consistent with the NASS crop production survey record. The P_M results showed the strongest correlation with P_N for durum wheat (r ≈ 1.00, p < 0.05), with moderate model bias (−10.2%) and low relative RMSE (21.9%). The P_M relationship with P_N for alfalfa had the lowest correlation (r = 0.85), the second largest relative bias (−11.0%) and the largest relative RMSE (42.3%) difference due to larger annual fluctuations in model estimated production than the county survey record. Alfalfa for hay production is mixed with grasses and occurs on both irrigated and non-irrigated lands across nearly all MT counties, which may explain the relatively larger estimation errors for this crop. The P_M and P_N relationship for winter wheat showed the strongest correspondence (r = 0.99, p < 0.05), with relatively low bias (−1.5%) and RMSE (21.0%). Crops that are mostly rainfed in MT (barley, peas, durum wheat, winter wheat) showed overall lower average estimation uncertainties (Table 5; Figure 5). The model crop production results showed strong correlation (r = 0.96) and low relative bias (−7.8%) with the reported NASS county-level data when the model results were combined for all major MT crops during the 2008 to 2015 study period. However, the model results showed a lower RMSE (37.0%) performance, mainly due to the relatively large difference between P_M and P_N for alfalfa.

The spatial pattern of estimated 30-m resolution Yield_M results is illustrated in Figure 6 for winter wheat in 2015. The 30-m results reveal the heterogeneous distribution of field dimensions and crops within the region, as well as the large variation of Yield_M levels for winter wheat driven largely by spatial and temporal variations in NDVI. Winter wheat accounted for approximately 52.2% of the total cropland area in Liberty County and 24.1% of the total cropland area across MT in 2015. However, this crop type is distributed with other crops (shown in grey) in a complex spatial mosaic that is effectively delineated by the 30-m results.

Despite the enhanced delineation of field characteristics, the 30-m resolution Yield_M results showed relatively low correspondence with the field-scale crop yield data (Yield_F) obtained at four farms representing barley, peas, spring wheat and winter wheat. These results showed moderate correspondence between Yield_M and Yield_F values for winter wheat (r = 0.57) and barley (r = 0.68), but much lower correlations for spring wheat and peas (r < 0.10), including large negative model biases (−12.0–−71.4%) and degraded RMSE (32.9–75.3%) relative to the available observations. The small number of field-scale yield observations makes it difficult to draw conclusions about the performance of the algorithm to retrieve field-scale yields. For instance, the Yield_F observations for peas only represented one year of data from a single farm, generating the largest relative model bias (−71.4%) and RMSE (75.3%) among the four crop types represented; whereas, there were only three site-years of data for spring wheat from three different farms, inducing large relative model bias (−40.1%) and RMSE (50.5%). When all of the field-scale crop yield data were combined, the Yield_M performance showed improved correspondence (r = 0.44, p < 0.05), relative bias (−36.1%) and RMSE (48.8%). Although inconclusive, these results indicate meaningful, but lower model accuracy in delineating field level variations in crop yields relative to the county-level results.

The spatial distribution of the model uncertainties for alfalfa, barley, spring wheat and winter wheat at the county level are presented in Figure 7 because these crops were grown in most counties during 2008 to 2015 in Montana and permit the evaluation of the spatial variability of estimation errors. Estimation errors for alfalfa production (Figure 7a) tended to be larger in the more mountainous western portion of MT and in the eastern end of the state near the border with North Dakota. In these two regions, irrigated agriculture has increased significantly in the last few decades, and the mix of dryland and irrigated farming has the potential to increase the variance in production in these counties. Additionally, western Montana’s complex terrain produces highly variable climatic conditions, which also translates into high production variability. Variability in farming practices (dryland versus irrigated) also applies to spring and winter wheat (Figure 7b,d). For these crops, the most accurate model predictions are concentrated in the north-central part of the state, termed the ‘Golden Triangle’ and known for highly productive dryland grain farming; model accuracy is generally lower moving in any cardinal direction from this region. The spatial patterns of model uncertainty for barley production (Figure 7c) is moderate, with generally better model performance in the north-central portions of MT around Pondera, Teton, Choteau and Glacier counties. However, the model results exhibit somewhat larger uncertainties in the eastern part of the state, where counties are less specialized in barley production and farming practices are more variable.

The eight-year study period was too short for effectively quantifying and diagnosing regional trends in planted area and crop yield/production. However, the total planted area for all seven MT crop types showed a significant increasing trend from 2008 (4.36 × 10⁶ hectares) to 2015 (5.95 × 10⁶ hectares). The statewide averaged crop production for alfalfa was 3.01 × 10⁶ ton during the study period, with two annual production peaks in 2010 (3.83 × 10⁶ tons) and 2013 (3.21 × 10⁶ tons) (Figure 8). The amount of land allocated to alfalfa increased during the study period at a rate of 1.52 × 10⁵ hectares per year. Alfalfa is relatively productive (3.77 ± 1.06 tons/hectare) and has a larger planted area (3.63 to 18.05 × 10⁵ hectares) in MT during 2008 to 2015 than the other six crop types. Barley showed a strong increasing trend in crop production from 2011 to 2015, together with a significant increase in planted area (Figure 8). Maize showed an increasing trend in planted area from 2008 to 2014, but did not show a congruent increase in crop production; maize production also peaked in 2013 (2.00 × 10⁵ ton) relative to the other years examined (Figure 8). Though maize represented the lowest planted area of the major MT crops from 2008 to 2015, peas showed significant and consistent increasing trends in crop production and planted area during the study period (Figure 8). However, the crop yield for peas was 1.24 ± 0.28 tons/hectare, covering approximately 3.7% of the total cropland area in MT. There were no obvious trends in planted area or production for MT wheat crops (Figure 8). Spring wheat and winter wheat were planted over large areas in MT (>1.10 × 10⁶ hectare), while estimated annual production and yield for these two cereal crops averaged 1.27 ± 0.27 and 1.90 ± 0.39 ton/hectare during the study period, respectively.

4. Discussion

Data fusion approaches have been applied to synthesize images from multiple satellite remote sensing sources at different spatial and temporal scales [23,24,25,26,27,28]. The simple linear regression model used in this study blended overlapping 250-m eight-day MODIS NDVI with finer spatial resolution (30-m), but less frequent Landsat NDVI data to generate a continuous NDVI record with 30-m resolution and eight-day temporal fidelity extending across MT from 2008 to 2015. A primary assumption of the data fusion approach is that the spatial pattern of the relationships between MODIS and Landsat NDVI are consistent over the sensor records. However, variations of MODIS and Landsat NDVI relationships may vary due to different acquisition dates, cloud contamination, aerosols, viewing angles and temporal compositing [24], or from sub-grid-scale spatial heterogeneity that is not resolved by the coarser MODIS footprint observations [26]; these factors can introduce uncertainties into the NDVI fusion results, which can propagate into model-derived fPAR error. However, the fPAR uncertainty is only one of several potential sources of LUE model GPP error. Other potential error sources include daily surface meteorological inputs and model BPLUT parameterizations defining maximum productivity rates and environmental response characteristics for different crop types.

The LUE model BPLUT parameterization used in this study distinguishes two major cropland categories (cereal and broadleaf), while the model parameters were calibrated using a global network of tower eddy covariance measurement sites and are consistent with the SMAP L4C operational GPP product [16]. The SMAP L4C GPP product shows favorable accuracy and performance in global croplands [16], indicating that the BPLUT parameters are also suitable for MT croplands. Nevertheless, variability in plant characteristics can occur within these general cropland categories, which can lead to significant model GPP error [52]. However, both the MODIS MOD17 and SMAP L4C operational GPP products utilize a static land cover (LC) classification and 500-m resolution MODIS (MCD12Q1) fPAR record as primary inputs to derive GPP at respective 500-m and 1-km spatial scales. In this study, a finer (30-m) resolution NASS Cropland Data Layer (CDL) and fused NDVI record was used to derive GPP. The GPP results from this study are thus expected to provide enhanced regional accuracy for agricultural applications due to the finer spatial delineation of cropland heterogeneity and model representation of dynamic cropland areas from the NASS CDL record. For example, approximately 1.16 × 10⁵ hectares changed from cereal to broadleaf crop types, and 3.50 × 10⁵ hectares changed from broadleaf to cereal crops, while cropland area varied by 3.90 × 10⁵ hectares in MT from 2014 to 2015 based on the NASS CDL record; these changes are represented by our model results and influence the seasonal phenology, pattern and magnitude of estimated productivity and crop yields across the region. The resulting 30-m GPP simulations over MT croplands produced results that were generally similar, but more productive than the MODIS MOD17 (MOD17A2H, Collection 6) operational GPP product (Figure 4). However, the MODIS GPP product has been reported to have low productivity in croplands [18,20,53], indicating that the model results from this study benefit from a finer spatial resolution and a refined model cropland calibration. Though the Gridmet database shows favorable accuracy in relation to in situ weather station network measurements [49], the daily meteorological inputs may contribute uncertainties that propagate to cropland productivity model errors. Moreover, finer scale heterogeneity in surface meteorology not resolved by the 4-km Gridmet spatial resolution may contribute to model GPP errors derived at 30-m resolution. Management factors, including fertilization and irrigation practices, are not directly considered in our model and may also introduce uncertainties into the resulting GPP_M calculations.

Despite the potential model uncertainties, the modeled crop production (P_M) results were favorably correlated with the county-scale NASS annual crop production (P_N) records for all seven major crop types in MT during the 2008 to 2015 study period (r = 0.96, relative bias = −7.8%, relative RMSE = 37.0%). The county level model performance is robust, and the spatial fidelity achieved with the fusion of sensors allows delineation of production variability within a county, which is an important indicator of potential imbalances within a region. However, the estimated crop yields (Yield_M) showed lower correspondence with field-scale yield (Yield_F) records (r = 0.44). Potential reasons for the lower model performance in estimated field-scale yields are discussed below, but when interpreting the performance metrics, it is important to bear in mind that the state-wise harvest index (HI_GPP) used to derive Yield_M from the 30-m GPP calculations was calibrated using county-level NASS annual crop yield data (Yield_N), which represents a major source of model Yield_M uncertainty at finer spatial scales. A second consideration is that the Yield_F observations used for model validation were obtained from a very limited set of 10 m by 10 m plots at four farm locations, which may not adequately reflect the 30-m resolution Yield_M estimates derived across MT using coarser scale satellite NDVI observations and ancillary data. Uncertainty in the CDL files used to delineate annual crop types and planted areas in this study may also contribute to low Yield_M and Yield_F correspondence. Indeed, the reported CDL accuracy from 2008 to 2015 across MT croplands varies from 69.6% (2010) to 85.0% (2014), while the CDL accuracy levels may be inflated because they do not account for edge effects [54]. CDL accuracy may also be degraded in regions with sparse or complex agriculture [55], which is characteristic of MT croplands.

Alternatively, reported crop yields for spring wheat in Lethbridge, Canada, near the MT border, ranged from 1.5 to 1.8 tons/hectare [56]; reported spring wheat yields also ranged from 0.5 to less than 3 tons/hectare in arid and semi-arid areas of the Canadian Prairies [22]. These reported yields are similar to our model Yield_m results for MT spring wheat (1.27 ± 0.27 tons/hectare). The mean modeled crop yield during 2008 to 2015 for barley in this study is 2.48 ± 0.61 tons/hectare, which is also similar to reported barley yields on the Canadian Prairies (0.9 to 4.5 tons/hectare; [22]). For alfalfa, the resulting crop yield estimates range from 1.58 to 8.06 tons/hectare, while another study reported alfalfa yields from 2.97 to 14.72 tons/hectare [57], similar to our model results. Moreover, alfalfa yields reported by the Food and Agricultural Organization range from 5 to 17 tons/hectare under rainfed conditions with 500 to 800 mm annual rainfall [58]; these values are similar to our model yields for MT alfalfa and indicate that the calibrated HI_GPP parameter used in this investigation is reasonable, though it is higher than the HI values for this crop type reported from prior limited studies. Thus, our model estimates of crop yield and production may be suitable for regional cropland productivity and yield assessments that can inform agricultural decisions. These findings are consistent with a similar regional analysis of smallholder agricultural productivity across Africa indicating that high-resolution satellite imagery can be used to make crop yield predictions that are roughly as accurate as survey-based measurements [59].

A limitation of our current calibration strategy that contributes to modest estimation performance at the field scale is that sub-grid scale heterogeneity in environmental conditions, including soil moisture, nutrient availability and resource competition from weeds, is not accounted for by the regional model, and probably most importantly, calibration does not differentiate between irrigated and rainfed crops because spatially explicit information about this practice is currently not available for MT. More specific model calibrations that account for irrigated and non-irrigated crops are expected to increase the precision of production estimates at sub-county scales. Moreover, the growing season in this study is assumed to extend from April to September for all seven MT crop types examined, which may not be accurate for each crop type. For example, winter wheat is planted in September or October and harvested the following July or August [45]. All of these factors likely contribute to the lower Yield_M and Yield_F correspondence.

For the seven major crop types examined over the 2008 to 2015 study period, although the total planted area in MT is increasing (2.79 × 10⁵ hectares year⁻¹), the planted area of each crop type shows large annual variability (Figure 8), which contributes to variations in annual crop production in addition to other influential factors including climate variability. Spring wheat, winter wheat and alfalfa are the three largest components of the total cropland area in MT (34.4%, 26.8% and 19.7%, respectively) from 2008 to 2015. The land dedicated to alfalfa increased by approximately 1.52 × 10⁵ hectares year⁻¹, and was also the largest contributor to annual variability in MT total planted area according to the NASS CDL record. The planted areas and estimated annual crop production for spring wheat and winter wheat did not show consistent trends during the study period (Figure 8). Simultaneously, barley, maize and peas showed increasing trends in planted area, but at relatively low levels (2.86, 0.48 and 3.73 × 10⁴ hectares year⁻¹ respectively), although the combined extent of all three crops occupied only 13.5% of the total planted cropland area and made only a small contribution to changes in total planted area and crop production for the state.

Because Montana’s agriculture is mostly rainfed and has relatively low diversification, planted areas and crop production are very sensitive to climate variations and to changes in agricultural pricing. The MT crop production for all selected crops showed peaks in 2010 and 2013, which had the largest growing season precipitation volumes during the eight-year study period (351.6 mm in 2010 and 336.4 mm in 2013 and approximately 24.6% and 19.2% above mean precipitation from 2008 to 2015), as well as large ratios of precipitation to potential evapotranspiration (0.45 and 0.40 for the years 2010 and 2013, respectively), indicating wet climate anomalies for these two years. Further research is needed to understand the sensitivity of annual crop yield and planted land variations to external factors such as climate fluctuations or changes in crop prices or cost of agricultural inputs.

In recent decades, agriculture has continuously outpaced all other industries as a proportion of total gross domestic product (GDP) in MT, while agricultural production accounted for over 30% of the state’s basic industry employment, labor income and gross sales [44], suggesting the importance for obtaining consistent crop yield assessments across the entire state. This study provides a potential approach for continuous crop yield monitoring at 30-m resolution, which may better inform agricultural management and policy decisions. However, further studies are needed to improve crop yield estimation accuracy and performance and clarify the role of socioeconomic factors, management practices and climate variability on crop yield patterns and trends. Promising areas for model improvement include the use of additional satellite observations to construct more accurate estimates of key vegetation parameters and drivers [60] and the use of high temporal frequency estimates of crop yield at key decision points during the growing season to determine the projected value of alternative decisions. More detailed cropland inventory data capable of resolving field and farm level heterogeneity may improve the delineation of HI_GPP and yield variations for different crop conditions and consequently improve the utility of yield estimates. In situ monitoring of GPP from flux tower sites representing major MT crop types would provide an effective means to improve model calibration and performance. The results of this study and other ongoing efforts are being used to improve crop monitoring capabilities using global satellite observations to inform resource management and policy decisions and enhance national and global food security.

5. Conclusions

We applied a satellite data-driven LUE model framework to produce 30-m resolution daily GPP and annual crop yield estimates for the seven major crop types in Montana from 2008 to 2015. A fused 30-m eight-day NDVI record was developed using an empirical regression model combining similar spectral information from overlapping MODIS and Landsat imagery. The fused NDVI record was used with 4-km resolution gridded daily surface meteorology as primary inputs to the LUE model to derive fPAR and GPP at 30-m resolution and an eight-day time step over MT cropland areas. The fused NDVI record overcomes many of the limitations of the contributing Landsat and MODIS sensor records, including Landsat data gaps and relatively coarse (250-m) MODIS spatial resolution (Figure 2). The fused record also shows generally consistent NDVI magnitudes and seasonal variations as the original Landsat and MODIS data (Figure 3), while providing complete spatial coverage, eight-day temporal fidelity and 30-m spatial resolution. The resulting 30-m GPP simulations over MT croplands are produced benefiting from the finer resolution (30-m, eight-day) NDVI and the refined model cropland calibration. Estimated annual crop production and yields are derived from the 30-m eight-day GPP simulations over MT agricultural areas defined from the annual USDA NASS cropland data layer (CDL) and calibrated HI_GPP values for the seven major MT crop types. The model results corresponded favorably to NASS county-level crop production data (Table 5), indicating that the modeling approach captures regional patterns and annual variations in annual crop yield across MT and that obtaining regional-scale crop production assessments using existing global production models is possible.