Skip to content

Commit

Permalink
STACs DeltaDTM v1.1 (#10)
Browse files Browse the repository at this point in the history
* correct syntax datetime.utc

* Changes to incorporate geoserver in collection

Also put stac_io to default from pystac as in CoCliCo
Note, geoserver link is still a dummy.. Cannot fix it yet (will be soon)

* rm storage account

* add cog media type

* forgiveness when testing if collection exists

* example for reading deltadtm data

* add storage account name

* deltadtm v1.1

* fix

---------

Co-authored-by: floriscalkoen <[email protected]>
Co-authored-by: EtienneKras <[email protected]>
  • Loading branch information
3 people authored Nov 22, 2024
1 parent 477fc7e commit c001877
Show file tree
Hide file tree
Showing 5 changed files with 215 additions and 15 deletions.
3 changes: 1 addition & 2 deletions .env.example
Original file line number Diff line number Diff line change
@@ -1,2 +1 @@
AZURE_STORAGE_SAS_TOKEN="paste-your-sas-token-here"
CLIENT_AZURE_STORAGE_CONNECTION_STRING="paste-your-connection-string-here"
AZURE_STORAGE_SAS_TOKEN="paste-your-sas-token-here"
2 changes: 1 addition & 1 deletion analytics/hypsometry.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@
"load_dotenv(override=True)\n",
"\n",
"sas_token = os.getenv(\"AZURE_STORAGE_SAS_TOKEN\")\n",
"storage_account_name = os.getenv(\"AZURE_STORAGE_ACCOUNT\")\n",
"storage_account_name = \"coclico\"\n",
"storage_options = {\"account_name\": storage_account_name, \"sas_token\": sas_token}\n",
"\n",
"LOWER_THAN = (\n",
Expand Down
36 changes: 25 additions & 11 deletions scripts/python/stac_deltadtm.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

import fsspec
import pystac
import pystac.media_type
import rasterio
import shapely
import stac_geoparquet
Expand All @@ -20,11 +21,12 @@
# os.system(f"git clone https://github.com/openearth/coclicodata.git {dev_dir / 'coclicodata'}")
# # Install the package in development mode
# os.system(f"pip install -e {dev_dir / 'coclicodata'}")
from coclicodata.coclico_stac.io import CoCliCoStacIO
from coclicodata.coclico_stac.layouts import CoCliCoCOGLayout
from dotenv import load_dotenv
from pystac.extensions import raster
from pystac.stac_io import DefaultStacIO
from stactools.core.utils import antimeridian
from tqdm import tqdm

# Load the environment variables from the .env file
load_dotenv(override=True)
Expand All @@ -33,8 +35,8 @@

# Get the SAS token and storage account name from environment variables
sas_token = os.getenv("AZURE_STORAGE_SAS_TOKEN")
storage_account_name = os.getenv("AZURE_STORAGE_ACCOUNT_NAME")
storage_options = {"account_name": storage_account_name, "credential": sas_token}
STORAGE_ACCOUNT_NAME = "coclico"
storage_options = {"account_name": STORAGE_ACCOUNT_NAME, "credential": sas_token}

# CoCliCo STAC
STAC_DIR = pathlib.Path.home() / "dev" / "coclicodata" / "current"
Expand Down Expand Up @@ -195,7 +197,7 @@ def create_item(block, item_id, antimeridian_strategy=antimeridian.Strategy.SPLI

antimeridian.fix_item(item, antimeridian_strategy)

item.common_metadata.created = datetime.datetime.utcnow()
item.common_metadata.created = datetime.datetime.now(datetime.UTC)

ext = pystac.extensions.projection.ProjectionExtension.ext(
item, add_if_missing=True
Expand Down Expand Up @@ -245,7 +247,7 @@ def create_asset(

fps = fs.glob(f"{root}/**/*.tif")

stac_io = CoCliCoStacIO()
stac_io = DefaultStacIO() # CoCliCoStacIO()
layout = CoCliCoCOGLayout()

collection = create_collection()
Expand All @@ -255,7 +257,7 @@ def create_asset(
# filename=b'ed4e8723-89d8-4075-b946-1ae1a21cca03/ed4e8723-89d8-4075-b946-1ae1a21cca03.aux'
logging.getLogger("rasterio").setLevel(logging.WARNING)

for fp in fps:
for fp in tqdm(fps, desc="Processing file paths"):
href = "az://" + fp
pp = PathParser(href)

Expand Down Expand Up @@ -293,11 +295,25 @@ def create_asset(
),
)

collection.add_asset(
"geoserver_link",
pystac.Asset(
# https://coclico.avi.deltares.nl/geoserver/%s/wms?bbox={bbox-epsg-3857}&format=image/png&service=WMS&version=1.1.1&request=GetMap&srs=EPSG:3857&transparent=true&width=256&height=256&layers=%s"%(COLLECTION_ID, COLLECTION_ID + ":" + ASSET_TITLE),
"https://coclico.avi.deltares.nl/geoserver/cfhp/wms?bbox={bbox-epsg-3857}&format=image/png&service=WMS&version=1.1.1&request=GetMap&srs=EPSG:3857&transparent=true&width=256&height=256&layers=cfhp:HIGH_DEFENDED_MAPS_Mean_spring_tide", # test
title="Geoserver Mosaic link",
media_type=pystac.media_type.MediaType.COG,
),
)

catalog = pystac.Catalog.from_file(str(STAC_DIR / "catalog.json"))

if catalog.get_child(collection.id):
catalog.remove_child(collection.id)
print(f"Removed child: {collection.id}.")
# TODO: there should be a cleaner method to remove the previous stac catalog and its items
try:
if catalog.get_child(collection.id):
catalog.remove_child(collection.id)
print(f"Removed child: {collection.id}.")
except Exception:
pass

catalog.add_child(collection)

Expand All @@ -310,5 +326,3 @@ def create_asset(
dest_href=str(STAC_DIR),
stac_io=stac_io,
)

catalog.validate_all()
187 changes: 187 additions & 0 deletions tutorials/deltadtm.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,187 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# DeltaDTM\n",
"\n",
"A global coastal digital terrain model, based on CopernicusDEM, ESA WorldCover, ICESat-2 and GEDI data. For more information, see [Pronk et al. (2024)](https://www.nature.com/articles/s41597-024-03091-9) DeltaDTM: A global coastal digital terrain model. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1",
"metadata": {},
"outputs": [],
"source": [
"import geopandas as gpd\n",
"import hvplot.xarray\n",
"import pystac\n",
"import rioxarray\n",
"import shapely\n",
"import xarray as xr\n",
"from ipyleaflet import Map, basemaps\n",
"\n",
"storage_options = {\"account_name\": \"coclico\"}"
]
},
{
"cell_type": "markdown",
"id": "2",
"metadata": {},
"source": [
"## Read a snapshot with the spatial extents of all tiles\n",
"\n",
"Connect to the CoCliCo STAC and read the spatial extents using stac-geoparquet."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3",
"metadata": {},
"outputs": [],
"source": [
"from coastpy.io.utils import read_items_extent\n",
"\n",
"coclico_catalog = pystac.Catalog.from_file(\n",
" \"https://coclico.blob.core.windows.net/stac/v1/catalog.json\"\n",
")\n",
"\n",
"ddtm_collection = coclico_catalog.get_child(\"deltares-delta-dtm\")\n",
"ddtm_extents = read_items_extent(\n",
" ddtm_collection, columns=[\"geometry\", \"assets\", \"href\"]\n",
")\n",
"\n",
"ddtm_extents.head()"
]
},
{
"cell_type": "markdown",
"id": "4",
"metadata": {},
"source": [
"## Zoom to your area of interest"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5",
"metadata": {},
"outputs": [],
"source": [
"m = Map(basemap=basemaps.Esri.WorldImagery, scroll_wheel_zoom=True)\n",
"m.center = 15.827, -95.96\n",
"m.zoom = 15\n",
"m.layout.height = \"800px\"\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6",
"metadata": {},
"outputs": [],
"source": [
"west, south, east, north = m.west, m.south, m.east, m.north\n",
"# Note: small little hack to ensure the notebook also works when running all cells at once\n",
"if not west:\n",
" west, south, east, north = (\n",
" 30.28415679931641,\n",
" 31.276790311057272,\n",
" 30.630912780761722,\n",
" 31.51123970051334,\n",
" )\n",
"roi = gpd.GeoDataFrame(\n",
" geometry=[shapely.geometry.box(west, south, east, north)], crs=4326\n",
")"
]
},
{
"cell_type": "markdown",
"id": "7",
"metadata": {},
"source": [
"## Find the tiles for your region of interest"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"metadata": {},
"outputs": [],
"source": [
"from coastpy.io.cloud import to_https_url\n",
"\n",
"hrefs = gpd.sjoin(ddtm_extents, roi).href.to_list()\n",
"href = hrefs[0]\n",
"href = to_https_url(href, storage_options={\"account_name\": \"coclico\"})"
]
},
{
"cell_type": "markdown",
"id": "9",
"metadata": {},
"source": [
"## Read data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"metadata": {},
"outputs": [],
"source": [
"da = rioxarray.open_rasterio(href, chunks={}, lock=False).squeeze().drop_vars(\"band\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11",
"metadata": {},
"outputs": [],
"source": [
"da.where(lambda xx: xx != xx.attrs[\"_FillValue\"]).hvplot(\n",
" x=\"x\", y=\"y\", geo=True, tiles=\"ESRI\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:coastal-full] *",
"language": "python",
"name": "conda-env-coastal-full-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
2 changes: 1 addition & 1 deletion tutorials/global_coastal_transect_system.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:coastal] *",
"display_name": "Python [conda env:coastal]",
"language": "python",
"name": "conda-env-coastal-py"
},
Expand Down

0 comments on commit c001877

Please sign in to comment.