Skip to content

Commit 936da03

Browse files
rouaultrouault
committed
Zarr: add write support for 'spatial' and 'proj' conventions if GEOREFERENCING_CONVENTION creation option set to SPATIAL_PROJ
also fix read-support for how spatial:registration=node influences GDAL geotransform
1 parent b79430a commit 936da03

File tree

8 files changed

+698
-61
lines changed

8 files changed

+698
-61
lines changed

‎autotest/gdrivers/zarr_driver.py‎

Lines changed: 244 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6971,3 +6971,247 @@ def test_zarr_read_spatial_geotransform():
69716971
assert ds.GetGeoTransform() == (450000, 10, 0.1, 5000000, 0.15, -10)
69726972
assert ds.GetSpatialRef() is None
69736973
assert ds.GetMetadata() == {"AREA_OR_POINT": "Area"}
6974+
6975+
6976+
###############################################################################
6977+
#
6978+
6979+
6980+
@gdaltest.enable_exceptions()
6981+
def test_zarr_write_spatial_geotransform(tmp_vsimem):
6982+
6983+
src_ds = gdal.Open("data/byte.tif")
6984+
gdal.GetDriverByName("Zarr").CreateCopy(
6985+
tmp_vsimem / "out.zarr",
6986+
src_ds,
6987+
options={"FORMAT": "ZARR_V3", "GEOREFERENCING_CONVENTION": "SPATIAL_PROJ"},
6988+
)
6989+
6990+
with gdal.VSIFile(tmp_vsimem / "out.zarr" / "zarr.json", "rb") as f:
6991+
j = json.loads(f.read())
6992+
6993+
assert j == {
6994+
"zarr_format": 3,
6995+
"node_type": "group",
6996+
"attributes": {},
6997+
"consolidated_metadata": {
6998+
"kind": "inline",
6999+
"must_understand": False,
7000+
"metadata": {
7001+
"X": {
7002+
"zarr_format": 3,
7003+
"node_type": "array",
7004+
"shape": [20],
7005+
"data_type": "float64",
7006+
"chunk_grid": {
7007+
"name": "regular",
7008+
"configuration": {"chunk_shape": [20]},
7009+
},
7010+
"chunk_key_encoding": {
7011+
"name": "default",
7012+
"configuration": {"separator": "/"},
7013+
},
7014+
"fill_value": "NaN",
7015+
"codecs": [
7016+
{"name": "bytes", "configuration": {"endian": "little"}}
7017+
],
7018+
"attributes": {},
7019+
"dimension_names": ["X"],
7020+
},
7021+
"Y": {
7022+
"zarr_format": 3,
7023+
"node_type": "array",
7024+
"shape": [20],
7025+
"data_type": "float64",
7026+
"chunk_grid": {
7027+
"name": "regular",
7028+
"configuration": {"chunk_shape": [20]},
7029+
},
7030+
"chunk_key_encoding": {
7031+
"name": "default",
7032+
"configuration": {"separator": "/"},
7033+
},
7034+
"fill_value": "NaN",
7035+
"codecs": [
7036+
{"name": "bytes", "configuration": {"endian": "little"}}
7037+
],
7038+
"attributes": {},
7039+
"dimension_names": ["Y"],
7040+
},
7041+
"out": {
7042+
"zarr_format": 3,
7043+
"node_type": "array",
7044+
"shape": [20, 20],
7045+
"data_type": "uint8",
7046+
"chunk_grid": {
7047+
"name": "regular",
7048+
"configuration": {"chunk_shape": [20, 20]},
7049+
},
7050+
"chunk_key_encoding": {
7051+
"name": "default",
7052+
"configuration": {"separator": "/"},
7053+
},
7054+
"fill_value": None,
7055+
"codecs": [
7056+
{"name": "bytes", "configuration": {"endian": "little"}}
7057+
],
7058+
"attributes": {
7059+
"COLOR_INTERPRETATION": "Gray",
7060+
"proj:code": "EPSG:26711",
7061+
"spatial:bbox": [
7062+
440720.0,
7063+
3750120.0,
7064+
441920.0,
7065+
3751320.0,
7066+
],
7067+
"spatial:transform_type": "affine",
7068+
"spatial:transform": [
7069+
60.0,
7070+
0.0,
7071+
440720.0,
7072+
0.0,
7073+
-60.0,
7074+
3751320.0,
7075+
],
7076+
"spatial:registration": "pixel",
7077+
"spatial:dimensions": ["Y", "X"],
7078+
"zarr_conventions": [
7079+
{
7080+
"schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json",
7081+
"spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md",
7082+
"uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f",
7083+
"name": "proj:",
7084+
"description": "Coordinate reference system information for geospatial data",
7085+
},
7086+
{
7087+
"schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json",
7088+
"spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md",
7089+
"uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4",
7090+
"name": "spatial:",
7091+
"description": "Spatial coordinate information",
7092+
},
7093+
],
7094+
},
7095+
"dimension_names": ["Y", "X"],
7096+
},
7097+
},
7098+
},
7099+
}
7100+
7101+
# Remove CF-like X and Y coordinate variables, to be sure we read from
7102+
# spatial:transform
7103+
gdal.RmdirRecursive(tmp_vsimem / "out.zarr" / "X")
7104+
gdal.RmdirRecursive(tmp_vsimem / "out.zarr" / "Y")
7105+
7106+
ds = gdal.Open(tmp_vsimem / "out.zarr")
7107+
assert ds.GetSpatialRef().IsSame(src_ds.GetSpatialRef())
7108+
assert ds.GetGeoTransform() == src_ds.GetGeoTransform()
7109+
assert ds.GetRasterBand(1).Checksum() == src_ds.GetRasterBand(1).Checksum()
7110+
7111+
7112+
###############################################################################
7113+
#
7114+
7115+
7116+
@gdaltest.enable_exceptions()
7117+
def test_zarr_write_spatial_geotransform_no_epsg_code_rotated_gt_and_pixel_center(
7118+
tmp_vsimem,
7119+
):
7120+
7121+
src_ds = gdal.GetDriverByName("MEM").Create("", 2, 3)
7122+
src_ds.SetGeoTransform([1, 2, 3, 4, 5, 6])
7123+
src_ds.SetSpatialRef(
7124+
osr.SpatialReference("+proj=longlat +ellps=GRS80 +towgs84=0,0,0")
7125+
)
7126+
src_ds.SetMetadataItem("AREA_OR_POINT", "Point")
7127+
7128+
gdal.GetDriverByName("Zarr").CreateCopy(
7129+
tmp_vsimem / "out.zarr",
7130+
src_ds,
7131+
options={"FORMAT": "ZARR_V3", "GEOREFERENCING_CONVENTION": "SPATIAL_PROJ"},
7132+
)
7133+
7134+
with gdal.VSIFile(tmp_vsimem / "out.zarr" / "zarr.json", "rb") as f:
7135+
j = json.loads(f.read())
7136+
7137+
assert j["consolidated_metadata"]["metadata"]["out"]["attributes"][
7138+
"proj:wkt2"
7139+
].startswith("BOUNDCRS")
7140+
assert (
7141+
j["consolidated_metadata"]["metadata"]["out"]["attributes"]["proj:projjson"][
7142+
"type"
7143+
]
7144+
== "BoundCRS"
7145+
)
7146+
7147+
# Set those 2 elements to dummy value for comparison, as their content might be PROJ dependent
7148+
j["consolidated_metadata"]["metadata"]["out"]["attributes"][
7149+
"proj:wkt2"
7150+
] = "redacted"
7151+
j["consolidated_metadata"]["metadata"]["out"]["attributes"][
7152+
"proj:projjson"
7153+
] = "redacted"
7154+
7155+
assert j == {
7156+
"zarr_format": 3,
7157+
"node_type": "group",
7158+
"attributes": {},
7159+
"consolidated_metadata": {
7160+
"kind": "inline",
7161+
"must_understand": False,
7162+
"metadata": {
7163+
"out": {
7164+
"zarr_format": 3,
7165+
"node_type": "array",
7166+
"shape": [3, 2],
7167+
"data_type": "uint8",
7168+
"chunk_grid": {
7169+
"name": "regular",
7170+
"configuration": {"chunk_shape": [3, 2]},
7171+
},
7172+
"chunk_key_encoding": {
7173+
"name": "default",
7174+
"configuration": {"separator": "/"},
7175+
},
7176+
"fill_value": None,
7177+
"codecs": [
7178+
{"name": "bytes", "configuration": {"endian": "little"}}
7179+
],
7180+
"attributes": {
7181+
"proj:wkt2": "redacted",
7182+
"proj:projjson": "redacted",
7183+
"spatial:bbox": [3.5, 9.5, 11.5, 26.5],
7184+
"spatial:transform_type": "affine",
7185+
"spatial:transform": [2.0, 3.0, 3.5, 5.0, 6.0, 9.5],
7186+
"spatial:dimensions": ["Y", "X"],
7187+
"spatial:registration": "node",
7188+
"zarr_conventions": [
7189+
{
7190+
"schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json",
7191+
"spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md",
7192+
"uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f",
7193+
"name": "proj:",
7194+
"description": "Coordinate reference system information for geospatial data",
7195+
},
7196+
{
7197+
"schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json",
7198+
"spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md",
7199+
"uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4",
7200+
"name": "spatial:",
7201+
"description": "Spatial coordinate information",
7202+
},
7203+
],
7204+
},
7205+
"dimension_names": ["Y", "X"],
7206+
},
7207+
},
7208+
},
7209+
}
7210+
7211+
ds = gdal.Open(tmp_vsimem / "out.zarr")
7212+
assert ds.GetSpatialRef().IsSame(src_ds.GetSpatialRef())
7213+
assert ds.GetGeoTransform() == src_ds.GetGeoTransform()
7214+
assert ds.GetMetadataItem("AREA_OR_POINT") == src_ds.GetMetadataItem(
7215+
"AREA_OR_POINT"
7216+
)
7217+
assert ds.GetRasterBand(1).Checksum() == src_ds.GetRasterBand(1).Checksum()

‎doc/source/drivers/raster/zarr.rst‎

Lines changed: 93 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -13,20 +13,11 @@ Zarr
1313

1414
Zarr is a format for the storage of chunked, compressed, N-dimensional arrays.
1515
This format is supported for read and write access, and using the traditional
16-
2D raster API or the newer multidimensional API
16+
2D raster API or the multidimensional API
1717

18-
The driver supports the Zarr V2 specification, and has experimental support
19-
for the in-progress Zarr V3 specification. It also supports Kerchunk reference
18+
The driver supports the Zarr V2 an V3 specifications. It also supports Kerchunk reference
2019
files since GDAL 3.11.
2120

22-
.. warning::
23-
24-
The current implementation of Zarr V3 before GDAL 3.8 is incompatible with
25-
the latest evolutions of the Zarr V3 specification.
26-
GDAL 3.8 is compatible with the Zarr V3 specification at date 2023-May-7,
27-
and is not interoperable with Zarr V3 datasets produced by earlier GDAL
28-
versions.
29-
3021
Local and cloud storage (see :ref:`virtual_file_systems`) are supported in read and write.
3122

3223
Driver capabilities
@@ -196,16 +187,28 @@ The driver support the
196187
`NCZarr v2 <https://www.unidata.ucar.edu/software/netcdf/documentation/NUG/nczarr_head.html>`__
197188
extensions of storing the dimension names of an array (read-only)
198189

199-
SRS encoding
200-
------------
190+
Georeferencing encoding (CRS and geotransformation matrix)
191+
----------------------------------------------------------
201192

202193
The Zarr specification has no provision for spatial reference system encoding.
203-
GDAL uses a ``_CRS`` attribute that is a dictionary that may contain one or
194+
Several conventions
195+
196+
GDAL convention
197+
+++++++++++++++
198+
199+
Before GDAL 3.13, the only convention supported both in reading and writing was
200+
the ``GDAL`` one, using a ``_CRS`` attribute. The geotransformation matrix, when
201+
no rotation terms is present, is encoded as ``X`` and ``Y`` one-dimensional
202+
coordinate arrays.
203+
204+
The ``_CRS`̀ attribute is a dictionary that may contain one or
204205
several of the following keys: ``url`` (using a OGC CRS URL), ``wkt`` (WKT:2019
205206
used by default on writing, WKT1 also supported on reading.), ``projjson``.
206207
On reading, it will use ``url`` by default, if not found will fallback to ``wkt``
207208
and then ``projjson``.
208209

210+
Example:
211+
209212
.. code-block:: json
210213
211214
{
@@ -333,6 +336,70 @@ and then ``projjson``.
333336
}
334337
}
335338
339+
340+
SPATIAL_PROJ convention
341+
+++++++++++++++++++++++
342+
343+
.. versionadded:: 3.13
344+
345+
Since GDAL 3.13, the Zarr `spatial <https://github.com/zarr-conventions/spatial>`__
346+
and `geo-proj <https://github.com/zarr-conventions/geo-proj>`__ conventions
347+
are supported in reading, and in writing when the ``GEOREFERENCING_CONVENTION``
348+
creation option is set to ``SPATIAL_PROJ``. X and Y coordinate arrays are
349+
written only if the geotransformation matrix has no rotation terms.
350+
351+
Example:
352+
353+
.. code-block:: json
354+
355+
{
356+
"attributes": {
357+
"proj:code": "EPSG:26711",
358+
"spatial:bbox": [
359+
440720.0,
360+
3750120.0,
361+
441920.0,
362+
3751320.0,
363+
],
364+
"spatial:transform_type": "affine",
365+
"spatial:transform": [
366+
60.0,
367+
0.0,
368+
440720.0,
369+
0.0,
370+
-60.0,
371+
3751320.0,
372+
],
373+
"spatial:registration": "pixel",
374+
"spatial:dimensions": ["Y", "X"],
375+
"zarr_conventions": [
376+
{
377+
"schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json",
378+
"spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md",
379+
"uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f",
380+
"name": "proj:",
381+
"description": "Coordinate reference system information for geospatial data",
382+
},
383+
{
384+
"schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json",
385+
"spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md",
386+
"uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4",
387+
"name": "spatial:",
388+
"description": "Spatial coordinate information",
389+
},
390+
]
391+
}
392+
}
393+
394+
395+
netCDF CF conventions
396+
+++++++++++++++++++++
397+
398+
.. versionadded:: 3.9
399+
400+
The driver supports reading a CRS using the `CF conventions <https://cfconventions.org/>`__.
401+
402+
336403
Particularities of the classic raster API
337404
-----------------------------------------
338405

@@ -449,6 +516,18 @@ The following options are creation options of the classic raster API, or
449516
array-level creation options for the multidimensional API (must be prefixed
450517
with ``ARRAY:`` using :program:`gdalmdimtranslate`):
451518

519+
- .. co:: GEOREFERENCING_CONVENTION
520+
:choices: GDAL, SPATIAL_PROJ
521+
:default: GDAL
522+
523+
Which convention is used to write georeferencing information: geotransformation
524+
and CRS.
525+
526+
The ``GDAL`` convention uses a ``_CRS`` attribute described above. The
527+
``SPATIAL_PROJ`` convention, added both in read and write support in GDAL 3.13,
528+
uses the Zarr `spatial <https://github.com/zarr-conventions/spatial>`__
529+
and `geo-proj <https://github.com/zarr-conventions/geo-proj>`__ conventions.
530+
452531
- .. co:: COMPRESS
453532
:choices: NONE, BLOSC, ZLIB, GZIP, LZMA, ZSTD, LZ4
454533
:default: NONE

0 commit comments

Comments
 (0)