Skip to content

r.in.gdal: Compare only the horizontal CRS with the project #7936

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.

Sign up for GitHub

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

Merged
petrasovaa merged 1 commit into from
Sep 24, 2026
Merged

r.in.gdal: Compare only the horizontal CRS with the project #7936

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
43 changes: 37 additions & 6 deletions raster/r.external/proj.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,41 @@

#include <proj.h>

// 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)
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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"
Expand Down
43 changes: 37 additions & 6 deletions raster/r.in.gdal/proj.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,41 @@

#include <proj.h>

// 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)
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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"
Expand Down
92 changes: 92 additions & 0 deletions raster/r.in.gdal/tests/r_in_gdal_crs_check_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"""<VRTDataset rasterXSize="10" rasterYSize="10">
<SRS>{crs}</SRS>
<GeoTransform>-79, 0.01, 0, 36, 0, -0.01</GeoTransform>
<VRTRasterBand dataType="Float32" band="1"/>
</VRTDataset>
"""
)
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)
43 changes: 37 additions & 6 deletions vector/v.external/proj.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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"
Expand Down
Loading
Loading