Metadata

title: "E-TRAINEE Tutorial - Sentinel-2 data STAC access and time series processing"
description: "This is a tutorial within the second theme of Module 1 of the E-TRAINEE course."
lastUpdate: 2023-07-13
authors: Andreas Mayr

Sentinel-2 data STAC access and time series processing¶

This Notebook illustrates how to query and download Sentinel-2 satellite imagery from Amazon Web Service (AWS) using a SpatioTemporal Asset Catalog (STAC). Furthermore, it shows how to work with the data using xarray (which was introduced in a tutorial in theme 1), and how to perform tasks like subsetting, cloud masking, spectral index calculation, temporal aggregation, and trend analysis. A few simple processing steps let you explore the development of the Normalized Difference Vegetation Index (NDVI) in spring and summer, including e.g. the seasonal greening of mountain grasslands.

The Notebook heavily relies on the Python packages pystac-client, stackstac, and xarray. These packages are contained in the requirements file provided for the course. Please see the instructions on the software page for setting up a Conda environment based on this file.

Parts of the Notebook use code from a tutorial provided here.

Query data via STAC¶

Import packages.

In [3]:

Copied!





import pystac_client
import stackstac
import pyproj
import xarray
import rioxarray
import pystac_client
import stackstac
import pyproj
import xarray
import rioxarray

Define a point of interest to intersect our query.

In [4]:

Copied!

lon, lat = 11.39, 47.26     # Innsbruck (Austria)
lon, lat = 11.39, 47.26     # Innsbruck (Austria)

Use the pystac_client package and element84’s search endpoint to look for items from the sentinel-s2-l2a-cogs collection on Amazon Web Service (AWS). As search criteria we use a maximum cloud cover (per scene) and a time period.

In [5]:

Copied!

URL = "https://earth-search.aws.element84.com/v1"
catalog = pystac_client.Client.open(URL)
URL = "https://earth-search.aws.element84.com/v1"
catalog = pystac_client.Client.open(URL)

In [7]:

Copied!





%%time                                  
items = catalog.search(
    intersects=dict(type="Point", coordinates=[lon, lat]),
    collections=["sentinel-2-l2a"],
    query =  ['eo:cloud_cover<50'],
    datetime="2021-04-01/2021-08-31"
).item_collection()
print(f'{len(items)} items found')      # How many items in the catalog match these criteria?
%%time                                  
items = catalog.search(
    intersects=dict(type="Point", coordinates=[lon, lat]),
    collections=["sentinel-2-l2a"],
    query =  ['eo:cloud_cover<50'],
    datetime="2021-04-01/2021-08-31"
).item_collection()
print(f'{len(items)} items found')      # How many items in the catalog match these criteria?

42 items found
CPU times: total: 141 ms
Wall time: 2.67 s

Use stackstac to turn those STAC items into a lazy xarray. Without specifying anything, this makes use of dask and reads the data in chunks, with one chunk shaped (1, 1, 1024, 1024) in the time, band, y, x dimensions (Read more about chunk sizes here, here, and here.). Using all the defaults, our data will be in its native coordinate reference system, at the finest resolution of all the assets.

In [34]:

Copied!

%time stack = stackstac.stack(items)
%time stack = stackstac.stack(items)

CPU times: total: 78.1 ms
Wall time: 87.7 ms

Now print info about this xarray.DataArray called 'stack'. Note that everything has run quickly so far, despite the large data volume covered by our query. Up to now, we are only working with metadata, no download or processing of any pixel values.

In [35]:

Copied!

stack
stack

Out[35]:

<xarray.DataArray 'stackstac-d35de160ec768a0c51870cec97c80190' (time: 42,
                                                                band: 32,
                                                                y: 10980,
                                                                x: 10980)>
dask.array<fetch_raster_window, shape=(42, 32, 10980, 10980), dtype=float64, chunksize=(1, 1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/52)
  * time                                     (time) datetime64[ns] 2021-04-01...
    id                                       (time) <U24 'S2B_32TPT_20210401_...
  * band                                     (band) <U12 'aot' ... 'wvp-jp2'
  * x                                        (x) float64 6e+05 ... 7.098e+05
  * y                                        (y) float64 5.3e+06 ... 5.19e+06
    s2:medium_proba_clouds_percentage        (time) float64 1.463 16.5 ... 2.701
    ...                                       ...
    title                                    (band) <U31 'Aerosol optical thi...
    gsd                                      (band) object None 10 ... None None
    common_name                              (band) object None 'blue' ... None
    center_wavelength                        (band) object None 0.49 ... None
    full_width_half_max                      (band) object None 0.098 ... None
    epsg                                     int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

We are only interested in a relatively small area around our main point of interest (defined above for intersection). So we convert the coordinates of this point into the native CRS of the satellite data, add a buffer and then use this rectangle to slice ("clip") the stacked array spatially.

In [36]:

Copied!

x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)     # Re-project the geographic point coordinates to the sat data CRS (i.e. EPSG:32632)
buffer = 5000  # meters
x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)     # Re-project the geographic point coordinates to the sat data CRS (i.e. EPSG:32632)
buffer = 5000  # meters

Now we have a much smaller DataArray (with 1000 by 1000 cells in x/y) called 'stack_aoi'.

In [37]:

Copied!

stack_aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]
stack_aoi
stack_aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]
stack_aoi

Out[37]:

<xarray.DataArray 'stackstac-d35de160ec768a0c51870cec97c80190' (time: 42,
                                                                band: 32,
                                                                y: 1000, x: 1000)>
dask.array<getitem, shape=(42, 32, 1000, 1000), dtype=float64, chunksize=(1, 1, 678, 610), chunktype=numpy.ndarray>
Coordinates: (12/52)
  * time                                     (time) datetime64[ns] 2021-04-01...
    id                                       (time) <U24 'S2B_32TPT_20210401_...
  * band                                     (band) <U12 'aot' ... 'wvp-jp2'
  * x                                        (x) float64 6.758e+05 ... 6.858e+05
  * y                                        (y) float64 5.242e+06 ... 5.232e+06
    s2:medium_proba_clouds_percentage        (time) float64 1.463 16.5 ... 2.701
    ...                                       ...
    title                                    (band) <U31 'Aerosol optical thi...
    gsd                                      (band) object None 10 ... None None
    common_name                              (band) object None 'blue' ... None
    center_wavelength                        (band) object None 0.49 ... None
    full_width_half_max                      (band) object None 0.098 ... None
    epsg                                     int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

We want only the four 10-m bands and the Scene Classification (SCL) layer.

In [39]:

Copied!

stack_aoi = stack_aoi.sel(band = ['blue', 'green', 'red', 'nir', 'scl'])
stack_aoi
stack_aoi = stack_aoi.sel(band = ['blue', 'green', 'red', 'nir', 'scl'])
stack_aoi

Out[39]:

<xarray.DataArray 'stackstac-d35de160ec768a0c51870cec97c80190' (time: 42,
                                                                band: 5,
                                                                y: 1000, x: 1000)>
dask.array<getitem, shape=(42, 5, 1000, 1000), dtype=float64, chunksize=(1, 1, 678, 610), chunktype=numpy.ndarray>
Coordinates: (12/52)
  * time                                     (time) datetime64[ns] 2021-04-01...
    id                                       (time) <U24 'S2B_32TPT_20210401_...
  * band                                     (band) <U12 'blue' ... 'scl'
  * x                                        (x) float64 6.758e+05 ... 6.858e+05
  * y                                        (y) float64 5.242e+06 ... 5.232e+06
    s2:medium_proba_clouds_percentage        (time) float64 1.463 16.5 ... 2.701
    ...                                       ...
    title                                    (band) <U31 'Blue (band 2) - 10m...
    gsd                                      (band) object 10 10 10 10 None
    common_name                              (band) object 'blue' ... None
    center_wavelength                        (band) object 0.49 0.56 ... None
    full_width_half_max                      (band) object 0.098 0.045 ... None
    epsg                                     int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

Now we trigger the stacking of these selected bands for our area-of-interest by calling compute(), i.e. we load the data into memory. This may take some minutes but then enables us to analyse and plot the data faster. Alternatively (also for larger datasets), we could continue with the lazy array and in subsequent steps only load what is needed.

In [40]:

Copied!

%%time
data = stack_aoi.compute()
%%time
data = stack_aoi.compute()

CPU times: total: 56.2 s
Wall time: 5min 41s

Try plotting some data (July 2021 only):

In [43]:

Copied!

data.sel(band='nir').sel(time='2021-07').plot.imshow(col='time')
data.sel(band='nir').sel(time='2021-07').plot.imshow(col='time')

Out[43]:

<xarray.plot.facetgrid.FacetGrid at 0x16ede375840>

No description has been provided for this image

Cloud masking¶

We use the scene classification provided by ESA via the SCL band as a simple (but not the best) means to mask clouds and other problematic pixels.

In [44]:

Copied!





import dask.array as da

scl = data.sel(band=["scl"])
# Sentinel-2 Scene Classification Map: nodata, saturated/defective, dark, cloud shadow, cloud med. prob., cloud high prob., cirrus
invalid_data = da.isin(scl, [0, 1, 2, 3, 8, 9, 10])
valid_data = data.where(~invalid_data)
# Show some of the cloud masked data
valid_data.sel(band='nir').sel(time='2021-07').plot.imshow(col='time', robust=True)
import dask.array as da

scl = data.sel(band=["scl"])
# Sentinel-2 Scene Classification Map: nodata, saturated/defective, dark, cloud shadow, cloud med. prob., cloud high prob., cirrus
invalid_data = da.isin(scl, [0, 1, 2, 3, 8, 9, 10])
valid_data = data.where(~invalid_data)
# Show some of the cloud masked data
valid_data.sel(band='nir').sel(time='2021-07').plot.imshow(col='time', robust=True)

Out[44]:

<xarray.plot.facetgrid.FacetGrid at 0x16eddff1960>

Let's show cloud-masked RGB true color composites.

In [45]:

Copied!

valid_data.sel(time='2021-07').sel(band=['red', 'green', 'blue']).plot.imshow(col='time', robust=True)
valid_data.sel(time='2021-07').sel(band=['red', 'green', 'blue']).plot.imshow(col='time', robust=True)

Out[45]:

<xarray.plot.facetgrid.FacetGrid at 0x16e84e5bb20>

$No description has been provided for this image$

Spectral index calculation (NDVI)¶

We calculate the Normalized Difference Vegetation Index (NDVI). If you need more indices consider the Spyndex package.

In [46]:

Copied!

# Compute the NDVI from the NIR (B8) and red (B4) bands
nir, red = valid_data.sel(band='nir'), valid_data.sel(band='red')
NDVI = (nir - red) / (nir + red)
# Compute the NDVI from the NIR (B8) and red (B4) bands
nir, red = valid_data.sel(band='nir'), valid_data.sel(band='red')
NDVI = (nir - red) / (nir + red)

NDVI is a new DataArray. If you want to select the NDVI raster closest to a specific date and export it as a GeoTIFF, just run this line:

In [47]:

Copied!

NDVI.sel(time='2021-07-15', method='nearest').rio.to_raster("C:\work\etrainee\stac_exports/S2-NDVI_midjuly.tif")
NDVI.sel(time='2021-07-15', method='nearest').rio.to_raster("C:\work\etrainee\stac_exports/S2-NDVI_midjuly.tif")

Let's plot one month of NDVI rasters.

In [48]:

Copied!

NDVI.sel(time='2021-07').plot.imshow(col='time', cmap="PRGn", robust=True)
NDVI.sel(time='2021-07').plot.imshow(col='time', cmap="PRGn", robust=True)

Out[48]:

<xarray.plot.facetgrid.FacetGrid at 0x16e83c1f010>

Time series operations¶

Interpolate missing values¶

Let's fill missing values by (linear) interpolation between observations.

In [49]:

Copied!

NDVI.sel(time='2021-07').interpolate_na(dim='time').plot(col='time', cmap="PRGn", robust=True)
NDVI.sel(time='2021-07').interpolate_na(dim='time').plot(col='time', cmap="PRGn", robust=True)

Out[49]:

<xarray.plot.facetgrid.FacetGrid at 0x16e80ce0be0>

Resampling¶

We resample/aggregate to one month temporal resolution by computing the mean for each month. Other aggregations (e.g. maximum) can be achieved similarly. Note how we can absolutely define the range of values to be displayed (important for visual comparisons).

In [50]:

Copied!

NDVI.resample(time='1M').mean().plot(col='time', cmap="PRGn", vmin=-0.2, vmax=0.8)
NDVI.resample(time='1M').mean().plot(col='time', cmap="PRGn", vmin=-0.2, vmax=0.8)

Out[50]:

<xarray.plot.facetgrid.FacetGrid at 0x16e85458d30>

Trend analysis¶

The xarrayutils package contains a convenient linear_trend() function. We use this to fit a linear least-squares trendline to our NDVI time series.

In [51]:

Copied!

from xarrayutils import linear_trend
NDVI_regressed = linear_trend(NDVI, 'time')
from xarrayutils import linear_trend
NDVI_regressed = linear_trend(NDVI, 'time')

Plot maps (with default colormaps) for the resulting slope of the estimated linear trend along with r-value (Pearson correlation coefficient), p-value (for a hypothesis test whose null hypothesis is that the slope is zero), and standard error (of the estimated slope). We use matplotlib syntax and functions to configure the plot.

In [52]:

Copied!





from matplotlib import pyplot as plt
fig, axes = plt.subplots(2,2, figsize=(15,10), sharex=True, sharey=True)
NDVI_regressed.slope.plot(robust=True, ax=axes[0,0])
NDVI_regressed.r_value.plot(robust=True, ax=axes[0,1])
NDVI_regressed.p_value.plot(robust=True, ax=axes[1,0])
NDVI_regressed.std_err.plot(robust=True, ax=axes[1,1])
for ax, param in zip(axes.flat, ['slope', 'r-value', 'p-value', 'std error']):
    ax.set_title(f'NDVI trend 04-08/2021: {param}')
fig.tight_layout()
from matplotlib import pyplot as plt
fig, axes = plt.subplots(2,2, figsize=(15,10), sharex=True, sharey=True)
NDVI_regressed.slope.plot(robust=True, ax=axes[0,0])
NDVI_regressed.r_value.plot(robust=True, ax=axes[0,1])
NDVI_regressed.p_value.plot(robust=True, ax=axes[1,0])
NDVI_regressed.std_err.plot(robust=True, ax=axes[1,1])
for ax, param in zip(axes.flat, ['slope', 'r-value', 'p-value', 'std error']):
    ax.set_title(f'NDVI trend 04-08/2021: {param}')
fig.tight_layout()

c:\Users\Andi\mambaforge\envs\etrainee_m1\lib\site-packages\dask\core.py:121: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))
c:\Users\Andi\mambaforge\envs\etrainee_m1\lib\site-packages\dask\array\numpy_compat.py:43: RuntimeWarning: invalid value encountered in divide
  x = np.divide(x1, x2, out)
c:\Users\Andi\mambaforge\envs\etrainee_m1\lib\site-packages\dask\core.py:121: RuntimeWarning: divide by zero encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

Create an interactive map to inspect the slope of the trend in more detail.

In [53]:

Copied!

import hvplot.xarray
import holoviews as hv
NDVI_regressed.slope.hvplot(x='x', y='y', aspect=1, clim=(-0.02,0.02), cmap="PRGn", title='NDVI slope April - Aug')
import hvplot.xarray
import holoviews as hv
NDVI_regressed.slope.hvplot(x='x', y='y', aspect=1, clim=(-0.02,0.02), cmap="PRGn", title='NDVI slope April - Aug')

c:\Users\Andi\mambaforge\envs\etrainee_m1\lib\site-packages\dask\array\numpy_compat.py:43: RuntimeWarning: invalid value encountered in divide
  x = np.divide(x1, x2, out)
c:\Users\Andi\mambaforge\envs\etrainee_m1\lib\site-packages\dask\core.py:121: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

Out[53]:

Create an animation¶

Let's create an animation of our NDVI time series using the xmovie package. We resample to weekly means of NDVI and interpolate missing values.

In [54]:

Copied!





%%time
NDVI_weekly = NDVI.resample(time='1W').mean()
NDVI_weekly = NDVI_weekly.chunk(dict(time=-1))                        # rechunk into a single dask array chunk along the time dimension
NDVI_weekly = NDVI_weekly.interpolate_na(dim='time')
%%time
NDVI_weekly = NDVI.resample(time='1W').mean()
NDVI_weekly = NDVI_weekly.chunk(dict(time=-1))                        # rechunk into a single dask array chunk along the time dimension
NDVI_weekly = NDVI_weekly.interpolate_na(dim='time')

CPU times: total: 109 ms
Wall time: 261 ms

Now be patient while the movie is created.

In [55]:

Copied!

%%time
from xmovie import Movie
mov = Movie(NDVI_weekly.chunk({'time':1}), dpi=100, vmin=-0.8, vmax=0.8, cmap="PRGn")
%%time
from xmovie import Movie
mov = Movie(NDVI_weekly.chunk({'time':1}), dpi=100, vmin=-0.8, vmax=0.8, cmap="PRGn")

CPU times: total: 1min 12s
Wall time: 1min 16s

Use parallelization by Dask for saving the frames to an MP4 video clip (or a GIF). For this to work we have chunked the input DataArray with a single step along the dimension that represent our frames. Again be patient ... (Code cell doesn't stop running even after video clip is saved!?)

In [56]:

Copied!

%%time
mov.save(r'C:\work\etrainee\stac_exports\NDVI_animation.mp4', overwrite_existing=True, parallel=True, framerate=3)
%%time
mov.save(r'C:\work\etrainee\stac_exports\NDVI_animation.mp4', overwrite_existing=True, parallel=True, framerate=3)