Next Article in Journal
Structure Tensor-Based Algorithm for Hyperspectral and Panchromatic Images Fusion
Next Article in Special Issue
Assessing the Spatial and Occupation Dynamics of the Brazilian Pasturelands Based on the Automated Classification of MODIS Images from 2000 to 2016
Previous Article in Journal
Burned Area Mapping of an Escaped Fire into Tropical Dry Forest in Western Madagascar Using Multi-Season Landsat OLI Data
Previous Article in Special Issue
High Spatial Resolution Visual Band Imagery Outperforms Medium Resolution Spectral Imagery for Ecosystem Assessment in the Semi-Arid Brazilian Sertão
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Numerical Terradynamic Simulation Group, College of Forestry & Conservation, University of Montana, Missoula, MT 59812, USA
2
Department of Ecosystem and Conservation Sciences, College of Forestry & Conservation, University of Montana, Missoula, MT 59812, USA
3
Department of Geosciences, University of Montana, Missoula, MT 59812, USA
4
Department of Land Resources and Environmental Science, Montana State University, Bozeman, MT 59717, USA
5
Estación Experimental de Aula Dei, Consejo Superior de Investigaciones Científicas (EEAD-CSIC), 50059 Zaragoza, Spain
6
Department of Microbiology and Plant Biology, Center for Spatial Analysis, University of Oklahoma, Norman, OK 73019, USA
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(3), 372; https://doi.org/10.3390/rs10030372
Submission received: 16 January 2018 / Revised: 20 February 2018 / Accepted: 22 February 2018 / Published: 28 February 2018
(This article belongs to the Collection Google Earth Engine Applications)

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 (HIGPP). 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.

Graphical Abstract

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 (R2 value of 0.46 and 0.33) and North Dakota (R2 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 km2 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 × 106 hectares in 2009 to a maximum of 6.34 × 106 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. 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]:
GPP = L U E m a x × f ( V P D ) × f ( T ) × f P A R × S W r a d × 0.45
f ( T ) = { 0 T T m i n T T m i n T m a x T m i n T m i n < T < T m a x 1 T T m a x
f ( V P D ) = { 0 V P D V P D m a x V P D m a x V P D V P D m a x V P D m i n V P D m i n < V P D < V P D m a x 1 V P D V P D m i n
where LUEmax is a prescribed maximum light use efficiency (g C MJ−1), which is specific for different biome types (Table 2); 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); VPDmax, VPDmin, Tmax and Tmin 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 (VPDmax, VPDmin, Tmin and Tmax) 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 CO2 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]:
f P A R = N D V I N D V I m i n N D V I m a x N D V I m i n ( f P A R m a x f P A R min )
where fPARmax and fPARmin, 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]. NDVImax and NDVImin 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 CO2 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 (GPPM) and flux tower observations (GPPF). We also compared GPPM 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 GPPM 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 GPPM and MOD17H2 products were extracted as the mean values within 1 × 1 km2 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 (PM) 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; 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 (HIGPP). The Monte Carlo Markov chain (MCMC) method was used to calibrate HIGPP (Table 3) by minimizing the root mean square error between modeled crop production (PM) and two-thirds of the NASS reported county scale crop production (PN) records for each crop type in MT.
We also compared the 30-m model crop yield estimates (YieldM) against independent field-scale crop yield observations (YieldF) obtained from 10 m by 10 m plots at four different MT farms. The YieldF 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 YieldM results and YieldF 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 (GPPM) across MT from April to September for the seven major crop types and eight-year (2008 to 2015) study period. The mean GPPM (2008 to 2015) results across the MT domain showed distinct spatial patterns corresponding to the different crop types. The GPPM results showed higher productivity (>4 g C m−2 day−1) in the northwest and central areas, where alfalfa and barley were located. Higher GPPM values were also found in the southeast portion of MT associated with maize production. The central and northeast areas of MT showed slightly lower GPPM (<2.8 g C m−2 day−1) associated with durum wheat, spring wheat and winter wheat production.
The GPPM 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 GPPM and MODIS GPP were extracted as the spatial mean values within 1 × 1 km2 windows representing each major crop type. The GPPM results showed consistent seasonal variations, but different magnitudes among the different MT crop types relative to the MOD17A2H global GPP product (Figure 4). The GPPM 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 GPPM 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 GPPM mean eight-day climatology derived from the 2008 to 2015 record was compared with the GPPF during the growing season (r = 0.77), but with a relatively large bias (1.73 g C m−2 day−1) and RMSE (2.16 g C m−2 day−1, Figure 4); these results are consistent with the GPPM depiction of a cereal cropland for this site relative to the GPPF observations representing a less productive natural grassland [51]. The GPPM results also showed stronger correspondence with the MODIS GPP record (r = 0.97, p < 0.05) at the US-FPe site than with the GPPF, including smaller positive GPPM bias.
The eight-day 30-m GPPM 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. The eight-day GPPM 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 GPPM 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 GPPM for barley were relatively close to the MODIS GPP record (r = 0.87; bias = −0.14 g C m−2 day−1). Significant differences in magnitude between GPPM 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 GPPM and MODIS GPP during the 2008 to 2015 growing seasons (r ≥ 0.89, Table 4). Moreover, the GPPM results showed the lowest correlations (r = 0.38, p < 0.05) and largest positive bias (2.69 g C m−2 day−1) 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 HIGPP values for each crop type (Table 3), crop yield (YieldM) 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 (PM) with reported annual county-level NASS crop production (PN) values for the 2008 to 2015 record and MT domain are summarized in Table 5 and Figure 5; these results indicated strong correspondence between PM and PN at the MT county level (r ≥ 0.85, p < 0.05). The PM results also showed large spatial and annual variability in MT crop yields, which were generally consistent with the NASS crop production survey record. The PM results showed the strongest correlation with PN for durum wheat (r ≈ 1.00, p < 0.05), with moderate model bias (−10.2%) and low relative RMSE (21.9%). The PM relationship with PN 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 PM and PN 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 PM and PN for alfalfa.
The spatial pattern of estimated 30-m resolution YieldM 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 YieldM 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 YieldM results showed relatively low correspondence with the field-scale crop yield data (YieldF) obtained at four farms representing barley, peas, spring wheat and winter wheat. These results showed moderate correspondence between YieldM and YieldF 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 YieldF 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 YieldM 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 × 106 hectares) to 2015 (5.95 × 106 hectares). The statewide averaged crop production for alfalfa was 3.01 × 106 ton during the study period, with two annual production peaks in 2010 (3.83 × 106 tons) and 2013 (3.21 × 106 tons) (Figure 8). The amount of land allocated to alfalfa increased during the study period at a rate of 1.52 × 105 hectares per year. Alfalfa is relatively productive (3.77 ± 1.06 tons/hectare) and has a larger planted area (3.63 to 18.05 × 105 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 × 105 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 × 106 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 × 105 hectares changed from cereal to broadleaf crop types, and 3.50 × 105 hectares changed from broadleaf to cereal crops, while cropland area varied by 3.90 × 105 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 GPPM calculations.
Despite the potential model uncertainties, the modeled crop production (PM) results were favorably correlated with the county-scale NASS annual crop production (PN) 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 (YieldM) showed lower correspondence with field-scale yield (YieldF) 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 (HIGPP) used to derive YieldM from the 30-m GPP calculations was calibrated using county-level NASS annual crop yield data (YieldN), which represents a major source of model YieldM uncertainty at finer spatial scales. A second consideration is that the YieldF 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 YieldM 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 YieldM and YieldF 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 Yieldm 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 HIGPP 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 YieldM and YieldF 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 × 105 hectares year−1), 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 × 105 hectares year−1, 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 × 104 hectares year−1 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 HIGPP 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 HIGPP 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.

Acknowledgments

This work was conducted at the University of Montana with funding provided by the USDA NIFA (National Institute of Food and Agriculture) contract 658 2016-67026-25067, USDA 365063, NASA EPSCoR (Established Program to Stimulate Competitive Research) contract 80NSSC18M0025 and NASA (NNX14AI50G, NNX14A169G, NNX08AG87A). The field-level crop yield observations were funded by the Montana Research and Economic Development Initiative.

Author Contributions

M.H., J.S.K. and M.P.M. developed this study. M.H. wrote the manuscript, and all authors discussed and reviewed the manuscript. B.D.M. provided field scale crop yield data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lobell, D.B.; Asner, G.P.; Ortiz-Monasterio, J.I.; Benning, T.L. Remote sensing of regional crop production in the Yaqui Valley, Mexico: Estimates and uncertainties. Agric. Ecosyst. Environ. 2003, 94, 205–220. [Google Scholar] [CrossRef]
  2. Lobell, D.B.; Cassman, K.G.; Field, C.B. Crop yield gaps: Their importance, magnitudes and causes. Annu. Rev. Environ. Resour. 2009, 34, 179–204. [Google Scholar] [CrossRef]
  3. Moriondo, M.; Maselli, F.; Bindi, M. A simple model for regional wheat yield based on NDVI data. Eur. J. Agron. 2007, 26, 266–274. [Google Scholar] [CrossRef]
  4. Yuan, W.; Cai, W.; Nguy-Robertson, A.L.; Fang, H.; Suyker, A.E.; Chen, Y.; Dong, W.; Liu, S.; Zhang, H. Uncertainty in simulating gross primary production of cropland ecosystem from satellite-based models. Agric. For. Meteorol. 2015, 207, 48–57. [Google Scholar] [CrossRef]
  5. Ines, A.V.M.; Das, N.N.; Hansen, J.W.; Njoku, E.G. Assimilation of remotely sensed soil moisture and vegetation with a crop simulation model for maize yield prediction. Remote Sens. Environ. 2013, 138, 149–164. [Google Scholar] [CrossRef]
  6. Singh, R.; Semwal, D.P.; Rai, A.; Chhikara, R.S. Small area estimation of crop yield using remote sensing satellite data. Int. J. Remote Sens. 2002, 23, 49–56. [Google Scholar] [CrossRef]
  7. Sakamoto, T.; Gitelson, A.A.; Arkebauer, T.J. MODIS-based cron grain yield estimation model incorporating crop phenology information. Remote Sens. Environ. 2013, 131, 215–231. [Google Scholar] [CrossRef]
  8. Fang, H.; Liang, S.; Hoogenboom, G.; Teasdale, J.; Cavigelli, M. Corn-yield estimation through assimilation of remotely sensed data into the CSM-CERES-Maize model. Int. J. Remote Sens. 2008, 29, 3011–3032. [Google Scholar] [CrossRef]
  9. Ren, J.; Chen, Z.; Zhou, Q.; Tang, H. Regional yield estimation for winter wheat with MODIS-NDVI data in Shandong, China. Int. J. Appl. Earth Obs. Geoinf. 2008, 10, 403–413. [Google Scholar] [CrossRef]
  10. Lobell, D.B. The use of satellite data for crop yield gap analysis. Field Crops Res. 2013, 143, 56–64. [Google Scholar] [CrossRef]
  11. Prasad, A.K.; Chai, L.; Singh, R.P.; Kafatos, M. Crop yield estimation model for Iowa using remote sensing and surface parameters. Int. J. Appl. Obs. Geoinf. 2006, 8, 26–33. [Google Scholar] [CrossRef]
  12. Bolton, D.K.; Friedl, M.A. Forecasting crop yield using remotely sensed vegetation indices and crop phenology metrics. Agric. For. Meteorol. 2013, 173, 74–84. [Google Scholar] [CrossRef]
  13. Mkhabela, M.S.; Bullock, P.; Raj, S.; Wang, S.; Yang, Y. Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agric. For. Meteorol. 2011, 151, 385–393. [Google Scholar] [CrossRef]
  14. Xiao, X.M.; Hollinger, D.; Aber, J.; Goltz, M.; Davidson, E.A.; Zhang, Q.Y.; Moore, B., III. Satellite-based modeling of gross primary production in an evergreen needleleaf forest. Remote Sens. Environ. 2004, 89, 519–534. [Google Scholar] [CrossRef]
  15. Running, S.W.; Nemani, R.R.; Heinsch, F.A.; Zhao, M.S.; Reeves, M.; Hashimoto, H. A continuous satellite-drived measure of global terrestrial primary production. BioScience 2004, 54, 547–560. [Google Scholar] [CrossRef]
  16. Hilker, T.; Coops, N.C.; Wulder, M.A.; Black, T.A.; Guy, R.D. The use of remote sensing in light use efficiency based models of gross primary production: A review of current status and future requirements. Remote Sens. Environ. 2008, 404, 411–423. [Google Scholar] [CrossRef] [PubMed]
  17. Jones, L.A.; Kimball, J.S.; Reichle, R.H.; Madani, N.; Glassy, J.M.; Ardizzone, J.V.; Colliander, A.; Cleverly, J.; Desai, A.R.; Eamus, D.; et al. The SMAP Level 4 carbon product for monitoring ecosystem Land-Atmosphere CO2 exchange. IEEE Trans. Geosci. Remote Sens. 2017, 99, 1–16. [Google Scholar] [CrossRef]
  18. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Gower, S.T.; Running, S.W.; Zhao, M.; Costa, M.H.; Kirschbaum, A.A.; Ham, J.M.; Saleska, S.R.; et al. Evaluation of MODIS NPP and GPP products across multiple biomes. Remote Sens. Environ. 2006, 102, 282–292. [Google Scholar] [CrossRef]
  19. Zhang, Y.; Yu, Q.; Jiang, J.; Tang, Y. Calibration of Terra/MODIS gross primary production over an irrigated cropland on the North China Plain and an alpine meadow on the Tibetan Plateau. Glob. Chang. Biol. 2008, 14, 757–767. [Google Scholar] [CrossRef]
  20. Zhao, M.; Heinsch, F.A.; Nemani, R.R.; Running, S.W. Improvements of the MODIS terrestrial gross and net primary prodcution global data set. Remote Sens. Environ. 2005, 95, 164–176. [Google Scholar] [CrossRef]
  21. Zhang, F.; Chen, J.M.; Chen, J.; Gough, C.M.; Martin, T.A.; Dragoni, D. Evaluating spatial and temporal patterns of MODIS GPP over the conterminous U.S. against flux measurements and a process model. Remote Sens. Environ. 2012, 124, 717–729. [Google Scholar] [CrossRef]
  22. Gao, F.; Hilker, T.; Zhu, X.; Anderson, M.C.; Masek, J.G.; Wang, P.; Yang, Y. Fusing Landsat and MODIS data for vegetation monitoring. IEEE Geosci. Remote Sens. Mag. 2015, 3, 47–60. [Google Scholar] [CrossRef]
  23. Hilker, T.; Wulder, M.A.; Coops, N.C.; Seitz, N.; White, J.C.; Gao, F.; Masek, J.G.; Stenhouse, G. Generation of dense time series synthetic Landsat data through data blending with MODIS using a spatial and temporal adaptive reflectance fusion model. Remote Sens. Environ. 2009, 113, 1988–1999. [Google Scholar] [CrossRef]
  24. Ju, J.C.; Roy, D.P. The availability of cloud-free Landsat ETM+ data over the conterminous United States and globally. Remote Sens. Environ. 2007, 112, 1196–1211. [Google Scholar] [CrossRef]
  25. Gao, F.; Masek, J.; Schwaller, M.; Hall, F. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2207–2218. [Google Scholar]
  26. Gao, F.; Anderson, M.C.; Zhang, X.; Yang, Z.; Alfieri, J.G.; Kustas, W.P.; Mueller, R.; Johnson, D.M.; Prueger, J.H. Toward mapping crop progress at field scales through fusion of Landsat and MODIS imagery. Remote Sens. Environ. 2017, 188, 9–25. [Google Scholar] [CrossRef]
  27. Zhu, X.; Helmer, E.H.; Gao, F.; Liu, D.; Chen, J.; Lefsky, M.A. A flexible spatiotemporal method for fusing satellite images with different resolutions. Remote Sens. Environ. 2016, 172, 165–177. [Google Scholar] [CrossRef]
  28. Schmidt, M.; Udelhoven, T.; Gill, T.; Roder, A. Long term data fusion for a dense time series analysis with MODIS and Landsat imagery in an Australian savanna. J. Appl. Remote Sens. 2012, 6, 063512. [Google Scholar]
  29. Zhang, W.; Li, A.; Jin, H.; Bian, J.; Zhang, Z.; Lei, G.; Qin, Z.; Huang, C. An enhanced spatial and temporal data fusion model for fusing Landsat and MODIS surface reflectance to generate high temporal Landsat-like data. Remote Sens. 2013, 5, 5346–5368. [Google Scholar] [CrossRef]
  30. Doraiswamy, P.C.; Hatfield, J.L.; Jackson, T.J.; Akhmedov, B.; Prueger, J.; Stern, A. Crop condition and yield simulations using Landsat and MODIS. Remote Sens. Environ. 2004, 92, 548–559. [Google Scholar] [CrossRef]
  31. Gu, Y.; Wylie, B.K. Downscaling 250-m MODIS growing season NDVI based on multiple-date Landsat images and data mining approaches. Remote Sens. 2015, 7, 3489–3506. [Google Scholar] [CrossRef]
  32. Hwang, T.; Song, C.; Bolstad, P.V.; Band, L.E. Downscaling real-time vegetation dynamics by fusing multi-temporal MODIS and Landsat NDVI in topographically complex terrain. Remote Sens. Environ. 2011, 115, 2499–2512. [Google Scholar] [CrossRef]
  33. Prince, S.D.; Haskett, J.; Steininger, M.; Strand, H.; Wright, R. Net Primary production of U.S. Midwest croplands from agricultural harvest yield data. Ecol. Appl. 2001, 11, 1194–1205. [Google Scholar] [CrossRef]
  34. Lorenz, A.J.; Gustafson, T.J.; Coors, J.G.; de Leon, N. Breeding maize for a bioeconomy: A literature survey examining harvest index and stover yield and their relationship to a grain yield. Crop Sci. 2010, 50, 1–12. [Google Scholar] [CrossRef]
  35. Kemanian, A.R.; Stöckle, C.O.; Huggins, D.R.; Viega, L.M. A simple method to estimate harvest index in grain crops. Field Crops Res. 2007, 103, 208–216. [Google Scholar] [CrossRef]
  36. Peltonen-Sainio, P.; Muurinen, S.; Rajala, A.; Jauhiainen, L. Variation in harvest index of modern spring barley, oat and wheat cultivars adapted to northern growing conditions. J. Agric. Sci. 2008, 146, 35–47. [Google Scholar] [CrossRef]
  37. Iannucci, A.; Di Fonzo, N.; Martiniello, P. Alfalfa (Medicago sativa L.) seed yield and quality under different forage management systems and irrigation treatments in a Mediterranean environment. Field Crops Res. 2002, 78, 65–74. [Google Scholar] [CrossRef]
  38. Reeves, M.C.; Zhao, M.; Running, S.W. Usefulness and limits of MODIS GPP for estimating wheat yield. Int. J. Remote Sens. 2005, 26, 1403–1421. [Google Scholar] [CrossRef]
  39. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google earth engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  40. Lobell, D.B.; Thau, D.; Seifert, C.; Engle, E.; Little, B. A scalable satellite-based crop yield mapper. Remote Sens. Environ. 2015, 164, 324–333. [Google Scholar] [CrossRef]
  41. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-resolution global maps of 21st-century forest cover change. Science 2013, 342, 850–853. [Google Scholar] [CrossRef] [PubMed]
  42. Pekel, J.-F.; Cottam, A.; Gorelick, N.; Belward, A.S. High-resolution mapping of global surface water and its long-term changes. Nature 2016, 540, 418–422. [Google Scholar] [CrossRef] [PubMed]
  43. Robinson, N.P.; Allread, B.W.; Jones, M.O.; Moreno, A.; Kimball, J.S.; Naugle, D.E.; Erickson, T.A.; Richardson, A.D. A dynamic Landsat derived Normalized Difference Vegetation Index (NDVI) product for the Conterminous United States. Remote Sens. 2017, 9, 863. [Google Scholar] [CrossRef]
  44. Labus, M.P.; Nielsen, G.A.; Lawrence, R.L.; Engel, R.; Long, D. Wheat yield estimates using multi-temporal NDVI satellite imagery. Int. J. Remote Sens. 2002, 23, 4169–4180. [Google Scholar] [CrossRef]
  45. NASS Field Crops. Usual Planting and Harvesting Dates. 2010. Available online: http://usda.mannlib.cornell.edu/usda/current/planting/planting-10–29–2010.pdf (accessed on 24 February 2018).
  46. Masek, J.G.; Vermote, E.F.; Saleous, N.E.; Wolfe, R.; Hall, F.G.; Huemmrich, K.F.; Gao, F.; Kutler, J.; Lim, T. A Landsat surface reflectance dataset for North America, 1990–2000. IEEE Geosci. Remote Sens. Lett. 2006, 3, 68–72. [Google Scholar] [CrossRef]
  47. U.S. Geological Survey. Landsat 4–7 Surface Reflectance Product Guide; U.S. Geological Survey: Reston, WV, USA, 2017.
  48. Heinsch, F.A.; Zhao, M.; Running, S.W.; Kimball, J.S.; Nemani, R.R.; Davis, K.J.; Bolstad, P.V.; Cook, B.D.; Desai, A.R.; Ricciuto, D.M.; et al. Evaluation of Remote Sensing based terrestrial productivity from MODIS using regional tower eddy flux network observations. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1908–1925. [Google Scholar] [CrossRef]
  49. Abatzoglou, J.T. Development of gridded surface meteorological data for ecological applications and modelling. Int. J. Climatol. 2013, 33, 121–131. [Google Scholar] [CrossRef]
  50. Peng, D.; Zhang, B.; Liu, L.; Fang, H.; Chen, D.; Hu, Y.; Liu, L. Characteristics and drivers of global NDVI-based FPAR from 1982 to 2006. Glob. Biogeochem. Cycles 2012, 26. [Google Scholar] [CrossRef]
  51. Zhang, L.; Wylie, B.; Loveland, T.; Fosnight, E.; Tieszen, L.L.; Ji, L.; Gilmanov, T. Evaluation and comparison of gross primary production estimates for the Northern Great Plains grassland. Remote Sens. Environ. 2007, 106, 173–189. [Google Scholar] [CrossRef]
  52. Madani, N.; Kimball, J.S.; Affleck, D.L.R.; Kattge, J.; Graham, J.; van Bodegom, P.M.; Reich, P.B.; Running, S.W. Improving ecosystem productivity modeling through spatially explicit estimation of optimal light use efficiency. J. Geophys. Res. Biogeosci. 2014, 119, 1755–1769. [Google Scholar] [CrossRef]
  53. Xin, Q.; Broich, M.; Suyker, A.E.; Yu, L.; Gong, P. Multi-scale evaluation of light use efficiency in MODIS gross primary productivity for croplands in the Midwestern United States. Agric. For. Meteorol. 2015, 201, 111–119. [Google Scholar] [CrossRef]
  54. Lark, T.J.; Mueller, R.M.; Johnson, D.M.; Gibbs, H.K. Measuring land-use and land-cover change using the U.S. department of agriculture’s cropland data layer: Cautions and recommendations. Int. J. Appl. Earth Obs. Geoinf. 2017, 62, 224–235. [Google Scholar] [CrossRef]
  55. Larsen, A.E.; Hendrickson, B.T.; Dedeic, N.; MacDonald, A.J. Taken as a given: Evaluating the accuracy of remotely sensed crop data in the USA. Agric. Syst. 2015, 141, 121–125. [Google Scholar] [CrossRef]
  56. Smith, W.N.; Grant, B.B.; Desjardins, R.L.; Kroebel, R.; Li, C.; Qian, B.; Worth, D.E.; McConkey, B.G.; Drury, C.F. Assessing the effects of climate change on crop production and GHG emissions in Canada. Agric. Ecosyst. Environ. 2013, 179, 139–150. [Google Scholar] [CrossRef]
  57. Li, Y.; Hunag, M. Pasture yield and soil water depletion of continuous growing alfalfa in the Loess Plateau of China. Agric. Ecosyst. Environ. 2008, 124, 24–32. [Google Scholar] [CrossRef]
  58. Steduto, P.; Hsiao, T.C.; Fereres, E.; Raes, D. Crop Yield Response to Water; Food and Agriculture Organization of the United Nations: Rome, Italy, 2012. [Google Scholar]
  59. Burke, M.; Lobell, D.B. Satellite-based assessment of yield variation and its determinants in smallholder African systems. Proc. Natl. Acad. Sci. USA 2017, 114, 2189–2194. [Google Scholar] [CrossRef] [PubMed]
  60. Guan, K.; Wu, J.; Kimball, J.S.; Anderson, M.C.; Frolking, S.; Li, B.; Hain, C.R.; Lobell, D. The shared and unique values of optical, fluorescence, thermal and microwave satellite data for estimating large-scale crop yields. Remote Sens. Environ. 2017, 199, 333–349. [Google Scholar] [CrossRef]
Figure 1. (a) Cropland distributions for the seven major crop types in Montana in 2015 based on the 30-m Cropland Data Layer from National Agricultural Statistics Service(NASS), including alfalfa, barley, maize, durum wheat, spring wheat, winter wheat and peas; areas in white denote other crop types and land cover areas; (b) planted areas for the different crop types per year from 2008 to 2015 are represented as a proportion (%) of the total planted area for all seven MT crop types.
Figure 1. (a) Cropland distributions for the seven major crop types in Montana in 2015 based on the 30-m Cropland Data Layer from National Agricultural Statistics Service(NASS), including alfalfa, barley, maize, durum wheat, spring wheat, winter wheat and peas; areas in white denote other crop types and land cover areas; (b) planted areas for the different crop types per year from 2008 to 2015 are represented as a proportion (%) of the total planted area for all seven MT crop types.
Remotesensing 10 00372 g001
Figure 2. Examples of (a) 30-m NDVI derived from Landsat 5 surface reflectance data after removing cloud-contaminated pixels; (b) 250-m NDVI from the MODIS MOD09Q1 eight-day temporal composite product; and (c) 30-m NDVI created by fusing Landsat and MODIS data from (a,b). The NDVI metric ranges from zero to one between respective low to high levels of vegetation greenness. The grey area represents no NDVI data. All three images are obtained for the selected period during 1 to 8 July 2008 over Valley County, Montana.
Figure 2. Examples of (a) 30-m NDVI derived from Landsat 5 surface reflectance data after removing cloud-contaminated pixels; (b) 250-m NDVI from the MODIS MOD09Q1 eight-day temporal composite product; and (c) 30-m NDVI created by fusing Landsat and MODIS data from (a,b). The NDVI metric ranges from zero to one between respective low to high levels of vegetation greenness. The grey area represents no NDVI data. All three images are obtained for the selected period during 1 to 8 July 2008 over Valley County, Montana.
Remotesensing 10 00372 g002
Figure 3. Seasonal NDVI variations for randomly-selected pixels representing the major MT crop types as derived from the different sensor records for the 2008 to 2015 growing season (April to September), including: Landsat 5 (LS5) and Landsat 7 (LS7) after removing cloud contamination; the MODIS MOD09Q1 product (MODIS); fusion of MODIS MOD09Q1 with Landsat 5 (LS5_fused) and Landsat 7 (LS7_fused). The MT crop types represented include alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat, while the locations of the selected pixels representing each crop type are indicated in the upper left, with the background image representing the seven major MT crop types depicted from the 2015 Cropland Data Layer (CDL).
Figure 3. Seasonal NDVI variations for randomly-selected pixels representing the major MT crop types as derived from the different sensor records for the 2008 to 2015 growing season (April to September), including: Landsat 5 (LS5) and Landsat 7 (LS7) after removing cloud contamination; the MODIS MOD09Q1 product (MODIS); fusion of MODIS MOD09Q1 with Landsat 5 (LS5_fused) and Landsat 7 (LS7_fused). The MT crop types represented include alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat, while the locations of the selected pixels representing each crop type are indicated in the upper left, with the background image representing the seven major MT crop types depicted from the 2015 Cropland Data Layer (CDL).
Remotesensing 10 00372 g003aRemotesensing 10 00372 g003b
Figure 4. The mean (2008 to 2015) eight-day climatologies of modeled GPP (GPPM) and the MOD17A2H v006 GPP (MODIS) for the selected major crop type locations and Fort Peck (US-FPe) natural grassland tower site in Montana. The eight-day climatology (2000 to 2006) of GPP estimated from the US-FPe flux tower measurements (GPPF) is also shown. The grey shading denotes the temporal standard deviation of GPPM (SD_GPPM), and the light orange shading indicates the standard deviation of MODIS GPP (SD_MODIS) for the different crop types, while the light blue shading represents the standard deviation of GPP estimated from flux tower measurements (SD_GPPF) at US-FPe. The locations of the randomly-selected pixels representing the different crop types are indicated in the upper left map, which also shows the CDL cropland distribution map in 2015 as the background.
Figure 4. The mean (2008 to 2015) eight-day climatologies of modeled GPP (GPPM) and the MOD17A2H v006 GPP (MODIS) for the selected major crop type locations and Fort Peck (US-FPe) natural grassland tower site in Montana. The eight-day climatology (2000 to 2006) of GPP estimated from the US-FPe flux tower measurements (GPPF) is also shown. The grey shading denotes the temporal standard deviation of GPPM (SD_GPPM), and the light orange shading indicates the standard deviation of MODIS GPP (SD_MODIS) for the different crop types, while the light blue shading represents the standard deviation of GPP estimated from flux tower measurements (SD_GPPF) at US-FPe. The locations of the randomly-selected pixels representing the different crop types are indicated in the upper left map, which also shows the CDL cropland distribution map in 2015 as the background.
Remotesensing 10 00372 g004
Figure 5. The comparisons between modeled annual crop production and USDA NASS reported county-level crop production data during the study period of 2008 to 2015 for the major crop types across Montana, including alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat. The red dotted line is the 1:1 line, while the black dashed line is the linear regression relationship.
Figure 5. The comparisons between modeled annual crop production and USDA NASS reported county-level crop production data during the study period of 2008 to 2015 for the major crop types across Montana, including alfalfa, barley, maize, peas, durum wheat, spring wheat and winter wheat. The red dotted line is the 1:1 line, while the black dashed line is the linear regression relationship.
Remotesensing 10 00372 g005
Figure 6. Example spatial distribution of estimated crop yield for winter wheat across MT and Liberty County in 2015; grey areas denote other CDL-defined crop types represented in this study, while white areas denote other land cover types excluded from the analysis.
Figure 6. Example spatial distribution of estimated crop yield for winter wheat across MT and Liberty County in 2015; grey areas denote other CDL-defined crop types represented in this study, while white areas denote other land cover types excluded from the analysis.
Remotesensing 10 00372 g006
Figure 7. Standard error of model estimated annual crop production compared with USDA NASS reported crop production records per county for (a) alfalfa, (b) spring wheat, (c) barley and (d) winter wheat across MT. The NA term and associated white areas denote MT counties with no reported NASS crop production data for the crop types examined in this study.
Figure 7. Standard error of model estimated annual crop production compared with USDA NASS reported crop production records per county for (a) alfalfa, (b) spring wheat, (c) barley and (d) winter wheat across MT. The NA term and associated white areas denote MT counties with no reported NASS crop production data for the crop types examined in this study.
Remotesensing 10 00372 g007
Figure 8. The variations of estimated annual crop yields from this study (bars) and USDA NASS reported planted areas (lines) for each major crop type from 2008 to 2015 in Montana; vertical error bars denote spatial variability (standard deviation) in estimated yields.
Figure 8. The variations of estimated annual crop yields from this study (bars) and USDA NASS reported planted areas (lines) for each major crop type from 2008 to 2015 in Montana; vertical error bars denote spatial variability (standard deviation) in estimated yields.
Remotesensing 10 00372 g008
Table 1. Descriptions of data used in this study.
Table 1. Descriptions of data used in this study.
DataTime PeriodSpatial ResolutionTemporal Resolution
Landsat 5 surface reflectance2008 to 201130-m16 day
Landsat 7 surface reflectance2012 to 201530-m16 day
MOD09Q12008 to 2015250-m8 day
MOD17A2H GPP2008 to 2015500-m8 day
Flux tower based GPP2000 to 20061-kmdaily
Cropland Date Layer2008 to 201530-mannual
Gridded Surface Meteorological Dataset2008 to 20154-kmdaily
USDA NASS crop yield/production data2008 to 2015Countyannual
Crop yield field measurements2008 to 201510-mannual
Table 2. LUE model Biome Properties Look-Up Table (BPLUT) used for the MT cropland GPP calculations.
Table 2. LUE model Biome Properties Look-Up Table (BPLUT) used for the MT cropland GPP calculations.
Crop TypeLUEmax (g C MJ−1)VPDmax (Pa)VPDmin (Pa)Tmax (°C)Tmin (°C)
Cereal crop2.556940145.85−23.15
Broadleaf crop2.57000150027.85−2.15
Table 3. Estimated annual harvest index (HI) values for the seven major MT crop types.
Table 3. Estimated annual harvest index (HI) values for the seven major MT crop types.
Crop TypeHI from Literature *Calibrated HIGPP
Alfalfa0.07 to 0.180.55
Barley0.30 to 0.620.42
Maize0.25 to 0.580.44
Durum Wheat0.31 to 0.430.22
Peas0.33 to 0.590.28
Spring Wheat0.31 to 0.530.24
Winter Wheat0.33 to 0.530.35
* HI from the literature includes HI values derived for converting biomass, gross and net primary production to crop yield.
Table 4. Summary of statistical metrics describing relations between the modeled GPP (GPPM) 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).
Table 4. Summary of statistical metrics describing relations between the modeled GPP (GPPM) 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 TyperBiasRMSE
(g C m−2 day−1)(g C m−2 day−1)
Alfalfa0.810.451.48
Barley0.87−0.141.51
Maize0.382.693.92
Durum wheat0.890.801.47
Peas0.910.741.44
Spring wheat0.920.821.39
Winter wheat0.901.201.66
Table 5. Comparisons between mean county-level crop production from the USDA NASS survey (PN) and mean modeled crop production (PM) for each county in Montana from 2008 to 2015 (all of the correlations are significant with p < 0.05).
Table 5. Comparisons between mean county-level crop production from the USDA NASS survey (PN) and mean modeled crop production (PM) for each county in Montana from 2008 to 2015 (all of the correlations are significant with p < 0.05).
Crop TypeMean of PN (103 Ton)Mean of PM (103 Ton)rBias (103 Ton)Relative BiasRMSE (103 Ton)Relative RMSE
Alfalfa67.22 ± 48.9759.81 ± 50.250.85−7.41−11.0%28.4442.3%
Barley33.19 ± 44.8131.50 ± 42.770.96−1.69−5.1%12.2236.8%
Maize10.27 ± 8.199.91 ± 7.840.91−0.36−3.5%3.4633.6%
Peas10.70 ± 13.8710.94 ± 14.390.980.242.2%3.1029.0%
Durum Wheat39.55 ± 68.9135.50 ± 72.211.00−4.05−10.2%8.6721.9%
Spring Wheat59.41 ± 70.9452.49 ± 70.450.97−6.91−11.6%19.0532.1%
Winter Wheat70.84 ± 110.1169.87 ± 107.930.99−0.97−1.4%14.8721.0%
All Crops51.08 ± 69.2247.10 ± 68.200.96−3.98−7.8%18.8937.0%

Share and Cite

MDPI and ACS Style

He, M.; Kimball, J.S.; Maneta, M.P.; Maxwell, B.D.; Moreno, A.; Beguería, S.; Wu, X. Regional Crop Gross Primary Productivity and Yield Estimation Using Fused Landsat-MODIS Data. Remote Sens. 2018, 10, 372. https://doi.org/10.3390/rs10030372

AMA Style

He M, Kimball JS, Maneta MP, Maxwell BD, Moreno A, Beguería S, Wu X. Regional Crop Gross Primary Productivity and Yield Estimation Using Fused Landsat-MODIS Data. Remote Sensing. 2018; 10(3):372. https://doi.org/10.3390/rs10030372

Chicago/Turabian Style

He, Mingzhu, John S. Kimball, Marco P. Maneta, Bruce D. Maxwell, Alvaro Moreno, Santiago Beguería, and Xiaocui Wu. 2018. "Regional Crop Gross Primary Productivity and Yield Estimation Using Fused Landsat-MODIS Data" Remote Sensing 10, no. 3: 372. https://doi.org/10.3390/rs10030372

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop