Time-Series Modeling

Create an average NDVI time series within a user-defined AOI

When using remote sensing to monitor activity on the ground, it is often useful to look at the change in the average within an area of interest (AOI) over time. In this example, we will be using the Normalized Difference Vegetation Index (NDVI) to monitor the growing season in an agricultural region of the Pacific Northwest.

First, we setup a Workflows proxy object that will pull Sentinel-2 imagery from January 1, 2017 through June 1, 2020, and use this imagery to calculate NDVI.

import descarteslabs.workflows as wf
import datetime

# Use the dataset id for Sentinel 2
dataset_id = "sentinel-2:L1C"

# Pull imagery from January 1, 2017 to June 1, 2020
img_stack = wf.ImageCollection.from_id(dataset_id,
                                       start_datetime=datetime.date(2017, 1, 1),
                                       end_datetime=datetime.date(2020, 6, 1)
                                       )

# Extract the Red, Near-Infrared, and two data quality bands
red, nir  = img_stack.unpack_bands('red nir')

# Calculate NDVI
ndvi = (nir - red) / (red + nir)

Using the cloud mask bands provided by the ESA, we mask out cloudy pixels and filter out any scenes with a high incidence of clouds.

cirrus_cloud_mask, cloud_mask, alpha = img_stack.unpack_bands('cirrus-cloud-mask cloud-mask alpha')
cloud_mask = ((cirrus_cloud_mask==1) & (cloud_mask==1) & (alpha==1))
ndvi = ndvi.mask(cloud_mask)
ndvi = ndvi.filter(lambda img: img.properties['cloud_fraction'] < 0.3)

With the stack of cloud free images, we group the images and create a weekly maximum NDVI mosaic. This gives us a regular time series of images, where each pixel is represented by the maximum NDVI value observed within a given week. 

Finally, we take the weekly average of all the NDVI values within the AOI. 

# As Sentinel 2 generally has multiple collection opportunities a week, we group 
# the images into 7 day bins, then mosaic the images 

ndvi_stack = ndvi.groupby(
                 lambda img: img.properties['date'] // wf.Timedelta(days=7)
                          ).max(axis="images")

ndvi_stack = ndvi_stack.rename_bands('ndvi')
ndvi_mean = ndvi_stack.mean(axis='pixels')

Using the Workflow outlined above, we can generate a NDVI time series for any AOI we are interested in. To actually run the computation, though, we need to define an AOI.

Irrigated agriculture is common in many parts of the Pacific Northwest, this area outside of Madras, OR produces a wide variety of specialty crops.

Using a GeoJSON of the AOI, we define a Workflows GeoContext. This is the area we will pass to the Workflows proxy object when we run the analysis.

# Define a function to read in a geojson, return a workflows GeoContext
def read_geo(input_file_name):
    with open(input_file_name, 'r') as f:
        input_feature = json.load(f)

    return wf.GeoContext(
        geometry=input_feature['features'][0]['geometry'],
        resolution=30.,
        crs='EPSG:3857',
        bounds_crs='EPSG:4326')

# Load a predefined AOI
aoi = read_geo("Madras_ag.geojson")

Calling ".compute(aoi)" on a Workflows proxy object runs the code we previously defined on the specific area of interest. 

# Calculate the mean NDVI in the AOI, returns a list of dictionaries
ndvi_ts = ndvi_mean.compute(aoi)

# Retrieve the properties for each image in the weekly mosaic 
ndvi_properties = ndvi_stack.properties.compute(aoi)

Once the weekly average NDVI has been calculated for the AOI, we need to add a date indicating the week of each average value. Finally, we create a Pandas DataFrame to hold the time series data and create a plot.

# Add the week to each 
ndvi_ts = [dict(img, 
                **{'date':pd.to_datetime(props['group']).date()}) 
           for img,props in zip(ndvi_ts, ndvi_properties)]


ndvi_data = pd.DataFrame.from_records(ndvi_ts)
ndvi_data = ndvi_data.set_index('date')

# Plot the data
fig, (ax0) = plt.subplots(nrows=1,figsize=(12,6))
plt.plot(ndvi_dat.index, ndvi_dat.ndvi)