STAC Item Creation

Once you have created a collection (see stac-collection-creation.ipynb) use this notebook to create items for that collection.

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 a VEDA JupyterHub instance

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

You are welcome to run this anywhere you like (Note: alternatively you can run this on https://daskhub.veda.smce.nasa.gov/, MAAP, locally, …), just make sure that the data is accessible, or get in contact with the VEDA team to enable access.

Install extra packages

!pip install -U rio_stac parse xpystac pystac nbss-upload --quiet

from parse import parse
from datetime import datetime
import rio_stac
import xarray as xr
import hvplot.xarray

Create pystac.Item

In this section we will be creating a pystac.Item object. This is the part of that notebook that you should update.

Declare constants

Start by declaring some string fields.

COLLECTION_ID = "no2-monthly-diff"
ITEM_ID = "OMI_trno2_0.10x0.10_202212_Col3_V4.nc"
SOURCE = "s3://veda-data-store-staging/no2-monthly-diff/OMI_trno2_0.10x0.10_202212_Col3_V4.nc.tif"

Calculate datetime

Create a function that calculates datetime when given an item_id. You can change this to depend on the source instead if that works better.

def datetime_func(item_id: str) -> datetime:
    """Given the item_id, figure out the datetime"""
    
    fields = parse("OMI_trno2_0.10x0.10_{year:4}{month:2}_Col3_V4.nc", item_id)
    year = int(fields["year"])
    month = int(fields["month"])
    day = 1
    return datetime(year, month, day)

Test out the datetime function:

datetime_func(ITEM_ID)

datetime.datetime(2022, 12, 1, 0, 0)

Put it together

Now take your constants and datetime function and create the STAC Item using rio_stac.

item = rio_stac.stac.create_stac_item(
  id=ITEM_ID,
  source=SOURCE,
  collection=COLLECTION_ID, 
  input_datetime=datetime_func(ITEM_ID),
  with_proj=True,
  with_raster=True,
  asset_name="cog_default",
  asset_roles=["data", "layer"],
  asset_media_type="image/tiff; application=geotiff; profile=cloud-optimized",
)

Try it out!

Now that you have an item you can try it out and make sure it looks good and passes validation checks.

item.validate()

['https://schemas.stacspec.org/v1.0.0/item-spec/json-schema/item.json',
 'https://stac-extensions.github.io/projection/v1.1.0/schema.json',
 'https://stac-extensions.github.io/raster/v1.1.0/schema.json']

item.to_dict()

{'type': 'Feature',
 'stac_version': '1.0.0',
 'id': 'OMI_trno2_0.10x0.10_202212_Col3_V4.nc',
 'properties': {'proj:epsg': 4326,
  '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:bbox': [-180.0, -90.0, 180.0, 90.0],
  'proj:shape': [1800, 3600],
  'proj:transform': [0.1, 0.0, -180.0, 0.0, -0.1, 90.0, 0.0, 0.0, 1.0],
  'proj:projjson': {'$schema': 'https://proj.org/schemas/v0.5/projjson.schema.json',
   'type': 'GeographicCRS',
   'name': 'WGS 84',
   'datum': {'type': 'GeodeticReferenceFrame',
    'name': 'World Geodetic System 1984',
    'ellipsoid': {'name': 'WGS 84',
     'semi_major_axis': 6378137,
     'inverse_flattening': 298.257223563}},
   'coordinate_system': {'subtype': 'ellipsoidal',
    'axis': [{'name': 'Geodetic latitude',
      'abbreviation': 'Lat',
      'direction': 'north',
      'unit': 'degree'},
     {'name': 'Geodetic longitude',
      'abbreviation': 'Lon',
      'direction': 'east',
      'unit': 'degree'}]},
   'id': {'authority': 'EPSG', 'code': 4326}},
  'datetime': '2022-12-01T00:00:00Z'},
 'geometry': {'type': 'Polygon',
  'coordinates': [[(-180.0, -90.0),
    (180.0, -90.0),
    (180.0, 90.0),
    (-180.0, 90.0),
    (-180.0, -90.0)]]},
 'links': [{'rel': <RelType.COLLECTION: 'collection'>,
   'href': 'no2-monthly-diff',
   'type': <MediaType.JSON: 'application/json'>}],
 'assets': {'cog_default': {'href': 's3://veda-data-store-staging/no2-monthly-diff/OMI_trno2_0.10x0.10_202212_Col3_V4.nc.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'raster:bands': [{'data_type': 'float32',
     'scale': 1.0,
     'offset': 0.0,
     'sampling': 'area',
     'nodata': -1.2676506002282294e+30,
     'statistics': {'mean': 12233282717799.0,
      'minimum': -1.30282195779584e+16,
      'maximum': 2.082349180465971e+16,
      'stddev': 416857512760678.5,
      'valid_percent': 82.7056884765625},
     'histogram': {'count': 11,
      'min': -1.30282195779584e+16,
      'max': 2.082349180465971e+16,
      'buckets': [20, 138, 881, 421049, 11300, 203, 20, 3, 1, 1]}}],
   'roles': ['data', 'layer']}},
 'bbox': [-180.0, -90.0, 180.0, 90.0],
 'stac_extensions': ['https://stac-extensions.github.io/projection/v1.1.0/schema.json',
  'https://stac-extensions.github.io/raster/v1.1.0/schema.json'],
 'collection': 'no2-monthly-diff'}

Plot it (optional)

Create a quick visual to make sure that data loads and visualizes properly.

data = xr.open_dataset(item).cog_default.isel(time=0)
data

<xarray.DataArray 'cog_default' (y: 1800, x: 3600)>
[6480000 values with dtype=float64]
Coordinates:
    time            datetime64[ns] 2022-12-01
    id              <U37 ...
  * x               (x) float64 -180.0 -179.9 -179.8 ... 179.7 179.8 179.9
  * y               (y) float64 90.0 89.9 89.8 89.7 ... -89.6 -89.7 -89.8 -89.9
    proj:transform  object ...
    proj:projjson   object ...
    proj:epsg       int64 ...
    proj:bbox       object ...
    proj:shape      object ...
    proj:geometry   object ...
    raster:bands    object ...
    epsg            int64 ...

xarray.DataArray

'cog_default'

• y: 1800
• x: 3600

• ...

  [6480000 values with dtype=float64]

• Coordinates: (12)
  • time
    ()
    datetime64[ns]
    2022-12-01

    array('2022-12-01T00:00:00.000000000', dtype='datetime64[ns]')

  • id
    ()
    <U37
    ...

    [1 values with dtype=<U37]

  • x
    (x)
    float64
    -180.0 -179.9 ... 179.8 179.9

    array([-180. , -179.9, -179.8, ...,  179.7,  179.8,  179.9])

  • y
    (y)
    float64
    90.0 89.9 89.8 ... -89.8 -89.9

    array([ 90. ,  89.9,  89.8, ..., -89.7, -89.8, -89.9])

  • proj:transform
    ()
    object
    ...

    [1 values with dtype=object]

  • proj:projjson
    ()
    object
    ...

    [1 values with dtype=object]

  • proj:epsg
    ()
    int64
    ...

    [1 values with dtype=int64]

  • proj:bbox
    ()
    object
    ...

    [1 values with dtype=object]

  • proj:shape
    ()
    object
    ...

    [1 values with dtype=object]

  • proj:geometry
    ()
    object
    ...

    [1 values with dtype=object]

  • raster:bands
    ()
    object
    ...

    [1 values with dtype=object]

  • epsg
    ()
    int64
    ...

    [1 values with dtype=int64]

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Float64Index([            -180.0,             -179.9,             -179.8,
                              -179.7,             -179.6,             -179.5,
                              -179.4,             -179.3,             -179.2,
                              -179.1,
                  ...
                               179.0, 179.10000000000002, 179.20000000000005,
                               179.3, 179.40000000000003,              179.5,
                  179.60000000000002, 179.70000000000005,              179.8,
                  179.90000000000003],
                 dtype='float64', name='x', length=3600))

  • y
    PandasIndex

    PandasIndex(Float64Index([              90.0,               89.9,               89.8,
                                89.7,               89.6,               89.5,
                                89.4,               89.3,               89.2,
                                89.1,
                  ...
                               -89.0, -89.10000000000002, -89.20000000000002,
                  -89.30000000000001,              -89.4,              -89.5,
                  -89.60000000000002, -89.70000000000002, -89.80000000000001,
                               -89.9],
                 dtype='float64', name='y', length=1800))

• Attributes: (0)

color_range = tuple(data.quantile([.02, .98]).values)

data.hvplot("x", "y", clim=color_range, cmap="jet", rasterize=True)

NOTE: Jet is a bad colormap because it overemphasizes certain values, but it has a very large number of colors so it is good for spotting odd patterns in the data.

Upload this notebook

You can upload the notebook to anyplace you like, but one of the easiest ones is notebook sharing space. Just change the following cell from “Raw” to “Code”, run it and copy the output link.

Before uploading make sure: 1) you have not hard-coded any secrets or access keys. 2) you have saved this notebook. Hint (ctrl+s) will do it

!nbss-upload new-item.ipynb