OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.67k stars 752 forks source link

Wrong transformation from vertical with geoid model to "vertical offset" #4175

Closed jjimenezshaw closed 1 week ago

jjimenezshaw commented 2 weeks ago

Example of problem

If I define this compound system with a vertical offset, sample.txt

COMPOUNDCRS["UTM30 with vertical offset",
    PROJCRS["WGS 84 / UTM zone 30N",
        BASEGEOGCRS["WGS 84",
            ENSEMBLE["World Geodetic System 1984 ensemble",
                MEMBER["World Geodetic System 1984 (Transit)"],
                MEMBER["World Geodetic System 1984 (G730)"],
                MEMBER["World Geodetic System 1984 (G873)"],
                MEMBER["World Geodetic System 1984 (G1150)"],
                MEMBER["World Geodetic System 1984 (G1674)"],
                MEMBER["World Geodetic System 1984 (G1762)"],
                MEMBER["World Geodetic System 1984 (G2139)"],
                MEMBER["World Geodetic System 1984 (G2296)"],
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ENSEMBLEACCURACY[2.0]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]],
        CONVERSION["UTM zone 30N",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",-3,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9996,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",500000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["(E)",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["(N)",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        USAGE[
            SCOPE["Navigation and medium accuracy spatial referencing."],
            AREA["Between 6°W and 0°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Burkina Faso. Côte' Ivoire (Ivory Coast). Faroe Islands - offshore. France. Ghana. Gibraltar. Ireland - offshore Irish Sea. Mali. Mauritania. Morocco. Spain. United Kingdom (UK)."],
            BBOX[0,-6,84,0]],
        ID["EPSG",32630]],
    VERTCRS["Derived verticalCRS",
        BASEVERTCRS["Ellipsoid (metre)",
            VDATUM["Ellipsoid"]],
        DERIVINGCONVERSION["Conv Vertical Offset",
            METHOD["Vertical Offset",
                ID["EPSG",9616]],
            PARAMETER["Vertical Offset",42.3,
                LENGTHUNIT["metre",1],
                ID["EPSG",8603]]],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]]]

And I make a conversion from WGS84+EGM96

$ PROJ_NETWORK=ON projinfo -s EPSG:4326+5773 -t @sample.txt -o proj
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of Transformation from Ellipsoid (metre) to EGM96 height (ballpark vertical transformation) + Conv Vertical Offset + UTM zone 30N, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=geogoffset +dh=42.3
  +step +proj=utm +zone=30 +ellps=WGS84

I would expect that the gridfile us_nga_egm96_15.tif were used. It is not the case.

Problem description

Apparently using a "vertical offset" over the ellipsoid is removing the vertical grid file from the source while doing the transformation.

Expected Output

Usage of both +step +proj=geogoffset +dh=42.3 and +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1

Environment Information

Installation method

jjimenezshaw commented 2 weeks ago

I have the impression that the reason is the "ellipsoid" vertical reference (with a "normal" vertical reference like egm2008 it works). But I need the vertical offset over the ellipsoid exactly to "mimic" a geoid height.