Get timeseries from COGs

Demonstrates gathering a timeseries of statistics for a given area using the Raster API

Author

Leo Thomas, Julia Signell

Published

February 1, 2023

Run this notebook

You can launch this notbook using mybinder, by clicking the button below.

Approach

Using a list of STAC items and a bouding box fetch stats from /cog/statistics endpoint
Generate a timeseries plot using statistics from each time step
Speed up workflow using Dask

import requests
import folium
import matplotlib.pyplot as plt

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://staging-stac.delta-backend.com/collections
STAC Browser: http://veda-staging-stac-browser.s3-website-us-west-2.amazonaws.com

STAC_API_URL = "https://staging-stac.delta-backend.com"
RASTER_API_URL = "https://staging-raster.delta-backend.com"

collection_name = "no2-monthly"

Fetch STAC collection

We will use requests to fetch all the metadata about the collection of interest from STAC.

collection = requests.get(f"{STAC_API_URL}/collections/{collection_name}").json()
collection

{'id': 'no2-monthly',
 'type': 'Collection',
 'links': [{'rel': 'items',
   'type': 'application/geo+json',
   'href': 'https://staging-stac.delta-backend.com/collections/no2-monthly/items'},
  {'rel': 'parent',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/'},
  {'rel': 'root',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/'},
  {'rel': 'self',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/collections/no2-monthly'}],
 'title': 'NO₂',
 'assets': None,
 'extent': {'spatial': {'bbox': [[-180, -90, 180, 90]]},
  'temporal': {'interval': [['2016-01-01 00:00:00+00',
     '2023-09-30 00:00:00+00']]}},
 'license': 'MIT',
 'keywords': None,
 'providers': None,
 'summaries': {'datetime': ['2016-01-01T00:00:00Z', '2023-09-30T00:00:00Z']},
 'description': 'Darker colors indicate higher nitrogen dioxide (NO₂) levels and more activity. Lighter colors indicate lower levels of NO₂ and less activity. Missing pixels indicate areas of no data most likely associated with cloud cover or snow.',
 'item_assets': {'cog_default': {'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'roles': ['data', 'layer'],
   'title': 'Default COG Layer',
   'description': 'Cloud optimized default layer to display on map'}},
 'stac_version': '1.0.0',
 'stac_extensions': None,
 'dashboard:is_periodic': True,
 'dashboard:time_density': 'month'}

Describe the periodic nature of the data

In the collection above we will pay particular attention to the fields that define the periodicity of the data.

collection["dashboard:is_periodic"]

True

collection["dashboard:time_density"]

'month'

collection["summaries"]

{'datetime': ['2016-01-01T00:00:00Z', '2023-09-30T00:00:00Z']}

Fetch STAC items

Get the list of all the STAC items within this collection.

items = requests.get(
    f"{STAC_API_URL}/collections/{collection_name}/items?limit=100"
).json()["features"]
len(items)

We can inspect one of these items to get a sense of what metadata is available.

items[0]

{'id': 'OMI_trno2_0.10x0.10_202309_Col3_V4.nc',
 'bbox': [-180.0, -90.0, 180.0, 90.0],
 'type': 'Feature',
 'links': [{'rel': 'collection',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/collections/no2-monthly'},
  {'rel': 'parent',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/collections/no2-monthly'},
  {'rel': 'root',
   'type': 'application/json',
   'href': 'https://staging-stac.delta-backend.com/'},
  {'rel': 'self',
   'type': 'application/geo+json',
   'href': 'https://staging-stac.delta-backend.com/collections/no2-monthly/items/OMI_trno2_0.10x0.10_202309_Col3_V4.nc'}],
 'assets': {'cog_default': {'href': 's3://veda-data-store-staging/no2-monthly/OMI_trno2_0.10x0.10_202309_Col3_V4.nc.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'roles': ['data', 'layer'],
   'title': 'Default COG Layer',
   'description': 'Cloud optimized default layer to display on map',
   'raster:bands': [{'scale': 1.0,
     'nodata': -1.2676506002282294e+30,
     'offset': 0.0,
     'sampling': 'area',
     'data_type': 'float32',
     'histogram': {'max': 11639639471292416,
      'min': -4024247454269440.0,
      'count': 11.0,
      'buckets': [24.0,
       1128.0,
       403168.0,
       64788.0,
       5583.0,
       1096.0,
       252.0,
       64.0,
       19.0,
       4.0]},
     'statistics': {'mean': 343251029873510.25,
      'stddev': 563904505347325.8,
      'maximum': 11639639471292416,
      'minimum': -4024247454269440.0,
      'valid_percent': 90.81382751464844}}]}},
 'geometry': {'type': 'Polygon',
  'coordinates': [[[-180, -90],
    [180, -90],
    [180, 90],
    [-180, 90],
    [-180, -90]]]},
 'collection': 'no2-monthly',
 'properties': {'proj:bbox': [-180.0, -90.0, 180.0, 90.0],
  'proj:epsg': 4326.0,
  'proj:shape': [1800.0, 3600.0],
  'end_datetime': '2023-09-30T00:00:00',
  'proj:geometry': {'type': 'Polygon',
   'coordinates': [[[-180.0, -90.0],
     [180.0, -90.0],
     [180.0, 90.0],
     [-180.0, 90.0],
     [-180.0, -90.0]]]},
  'proj:transform': [0.1, 0.0, -180.0, 0.0, -0.1, 90.0, 0.0, 0.0, 1.0],
  'start_datetime': '2023-09-01T00:00:00'},
 'stac_version': '1.0.0',
 'stac_extensions': ['https://stac-extensions.github.io/projection/v1.0.0/schema.json',
  'https://stac-extensions.github.io/raster/v1.1.0/schema.json']}

Define an area of interest

We will be using a bounding box over metropolitan france. We’ll use that bounding box to subset the data when calculating the timeseries.

france_bounding_box = {
    "type": "Feature",
    "properties": {"ADMIN": "France", "ISO_A3": "FRA"},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [-5.183429, 42.332925],
                [8.233998, 42.332925],
                [8.233998, 51.066135],
                [-5.183429, 51.066135],
                [-5.183429, 42.332925],
            ]
        ],
    },
}

Let’s take a look at that box.

m = folium.Map(
    location=[40, 0],
    zoom_start=3,
)

folium.GeoJson(france_bounding_box, name="France").add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Use `/cog/statistics` to get data for the AOI

First, we create a generate_stats function and then we call it with the bounding box defined for France.

# the bounding box should be passed to the geojson param as a geojson Feature or FeatureCollection
def generate_stats(item, geojson):
    result = requests.post(
        f"{RASTER_API_URL}/cog/statistics",
        params={"url": item["assets"]["cog_default"]["href"]},
        json=geojson,
    ).json()
    return {
        **result["properties"],
        "start_datetime": item["properties"]["start_datetime"],
    }

Generate stats

This may take some time due to the complexity of the shape we’re requesting. See the end of this notebook for tips on how to speed this up.

%%time
stats = [generate_stats(item, france_bounding_box) for item in items]

CPU times: user 1.21 s, sys: 50.7 ms, total: 1.26 s
Wall time: 52.7 s

Inspect one result

stats[0]

{'statistics': {'b1': {'min': -387779578036224.0,
   'max': 5971174383157248.0,
   'mean': 1633500418513795.8,
   'count': 11658.0,
   'sum': 1.904334787903383e+19,
   'std': 807808763419533.5,
   'median': 1504613211570176.0,
   'majority': 1143241374171136.0,
   'minority': -387779578036224.0,
   'unique': 11653.0,
   'histogram': [[118.0,
     1594.0,
     4202.0,
     3374.0,
     1350.0,
     613.0,
     233.0,
     122.0,
     40.0,
     12.0],
    [-387779578036224.0,
     248115814727680.0,
     884011241046016.0,
     1519906650587136.0,
     2155802060128256.0,
     2791697335451648.0,
     3427592744992768.0,
     4063488154533888.0,
     4699383564075008.0,
     5335279242051584.0,
     5971174383157248.0]],
   'valid_percent': 100.0,
   'masked_pixels': 0.0,
   'valid_pixels': 11658.0,
   'percentile_2': 395805562618511.4,
   'percentile_98': 3845241758988372.5}},
 'ADMIN': 'France',
 'ISO_A3': 'FRA',
 'start_datetime': '2023-09-01T00:00:00'}

Plot timeseries

It is easier to interact with these results as a pandas dataframe. The following function takes the json, passes it to pandas, cleans up the column names and parses the date column.

import pandas as pd


def clean_stats(stats_json) -> pd.DataFrame:
    df = pd.json_normalize(stats_json)
    df.columns = [col.replace("statistics.b1.", "") for col in df.columns]
    df["date"] = pd.to_datetime(df["start_datetime"])
    return df


df = clean_stats(stats)

Construct the plot

fig = plt.figure(figsize=(20, 10))

plt.plot(df["date"], df["mean"], "black", label="Mean monthly NO2 values")

plt.fill_between(
    df["date"],
    df["mean"] + df["std"],
    df["mean"] - df["std"],
    facecolor="lightgray",
    interpolate=False,
    label="+/- one standard devation",
)

plt.plot(
    df["date"],
    df["min"],
    color="blue",
    linestyle="-",
    linewidth=0.5,
    label="Min monthly NO2 values",
)
plt.plot(
    df["date"],
    df["max"],
    color="red",
    linestyle="-",
    linewidth=0.5,
    label="Max monhtly NO2 values",
)

plt.legend()
plt.title("NO2 Values in France (2016-2022)")

Text(0.5, 1.0, 'NO2 Values in France (2016-2022)')

In this graph we can see the yearly cycles in NO2 values due to seasonal variations, as well as a slight downward slope in maximum NO2 values

Complex AOI

The values plotted above don’t correspond exactly to Fance, since the bounding box excludes Corsica and overseas territories such as Mayotte and French Polynesia, and covers parts of neighboring countries including Spain, Italy, Germany and the entirety of Luxembourg. We can fetch GeoJSON from an authoritative online source (https://gadm.org/download_country.html).

While the NO2 values above correspond more or less to those of in France, we can be much more precise by using a complex geojson that represents the boundaries of France exactly, including overseas territories in the Carribean and Indian Oceans, and South America.

Note: In this notebook we write out the whole perimeter as a MultiPolygon in geojson. In practice you will often be reading this kind of shape from a file (usually with the help of geopandas).

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]

Let’s take a look at this AOI on a map

m = folium.Map(
    location=[40, 0],
    zoom_start=3,
)

folium.GeoJson(france_aoi, name="France").add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

We can now request the NO2 values for this complex AOI the same way as for the bounding box.

Notice, however, that due to the complexity of the shape, it takes much longer to gather the requested data about 4 times as long as for the bounding box example above.

%%time
aoi_df = clean_stats([generate_stats(item, france_aoi) for item in items])

CPU times: user 4.08 s, sys: 97.2 ms, total: 4.17 s
Wall time: 1min 36s

We can compare the mean monthly NO2 values calculated when using the bounding box and when using the country’s exact borders

fig = plt.figure(figsize=(20, 10))

plt.plot(
    df["date"],
    df["mean"],
    color="blue",
    label="Mean monthly NO2 values using bounding box",
)
plt.plot(
    aoi_df["date"],
    aoi_df["mean"],
    color="red",
    label="Mean monthly NO2 values using complex AOI",
)

plt.legend()
plt.title("NO2 Values in France (2016-2022)")

Text(0.5, 1.0, 'NO2 Values in France (2016-2022)')

While the difference is small, it is very interesting to note that the NO2 values calculated using the exact borders are systematically less than when using the bounding box. This may be due to the fact that the bounding box includes parts of western Germany and northern Italy that have a lot industrial activity, whereas the areas included when using the exact borders that are not included in the bounding box case, are overseas territories much less industrial activity.

Speed things up: parallelize computation with Dask

We can drastically reduce the time it takes to generate the timeseries, even with the complex AOI above, by parallelizing our code. The cogs/statistics API is powered by AWS Lambda which executes each request in a separate instance. This means the requests are highly scalable. Since each statistics request is for a single timestamp, we can request statistics for multiple timesteps concurrently, and greatly reduce the amount of time needed. We will demonstrate this by using the Dask.

Create a Dask client

First we will create a Dask client. In this case we will use the threads on the same server that is running this jupyter notebook.

import dask.distributed

client = dask.distributed.Client()

Submit work

We will submit the generate_stats function for each item in our list and collect a list of futures. This will immediately kick off work in dask. We can then gather all the results.

%%time
futures = [client.submit(generate_stats, item, france_aoi) for item in items]
stats = client.gather(futures)

CPU times: user 19.4 s, sys: 338 ms, total: 19.7 s
Wall time: 34 s

Close the Dask client

It is good practice to close the client when you are done.

client.shutdown()

Alternate approach

If you are familiar with the concurrent.futures library you can use that instead of Dask.