IMPUTATION OF REMOTE SENSING TIME SERIES FOR LOW-LATENCY AGRICULTURAL APPLICATIONS

Imputation of remote sensing time series for low-latency agricultural applications is provided. In various embodiments, a first time series of raster data is read. The first time series spans a geographic region and has a first resolution and a first frequency. A second time series of raster data is read. The second time series spans the geographic region and has a second resolution and a second frequency. The second resolution is lower than the first resolution. The second frequency is higher than the first frequency. A mean time series is determined from the first time series of raster data. A predicted time series of values for a location within the geographic region is determined at the first resolution by determining a first time series of values for the location from the first time series of raster data, determining a second time series of values of the location from the second time series of raster data, and determining the predicted time series by multiple linear regression with the first time series dependent on the mean time series and the second time series.

BACKGROUND

Embodiments of the present disclosure relate to remote sensing data analysis, and more specifically, to imputation of remote sensing time series for low-latency agricultural applications.

BRIEF SUMMARY

According to embodiments of the present disclosure, methods of and computer program products for agricultural index prediction are provided. A first time series of raster data is read. The first time series spans a geographic region and has a first resolution and a first frequency. A second time series of raster data is read. The second time series spans the geographic region and has a second resolution and a second frequency. The second resolution is lower than the first resolution. The second frequency is higher than the first frequency. A mean time series is determined from the first time series of raster data. A predicted time series of values for a location within the geographic region is determined at the first resolution by determining a first time series of values for the location from the first time series of raster data, determining a second time series of values of the location from the second time series of raster data, and determining the predicted time series by multiple linear regression with the first time series dependent on the mean time series and the second time series.

In various embodiments, the mean time series is smoothed. In various embodiments, the predicted time series of values is smoothed.

In various embodiments, each of the first time series of raster data, the second time series of raster data, the mean time series, and the predicted time series correspond to an agricultural index. In some embodiments, the agricultural index comprises normalized difference vegetation index, land surface water index, or and mean brightness.

In various embodiments, a crop type mask is read and determining the mean time series comprises masking the first time series of raster data according to the crop type mask. In various embodiments, a plurality of crop type masks is read, a plurality of mean time series is determined from the first time series of raster data, each mean time series corresponding to one of the plurality of crop type masks, and determining each mean time series comprises masking the first time series of raster data according to the respective crop type mask.

In various embodiments, determining the mean time series comprises selecting the mean time series from the plurality of mean time series. In various embodiments, determining the mean time series comprises selecting one of the plurality of mean time series most similar to the first time series of values. In some embodiments, a crop type associated with the selected one of the plurality of mean time series is determined.

In various embodiments, a time shift is applied to the mean time series based on the first time series. In some embodiments, the time shift is determined by cross-correlation between the mean time series and the first time series.

In various embodiments, the second time series is rescaled. In some embodiments, rescaling the second time series comprises CDF matching.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

FIG. 1illustrates a general process for imputation of remote sensing time series according to embodiments of the present disclosure.

FIG. 2illustrates a method of generating an imputed series according to embodiments of the present disclosure.

FIG. 3illustrates an example of CDF matching for a given pixel according to embodiments of the present disclosure.

FIGS. 4A-Bis a graph of NDVI over time as observed and as derived for a given pixel according to embodiments of the present disclosure.

FIG. 5Ashows raw MODIS NDVI for an exemplary region according to embodiments of the present disclosure.

FIG. 5Bshows the corresponding corrected MODIS NDVI according to embodiments of the present disclosure.

FIG. 6Ashows raw MODIS NDVI for the same exemplary region according to embodiments of the present disclosure.

FIG. 6Bshows corresponding raw HLS NDVI according to embodiments of the present disclosure.

FIG. 6Cshows the ECDF-matched MODIS NDVI according to embodiments of the present disclosure.

FIG. 7Ashows the r2values for HLS versus CDF-matched MODIS NDVI according to embodiments of the present disclosure.

FIG. 7Bshows the corresponding mean error according to embodiments of the present disclosure.

FIG. 8illustrates a method of generating an imputed series according to embodiments of the present disclosure.

FIGS. 9-12show exemplary mean annual time series for a plurality of exemplary crop types according to embodiments of the present disclosure.

FIG. 13shows exemplary predicted HLS values according to embodiments of the present disclosure.

FIG. 14is a graph of NDVI over time as observed and as derived for a given pixel and a given crop according to embodiments of the present disclosure.

FIG. 15Ais a graph illustrating the correlation between HLS observations of NDVI and mean crop time series according to embodiments of the present disclosure.

FIG. 15Bis a graph illustrating the correlation between HLS observations of NDVI and co-located MODIS observations according to embodiments of the present disclosure.

FIG. 16is a graph illustrating the observed versus predicted HLS NDVI based on multiple regression coefficients according to embodiments of the present disclosure.

FIG. 17is a graph illustrating the observed versus HLS NDVI time series according to embodiments of the present disclosure.

FIG. 18Ashows raw MODIS NDVI for an exemplary region according to embodiments of the present disclosure.

FIG. 18Bshows the corresponding corrected MODIS NDVI according to embodiments of the present disclosure.

FIG. 19shows raw corrected NDVI for the same exemplary region according to embodiments of the present disclosure.

FIG. 20Ashows the r2values for HLS versus linear unmixing NDVI according to embodiments of the present disclosure.

FIG. 20Bshows the corresponding mean error according to embodiments of the present disclosure.

FIG. 21shows error bars for r2over multiple exemplary crop types according to embodiments of the present disclosure.

FIG. 22shows error bars for mean error over multiple exemplary crop types according to embodiments of the present disclosure.

FIG. 23illustrates a method of generating an imputed series according to embodiments of the present disclosure.

FIG. 24illustrates an exemplary mean annual time series for soybeans according to embodiments of the present disclosure.

FIGS. 25A-Iillustrate exemplary mean annual time series for a plurality of exemplary crop types and exemplary observed data according to embodiments of the present disclosure.

FIGS. 26A-Iillustrate the Manhattan distance for various crop types and exemplary observed data according to embodiments of the present disclosure.

FIGS. 27-28illustrate the results of cross-validation according to embodiments of the present disclosure.

FIG. 29shows an exemplary analysis tile according to embodiments of the present disclosure.

FIG. 30illustrates exemplary results according to embodiments of the present disclosure for the tile ofFIG. 29.

FIG. 31illustrates exemplary results according embodiments of to the present disclosure for the tile ofFIG. 29.

FIGS. 32-33illustrate cross validation results according embodiments of to the present disclosure for the tile ofFIG. 29.

FIG. 34illustrates a method for agricultural index prediction according to embodiments of the present disclosure.

FIG. 35depicts a computing node according to an embodiment of the present disclosure.

DETAILED DESCRIPTION

Growers rely on frequently updated data layers in order to be proactive and informed with their daily decision-making, especially during the growing season. In addition to weather, satellite-derived vegetation indices such as NDVI can be a very useful source of information about crop health and growing stages.

The Harmonized Landsat Sentinel-2 (HLS) product provides observations at spatial resolutions high enough to resolve sub-field processes, but the latency remains on a weekly time scale given the overpass times of Landsat 8 and Sentinel-2 sensors, as well as cloud and snow coverage. On the other hand, MODIS-based imagery can be available daily, but the spatial resolution of 250+ m is too coarse to work on a field scale.

More generally, there exists a need to reconcile multimodal data where there is a tradeoff between frequency and resolution in order to arrive at a consistent and timely high-resolution combined product suitable for downstream processing such as agricultural index computation.

To address this and other needs in the art, the present disclosure provides statistical solutions based on the synthesis of MODIS and HLS data to generate vegetation indices with low latency and high spatial resolution. These approaches are suitable for supporting analysis and decision making for agriculture in near-real time. It will be appreciated that while various examples provided herein focus on MODIS and HLS, the present disclosure is applicable to other data sources where disparities are present between frequency and resolution.

In various embodiments, the statistical problem is isolated from the interpretation problem by focusing on the creation of an optimal but sensor-specific data set, upstream of analysis and modeling solutions. This allows creation of a daily, spatially comprehensive HLS time series that fills gaps due to missing imagery, or clouds and snow within available imagery.

In a first approach, interpolation is applied. In such embodiments, the value at a specific time is predicted based on the other values in a time series. In a second approach, statistical modeling is applied. Such embodiments utilize covariates available at a specific time, for example using the interpolated value plus CDF-transformed MODIS. In a third approach, spatial unmixing is applied. In such embodiments, mean annual time series for each crop type and co-located MODIS pixel are used. Combinations of these approaches can be used, e.g., the prediction from interpolation, then transformation of MODIS indices based on the co-located HLS time series and/or the mean HLS crop signal can be used in a statistical approach.

As set out below, a goodness of fit and out of sample error statistics are calculated to quantify the added value of each new approach through a cross-validation strategy. Predicted values are evaluated at a selected set of times and places for which HLS data are held out.

In addition to generating a comprehensive time-series of data based on prior years, the present disclosure is also applicable to generate in-season predictions, as set out below.

With reference now toFIG. 1, a general process for imputation of remote sensing time series is illustrated according to embodiments of the present disclosure. A first data source101, containing higher resolution, lower frequency data102, is read. In various embodiments, these data comprise one or more agricultural indices, such as the normalized difference vegetation index (NDVI) computed from satellite imagery, such as the Harmonized Landsat Sentinel-2 (HLS). HLS provides medium resolution imagery (approximately 30 m resolution) with low spatial mixing. HLS is available on an approximately weekly basis, and contains data gaps such as clouds and shadows.

A second data source103, containing lower resolution, higher frequency data104, is read. In various embodiments, these data comprise one or more agricultural indices, such as the normalized difference vegetation index (NDVI) computed from satellite imagery, such as the Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS provides low resolution imagery (approximately 250 m resolution) with high spatial mixing. MODIS is available on an approximately daily basis.

In some embodiments, data source101,103contain raw data from which indices may be computed. In such cases, it will be appreciated that indices may be computed before further computation.

All available data for a period of interest is retrieved from each available data source, resulting in temporally overlapping series105,106. Since it comes from the lower frequency data source, series105will contain fewer snapshots. Based on the higher-frequency, lower-resolution series106, and the lower-frequency, higher-resolution series105, a combined series107is generated having both the higher frequency and the higher resolution.

In various embodiments, each series comprises one or more indices for each point in time for which data is available at each pixel in a region of interest. Various indices, for example, normalized difference vegetation index (NDVI), land surface water index (LSWI), and mean brightness (BRT), may be used. Each snapshot in a series is a raster, or image, whose pixel intensity indicates the index value.

Referring toFIG. 2, a method of generating an imputed series is illustrated according to the present disclosure. At201, for each pixel in a region of interest, all available HLS observations for a predetermined period (e.g., a year) are loaded. At202, for each pixel in a region of interest, all available MODIS observations for a predetermined period (e.g., a year) are loaded. In various embodiments, the raw satellite data are loaded and indices are computed, while in some embodiments, precomputed indices are loaded. At203, the MODIS data is gap filled using linear interpolation. At204, the MODIS data are extracted that correspond to the available days of HLS data. At205, an empirical cumulative distribution function (ECDF) is built for HLS and MODIS. At206, the MODIS ECDF is mapped to the HLS ECDF using a linear interpolation. At207, the value of the MODIS observation is scaled up or down using that mapping function. At208, the resulting time series is smoothed to reduce noise.

Referring toFIG. 3, an example of cumulative distribution function (CDF) matching for a given pixel is illustrated. For this pixel, a MODIS NDVI value of0.68will be shifted to 0.82 based on the HLS ECDF. It will be appreciated that while various examples provided herein use an empirical cumulative distribution function (ECDF), alternative methods of determining a cumulative distribution function (CDF) may be used in accordance with the present disclosure.

Referring toFIGS. 4A-B, a graph is provided of NDVI over time as observed and as derived for a given pixel. In this example, the pixel is 10000 in tile 599_361 and the graph is split over two sheets. In particular, the graph shows each observed value as a dot, with HLS corrected MODIS, both smoothed and unsmoothed, shown as solid lines. Running the CDF-matching algorithm took about 500 seconds for one tile/year without smoothing, and about 700 s with smoothing. In this example, every pixel in the tile is computed separately, which is computationally expensive. In alternative embodiments, the data is vectorized, allowing greater efficiency in computation.

Referring toFIG. 5A, raw MODIS NDVI for an exemplary region is shown. InFIG. 5B, the corresponding corrected MODIS NDVI is shown.

Referring toFIG. 6A, raw MODIS NDVI for the same exemplary region is shown for a different date. InFIG. 6B, corresponding raw HLS NDVI is shown. InFIG. 6C, the ECDF-matched MODIS NDVI is shown.

Referring toFIG. 7A, the r2values for HLS versus CDF-matched MODIS NDVI are shown. InFIG. 7B, the corresponding mean error is shown.

Referring now toFIG. 8, a method of generating an imputed series is illustrated according to the present disclosure. At801, for each pixel in a region of interest, all available HLS observations for a predetermined period (e.g., a year) are loaded. At802, for each pixel in a region of interest, all available MODIS observations for a predetermined period (e.g., a year) are loaded. In various embodiments, the raw satellite data are loaded and indices are computed, while in some embodiments, precomputed indices are loaded.

At803, crop labels are loaded. In some embodiments, crop type labels are drawn from the NASS Cropland Data Layer (CDL). In some embodiments, the crop type labels are associated with the same predetermined period (e.g., a year) as the observed data.

At804, a mean annual time series is computed from the HLS data for each crop type in the region of interest. In some embodiments, the HLS data is masked using the crop type labels as masks in order to determine the subsets of the HLS data containing each crop type. At805, the mean annual time series are smoothed, for example using a cubic spline. At806, multiple linear regression is performed using the mean crop type with the MODIS pixel time series as independent variables and observed HLS time series as the dependent variable. At807, the HLS values are predicted using the regression coefficients resulting from the multiple linear regression. At808, a cubic smoothing spline is fit through the predicted HLS observations.

In various embodiments, a library of mean annual time series is generated. Such a library contains a mean annual curve for each of a plurality of crop types for each of a plurality of regions. For example, a curve may be generated for each crop type in each 0.5 degree region of the Earth's surface. The library of curves may be used as described above for determining the predicted HLS value when a crop type is known.

In addition, predicted HLS values may be determined when a crop type is not known. Crop-specific priors are computed for an arbitrary time window using all available HLS NDVI data and associated CDL maps. The observed HLS time series is compared with each crop-specific prior using Manhattan Distance (MD). The crop with the minimum MD is chosen and this time series and MODIS time series are used as independent variables in multivariate linear regression according to Equation 1. In addition to allowing determination of predicted HLS values this allows enables in-season crop detection, where the detected crop type is inferred from the crop-specific prior having the least Manhattan distance, as described above.

Referring toFIGS. 9-12, exemplary mean annual time series are shown for a plurality of exemplary crop types. In these figures, the y-axis corresponds to NDVI value, and the x-axis corresponds to days.

Referring toFIG. 13, exemplary predicted HLS values are shown. Actual HLS observations are shown as large dots, while predicted values are shown as small dots.

Referring toFIG. 14, graph is provided of NDVI over time as observed and as derived for a given pixel and a given crop. In particular, the graph shows observed HLS values as dots, HLS crop mean values as xs, and MODIS as a solid line.

Referring toFIG. 15A, a graph is provided illustrating the correlation between HLS observations of NDVI and mean crop time series. InFIG. 15B, a graph is provided illustrating the correlation between HLS observations of NDVI and co-located MODIS observations.

Referring toFIG. 16, a graph is provided illustrating the observed versus predicted HLS NDVI based on multiple regression coefficients.

Referring toFIG. 17, a graph is provided illustrating the observed versus HLS NDVI time series.

Running the spatial unmixing algorithm on a single analysis tile takes approximately 600-700 seconds. In this example, every pixel in the tile is computed separately, which is computationally expensive. In alternative embodiments, the data is vectorized, allowing greater efficiency in computation.

Referring toFIG. 18A, raw MODIS NDVI for an exemplary region is shown. InFIG. 18B, the corresponding corrected MODIS NDVI is shown.

Referring toFIG. 19, raw corrected NDVI for the same exemplary region is shown for a different date.

Referring toFIG. 20A, the r2values for HLS versus linear unmixing NDVI are shown. InFIG. 20B, the corresponding mean error is shown. Referring toFIG. 21, error bars are shown for r2over multiple exemplary crop types. Referring toFIG. 22, error bars are shown for mean error over multiple exemplary crop types.

In various embodiments, a hybrid approach between the linear unmixing and CDF matching is adopted. The linear unmixing approach described above cleans the crop signal, and provides a high quality approximation of the HLS, yielding both a good visual match and high r2. The approach relies on a CDL mask, and provides gap filling with a spline smoother. CDF-matching, by comparison, does not clean the crop signal. While good at approximating HLS, linear unmixing provides a higher quality output. It does not rely on any masks, and provides only limited gap-filling through the linear gap filling done to the raw MODIS time series.

In an exemplary hybrid approach, CDF matching is applied on MODIS, and then linear unmixing is performed. This approach cleans the original MODIS, improves over the linear unmixing performance, relies on CDL, and provides only the limited gap filling done to the raw MODIS time series. This may be advantageous for both historical products, and for in-season estimates before any crop mask becomes available.

Referring now toFIG. 23, a method of generating an imputed series is illustrated according to the present disclosure. At2301, for each pixel in a region of interest, all available HLS observations for a predetermined period (e.g., a year) are loaded. At2302, for each pixel in a region of interest, all available MODIS observations for a predetermined period (e.g., a year) are loaded. In various embodiments, the raw satellite data are loaded and indices are computed, while in some embodiments, precomputed indices are loaded.

At2303, crop labels are loaded. In some embodiments, crop type labels are drawn from the NASS Cropland Data Layer (CDL). In some embodiments, the crop type labels are associated with the same predetermined period (e.g., a year) as the observed data.

At2304, a mean annual time series is computed from the HLS data for each crop type in the region of interest. In some embodiments, the HLS data is masked using the crop type labels as masks in order to determine the subsets of the HLS data containing each crop type. At2305, the mean annual time series are smoothed, for example using a cubic spline.

At2306, a cross-correlation is computed between an observed HLS time series and the mean crop-specific time series to estimate any time shift resulting from weather or agricultural management practices. At2307, the estimated time shift is applied to the mean crop time series in order to adjust it for weather or agricultural management practices.

At2308, CDF matching is applied on the MODIS data, for example as described above, to rescale the observed MODIS data.

At2309, multiple linear regression is performed using the mean crop type (as time-shifted) with the MODIS pixel time series (as rescaled by CDF matching) as independent variables and observed HLS time series as the dependent variable. At2310, the HLS values are predicted using the regression coefficients resulting from the multiple linear regression. At2311, a cubic smoothing spline is fit through the predicted HLS observations.

It will be appreciated that while the above example includes computation and application of a time-shift, in additional embodiments these steps may be omitted. Likewise, while the above example includes CDF-matching of MODIS data, in additional embodiments, this step may be omitted.

In various embodiments, a library of mean annual time series is generated. Such a library contains a mean annual curve for each of a plurality of crop types for each of a plurality of regions. For example, a curve may be generated for each crop type in each 0.5 degree region of the Earth's surface. The library of curves may be used as described above for determining the predicted HLS value when a crop type is known.

In addition, predicted HLS values may be determined when a crop type is not known. Crop-specific priors are computed for an arbitrary time window using all available HLS NDVI data and associated CDL maps. The observed HLS time series is compared with each crop-specific prior using Manhattan Distance (MD). The crop with the minimum MD is chosen and this time series and MODIS time series are used as independent variables in multivariate linear regression according to Equation 1. In addition to allowing determination of predicted HLS values this allows enables in-season crop detection, where the detected crop type is inferred from the crop-specific prior having the least Manhattan distance, as described above.

Referring toFIG. 24, an exemplary mean annual time series is provided for soybeans, computed as described above.

Referring toFIGS. 25A-I, exemplary mean annual time series are shown for a plurality of exemplary crop types and exemplary observed data. In this example, both the time adjusted (solid) and original (dotted) mean annual time series are shown, which are determined as set forth above. In these figures, the y-axis corresponds to NDVI value, and the x-axis corresponds to days.

Referring toFIGS. 26A-I, the Manhattan distance for each crop type and the exemplary observed data are shown. In this example, the minimum Manhattan distance correspond to soybeans (FIG. 26G), indicating that this is crop reflected in the exemplary data illustrated. In these figures, the y-axis corresponds to NDVI value, and the x-axis corresponds to days.

Referring toFIGS. 27-28, the results of cross-validation are shown. InFIG. 27, the fitting observations of the exemplary data are shown as filled circles. The random sample of HLS observations withheld for accuracy testing are shown as exes. InFIG. 28, the observed values and predicted values (computed as set forth above) are compared to the results of a Savitzky-Golay filter.

Referring now toFIG. 29, an exemplary analysis tile is provided. This tile measures 0.15°×0.15°, and is located in south central Iowa. Scanline artifacts are present in the 2019 monthly NDVI composites. The area is largely corn and soy with neighboring pasture and forest.

FIG. 30illustrates results according to the present disclosure for the tile ofFIG. 29on May 1, 2019. In particular,FIG. 30juxtaposes the raw HLS NDVI values with the values resulting from a Savitzky-Golay filter and those obtained applying the methods described herein.

FIG. 31illustrates results according to the present disclosure for the tile ofFIG. 29on Jul. 8, 2019. In particular,FIG. 31juxtaposes the raw HLS NDVI values with the values resulting from a Savitzky-Golay filter and those obtained applying the methods described herein.

Referring toFIGS. 32-33, cross validation results are shown for the tile ofFIG. 29. In particular,FIG. 32juxtaposes the mean absolute error (MAE) for a Savitzky-Golay filter, the methods described herein, and raw MODIS. The MAE is charted inFIG. 33, where3301corresponds to the methods described herein (denoted Fusion),3302corresponds to a Savitzky-Golay filter (denoted SavGol), and3303corresponds to raw MODIS (denoted MODIS).

As set out above, the present disclosure provides a temporally-driven, pixel-wise fusion approach using a multi-year HLS record. These methods exhibit improved cross-validation results relative to a Savitzky-Golay filter. These improved results enhance data-driven agricultural products such as crop type mapping, compositing, and field boundary delineation.

Referring toFIG. 34, a method agricultural index prediction is illustrated according to embodiments of the present disclosure. At3401, a first time series of raster data is read. The first time series spans a geographic region and has a first resolution and a first frequency. At3402, a second time series of raster data is read. The second time series spans the geographic region and has a second resolution and a second frequency. The second resolution is lower than the first resolution. The second frequency is higher than the first frequency. At3403, a mean time series is determined from the first time series of raster data. At3404, a predicted time series of values for a location within the geographic region is determined at the first resolution by determining a first time series of values for the location from the first time series of raster data, determining a second time series of values of the location from the second time series of raster data, and determining the predicted time series by multiple linear regression with the first time series dependent on the mean time series and the second time series.

Referring now toFIG. 35, a schematic of an example of a computing node is shown. Computing node10is only one example of a suitable computing node and is not intended to suggest any limitation as to the scope of use or functionality of embodiments described herein. Regardless, computing node10is capable of being implemented and/or performing any of the functionality set forth hereinabove.