From 87e40b8bd9f9a61cf152596ae39e37ef50ff6adb Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Tue, 22 Sep 2026 16:58:50 -0400 Subject: [PATCH] r.in.gdal: Compare only the horizontal CRS with the project (#XXXX) Data with a compound or geographic 3D CRS was rejected even by a project created from the same file. Same fix in r.external, v.in.ogr and v.external, plus a message that the vertical component is ignored. --- raster/r.external/proj.c | 43 ++++++-- raster/r.in.gdal/proj.c | 43 ++++++-- .../tests/r_in_gdal_crs_check_test.py | 92 +++++++++++++++++ vector/v.external/proj.c | 43 ++++++-- vector/v.in.ogr/proj.c | 43 ++++++-- .../v.in.ogr/tests/v_in_ogr_crs_check_test.py | 99 +++++++++++++++++++ 6 files changed, 339 insertions(+), 24 deletions(-) create mode 100644 raster/r.in.gdal/tests/r_in_gdal_crs_check_test.py create mode 100644 vector/v.in.ogr/tests/v_in_ogr_crs_check_test.py diff --git a/raster/r.external/proj.c b/raster/r.external/proj.c index 670d1172b1f..534e330a9c4 100644 --- a/raster/r.external/proj.c +++ b/raster/r.external/proj.c @@ -8,6 +8,41 @@ #include +// GRASS has no vertical CRS, so let the user know when one is dropped. +static void message_ignored_vertical_crs(OGRSpatialReferenceH hSRS) +{ + if (OSRGetAxesCount(hSRS) != 3) + return; + + const char *name = OSRGetAttrValue(hSRS, "VERT_CS", 0); + if (name) + G_message(_("Vertical CRS <%s> of the dataset is ignored"), name); + else + G_message(_("The vertical component of the dataset CRS is ignored")); +} + +// GRASS projects have no vertical CRS, so only the horizontal component +// of a compound or geographic 3D CRS can be compared. Data axis mapping +// is handled separately and is ignored here. +static int is_same_horizontal_srs(OGRSpatialReferenceH hSRS1, + OGRSpatialReferenceH hSRS2) +{ + if (!hSRS1 || !hSRS2) + return 0; + + OGRSpatialReferenceH hSRS1_2d = OSRClone(hSRS1); + OGRSpatialReferenceH hSRS2_2d = OSRClone(hSRS2); + OSRDemoteTo2D(hSRS1_2d, NULL); + OSRDemoteTo2D(hSRS2_2d, NULL); + + const char *options[] = {"IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", NULL}; + int same = OSRIsSameEx(hSRS1_2d, hSRS2_2d, options); + OSRDestroySpatialReference(hSRS1_2d); + OSRDestroySpatialReference(hSRS2_2d); + + return same; +} + /* keep in sync with r.in.gdal, v.in.ogr, v.external */ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, int create_only, int override, int check_only) @@ -78,6 +113,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, G_asprintf(&srid, "%s:%s", authname, authcode); } } + message_ignored_vertical_crs(hSRS); } GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0); @@ -130,7 +166,6 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, char *loc_wkt = NULL, *loc_srid = NULL; void (*msg_fn)(const char *, ...); OGRSpatialReferenceH hSRS_loc = NULL; - char *papszOptions[2]; if (check_only && override) { /* can't check when over-riding check */ @@ -209,17 +244,13 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, GPJ_grass_to_osr2(loc_proj_info, loc_proj_units, loc_epsg); } - /* ignore data axis mapping, this is handled separately */ - papszOptions[0] = G_store("IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES"); - papszOptions[1] = NULL; - if (override) { cellhd->proj = loc_wind.proj; cellhd->zone = loc_wind.zone; G_message(_("Over-riding projection check")); } else if (loc_wind.proj != cellhd->proj || - !OSRIsSameEx(hSRS, hSRS_loc, (const char **)papszOptions)) { + !is_same_horizontal_srs(hSRS, hSRS_loc)) { strcpy(error_msg, _("Coordinate reference system of dataset does not" diff --git a/raster/r.in.gdal/proj.c b/raster/r.in.gdal/proj.c index 3d6aad58aae..9d70418f020 100644 --- a/raster/r.in.gdal/proj.c +++ b/raster/r.in.gdal/proj.c @@ -8,6 +8,41 @@ #include +// GRASS has no vertical CRS, so let the user know when one is dropped. +static void message_ignored_vertical_crs(OGRSpatialReferenceH hSRS) +{ + if (OSRGetAxesCount(hSRS) != 3) + return; + + const char *name = OSRGetAttrValue(hSRS, "VERT_CS", 0); + if (name) + G_message(_("Vertical CRS <%s> of the dataset is ignored"), name); + else + G_message(_("The vertical component of the dataset CRS is ignored")); +} + +// GRASS projects have no vertical CRS, so only the horizontal component +// of a compound or geographic 3D CRS can be compared. Data axis mapping +// is handled separately and is ignored here. +static int is_same_horizontal_srs(OGRSpatialReferenceH hSRS1, + OGRSpatialReferenceH hSRS2) +{ + if (!hSRS1 || !hSRS2) + return 0; + + OGRSpatialReferenceH hSRS1_2d = OSRClone(hSRS1); + OGRSpatialReferenceH hSRS2_2d = OSRClone(hSRS2); + OSRDemoteTo2D(hSRS1_2d, NULL); + OSRDemoteTo2D(hSRS2_2d, NULL); + + const char *options[] = {"IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", NULL}; + int same = OSRIsSameEx(hSRS1_2d, hSRS2_2d, options); + OSRDestroySpatialReference(hSRS1_2d); + OSRDestroySpatialReference(hSRS2_2d); + + return same; +} + /* keep in sync with r.external, v.in.ogr, v.external */ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, int create_only, int override, int check_only) @@ -78,6 +113,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, G_asprintf(&srid, "%s:%s", authname, authcode); } } + message_ignored_vertical_crs(hSRS); } GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0); @@ -130,7 +166,6 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, char *loc_wkt = NULL, *loc_srid = NULL; void (*msg_fn)(const char *, ...); OGRSpatialReferenceH hSRS_loc = NULL; - char *papszOptions[2]; if (check_only && override) { /* can't check when over-riding check */ @@ -209,17 +244,13 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, char *outloc, GPJ_grass_to_osr2(loc_proj_info, loc_proj_units, loc_epsg); } - /* ignore data axis mapping, this is handled separately */ - papszOptions[0] = G_store("IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES"); - papszOptions[1] = NULL; - if (override) { cellhd->proj = loc_wind.proj; cellhd->zone = loc_wind.zone; G_message(_("Over-riding projection check")); } else if (loc_wind.proj != cellhd->proj || - !OSRIsSameEx(hSRS, hSRS_loc, (const char **)papszOptions)) { + !is_same_horizontal_srs(hSRS, hSRS_loc)) { strcpy(error_msg, _("Coordinate reference system of dataset does not" diff --git a/raster/r.in.gdal/tests/r_in_gdal_crs_check_test.py b/raster/r.in.gdal/tests/r_in_gdal_crs_check_test.py new file mode 100644 index 00000000000..c6152e580af --- /dev/null +++ b/raster/r.in.gdal/tests/r_in_gdal_crs_check_test.py @@ -0,0 +1,92 @@ +# SPDX-License-Identifier: GPL-2.0-or-later + +"""Tests of the CRS check with compound and geographic 3D CRSs + +r.external shares the CRS check code with r.in.gdal, so both are tested. +""" + +import os + +import pytest + +import grass.script as gs +from grass.tools import Tools, ToolError + +# NAD83(2011) / North Carolina (ftUS), also with NAVD88 height (ftUS) +HORIZONTAL_CRS = "EPSG:6543" +COMPOUND_CRS = "EPSG:6543+6360" +# WGS 84 in 2D and 3D +GEOGRAPHIC_2D_CRS = "EPSG:4326" +GEOGRAPHIC_3D_CRS = "EPSG:4979" +CRSS_WITH_VERTICAL = {COMPOUND_CRS, GEOGRAPHIC_3D_CRS} +UNRELATED_CRS = "EPSG:32617" + + +def write_raster(path, crs): + """Write a 10x10 VRT raster with the given CRS (GDAL fills it with zeros) + + The extent is valid in both projected and geographic CRSs. + """ + path.write_text( + f""" + {crs} + -79, 0.01, 0, 36, 0, -0.01 + + +""" + ) + return path + + +@pytest.fixture(params=["r_in_gdal", "r_external"]) +def import_raster(request): + """Function which imports (or links) a raster file into a session + + Returns the r.info of the imported raster and the messages of the import. + """ + + def run(session, path): + tools = Tools(session=session, consistent_return_value=True) + result = getattr(tools, request.param)(input=str(path), output="imported") + return tools.r_info(map="imported", format="json"), result.stderr + + return run + + +@pytest.mark.parametrize( + ("project_crs", "file_crs"), + [ + (HORIZONTAL_CRS, COMPOUND_CRS), + (COMPOUND_CRS, COMPOUND_CRS), + (COMPOUND_CRS, HORIZONTAL_CRS), + (GEOGRAPHIC_2D_CRS, GEOGRAPHIC_3D_CRS), + (GEOGRAPHIC_3D_CRS, GEOGRAPHIC_3D_CRS), + (GEOGRAPHIC_3D_CRS, GEOGRAPHIC_2D_CRS), + ], +) +def test_same_horizontal_crs_matches(tmp_path, project_crs, file_crs, import_raster): + """A file imports when its horizontal CRS matches the project + + The user is told when the vertical component of the file CRS is ignored. + """ + file = write_raster(tmp_path / "file.vrt", file_crs) + project = tmp_path / "project" + gs.create_project(project, filename=write_raster(tmp_path / "crs.vrt", project_crs)) + with gs.setup.init(project, env=os.environ.copy()) as session: + info, messages = import_raster(session, file) + assert info["rows"] == 10 + assert ("vertical" in messages.lower()) == (file_crs in CRSS_WITH_VERTICAL) + + +def test_different_horizontal_crs_is_rejected(tmp_path, import_raster): + """Ignoring the vertical CRS does not make different horizontal CRSs match""" + unrelated_file = write_raster(tmp_path / "unrelated.vrt", UNRELATED_CRS) + project = tmp_path / "project" + gs.create_project( + project, filename=write_raster(tmp_path / "crs.vrt", COMPOUND_CRS) + ) + with ( + gs.setup.init(project, env=os.environ.copy()) as session, + pytest.raises(ToolError, match="does not appear to match"), + ): + import_raster(session, unrelated_file) diff --git a/vector/v.external/proj.c b/vector/v.external/proj.c index 5969ca34dce..572890136e0 100644 --- a/vector/v.external/proj.c +++ b/vector/v.external/proj.c @@ -47,6 +47,41 @@ OGRSpatialReferenceH get_layer_proj(OGRLayerH Ogr_layer, char *geom_col, return hSRS; } +// GRASS has no vertical CRS, so let the user know when one is dropped. +static void message_ignored_vertical_crs(OGRSpatialReferenceH hSRS) +{ + if (OSRGetAxesCount(hSRS) != 3) + return; + + const char *name = OSRGetAttrValue(hSRS, "VERT_CS", 0); + if (name) + G_message(_("Vertical CRS <%s> of the dataset is ignored"), name); + else + G_message(_("The vertical component of the dataset CRS is ignored")); +} + +// GRASS projects have no vertical CRS, so only the horizontal component +// of a compound or geographic 3D CRS can be compared. Data axis mapping +// is handled separately and is ignored here. +static int is_same_horizontal_srs(OGRSpatialReferenceH hSRS1, + OGRSpatialReferenceH hSRS2) +{ + if (!hSRS1 || !hSRS2) + return 0; + + OGRSpatialReferenceH hSRS1_2d = OSRClone(hSRS1); + OGRSpatialReferenceH hSRS2_2d = OSRClone(hSRS2); + OSRDemoteTo2D(hSRS1_2d, NULL); + OSRDemoteTo2D(hSRS2_2d, NULL); + + const char *options[] = {"IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", NULL}; + int same = OSRIsSameEx(hSRS1_2d, hSRS2_2d, options); + OSRDestroySpatialReference(hSRS1_2d); + OSRDestroySpatialReference(hSRS2_2d); + + return same; +} + /* keep in sync with r.in.gdal, r.external, v.in.ogr */ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, char *geom_col, char *outloc, int create_only, @@ -114,6 +149,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, G_asprintf(&srid, "%s:%s", authname, authcode); } } + message_ignored_vertical_crs(hSRS); } GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0); @@ -166,7 +202,6 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, char *loc_wkt = NULL, *loc_srid = NULL; void (*msg_fn)(const char *, ...); OGRSpatialReferenceH hSRS_loc = NULL; - char *papszOptions[2]; if (check_only && override) { /* can't check when over-riding check */ @@ -245,17 +280,13 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, GPJ_grass_to_osr2(loc_proj_info, loc_proj_units, loc_epsg); } - /* ignore data axis mapping, this is handled separately */ - papszOptions[0] = G_store("IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES"); - papszOptions[1] = NULL; - if (override) { cellhd->proj = loc_wind.proj; cellhd->zone = loc_wind.zone; G_message(_("Over-riding projection check")); } else if (loc_wind.proj != cellhd->proj || - !OSRIsSameEx(hSRS, hSRS_loc, (const char **)papszOptions)) { + !is_same_horizontal_srs(hSRS, hSRS_loc)) { strcpy(error_msg, _("Coordinate reference system of dataset does not" diff --git a/vector/v.in.ogr/proj.c b/vector/v.in.ogr/proj.c index 8338165f2ce..2144ac3344f 100644 --- a/vector/v.in.ogr/proj.c +++ b/vector/v.in.ogr/proj.c @@ -117,6 +117,41 @@ int cmp_layer_srs(GDALDatasetH Ogr_ds, int nlayers, int *layers, return 0; } +// GRASS has no vertical CRS, so let the user know when one is dropped. +static void message_ignored_vertical_crs(OGRSpatialReferenceH hSRS) +{ + if (OSRGetAxesCount(hSRS) != 3) + return; + + const char *name = OSRGetAttrValue(hSRS, "VERT_CS", 0); + if (name) + G_message(_("Vertical CRS <%s> of the dataset is ignored"), name); + else + G_message(_("The vertical component of the dataset CRS is ignored")); +} + +// GRASS projects have no vertical CRS, so only the horizontal component +// of a compound or geographic 3D CRS can be compared. Data axis mapping +// is handled separately and is ignored here. +static int is_same_horizontal_srs(OGRSpatialReferenceH hSRS1, + OGRSpatialReferenceH hSRS2) +{ + if (!hSRS1 || !hSRS2) + return 0; + + OGRSpatialReferenceH hSRS1_2d = OSRClone(hSRS1); + OGRSpatialReferenceH hSRS2_2d = OSRClone(hSRS2); + OSRDemoteTo2D(hSRS1_2d, NULL); + OSRDemoteTo2D(hSRS2_2d, NULL); + + const char *options[] = {"IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", NULL}; + int same = OSRIsSameEx(hSRS1_2d, hSRS2_2d, options); + OSRDestroySpatialReference(hSRS1_2d); + OSRDestroySpatialReference(hSRS2_2d); + + return same; +} + /* keep in sync with r.in.gdal, r.external, v.external */ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, char *geom_col, char *outloc, int create_only, @@ -184,6 +219,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, G_asprintf(&srid, "%s:%s", authname, authcode); } } + message_ignored_vertical_crs(hSRS); } GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0); @@ -236,7 +272,6 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, char *loc_wkt = NULL, *loc_srid = NULL; void (*msg_fn)(const char *, ...); OGRSpatialReferenceH hSRS_loc = NULL; - char *papszOptions[2]; if (check_only && override) { /* can't check when over-riding check */ @@ -315,17 +350,13 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int layer, GPJ_grass_to_osr2(loc_proj_info, loc_proj_units, loc_epsg); } - /* ignore data axis mapping, this is handled separately */ - papszOptions[0] = G_store("IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES"); - papszOptions[1] = NULL; - if (override) { cellhd->proj = loc_wind.proj; cellhd->zone = loc_wind.zone; G_message(_("Over-riding projection check")); } else if (loc_wind.proj != cellhd->proj || - !OSRIsSameEx(hSRS, hSRS_loc, (const char **)papszOptions)) { + !is_same_horizontal_srs(hSRS, hSRS_loc)) { strcpy(error_msg, _("Coordinate reference system of dataset does not" diff --git a/vector/v.in.ogr/tests/v_in_ogr_crs_check_test.py b/vector/v.in.ogr/tests/v_in_ogr_crs_check_test.py new file mode 100644 index 00000000000..4bad4695522 --- /dev/null +++ b/vector/v.in.ogr/tests/v_in_ogr_crs_check_test.py @@ -0,0 +1,99 @@ +# SPDX-License-Identifier: GPL-2.0-or-later + +"""Tests of the CRS check with compound and geographic 3D CRSs + +v.external shares the CRS check code with v.in.ogr, so both are tested. +""" + +import os + +import pytest + +import grass.script as gs +from grass.tools import Tools, ToolError + +# NAD83(2011) / North Carolina (ftUS), also with NAVD88 height (ftUS) +HORIZONTAL_CRS = "EPSG:6543" +COMPOUND_CRS = "EPSG:6543+6360" +# WGS 84 in 2D and 3D +GEOGRAPHIC_2D_CRS = "EPSG:4326" +GEOGRAPHIC_3D_CRS = "EPSG:4979" +CRSS_WITH_VERTICAL = {COMPOUND_CRS, GEOGRAPHIC_3D_CRS} +UNRELATED_CRS = "EPSG:32617" + + +def write_points(path, crs): + """Write an OGR VRT with two points read from a CSV file next to it + + The coordinates are valid in both projected and geographic CRSs. + """ + csv_path = path.with_suffix(".csv") + csv_path.write_text("x,y\n-79,36\n-78.9,35.9\n") + path.write_text( + f""" + + {csv_path.name} + {csv_path.stem} + wkbPoint + {crs} + + + +""" + ) + return path + + +@pytest.fixture(params=["v_in_ogr", "v_external"]) +def import_vector(request): + """Function which imports (or links) a vector file into a session + + Returns the v.info of the imported vector and the messages of the import. + """ + + def run(session, path): + tools = Tools(session=session, consistent_return_value=True) + result = getattr(tools, request.param)(input=str(path), output="imported") + info = tools.v_info(map="imported", flags="t", format="json") + return info, result.stderr + + return run + + +@pytest.mark.parametrize( + ("project_crs", "file_crs"), + [ + (HORIZONTAL_CRS, COMPOUND_CRS), + (COMPOUND_CRS, COMPOUND_CRS), + (COMPOUND_CRS, HORIZONTAL_CRS), + (GEOGRAPHIC_2D_CRS, GEOGRAPHIC_3D_CRS), + (GEOGRAPHIC_3D_CRS, GEOGRAPHIC_3D_CRS), + (GEOGRAPHIC_3D_CRS, GEOGRAPHIC_2D_CRS), + ], +) +def test_same_horizontal_crs_matches(tmp_path, project_crs, file_crs, import_vector): + """A file imports when its horizontal CRS matches the project + + The user is told when the vertical component of the file CRS is ignored. + """ + file = write_points(tmp_path / "file.vrt", file_crs) + project = tmp_path / "project" + gs.create_project(project, filename=write_points(tmp_path / "crs.vrt", project_crs)) + with gs.setup.init(project, env=os.environ.copy()) as session: + info, messages = import_vector(session, file) + assert info["points"] == 2 + assert ("vertical" in messages.lower()) == (file_crs in CRSS_WITH_VERTICAL) + + +def test_different_horizontal_crs_is_rejected(tmp_path, import_vector): + """Ignoring the vertical CRS does not make different horizontal CRSs match""" + unrelated_file = write_points(tmp_path / "unrelated.vrt", UNRELATED_CRS) + project = tmp_path / "project" + gs.create_project( + project, filename=write_points(tmp_path / "crs.vrt", COMPOUND_CRS) + ) + with ( + gs.setup.init(project, env=os.environ.copy()) as session, + pytest.raises(ToolError, match="does not appear to match"), + ): + import_vector(session, unrelated_file)