Calculate timeseries from COGs

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

Aimee Barciauskas, Julia Signell

Published

February 1, 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

Inside the Hub

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.

Outside the Hub

The data is in a protected bucket. Please request access by emailng aimee@developmentseed.org or alexandra@developmentseed.org and providing your affiliation, interest in or expected use of the dataset and an AWS IAM role or user Amazon Resource Name (ARN). The team will help you configure the cognito client.

You should then run:

%run -i 'cognito_login.py'

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 cropped to AOI
  3. Calculate the mean for each timestep over the AOI
from pystac_client import Client
import pandas as pd
import odc.stac

import rioxarray  # noqa
import hvplot.xarray  # noqa

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 = "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.

china_bbox = [
    73.675,
    18.198,
    135.026,
    53.459,
]
datetime = "2000-01-01/2025-07-25"
catalog = Client.open(STAC_API_URL)
search = catalog.search(
    bbox=china_bbox, datetime=datetime, collections=[collection], limit=1000
)
item_collection = search.item_collection()
print(f"Found {len(item_collection)} items")
Found 93 items

Read data

Read in data using xarray using a combination of xpystac, odc-stac, and rasterio.

ds = odc.stac.load(item_collection, chunks={"latitude": "auto", "longitude": "auto"})
da = ds['cog_default']
da = da.where(da != da.nodata)

Clip the data to the bounding box for China

# Subset to Bounding Box for China
subset = da.rio.clip_box(*china_bbox)
subset
<xarray.DataArray 'cog_default' (time: 93, latitude: 354, longitude: 615)> Size: 81MB
dask.array<getitem, shape=(93, 354, 615), dtype=float32, chunksize=(1, 354, 615), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 744B 2016-01-01 2016-02-01 ... 2023-09-01
  * latitude     (latitude) float64 3kB 53.45 53.35 53.25 ... 18.35 18.25 18.15
  * longitude    (longitude) float64 5kB 73.65 73.75 73.85 ... 134.9 135.0 135.1
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  -1.2676506e+30

Aggregate the data

Calculate the mean at each time across regional data. Note this is the first time that the data is actually loaded.

means = subset.mean(dim=("longitude", "latitude")).compute()

Plot the mean monthly NO2 using hvplot

means.hvplot.line(x="time", ylabel="NO2", title="Mean Monthly NO2 in China")