Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Backport release/3.10] Zarr: fix incorrect DataAxisToSRSAxisMapping for EPSG geographic CRS #11415

Merged
merged 1 commit into from
Dec 2, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 57 additions & 15 deletions autotest/gdrivers/zarr_driver.py
Original file line number Diff line number Diff line change
Expand Up @@ -819,15 +819,21 @@ def test_zarr_read_crs(crs_member):
assert rg
ar = rg.OpenMDArray(rg.GetMDArrayNames()[0])
srs = ar.GetSpatialRef()
if (
not (osr.GetPROJVersionMajor() > 6 or osr.GetPROJVersionMinor() >= 2)
and crs_member == "projjson"
):
assert srs is None
else:
assert srs is not None
assert srs.GetAuthorityCode(None) == "4326"
assert len(ar.GetAttributes()) == 0
assert srs is not None
assert srs.GetAuthorityCode(None) == "4326"
# Mapping is 1, 2 since the slowest varying axis in multidim
# mode is the lines, which matches latitude as the first axis of the CRS.
assert srs.GetDataAxisToSRSAxisMapping() == [1, 2]
assert len(ar.GetAttributes()) == 0

# Open as classic CRS
ds = gdal.Open("/vsimem/test.zarr")
srs = ds.GetSpatialRef()
assert srs.GetAuthorityCode(None) == "4326"
# Inverted mapping in classic raster mode compared to multidim mode,
# because the first "axis" in our data model is columns.
assert srs.GetDataAxisToSRSAxisMapping() == [2, 1]

finally:
gdal.RmdirRecursive("/vsimem/test.zarr")

Expand Down Expand Up @@ -3470,7 +3476,8 @@ def set_crs():
assert ar
crs = osr.SpatialReference()
crs.ImportFromEPSG(4326)
crs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
# lat first
crs.SetDataAxisToSRSAxisMapping([1, 2])
crs.SetCoordinateEpoch(2021.2)
assert ar.SetSpatialRef(crs) == gdal.CE_None

Expand All @@ -3497,7 +3504,7 @@ def check_crs():
crs = ar.GetSpatialRef()
assert crs is not None
assert crs.GetAuthorityCode(None) == "4326"
assert crs.GetDataAxisToSRSAxisMapping() == [2, 1]
assert crs.GetDataAxisToSRSAxisMapping() == [1, 2]
assert crs.GetCoordinateEpoch() == 2021.2

check_crs()
Expand All @@ -3506,6 +3513,8 @@ def check_crs_classic_dataset():
ds = gdal.Open("/vsimem/test.zarr")
crs = ds.GetSpatialRef()
assert crs is not None
assert crs.GetAuthorityCode(None) == "4326"
assert crs.GetDataAxisToSRSAxisMapping() == [2, 1]

check_crs_classic_dataset()

Expand Down Expand Up @@ -5506,10 +5515,9 @@ def test_zarr_read_cf1():

ds = gdal.Open("data/zarr/byte_cf1.zarr")
assert ds
assert (
ds.GetSpatialRef().ExportToProj4()
== "+proj=utm +zone=11 +ellps=clrk66 +units=m +no_defs"
)
srs = ds.GetSpatialRef()
assert srs.ExportToProj4() == "+proj=utm +zone=11 +ellps=clrk66 +units=m +no_defs"
assert srs.GetDataAxisToSRSAxisMapping() == [1, 2]


###############################################################################
Expand All @@ -5526,6 +5534,40 @@ def test_zarr_read_cf1_zarrv3():
)


###############################################################################


@gdaltest.enable_exceptions()
@pytest.mark.require_proj(9)
def test_zarr_write_WGS84_and_EGM96_height(tmp_vsimem):

tmp_filename = str(tmp_vsimem / "out.zarr")
with gdal.GetDriverByName("Zarr").Create(tmp_filename, 1, 1) as ds:
srs = osr.SpatialReference()
srs.ImportFromEPSG(9707)
ds.SetSpatialRef(srs)
with gdal.Open(tmp_filename) as ds:
srs = ds.GetSpatialRef()
assert srs.GetAuthorityCode(None) == "9707"
assert srs.GetDataAxisToSRSAxisMapping() == [2, 1]


###############################################################################


@gdaltest.enable_exceptions()
def test_zarr_write_UTM31N_and_EGM96_height(tmp_vsimem):

tmp_filename = str(tmp_vsimem / "out.zarr")
with gdal.GetDriverByName("Zarr").Create(tmp_filename, 1, 1) as ds:
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:32631+5773")
ds.SetSpatialRef(srs)
with gdal.Open(tmp_filename) as ds:
srs = ds.GetSpatialRef()
assert srs.GetDataAxisToSRSAxisMapping() == [1, 2]


###############################################################################
# Test bug fix for https://github.com/OSGeo/gdal/issues/11016

Expand Down
10 changes: 7 additions & 3 deletions frmts/zarr/zarr_array.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2724,6 +2724,7 @@ void ZarrArray::ParseSpecialAttributes(
if (item.IsValid())
{
poSRS = std::make_shared<OGRSpatialReference>();
poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
if (poSRS->SetFromUserInput(
item.ToString().c_str(),
OGRSpatialReference::
Expand All @@ -2748,6 +2749,7 @@ void ZarrArray::ParseSpecialAttributes(
if (gridMappingArray)
{
poSRS = std::make_shared<OGRSpatialReference>();
poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
CPLStringList aosKeyValues;
for (const auto &poAttr : gridMappingArray->GetAttributes())
{
Expand Down Expand Up @@ -2798,10 +2800,12 @@ void ZarrArray::ParseSpecialAttributes(
}
if (iDimX > 0 && iDimY > 0)
{
if (poSRS->GetDataAxisToSRSAxisMapping() == std::vector<int>{2, 1})
const auto &oMapping = poSRS->GetDataAxisToSRSAxisMapping();
if (oMapping == std::vector<int>{2, 1} ||
oMapping == std::vector<int>{2, 1, 3})
poSRS->SetDataAxisToSRSAxisMapping({iDimY, iDimX});
else if (poSRS->GetDataAxisToSRSAxisMapping() ==
std::vector<int>{1, 2})
else if (oMapping == std::vector<int>{1, 2} ||
oMapping == std::vector<int>{1, 2, 3})
poSRS->SetDataAxisToSRSAxisMapping({iDimX, iDimY});
}

Expand Down
Loading