We propose to showcase the capabilities of the recently created data catalog at the Computer Science department, West University of Timișoara. The catalog exposes Sentinel-2 imagery (currently) using STAC Catalogs.

We store Sentinel-2 L1C products from 2017-01-01 till present, covering the whole territory of Romania. This post will illustrate the usage of stackstac in order to turn a STAC Collection into an xarray dataframe that is later processed using dask showing a timelapse of some big cities in Romania throughout the years. The work presented is based on the stackstac example.

The approach for generating the timelapse GIF’s is rather simple:

1. Generate a bounding box from national codes (aka. SIRUTA in Romania) using an OWS feature containing Romania’s administrative territorial bounds from geo-spatial.org.
2. Query the data intersecting with said bounding box using stackstac, selecting the RGB bands for Sentinel-2 (B02, B03 and B04) and convert the data to an xarray array.
3. Filter scenes based on cloud coverage (this metadata is available in the STAC catalog). For the examples listed here, I filtered the images so that all of them have less than 1% cloud coverage.
4. Resample the xarray array on the time dimension at each 10 days.
5. Crop the array at a given buffer (in meters) from the centroid of the bounding box.
6. Apply gamma correction in order to create True Color Image (TCI) composites with a gamma value of 2.2 (Results shown in a paper published by Sinergise show that this is a value that works well with Sentinel-2 data).
7. Use dask and geogif to render the timelapse GIF.

So, how do we do this with Python?

Simply start off by interrogating the STAC catalog:

1
2
3
4
5
6
7
8
9

url = 'https://registry.sage.uvt.ro/geoserver/ogc/stac?f=application%2Fgeo%2Bjson'

catalog = pystac_client.Client.open(url)

result = catalog.search(
    bbox=bounds_wgs84,
    collections=["S2L1C"],
    datetime='2017-01-12T00:00:00Z/2022-05-31T12:31:12Z',
)

Then, we ask stackstac.stack() to create our array:

1

stack = stackstac.stack(items, epsg=32635, resolution=10, assets=["B04", "B03", "B02"], errors_as_nodata=Exception(".*"))

This is what our datacube looks like so far:

We can then start filtering scenes based on cloud cover and resample the data

1
2
3

stack = stack[stack["eo:cloud_cover"] < 1]

composites = stack.resample(time="10D").median("time")

Afterwards, convert the centroid lat-lon coordinates to UTM using proj and crop the array to our desired Area of Interest:

1
2
3
4

x_utm, y_utm = pyproj.Proj(stack.crs)(lat, lon)
buffer = 5000 # meters

aoi = composites.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]

We can use dask.array.power() to apply gamma correction to the cube:

1
2

def gamma_correction(xarr, gamma):
    return da.power(xarr, 1/gamma)

And finally, use dask and geogif to produce our GIF’s

1
2
3

with dask.diagnostics.ProgressBar():
    gif_img = geogif.dgif(aoi, fps=7, date_color=(255,255,255), date_bg=(0,10,0)).compute()
    dask.utils.format_bytes(len(gif_img.data))

These are the results:

Special thanks to my colleague, Marian Neagul for the support and guidance, to geo-spatial.org for providing the territorial administrative bounds, and to the Romanian Space Agency for supplying Jupyter Hub environments.