STAC Tutorial – OpenPlains

In this tutorial we will learn about the t.stac suite of tools for GRASS.

t.stac.catalog - is a tool for exploring SpatioTemporal Asset Catalogs metadata from a STAC API.
t.stac.collection - is a tool for exploring SpatioTemporal Asset Catalog (STAC) collection metadata.
t.stac.item - is a tool for exploring and importing SpatioTemporal Asset Catalog item metadata and assets into GRASS.

Introduction

SpatioTemporal Asset Catalog (STAC)…

The tutorial assumes you already have GRASS installed on your machine. If not please download and install GRASS before continuing the tutorial.

Getting Started

Start GRASS Session

# Import standard python packages
import sys
import subprocess
from IPython.display import display, JSON, HTML
import json
import pandas as pd
from pathlib import Path
import seaborn as sns
# import polars as pl

# Ask GRASS where its Python packages are and add them to the path
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

# Import the GRASS python packages we need
import grass.script as gs
import grass.jupyter as gj

Create a GRASS project

gs.create_project("stac", epsg="32119", overwrite=True)

and start a GRASS session.

session = gj.init("stac")

Install Addon

The t.stac addons require the following Python packages.

You can install them using pip or the Python package manager of your choice.

pip install pystac pystac_client tqdm

Now you can install the t.stac extensions.

gs.run_command("g.extension", extension="t.stac")

Define computational region

gs.run_command(
    "g.region", n=236687, s=210391, e=616042, w=598921, nsres=10, ewres=10, flags="pa"
)

STAC API

stac_url = "https://earth-search.aws.element84.com/v1/"

Table 1

catalogs_json = gs.parse_command(
    "t.stac.catalog", url=stac_url, format="json"
)

df_catalogs = pd.json_normalize(catalogs_json)
df_catalogs.head()

Table 2

!t.stac.catalog url={stac_url} format=plain -b

Collection Metadata

Python
Terminal

Table 3

items_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    flags="m",
    format="json",
)

df_items = pd.json_normalize(items_json, max_level=0)
display(df_items.T.head(5))

Table 4

!t.stac.item url={stac_url} collection_id="sentinel-2-l2a" format=plain -m

Query Items by datetime and display table

Python
Terminal

Table 5: Sentinel 2 Items

items_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    flags="i",
    datetime="2024-04-01/2024-09-30",
    format="json",
)
num_items = len(df_items)
df_items = pd.json_normalize(items_json, max_level=0)
display(df_items.head(5))

Table 6

!t.stac.item url={stac_url} collection_id="sentinel-2-l2a" datetime="2024-04-01/2024-09-30" format=plain -i

Let’s look closer at the properties of item id S2A_17SPV_20240919_0_L2A.

df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'")
df_properties = df_item["properties"].apply(pd.Series)
df_item = pd.concat([df_item, df_properties], axis=1)
df_item = df_item.drop("properties", axis=1)
display(df_item.T)

Let’s look at the items properties

df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'")
df_properties = df_item["properties"].apply(pd.Series)
display(df_properties.T)

Now we let’s look at the items assets.

df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'")
df_assets = df_item["assets"].apply(pd.Series)
display(df_assets.T.reset_index())

Assets

Each STAC Item may contain multiple asset types (e.g. red, nir, thumbnail). The snippet below enumerates assets across items, shows per-item asset rows, and prints simple summaries.

Table 7: Sentinel 2 Items Assets Metadata

stac_query = {"eo:cloud_cover": {"lt": 10}}

assets_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    flags="a",
    datetime="2024-04-01/2024-09-30",
    asset_keys="red,nir",
    query=json.dumps(stac_query),
    format="json",
)

df_items_assets = pd.json_normalize(assets_json, max_level=0)
display(df_items_assets.head(10))

display(assets_json)

EO Parameters

STAC item properties expose a number of Earth Observation (EO) parameters that are useful for filtering and selecting scenes. Common properties you will see in Sentinel-2 and other EO STAC records include:

datetime — acquisition time (ISO 8601)
eo:cloud_cover — percent cloud cover (0–100)
eo:gsd — ground sample distance / native resolution (meters)
view:sun_azimuth, view:sun_elevation — sun geometry
proj:epsg — EPSG code (when the projection extension is present)
collection / platform / instrument fields — identify sensor and mission

How to use them: - Use the datetime argument on t.stac.item to constrain time ranges (e.g. 2024-04-01/2024-09-30). - Use the query argument to pass a JSON object of property filters. Comparison operators follow the STAC filtering convention (e.g. lt, lte, gt, gte).

Example query (Python): stac_query = {“eo:cloud_cover”: {“lt”: 10}, “eo:gsd”: {“lte”: 10}, “view:sun_elevation”: {“gte”: 30}}

Pass this to t.stac.item (serialize with json.dumps) to fetch only low-cloud, high-sun-elevation, high-resolution items. Combine asset_keys to limit returned assets (e.g. asset_keys="red,nir"), and use flags (-a, -m, -i) to control output detail (assets, metadata, item list).

Tip: Always inspect a few returned items (format="json" or flags=“m” / flags=“a”`) to see the exact property names used by the catalog you query — naming can vary between providers and extension support.

Create Metadata Vector

stac_query = {"eo:cloud_cover": {"lt": 10}}

assets_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    flags="a",
    datetime="2024-04-01/2024-09-30",
    asset_keys="red,nir",
    query=json.dumps(stac_query),
    items_vector="sentinel_2_items",
    format="json",
)


df_items_assets = pd.json_normalize(assets_json, max_level=0)
display(df_items_assets.head(5))

View the items with iPyleaflet

View STAC items with ipyleaflet (install with pip install ipyleaflet if needed) from ipyleaflet use gj.Interactive and add the sentinel_2_items bounding boxes.

import random

m = gj.InteractiveMap(use_region=False)
m.add_vector(
    "sentinel_2_items",
    style={"opacity": 1, "dashArray": "9", "fillOpacity": 0.1, "weight": 1},
    hover_style={"color": "white", "dashArray": "0", "fillOpacity": 0.5},
    style_callback=lambda f: {
        "color": "black",
        "fillColor": random.choice(["red", "yellow", "green", "orange"]),
    }
)
display(m.show())

Download and Import Rasters

Register STRDS

Python
Terminal

Table 8: Sentinel 2 Items Assets Metadata

import grass.script as gs
import grass.jupyter as gj
stac_query = {"eo:cloud_cover": {"lt": 10}}

items_assets_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    datetime="2024-04-01/2024-09-30",
    asset_keys="red,nir",
    format="json",
    query=json.dumps(stac_query)
)

Table 9: Sentinel 2 Items Assets Metadata

!t.stac.item url="https://earth-search.aws.element84.com/v1/" collection_id="sentinel-2-l2a" datetime="2024-04-01/2024-09-30" asset_keys="red,nir" format="json" query='{"eo:cloud_cover": {"lt": 10}}'

print(f"""
Download Estimate
 
Files:  {items_assets_json["count"]}
Total Download Size: {items_assets_json["bytes"] / 1e9:.2f} GB
""")

Table 10: Sentinel 2 Items Assets Metadata

stac_query = {"eo:cloud_cover": {"lt": 10}}

collection_items_assets_json = gs.parse_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    datetime="2024-04-01/2024-09-30",
    asset_keys="red,green,blue,nir",
    query=json.dumps(stac_query),
    strds_output="outputs/s2_red_nir.txt",
    format="json",
)

print(Path("outputs/s2_red.txt").read_text())

startdate = df_items_assets["datetime"].min()
enddate = df_items_assets["datetime"].max()
print(f"Start Date: {startdate}\nEnd Date: {enddate}")

Create STRDS (spatio-temporal raster dataset) to manage the Sentiel data you want to analyze with GRASS.

# Create the space time dataset
gs.run_command(
    "t.create",
    output="sentinel_2",
    type="strds",
    temporaltype="absolute",
    title="Sentinel 2 Red/NIR Bands",
    description="Sentinel 2 Red/NIR Bands from 2024-04-01/2024-09-30",
    overwrite=True,
)

Download the data

print(gs.run_command(
    "t.stac.item",
    url=stac_url,
    collection_id="sentinel-2-l2a",
    datetime="2024-04-01/2024-09-30",
    asset_keys="red,green,blue,nir",
    query=json.dumps(stac_query),
    strds_output="outputs/s2_red_nir.txt",
    format="json",
    method="nearest",
    resolution="value",
    resolution_value="10",
    extent="region",
    nprocs=12,
    memory=1024 * 7,
    flags="d",
))

# Register the output maps into a space time dataset
gs.run_command(
    "t.register",
    input="sentinel_2",
    file="outputs/s2_red_nir.txt",
    type="raster"
)

gs.run_command("t.info", type="strds", input="sentinel_2")

Now let’s view a RGB composite of one of our tiles

m = gj.Map(use_region=True)
m.d_rgb(
    red="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.red",
    blue="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.blue",
    green="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.green",
)
m.d_barscale()
m.show()

Temporal Map Calculations

Run a temporal map calculation on your entire STRDS to generate an NDVI layer.

gs.run_command(
    "t.rast.mapcalc",
    inputs=["sentinel_2.nir", "sentinel_2.red"],
    output="ndvi",
    basename="ndvi",
    method="follows",
    expression=(
        "float(if('sentinel_2.nir' >= 10000, 10000, 'sentinel_2.nir') - if('sentinel_2.red' >= 10000, 10000, 'sentinel_2.red')) / (if('sentinel_2.nir' >= 10000, 10000, 'sentinel_2.nir') + if('sentinel_2.red' >= 10000, 10000, 'sentinel_2.red'))"
    ),
    overwrite=True,
    flags="n",
)

Check to see it worked.

gs.run_command("t.list", type="strds")

and plot the NDVI

strds_df = pd.DataFrame(
    gs.parse_command(
        "t.rast.list",
        input="ndvi",
        order="start_time",
        columns=["name", "semantic_label", "start_time", "end_time", "min", "max"],
        format="csv",
    )
)

# Convert start_time to datetime for better formatting
strds_df["start_time"] = pd.to_datetime(strds_df["start_time"])

sns.lineplot(
    x="start_time", y="max", data=strds_df
)

# Format x-axis labels for better readability
import matplotlib.dates as mdates
import matplotlib.pyplot as plt

plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
plt.gcf().autofmt_xdate()  # Auto-rotate date labels
plt.show()

ndvi_list = gs.parse_command(
    "t.rast.list",
    input="ndvi",
    order="start_time",
    columns=["name", "semantic_label", "start_time", "end_time", "min", "max"],
    format="csv",
)
gs.run_command("t.rast.colors", input="ndvi", color="ndvi")

for ndvi in ndvi_list:
    map_name = ndvi.get("name")
    m = gj.Map(use_region=True, width=500)
    m.d_rast(map=map_name)
    m.d_barscale()
    m.d_legend(raster=map_name, flags="b")
    m.show()

ndvi_rasters = [
    i.get("name")
    for i in gs.parse_command(
        "t.rast.list",
        input="ndvi",
        order="start_time",
        columns=["name"],
        format="csv",
    )
]

# Create Time Series Map
m = gj.SeriesMap(use_region=True)
m.add_rasters(
    rasters=ndvi_rasters
)
m.d_barscale()
m.render()
m.show()

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from PIL import Image


def generate_map(map_name):
    m = gj.Map(use_region=True)
    m.d_rast(map=map_name)
    m.d_legend(raster=map_name, flags="b")
    m.d_barscale()
    return m


for raster in ndvi_rasters:
    ndvi_map = generate_map(raster)
    ndvi_map.show()

References

--- title: "STAC Tutorial" description: "Learn how to use the GRASS addon t.stac to search and import STAC API data into a GRASS project." image: images/stac-to-the-future.png author: Corey T. White date: '2025-03-01' categories: [stac, grass-addon, tutorial] date-modified: today page-layout: full title-block-banner: true draft: false format: html: toc: true toc-depth: 4 tabsets: true code-copy: true code-fold: false code-tools: true mermaid: theme: forest other-links: - text: grass href: https://grass.osgeo.org - text: STAC href: https://stacspec.org/en lightbox: true canonical-url: true notebook-links: true engine: jupyter jupyter: python3 execute: enabled: true keep-ipynb: true eval: false freeze: auto cache: true --- In this tutorial we will learn about the [t.stac](https://grass.osgeo.org/grass-stable/manuals/addons/t.stac.html) suite of tools for [GRASS](https://grass.osgeo.org). * [t.stac.catalog](https://grass.osgeo.org/grass-stable/manuals/addons/t.stac.catalog.html) - is a tool for exploring SpatioTemporal Asset Catalogs metadata from a STAC API. * [t.stac.collection](https://grass.osgeo.org/grass-stable/manuals/addons/t.stac.collection.html) - is a tool for exploring SpatioTemporal Asset Catalog (STAC) collection metadata. * [t.stac.item](https://grass.osgeo.org/grass-stable/manuals/addons/t.stac.item.html) - is a tool for exploring and importing SpatioTemporal Asset Catalog item metadata and assets into GRASS. ## Introduction SpatioTemporal Asset Catalog (STAC)... The tutorial assumes you already have GRASS installed on your machine. If not please [download and install GRASS](https://grass.osgeo.org/download/) before continuing the tutorial. ## Getting Started ### Start GRASS Session ```{python} # | label: config-grass # | eval: false # Import standard python packages import sys import subprocess from IPython.display import display, JSON, HTML import json import pandas as pd from pathlib import Path import seaborn as sns # import polars as pl # Ask GRASS where its Python packages are and add them to the path sys.path.append( subprocess.check_output(["grass", "--config", "python_path"], text=True).strip() ) # Import the GRASS python packages we need import grass.script as gs import grass.jupyter as gj ``` Create a GRASS project ```{python} # | label: create-project gs.create_project("stac", epsg="32119", overwrite=True) ``` and start a GRASS session. ```{python} # | label: start-grass-session # | eval: false session = gj.init("stac") ``` ### Install Addon The [t.stac](https://grass.osgeo.org/grass-stable/manuals/addons/t.stac.html) addons require the following Python packages. * [pystac (>=1.12)](https://pystac.readthedocs.io/en/latest/installation.html) * [pystac_client (>=0.8)](https://pystac-client.readthedocs.io/en/stable/) * [tqdm (>=4.67)](https://tqdm.github.io/) You can install them using pip or the Python package manager of your choice. ```bash pip install pystac pystac_client tqdm ``` Now you can install the *t.stac* extensions. ```{python} # | label: install-t-stac gs.run_command("g.extension", extension="t.stac") ``` ### Define computational region ```{python} # | label: grass-set-region # | echo: true # | eval: false # | include: true # | output: true gs.run_command( "g.region", n=236687, s=210391, e=616042, w=598921, nsres=10, ewres=10, flags="pa" ) ``` ## STAC API ```{python} # | label: set-stac-catalog-url # | eval: false stac_url = "https://earth-search.aws.element84.com/v1/" ``` ### Searching STAC Catalogs ::: {.panel-tabset group="language"} ## Python ```{python} # | label: tbl-stac-catalog # | echo: true # | eval: false # | include: true # | output: true catalogs_json = gs.parse_command( "t.stac.catalog", url=stac_url, format="json" ) df_catalogs = pd.json_normalize(catalogs_json) df_catalogs.head() ``` ## Terminal ```{python} # | label: tbl-stac-catalog-shell # | echo: true # | eval: false # | include: true # | output: true !t.stac.catalog url={stac_url} format=plain -b ``` ::: #### Collection Metadata ::: {.panel-tabset group="language"} ## Python ```{python} # | label: tbl-stac-items # | echo: true # | eval: false # | include: true # | output: true items_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", flags="m", format="json", ) df_items = pd.json_normalize(items_json, max_level=0) display(df_items.T.head(5)) ``` ## Terminal ```{python} # | label: tbl-stac-items-m-shell # | echo: true # | eval: false # | include: true # | output: true !t.stac.item url={stac_url} collection_id="sentinel-2-l2a" format=plain -m ``` ::: #### Query Items by datetime and display table ::: {.panel-tabset group="language"} ## Python ```{python} # | label: tbl-items # | tbl-cap: Sentinel 2 Items # | echo: true # | eval: false # | include: true # | output: true items_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", flags="i", datetime="2024-04-01/2024-09-30", format="json", ) num_items = len(df_items) df_items = pd.json_normalize(items_json, max_level=0) display(df_items.head(5)) ``` ## Terminal ```{python} # | label: tbl-stac-items-i-shell # | echo: true # | eval: false # | include: true # | output: true !t.stac.item url={stac_url} collection_id="sentinel-2-l2a" datetime="2024-04-01/2024-09-30" format=plain -i ``` ::: ```{python} # | echo: false # | eval: false # | include: true # | output: true num_items = len(df_items) # Can't use inline python with using Jupyter cache print(f"There are {num_items} Sentinel 2 items found in our computational region.") ``` Let's look closer at the properties of item id `S2A_17SPV_20240919_0_L2A`. ```{python} # | echo: true # | eval: false # | include: true # | output: true df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'") df_properties = df_item["properties"].apply(pd.Series) df_item = pd.concat([df_item, df_properties], axis=1) df_item = df_item.drop("properties", axis=1) display(df_item.T) ``` Let's look at the items properties ```{python} # | echo: true # | eval: false # | include: true # | output: true df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'") df_properties = df_item["properties"].apply(pd.Series) display(df_properties.T) ``` Now we let's look at the items assets. ```{python} # | echo: true # | eval: false # | include: true # | output: true df_item = df_items.query("id == 'S2A_17SPV_20240919_0_L2A'") df_assets = df_item["assets"].apply(pd.Series) display(df_assets.T.reset_index()) ``` #### Assets Each STAC Item may contain multiple asset types (e.g. red, nir, thumbnail). The snippet below enumerates assets across items, shows per-item asset rows, and prints simple summaries. ```{python} # | label: tbl-items-assets # | tbl-cap: Sentinel 2 Items Assets Metadata # | echo: true # | eval: false # | include: true # | output: true stac_query = {"eo:cloud_cover": {"lt": 10}} assets_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", flags="a", datetime="2024-04-01/2024-09-30", asset_keys="red,nir", query=json.dumps(stac_query), format="json", ) df_items_assets = pd.json_normalize(assets_json, max_level=0) display(df_items_assets.head(10)) ``` ```{python} display(assets_json) ``` ```{python} # | echo: false # | eval: false # | include: true # | output: true num_assets = len(assets_json) # Can't use inline python with using Jupyter cache print(f"There are {num_assets} Sentinel 2 assets found in our computational region.") ``` #### EO Parameters #### EO Parameters STAC item properties expose a number of Earth Observation (EO) parameters that are useful for filtering and selecting scenes. Common properties you will see in Sentinel-2 and other EO STAC records include: - `datetime` — acquisition time (ISO 8601) - `eo:cloud_cover` — percent cloud cover (0–100) - `eo:gsd` — ground sample distance / native resolution (meters) - `view:sun_azimuth`, `view:sun_elevation` — sun geometry - `proj:epsg` — EPSG code (when the projection extension is present) - collection / platform / instrument fields — identify sensor and mission How to use them: - Use the `datetime` argument on t.stac.item to constrain time ranges (e.g. `2024-04-01/2024-09-30`). - Use the `query` argument to pass a JSON object of property filters. Comparison operators follow the STAC filtering convention (e.g. `lt`, `lte`, `gt`, `gte`). Example query (Python): stac_query = {"eo:cloud_cover": {"lt": 10}, "eo:gsd": {"lte": 10}, "view:sun_elevation": {"gte": 30}} Pass this to t.stac.item (serialize with json.dumps) to fetch only low-cloud, high-sun-elevation, high-resolution items. Combine `asset_keys` to limit returned assets (e.g. `asset_keys="red,nir"`), and use flags (`-a`, `-m`, `-i`) to control output detail (assets, metadata, item list). Tip: Always inspect a few returned items (`format="json"` or flags="m" / flags="a"`) to see the exact property names used by the catalog you query — naming can vary between providers and extension support. ### Create Metadata Vector ```{python} # | label: assets-metadata-vector # | tbl-cap: Sentinel 2 Items Assets Metadata # | echo: true # | eval: false # | include: true # | output: true stac_query = {"eo:cloud_cover": {"lt": 10}} assets_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", flags="a", datetime="2024-04-01/2024-09-30", asset_keys="red,nir", query=json.dumps(stac_query), items_vector="sentinel_2_items", format="json", ) df_items_assets = pd.json_normalize(assets_json, max_level=0) display(df_items_assets.head(5)) ``` #### View the items with iPyleaflet View STAC items with ipyleaflet (install with `pip install ipyleaflet` if needed) from ipyleaflet use `gj.Interactive` and add the sentinel_2_items bounding boxes. ```{python} # | echo: true # | eval: false # | include: true # | output: true import random m = gj.InteractiveMap(use_region=False) m.add_vector( "sentinel_2_items", style={"opacity": 1, "dashArray": "9", "fillOpacity": 0.1, "weight": 1}, hover_style={"color": "white", "dashArray": "0", "fillOpacity": 0.5}, style_callback=lambda f: { "color": "black", "fillColor": random.choice(["red", "yellow", "green", "orange"]), } ) display(m.show()) ``` ### Download and Import Rasters ### Register STRDS ::: {.panel-tabset} ## Python ```{python} # | label: tbl-collection-items-assets # | tbl-cap: Sentinel 2 Items Assets Metadata # | echo: true # | eval: false # | include: true # | output: true import grass.script as gs import grass.jupyter as gj stac_query = {"eo:cloud_cover": {"lt": 10}} items_assets_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", datetime="2024-04-01/2024-09-30", asset_keys="red,nir", format="json", query=json.dumps(stac_query) ) ``` ## Terminal ```{python} # | label: tbl-collection-items-assets-terminal # | tbl-cap: Sentinel 2 Items Assets Metadata # | echo: true # | eval: false # | include: true # | output: true !t.stac.item url="https://earth-search.aws.element84.com/v1/" collection_id="sentinel-2-l2a" datetime="2024-04-01/2024-09-30" asset_keys="red,nir" format="json" query='{"eo:cloud_cover": {"lt": 10}}' ``` ::: ```{python} # | label: download-estimate # | echo: true # | eval: false # | include: true # | output: true print(f""" Download Estimate Files: {items_assets_json["count"]} Total Download Size: {items_assets_json["bytes"] / 1e9:.2f} GB """) ``` ```{python} # | label: tbl-collection-items-assets-strds # | tbl-cap: Sentinel 2 Items Assets Metadata # | echo: true # | eval: false # | include: true # | output: true stac_query = {"eo:cloud_cover": {"lt": 10}} collection_items_assets_json = gs.parse_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", datetime="2024-04-01/2024-09-30", asset_keys="red,green,blue,nir", query=json.dumps(stac_query), strds_output="outputs/s2_red_nir.txt", format="json", ) ``` ```{python} print(Path("outputs/s2_red.txt").read_text()) ``` ```{python} startdate = df_items_assets["datetime"].min() enddate = df_items_assets["datetime"].max() print(f"Start Date: {startdate}\nEnd Date: {enddate}") ``` Create STRDS (spatio-temporal raster dataset) to manage the Sentiel data you want to analyze with GRASS. ```{python} # | label: create-strds # | echo: true # | eval: false # | include: true # | output: true # Create the space time dataset gs.run_command( "t.create", output="sentinel_2", type="strds", temporaltype="absolute", title="Sentinel 2 Red/NIR Bands", description="Sentinel 2 Red/NIR Bands from 2024-04-01/2024-09-30", overwrite=True, ) ``` Download the data ```{python} # | label: download-assets # | echo: true # | eval: false # | include: true # | output: true print(gs.run_command( "t.stac.item", url=stac_url, collection_id="sentinel-2-l2a", datetime="2024-04-01/2024-09-30", asset_keys="red,green,blue,nir", query=json.dumps(stac_query), strds_output="outputs/s2_red_nir.txt", format="json", method="nearest", resolution="value", resolution_value="10", extent="region", nprocs=12, memory=1024 * 7, flags="d", )) ``` Register the data ```{python} # | label: register-strds # | echo: true # | eval: false # | include: true # | output: true # Register the output maps into a space time dataset gs.run_command( "t.register", input="sentinel_2", file="outputs/s2_red_nir.txt", type="raster" ) ``` ```{python} gs.run_command("t.info", type="strds", input="sentinel_2") ``` Now let's view a RGB composite of one of our tiles ```{python} # | label: display-composite # | echo: true # | eval: false # | include: true # | output: true m = gj.Map(use_region=True) m.d_rgb( red="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.red", blue="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.blue", green="sentinel-2-l2a.S2A_17SPV_20240909_0_L2A.green", ) m.d_barscale() m.show() ``` #### Temporal Map Calculations Run a temporal map calculation on your entire STRDS to generate an NDVI layer. ```{python} gs.run_command( "t.rast.mapcalc", inputs=["sentinel_2.nir", "sentinel_2.red"], output="ndvi", basename="ndvi", method="follows", expression=( "float(if('sentinel_2.nir' >= 10000, 10000, 'sentinel_2.nir') - if('sentinel_2.red' >= 10000, 10000, 'sentinel_2.red')) / (if('sentinel_2.nir' >= 10000, 10000, 'sentinel_2.nir') + if('sentinel_2.red' >= 10000, 10000, 'sentinel_2.red'))" ), overwrite=True, flags="n", ) ``` Check to see it worked. ```{python} gs.run_command("t.list", type="strds") ``` and plot the NDVI ```{python} strds_df = pd.DataFrame( gs.parse_command( "t.rast.list", input="ndvi", order="start_time", columns=["name", "semantic_label", "start_time", "end_time", "min", "max"], format="csv", ) ) # Convert start_time to datetime for better formatting strds_df["start_time"] = pd.to_datetime(strds_df["start_time"]) sns.lineplot( x="start_time", y="max", data=strds_df ) # Format x-axis labels for better readability import matplotlib.dates as mdates import matplotlib.pyplot as plt plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d')) plt.gcf().autofmt_xdate() # Auto-rotate date labels plt.show() ``` ```{python} # | layout-ncol: 3 ndvi_list = gs.parse_command( "t.rast.list", input="ndvi", order="start_time", columns=["name", "semantic_label", "start_time", "end_time", "min", "max"], format="csv", ) gs.run_command("t.rast.colors", input="ndvi", color="ndvi") for ndvi in ndvi_list: map_name = ndvi.get("name") m = gj.Map(use_region=True, width=500) m.d_rast(map=map_name) m.d_barscale() m.d_legend(raster=map_name, flags="b") m.show() ``` ```{python} # | label: display-time-series # | echo: true # | eval: false # | include: true # | output: true ndvi_rasters = [ i.get("name") for i in gs.parse_command( "t.rast.list", input="ndvi", order="start_time", columns=["name"], format="csv", ) ] # Create Time Series Map m = gj.SeriesMap(use_region=True) m.add_rasters( rasters=ndvi_rasters ) m.d_barscale() m.render() m.show() ``` ```{python} # | label: display-ndvi-maps # | echo: true # | eval: false # | include: true # | output: true # | layout-ncol: 3 # | fig-cap: # | - Day 1 # | - Day 2 # | - Day 3 # | - Day 4 # | - Day 5 # | - Day 6 # | - Day 7 import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from PIL import Image def generate_map(map_name): m = gj.Map(use_region=True) m.d_rast(map=map_name) m.d_legend(raster=map_name, flags="b") m.d_barscale() return m for raster in ndvi_rasters: ndvi_map = generate_map(raster) ndvi_map.show() ``` ## References - [pystac Documentation](https://pystac.readthedocs.io/) - [pystac_client Documentation](https://pystac-client.readthedocs.io/)