Open and visualize COGs

Demonstrates how to visualize COGs using pystac_client, rioxarray, odc-stac, and hvplot
Author

Julia Signell

Published

February 2, 2023

Run this notebook

You can launch this notebook in VEDA JupyterHub by clicking the link below.

Launch in VEDA JupyterHub (requires access)

Learn more

This notebook was written on the VEDA JupyterHub and as such is designed to be run on a jupyterhub which is associated with an AWS IAM role which has been granted permissions to the VEDA data store via its bucket policy. The instance used provided 16GB of RAM.

See (VEDA Analytics JupyterHub Access)[https://nasa-impact.github.io/veda-docs/veda-jh-access.html] for information about how to gain access.

Approach

  1. Use pystac_client to open the STAC catalog and retrieve the items in the collection
  2. Use odc-stac to create an xarray dataset containing all the items
  3. Use rioxarray to crop data to AOI
  4. Use hvplot to render the COG at every timestep
import hvplot.xarray  # noqa
import numpy as np
import odc.stac
import requests
import rioxarray  # noqa
from pystac_client import Client

Declare your collection of interest

You can discover available collections the following ways:

  • Programmatically: see example in the list-collections.ipynb notebook
  • JSON API: https://openveda.cloud/api/stac/collections
  • STAC Browser: https://openveda.cloud
STAC_API_URL = "https://openveda.cloud/api/stac"
collection_id = "no2-monthly"

Discover items in collection for region and time of interest

Use pystac_client to search the STAC collection for a particular area of interest within specified datetime bounds.

catalog = Client.open(STAC_API_URL)
search = catalog.search(collections=[collection_id], sortby="start_datetime")

item_collection = search.item_collection()
print(f"Found {len(item_collection)} items")
Found 93 items

Define an AOI

We can fetch GeoJSON for metropolitan France and Corsica (excluding overseas territories) from an authoritative online source (https://gadm.org/download_country.html).

response = requests.get("https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_FRA_0.json")

# If anything goes wrong with this request output error contents
assert response.ok, response.text

result = response.json()
print(f"There are {len(result['features'])} features in this collection")
There are 1 features in this collection

That is the geojson for a feature collection, but since there is only one feature in it we can grab just that.

france_aoi = result["features"][0]

Read data

Create an xarray.DataArray using odc-stac. Geopolygon being set will result in the dataset / data array only being loaded bounded by the bbox of the geopolygon.

ds = odc.stac.load(
    item_collection,
    chunks={"latitude": "auto", "longitude": "auto", "time": "auto"},
    nodata=np.nan,
    geopolygon=france_aoi,
)
da = ds["cog_default"]
da
<xarray.DataArray 'cog_default' (time: 93, latitude: 98, longitude: 148)> Size: 5MB
dask.array<cog_default, shape=(93, 98, 148), dtype=float32, chunksize=(93, 98, 148), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 744B 2016-01-01 2016-02-01 ... 2023-09-01
  * latitude     (latitude) float64 784B 51.05 50.95 50.85 ... 41.55 41.45 41.35
  * longitude    (longitude) float64 1kB -5.15 -5.05 -4.95 ... 9.35 9.45 9.55
    spatial_ref  int32 4B 4326
Attributes:
    nodata:   nan

Clip the data to AOI

This will mask the raster by the geometry shape itself

subset = da.rio.clip([france_aoi["geometry"]])
subset
<xarray.DataArray 'cog_default' (time: 93, latitude: 97, longitude: 144)> Size: 5MB
dask.array<getitem, shape=(93, 97, 144), dtype=float32, chunksize=(93, 97, 144), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 744B 2016-01-01 2016-02-01 ... 2023-09-01
  * latitude     (latitude) float64 776B 51.05 50.95 50.85 ... 41.65 41.55 41.45
  * longitude    (longitude) float64 1kB -4.75 -4.65 -4.55 ... 9.35 9.45 9.55
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

Compute and plot

So far we have just been setting up a calculation lazily in Dask. Now we can trigger computation using .compute().

%%time

image_stack = subset.compute()
CPU times: user 7.83 s, sys: 1.25 s, total: 9.08 s
Wall time: 25.9 s
# get the 2% and 98% percentiles for min and max bounds of color
vmin, vmax = (
    image_stack.quantile(0.02).item(),
    image_stack.quantile(0.98).item(),
)

image_stack.hvplot.quadmesh(
    groupby="time",
    tiles=True,
    colorbar=False,
    clim=(vmin, vmax),
    cmap="viridis",
    alpha=0.8,
    frame_height=512,
    widget_location="bottom",
    aspect=1,
)