From 2aeab5733382871beafe88a7d40e5d8f5863d14f Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 7 Apr 2026 13:53:17 -0400 Subject: [PATCH 1/4] scratch work to assign spatial:transform when creating geozarr metadata --- docs/examples/sentinel2-tci-to-geozarr.py | 84 +++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 docs/examples/sentinel2-tci-to-geozarr.py diff --git a/docs/examples/sentinel2-tci-to-geozarr.py b/docs/examples/sentinel2-tci-to-geozarr.py new file mode 100644 index 0000000..301e571 --- /dev/null +++ b/docs/examples/sentinel2-tci-to-geozarr.py @@ -0,0 +1,84 @@ +import zarr +from async_geotiff import GeoTIFF +from geozarr_toolkit import ( + MultiscalesConventionMetadata, + ProjConventionMetadata, + SpatialConventionMetadata, + create_multiscales_layout, + create_proj_attrs, + create_spatial_attrs, + create_zarr_conventions, +) +from obstore.store import S3Store +from zarr.storage import LocalStore + +store = S3Store("sentinel-cogs", region="us-west-2", skip_signature=True) +path = "sentinel-s2-l2a-cogs/12/S/UF/2022/6/S2B_12SUF_20220609_0_L2A/TCI.tif" + +geotiff = await GeoTIFF.open(path, store=store) + +proj_attrs = create_proj_attrs(code=f"EPSG:{geotiff.crs.to_epsg()}") +spatial_attrs = create_spatial_attrs( + dimensions=["Y", "X"], + bbox=geotiff.bounds, +) + +# Build multiscales layout from the COG's overviews +# The base (full-resolution) image is level 0; each overview is a coarser level. +levels = [ + {"asset": "0", "transform": {"scale": [1.0, 1.0], "translation": [0.0, 0.0]}}, +] +for i, overview in enumerate(geotiff.overviews): + ov_res = overview.transform.a + scale_factor_x = overview.transform.a / geotiff.transform.a + scale_factor_y = overview.transform.e / geotiff.transform.e + levels.append( + { + "asset": str(i + 1), + "derived_from": "0", + "transform": { + "scale": [scale_factor_x, scale_factor_y], + "translation": [0.0, 0.0], + }, + } + ) + +multiscales_attrs = create_multiscales_layout(levels) +multiscales_attrs["multiscales"]["layout"][0]["spatial:transform"] = list(geotiff.transform) +for item, overview in zip( + multiscales_attrs["multiscales"]["layout"][1:], geotiff.overviews, strict=True +): + item["spatial:transform"] = list(overview.transform) + + +zarr_conventions = create_zarr_conventions( + MultiscalesConventionMetadata(), + ProjConventionMetadata(), + SpatialConventionMetadata(), +) + +geozarr_attrs = { + **proj_attrs, + **spatial_attrs, + **multiscales_attrs, + "zarr_conventions": zarr_conventions, +} + +local_path = "data/TCI.zarr" +zarr_store = LocalStore(local_path) + +root: zarr.Group = zarr.open_group(zarr_store, mode="w", zarr_format=3) + +# Set convention attributes on the group +root.attrs.update(geozarr_attrs) + +# Write the full-resolution image as level "0" +base_array = await geotiff.read() +root.create_array("0", data=base_array.data, chunks=(3, 512, 512)) +print(f"Level 0 (base): shape={base_array.data.shape}, dtype={base_array.data.dtype}") + +# Write each overview as a separate level +for i, overview in enumerate(geotiff.overviews): + ov_array = await overview.read() + root.create_array(str(i + 1), data=ov_array.data, chunks=(3, 512, 512)) + print(f"Level {i + 1} (overview): shape={ov_array.data.shape}") From 5c2ee782dccc2851c312413f964a53b726a5f052 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 7 Apr 2026 14:07:04 -0400 Subject: [PATCH 2/4] truncate transform to 6 elements --- docs/examples/sentinel2-tci-to-geozarr.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/examples/sentinel2-tci-to-geozarr.py b/docs/examples/sentinel2-tci-to-geozarr.py index 301e571..ddd3daf 100644 --- a/docs/examples/sentinel2-tci-to-geozarr.py +++ b/docs/examples/sentinel2-tci-to-geozarr.py @@ -44,11 +44,11 @@ ) multiscales_attrs = create_multiscales_layout(levels) -multiscales_attrs["multiscales"]["layout"][0]["spatial:transform"] = list(geotiff.transform) +multiscales_attrs["multiscales"]["layout"][0]["spatial:transform"] = geotiff.transform[:6] for item, overview in zip( multiscales_attrs["multiscales"]["layout"][1:], geotiff.overviews, strict=True ): - item["spatial:transform"] = list(overview.transform) + item["spatial:transform"] = overview.transform[:6] zarr_conventions = create_zarr_conventions( From 38bbff64949badfed4db6ba5a2eb81bd4ebe4681 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 7 Apr 2026 14:13:35 -0400 Subject: [PATCH 3/4] set `spatial:shape` on each level --- docs/examples/sentinel2-tci-to-geozarr.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/examples/sentinel2-tci-to-geozarr.py b/docs/examples/sentinel2-tci-to-geozarr.py index ddd3daf..07eded8 100644 --- a/docs/examples/sentinel2-tci-to-geozarr.py +++ b/docs/examples/sentinel2-tci-to-geozarr.py @@ -45,10 +45,12 @@ multiscales_attrs = create_multiscales_layout(levels) multiscales_attrs["multiscales"]["layout"][0]["spatial:transform"] = geotiff.transform[:6] +multiscales_attrs["multiscales"]["layout"][0]["spatial:shape"] = geotiff.shape for item, overview in zip( multiscales_attrs["multiscales"]["layout"][1:], geotiff.overviews, strict=True ): item["spatial:transform"] = overview.transform[:6] + item["spatial:shape"] = overview.height, overview.width zarr_conventions = create_zarr_conventions( From 9846dedecfce30ae61da63f2ce9cc9d2c61df424 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 20 Apr 2026 13:46:09 -0400 Subject: [PATCH 4/4] Define band names --- docs/examples/sentinel2-tci-to-geozarr.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/docs/examples/sentinel2-tci-to-geozarr.py b/docs/examples/sentinel2-tci-to-geozarr.py index 07eded8..9f2ea44 100644 --- a/docs/examples/sentinel2-tci-to-geozarr.py +++ b/docs/examples/sentinel2-tci-to-geozarr.py @@ -12,8 +12,12 @@ from obstore.store import S3Store from zarr.storage import LocalStore +DIMENSIONS = ("band", "Y", "X") + + store = S3Store("sentinel-cogs", region="us-west-2", skip_signature=True) -path = "sentinel-s2-l2a-cogs/12/S/UF/2022/6/S2B_12SUF_20220609_0_L2A/TCI.tif" +# path = "sentinel-s2-l2a-cogs/12/S/UF/2022/6/S2B_12SUF_20220609_0_L2A/TCI.tif" +path = "sentinel-s2-l2a-cogs/18/T/WL/2026/1/S2B_18TWL_20260101_0_L2A/TCI.tif" geotiff = await GeoTIFF.open(path, store=store) @@ -76,11 +80,21 @@ # Write the full-resolution image as level "0" base_array = await geotiff.read() -root.create_array("0", data=base_array.data, chunks=(3, 512, 512)) +root.create_array( + "0", + data=base_array.data, + chunks=(3, 512, 512), + dimension_names=DIMENSIONS, +) print(f"Level 0 (base): shape={base_array.data.shape}, dtype={base_array.data.dtype}") # Write each overview as a separate level for i, overview in enumerate(geotiff.overviews): ov_array = await overview.read() - root.create_array(str(i + 1), data=ov_array.data, chunks=(3, 512, 512)) + root.create_array( + str(i + 1), + data=ov_array.data, + chunks=(3, 512, 512), + dimension_names=DIMENSIONS, + ) print(f"Level {i + 1} (overview): shape={ov_array.data.shape}")