Skip to content
Snippets Groups Projects
Commit 5c4b37af authored by Snyder, Amelia Marie's avatar Snyder, Amelia Marie
Browse files

add archival thredds asset

parent 15482ba1
No related branches found
No related tags found
1 merge request!85add archival thredds asset
Hide whitespace changes
Inline Side-by-side
catalog/alaska_et_2020_ccsm4_historical_simulation/collection.json
+ 37
− 2
View file @ 5c4b37af
......@@ -112,7 +112,8 @@
"zarr-s3-osn": {
"href": "s3://mdmf/gdp/alaska_et_2020_ccsm4_historical_simulation.zarr/",
"type": "application/vnd+zarr",
"description": "Open Storage Network Pod S3 API access to collection zarr group",
"title": "Requester pays access to zarr via S3 API",
"description": "Free, public access to zarr data store via the S3 API. This data is stored on an Open Storage Network Pod.",
"xarray:open_kwargs": {
"chunks": {},
"engine": "zarr",
......@@ -127,7 +128,8 @@
"roles": [
"data",
"zarr",
"s3"
"s3",
"osn"
]
},
"zarr-s3": {
......@@ -147,6 +149,39 @@
"zarr",
"s3"
]
},
"zarr-s3-aws": {
"href": "s3://nhgf-development/workspace/DataConversion/alaska_et_2020_ccsm4_historical_simulation.zarr/",
"type": "application/vnd+zarr",
"title": "Free access to zarr via S3 API",
"description": "Requester pays, public access to zarr data store via the S3 API. This data is stored in an AWS S3 bucket.",
"xarray:open_kwargs": {
"chunks": {},
"engine": "zarr",
"consolidated": true
},
"xarray:storage_options": {
"requester_pays": true
},
"roles": [
"data",
"zarr",
"s3",
"aws",
"requester-pays"
]
},
"legacy-s3-aws": {
"href": "s3://nhgf-development/thredds/alaska_et_2020/",
"title": "Requester pays access to archival legacy files via S3 API",
"description": "Requester pays, public access (via the S3 API) to archival legacy files from WMA THREDDS server that were used to create this zarr store. This data is stored in an AWS S3 bucket.",
"roles": [
"data",
"directory",
"s3",
"aws",
"requester-pays"
]
}
}
}
\ No newline at end of file
workflows/archive/alaska_et_2020_ccsm4_historical_simulation_create_collection_from_zarr.ipynb
+ 36
− 5
View file @ 5c4b37af
%% Cell type:markdown id:6c10e07b-1e60-4926-af1d-fa75dc78e5d4 tags:
# alaska_et_2020_ccsm4_historical_simulation Zarr -> Collection Workflow
This is a workflow to build a [STAC collection](https://github.com/radiantearth/stac-spec/blob/master/collection-spec/collection-spec.md) from the zarr asset for the dataset named above. We use the [datacube extension](https://github.com/stac-extensions/datacube) to define the spatial and temporal dimensions of the zarr store, as well as the variables it contains.
To simplify this workflow so that it can scale to many datasets, a few simplifying suggestions and assumptions are made:
1. For USGS data, we can use the CC0-1.0 license. For all other data we can use Unlicense. Ref: https://spdx.org/licenses/
2. I am assuming all coordinates are from the WGS84 datum if not specified.
%% Cell type:code id:201e0945-de55-45ff-b095-c2af009a4e62 tags:
``` python
import pystac
from pystac.extensions.datacube import CollectionDatacubeExtension, AssetDatacubeExtension, AdditionalDimension, DatacubeExtension
import xarray as xr
import cf_xarray
import os
import fsspec
import cf_xarray
import hvplot.xarray
import pandas as pd
import json
import numpy as np
import pyproj
from pyproj import Transformer
import cartopy.crs as ccrs
import cfunits
import json
import sys
sys.path.insert(1, '..')
import stac_helpers
```
%% Cell type:markdown id:e19630b2-ffe5-4b8e-a903-07dcf7e73d8c tags:
## Collection ID
%% Cell type:code id:b9294785-fd54-4f8e-8002-ffacf1a9e6ff tags:
``` python
# name for STAC collection
collection_id = 'alaska_et_2020_ccsm4_historical_simulation'
```
%% Cell type:markdown id:116b5837-8e85-4ae7-964a-803533ded714 tags:
## Asset Metadata Input
%% Cell type:code id:dd6fa323-132a-4794-8c80-576933f547a0 tags:
``` python
# url to zarr store that you want to create a collection for
zarr_url = f's3://mdmf/gdp/{collection_id}.zarr/'
# define keyword arguments needed for opening the dataset with xarray
# ref: https://github.com/stac-extensions/xarray-assets
xarray_opendataset_kwargs = {"xarray:open_kwargs":{"chunks":{},"engine":"zarr","consolidated":True},
"xarray:storage_options": {"anon": True, "client_kwargs": {"endpoint_url":"https://usgs.osn.mghpcc.org/"}}}
# title - The displayed title for clients and users.
asset_title = 'Requester pays access to zarr via S3 API'
# description for zarr url asset attached to collection (zarr_url)
asset_description = "Open Storage Network Pod S3 API access to collection zarr group"
asset_description = "Free, public access to zarr data store via the S3 API. This data is stored on an Open Storage Network Pod."
# roles to tag zarr url asset with
asset_roles = ["data","zarr","s3"]
asset_roles = ["data","zarr","s3","osn"]
```
%% Cell type:code id:e1441cd4-e94c-4902-af46-8f1af470eb6b tags:
``` python
# url to zarr store that you want to create a collection for
zarr_url2 = f's3://nhgf-development/workspace/DataConversion/{collection_id}.zarr/'
# define keyword arguments needed for opening the dataset with xarray
# ref: https://github.com/stac-extensions/xarray-assets
xarray_opendataset_kwargs2 = {"xarray:open_kwargs":{"chunks":{},"engine":"zarr","consolidated":True},
"xarray:storage_options":{"requester_pays":True}}
# title - The displayed title for clients and users.
asset_title2 = 'Free access to zarr via S3 API'
# description for zarr url asset attached to collection (zarr_url)
asset_description2 = "S3 access to collection zarr group"
asset_description2 = "Requester pays, public access to zarr data store via the S3 API. This data is stored in an AWS S3 bucket."
# roles to tag zarr url asset with
asset_roles2 = ["data","zarr","s3"]
asset_roles2 = ["data","zarr","s3","aws","requester-pays"]
```
%% Cell type:markdown id:b213b74f-ad17-4774-93b6-3b62be616b45 tags:
## Data Exploration
%% Cell type:code id:708f2cf5-79ab-49af-8067-de31d0d13ee6 tags:
``` python
# open and view zarr dataset
fs2 = fsspec.filesystem('s3', anon=True, endpoint_url='https://usgs.osn.mghpcc.org/')
ds = xr.open_dataset(fs2.get_mapper(zarr_url), engine='zarr',
backend_kwargs={'consolidated':True}, chunks={})
ds
```
%% Cell type:markdown id:996e60ba-13e4-453a-8534-e62ce747f0fa tags:
## Collection Metadata Input
%% Cell type:code id:482d204d-b5b6-40e5-ac42-55b459be1097 tags:
``` python
# description of STAC collection
collection_description = ds.attrs['title']
print(f'collection description: {collection_description}')
```
%% Cell type:code id:d4656abd-f435-4ff6-b173-eb7bd3e1c56c tags:
``` python
# license for dataset
collection_license = stac_helpers.license_picker(ds.attrs['license'])
```
%% Cell type:markdown id:0bc7e9b3-ad62-4b10-a18e-66b7ed2d35dc tags:
## Identify x, y, t dimensions of dataset
May require user input if dimensions cannot be auto-detected.
%% Cell type:code id:ab91268f-7200-4cb1-979a-c7d75531d2c0 tags:
``` python
dims_auto_extract = ['X', 'Y', 'T']
dim_names_dict = {}
for d in dims_auto_extract:
dim_names_dict[d] = stac_helpers.extract_dim(ds, d)
print(f"Dimension dictionary: {dim_names_dict}")
```
%% Cell type:markdown id:bad1118b-32ce-4191-939a-21fd58167ea6 tags:
## Get crs info
%% Cell type:code id:f168c59f-9326-42a7-8db3-28dc4b50a416 tags:
``` python
crs_var = 'polar_stereographic'
```
%% Cell type:code id:8307fb3a-ca11-4c19-bcdc-c95ac7275b1c tags:
``` python
# use pyproj to automatically extract crs info
crs = pyproj.CRS.from_cf(ds[crs_var].attrs)
# alternatively, create the appropriate cartopy projection
# crs = ccrs.LambertConformal(central_longitude=crs_info.longitude_of_central_meridian,
# central_latitude=crs_info.latitude_of_projection_origin,
# standard_parallels=crs_info.standard_parallel)
```
%% Cell type:markdown id:9776a7aa-afd1-4de0-b180-1c978f101088 tags:
### Compare dataset crs var to generated proj4 string to make sure it looks ok
%% Cell type:code id:03b05ed7-dfd7-497e-abb2-0711691502e3 tags:
``` python
ds[crs_var]
```
%% Cell type:code id:700f39b4-bb45-40d2-ac0b-dfbf813007e3 tags:
``` python
crs.to_proj4()
```
%% Cell type:markdown id:a8c3ed37-8564-400b-a7fb-25bd5e43d21c tags:
## Create Collection Extent
%% Cell type:markdown id:69f0d837-68a5-4fed-9a14-5d75cfbb0da4 tags:
### Spatial Extent
%% Cell type:code id:aff4370e-205c-4a35-9c14-1d74be8d1f11 tags:
``` python
# pull out lat/lon bbox for data
# coordinates must be from WGS 84 datum
# left, bottom, right, top
# Note: try changing around the commented out lines below to get type float ratherthan a numpy float
#spatial_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float), ds[dim_names_dict['Y']].data.min().compute().astype(float), ds[dim_names_dict['X']].data.max().compute().astype(float), ds[dim_names_dict['Y']].data.max().compute().astype(float)]
#spatial_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['X']].data.max().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.max().compute().astype(float).tolist()]
spatial_bounds = [ds[dim_names_dict['X']].data.min().astype(float).item(), ds[dim_names_dict['Y']].data.min().astype(float).item(), ds[dim_names_dict['X']].data.max().astype(float).item(), ds[dim_names_dict['Y']].data.max().astype(float).item()]
print(spatial_bounds)
print(f'\nspatial_bounds data type: {type(spatial_bounds[0])}')
```
%% Cell type:code id:1a958e87-8dc0-4a4d-985f-302fa09e13fb tags:
``` python
XX, YY = np.meshgrid(ds[dim_names_dict['X']].data, ds[dim_names_dict['Y']].data)
```
%% Cell type:code id:d14bd8b0-d11c-4927-a6b1-7a8fcb35f5b8 tags:
``` python
transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
lon, lat = transformer.transform(XX.ravel(), YY.ravel())
```
%% Cell type:code id:5bd4dc6e-45a0-4af3-a29a-336b4f2d2532 tags:
``` python
print(f'lower left coordinates (WGS84): {min(lon)}, {min(lat)}')
print(f'upper right coordinates (WGS84): {max(lon)}, {max(lat)}')
```
%% Cell type:code id:dd752c5a-778a-4821-9478-55905d7ec121 tags:
``` python
# create a spatial extent object
spatial_extent = pystac.SpatialExtent(bboxes=[[min(lon).item(), min(lat).item(), max(lon).item(), max(lat).item()]])
```
%% Cell type:markdown id:a04c8fca-1d33-43ac-9e2b-62d7be2887f7 tags:
### Temporal Extent
%% Cell type:code id:41a84995-867c-4152-8c57-85e3758bbb77 tags:
``` python
# pull out first and last timestamps
temporal_extent_lower = pd.Timestamp(ds[dim_names_dict['T']].data.min())
temporal_extent_upper = pd.Timestamp(ds[dim_names_dict['T']].data.max())
print(f'min: {temporal_extent_lower} \nmax: {temporal_extent_upper}')
# create a temporal extent object
temporal_extent = pystac.TemporalExtent(intervals=[[temporal_extent_lower, temporal_extent_upper]])
```
%% Cell type:code id:1b1e37c4-5348-46ad-abc9-e005b5d6c02b tags:
``` python
collection_extent = pystac.Extent(spatial=spatial_extent, temporal=temporal_extent)
```
%% Cell type:markdown id:20b00e88-5a13-46b3-9787-d9ac2d4e7bd6 tags:
## Open up NHGF STAC Catalog and create a collection
%% Cell type:code id:adf6c59d-58cd-48b1-a5fd-3bb205a3ef56 tags:
``` python
# define folder location where your STAC catalog json file is
catalog_path = os.path.join('..', '..', 'catalog')
# open catalog
catalog = pystac.Catalog.from_file(os.path.join(catalog_path, 'catalog.json'))
```
%% Cell type:code id:7e96811b-95ae-406a-9728-55fc429d4e1f tags:
``` python
if catalog.get_child(collection_id):
collection = catalog.get_child(collection_id)
print("existing collection opened")
collection.extent=collection_extent
collection.description=collection_description
collection.license=collection_license
else:
collection = pystac.Collection(id=collection_id,
description=collection_description,
extent=collection_extent,
license=collection_license)
print("new collection created")
```
%% Cell type:markdown id:a21c76e8-cd57-4eb5-a33f-7c668a3b3205 tags:
## Add zarr url asset to collection
%% Cell type:code id:094832af-d22b-4359-b0f6-cf687acce5cc tags:
``` python
asset_id = "zarr-s3-osn"
asset = pystac.Asset(href=zarr_url,
title=asset_title,
description=asset_description,
media_type="application/vnd+zarr",
roles=asset_roles,
extra_fields = xarray_opendataset_kwargs)
collection.add_asset(asset_id, asset)
```
%% Cell type:code id:0c298d07-f234-4a08-986d-87f4a39e9ae6 tags:
``` python
asset_id2 = "zarr-s3"
asset_id2 = "zarr-s3-aws"
asset2 = pystac.Asset(href=zarr_url2,
title=asset_title2,
description=asset_description2,
media_type="application/vnd+zarr",
roles=asset_roles2,
extra_fields = xarray_opendataset_kwargs2)
collection.add_asset(asset_id2, asset2)
```
%% Cell type:markdown id:0428bc79-153a-4185-9b69-87d9b4e994fc tags:
### Add Collection Link to Archival Assets on S3
%% Cell type:code id:92ae6437-0404-4078-8104-33b38cfad7a4 tags:
``` python
archival_id = "legacy-s3-aws"
archival_asset = pystac.Asset(href='s3://nhgf-development/thredds/alaska_et_2020/',
title='Requester pays access to archival legacy files via S3 API',
description='Requester pays, public access (via the S3 API) to archival legacy files from WMA THREDDS server that were used to create this zarr store. This data is stored in an AWS S3 bucket.',
roles=["data","directory","s3","aws","requester-pays"])
collection.add_asset(archival_id, archival_asset)
```
%% Cell type:markdown id:f67cd5c9-db33-45c2-bc21-480cd67354f4 tags:
## Add datacube extension to collection
%% Cell type:code id:fc00946d-2880-491d-9b3b-3aeeb4414d6c tags:
``` python
# instantiate extention on collection
dc = DatacubeExtension.ext(collection, add_if_missing=True)
```
%% Cell type:markdown id:8bdd77a2-7587-485e-afb7-42af3a822241 tags:
### Add cube dimensions (required field for extension)
%% Cell type:code id:120a4914-3302-44a5-a282-0308ac84f040 tags:
``` python
# list out dataset dimensions
# When writing data to Zarr, Xarray sets this attribute on all variables based on the variable dimensions. When reading a Zarr group, Xarray looks for this attribute on all arrays,
# raising an error if it can’t be found.
dims = list(ds.dims)
print(dims)
```
%% Cell type:code id:00a18a29-fb9a-4b56-8009-493122997b16 tags:
``` python
# get x, y bounds for extent of those dimensions (required)
xy_bounds = [ds[dim_names_dict['X']].data.min().astype(float).item(), ds[dim_names_dict['Y']].data.min().astype(float).item(), ds[dim_names_dict['X']].data.max().astype(float).item(), ds[dim_names_dict['Y']].data.max().astype(float).item()]
print(xy_bounds)
```
%% Cell type:markdown id:e7dc357c-91ec-49ae-83e5-400f791f9792 tags:
#### user input needed - you will need to look at the crs information and create a cartopy crs object after identifying the projection type:
reference list of cartopy projections: https://scitools.org.uk/cartopy/docs/latest/reference/projections.html
%% Cell type:code id:ea452f62-5644-49b6-8a4e-7dc4f649fd1a tags:
``` python
# print ot crs information in dataset
print(crs)
```
%% Cell type:code id:686ce5f8-f10f-40a5-9780-0bbf85507f26 tags:
``` python
# the datacube extension can accept reference_system information as a numerical EPSG code,
# WKT2 (ISO 19162) string or PROJJSON object.
# we will use a projjson, as was done by Microsoft Planetary Computer here:
# https://planetarycomputer.microsoft.com/dataset/daymet-annual-na
# https://planetarycomputer.microsoft.com/api/stac/v1/collections/daymet-annual-na
projjson = crs.to_json()
print(projjson)
```
%% Cell type:markdown id:205a88b5-5d1b-4a0f-9737-ab308d9f4553 tags:
#### user review needed - looks at the steps pulled out and make sure they make sense
%% Cell type:code id:1a7f6872-8c5c-4bf1-819c-b0b6f5103741 tags:
``` python
time_step = stac_helpers.get_step(ds, dim_names_dict['T'])
print(f'time step: {time_step}')
```
%% Cell type:code id:50d891d3-b269-4ac0-9ae0-61f1c9bed3de tags:
``` python
# debugging for time steps: get all step values and locations
time_step = stac_helpers.get_step(ds, dim_names_dict['T'], debug=True, step_ix=2)
```
%% Cell type:code id:4c56aca3-79a8-4b0b-873f-c6539805f241 tags:
``` python
# debugging for time steps, cont:
# please choose one of the index locations printed above
# this will print the time steps adjacent to it
ix = 182
ds.isel(time=slice(ix-1,ix+3)).time
```
%% Cell type:markdown id:f6ab5e50-a6b8-42a4-a34c-733c5ee6d0c1 tags:
**Irregular time step**: This dataset covers April 1 to Septemper 30 of each year (daily time step), with no data the rest of the time. Because this is an irregular time step, we must specify null for the time step.
%% Cell type:code id:88e06715-e6ba-45f4-a94f-e3414ecb53aa tags:
``` python
x_step = stac_helpers.get_step(ds, dim_names_dict['X'])
print(f'x step: {x_step}')
```
%% Cell type:code id:d5bfded3-5da7-4bbb-9521-4ac0dee51fca tags:
``` python
y_step = stac_helpers.get_step(ds, dim_names_dict['Y'])
print(f'y step: {y_step}')
```
%% Cell type:markdown id:00a5e041-081d-428d-ac2e-75d16de205e6 tags:
#### user input needed - you will need to copy all of the dimensions from above into the dict and fill in the appropriate attributes(type, axis, extent):
%% Cell type:code id:3e86d82d-9e58-4bc8-a826-2191ed2b7477 tags:
``` python
print(dims)
```
%% Cell type:code id:5a443497-67a9-4dce-a8e9-b08d31a88223 tags:
``` python
# create a dictionary of datacube dimensions you would like to assign to this dataset
# dimension name should come from the coordinates printed above
# we do not recommend including redundant dimensions (do not include x,y if you have lon,lat)
# note that the extent of each dimension should be pulled from the dataset
dims_dict = {dim_names_dict['T']: pystac.extensions.datacube.Dimension({'type': 'temporal', 'description': stac_helpers.get_long_name(ds, dim_names_dict['T']), 'extent': [temporal_extent_lower.strftime('%Y-%m-%dT%XZ'), temporal_extent_upper.strftime('%Y-%m-%dT%XZ')], 'step':time_step}),
dim_names_dict['X']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'x', 'description': stac_helpers.get_long_name(ds, dim_names_dict['X']), 'extent': [xy_bounds[0], xy_bounds[2]], 'step': x_step, 'reference_system': projjson}),
dim_names_dict['Y']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'y', 'description': stac_helpers.get_long_name(ds, dim_names_dict['Y']), 'extent': [xy_bounds[1], xy_bounds[3]], 'step': y_step, 'reference_system': projjson}),
}
```
%% Cell type:markdown id:0f277883-a3fd-425f-966a-ca2140d0ef2f tags:
### Add cube variables (optional field for extension)
%% Cell type:code id:e9272931-fc0b-4f2a-9546-283033e9cde8 tags:
``` python
# drop metpy_crs coordinate we have added
if 'metpy_crs' in ds.coords:
ds = ds.drop_vars('metpy_crs')
# pull list of vars from dataset
vars = list(ds.variables)
# spec says that the keys of cube:dimensions and cube:variables should be unique together; a key like lat should not be both a dimension and a variable.
# we will drop all values in dims from vars
vars = [v for v in vars if v not in dims]
# Microsoft Planetary Computer includes coordinates and crs as variables here:
# https://planetarycomputer.microsoft.com/dataset/daymet-annual-na
# https://planetarycomputer.microsoft.com/api/stac/v1/collections/daymet-annual-na
# we will keep those in the var list
# create dictionary of dataset variables and associated dimensions
vars_dict={}
for v in vars:
unit = stac_helpers.get_unit(ds, v)
var_type = stac_helpers.get_var_type(ds, v, crs_var)
long_name = stac_helpers.get_long_name(ds, v)
vars_dict[v] = pystac.extensions.datacube.Variable({'dimensions':list(ds[v].dims), 'type': var_type, 'description': long_name, 'unit': unit})
```
%% Cell type:markdown id:11ad5352-884c-4472-8864-4570a96f66e5 tags:
### Finalize extension
%% Cell type:code id:10141fd4-91d6-491d-878b-02653720891d tags:
``` python
# add dimesions and variables to collection extension
dc.apply(dimensions=dims_dict, variables=vars_dict)
```
%% Cell type:markdown id:615ca168-75fb-4135-9941-0ef5fe4fd1cb tags:
## Add STAC Collection to Catalog and Save
%% Cell type:code id:e2120a55-3d04-4122-a93f-29afcdb8cb1b tags:
``` python
# # helper to find items of wrong type
# d = collection.to_dict()
# print(*stac_helpers.find_paths(d))
```
%% Cell type:code id:4b75791b-6b2d-40be-b7c6-330a60888fb5 tags:
``` python
if catalog.get_child(collection_id):
collection.normalize_and_save(root_href=os.path.join(catalog_path, collection_id), catalog_type=pystac.CatalogType.SELF_CONTAINED)
else:
catalog.add_child(collection)
catalog.normalize_and_save(root_href=catalog_path, catalog_type=pystac.CatalogType.SELF_CONTAINED)
```
%% Cell type:code id:34537e6e-0ed9-4be5-b55f-8139eed55392 tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment