How to deal with a vertical CRS in data to be imported?

[metzm]

For a given GRASS project with a regular horizontal CRS, e.g. EPSG:25833 (ETRS89 / UTM zone 33N), import of data with a compound CRS consisting of the same horizontal CRS and an additional vertical CRS, e.g. ETRS89 / UTM zone 33N + ETRS89 / UTM zone 33N, EPSG:25833 + EPSG:7837, fails with

Coordinate reference system of dataset does not appear to match current
project.

Is this a desired behaviour?

A vertical CRS could be stripped from the compound CRS with GDAL fns OSRIsCompound() + OSRStripVertical() if the GRASS project has only a horizontal CRS and external data have a compound CRS consisting of a horizontal and a vertical CRS. Would this be an acceptable solution?

[metzm]

Ignoring any vertical CRS is not an option. GRASS projection info must be extended to support also vertical CRSs.

A test with data with a compound CRS, ETRS89 / UTM zone 33N + ETRS89 / UTM zone 33N, EPSG:25833 + EPSG:7837, converted to EPSG:25833 without a vertical CRS showed substantial differences in height values:
original with vertical CRS

    STATISTICS_MINIMUM=34.669998168945
    STATISTICS_MAXIMUM=75.669998168945
    STATISTICS_MEAN=40.936102139794

warped, now without vertical CRS

    STATISTICS_MINIMUM=75.505447387695
    STATISTICS_MAXIMUM=116.48053741455
    STATISTICS_MEAN=81.772083701518

[nilason]

Maybe https://proj.org/en/stable/operations/operations_computation.html#compound-crs-to-a-geographic-crs can provide some hints?

[nilason]

See also https://proj.org/en/stable/development/reference/functions.html#c.proj_create_operations:

Starting with PROJ 9.1, vertical transformations are only done if both source CRS and target CRS are 3D CRS or Compound CRS with a vertical component. You may need to use proj_crs_promote_to_3D().

[metzm]

AFAIU, a vertical CRS applies only to elevation data because the vertical CRS provides a reference to these elevation data. Terefore the CRS of a GRASS project must not contain a vertical CRS because it can contain any kind of data.

If some data to be imported come with a vertical CRS, the information of this vertical CRS needs to be added to the data such that we can get rid off the vertical CRS. It seems that gdalwarp or now gdal raster reproject can do that. See also https://proj.org/en/stable/operations/operations_computation.html#vertical-crs-to-a-geographic-crs

[metzm]

@nilason you are right, this https://proj.org/en/stable/operations/operations_computation.html#compound-crs-to-a-geographic-crs would be the way to go.

But there are some mysteries with the test data, for example https://data.geobasis-bb.de/geobasis/daten/bdom/tif/bdom_33366-5784.zip

Using a simple gdalwarp -t_srs EPSG:25833 ... produces different results depending on the GDAL/PROJ versions.

• with GDAL 3.8.4, PROJ 9.4.0, the vertical CRS is removed, data are unchanged
• with GDAL 3.12.2, PROJ 9.7.1, the vertical CRS is removed, data are lifted up by about 40 meters

Compared to other DEMs without a vertical CRS, the original raster values in the example, ignoring the vertical CRS, seem to be correct. Lifting them up (like adding the vertical grid) seems to produce incorrect results.

projinfo -s -t from the compound CRS to the same horizontal CRS EPSG:25833 without the vertical CRS produces a noop pipeline, nothing changes.

projinfo -s -t from the compound CRS to the included BASEGEOGCRS ETRS89 produces different possible pipelines, all with a vertical shift.

Therefore the current behaviour that even r.import fails with an error about CRS mismatch is not the worst solution: the data need to be carefully investigated before working with them and maybe need modification with GDAL tools first, before import into GRASS.