8000 Bad transformation using a derived vertical CRS with base vertical that has a geoid · Issue #3407 · OSGeo/PROJ · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Bad transformation using a derived vertical CRS with base vertical that has a geoid #3407

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

Jump to bottom
Closed
hernando opened this issue Oct 21, 2022 · 3 comments · Fixed by #3408 or #3411
Closed

Bad transformation using a derived vertical CRS with base vertical that has a geoid #3407

hernando opened this issue Oct 21, 2022 · 3 comments · Fixed by #3408 or #3411
Assignees
@rouault
Labels
bug

Comments

@hernando
Copy link
Contributor
hernando commented Oct 21, 2022 •
edited
Loading

If I try to declare a vertical CRS as a derived on top of a vertical for which a geoid is available, the transformation from ellipsoidal heights to my vertical CRS will include the geoid grid, but ignore the deriving conversion.

CRS='COMPOUNDCRS["CH1903+ + custom height",
GEOGCRS["CH1903+",
    DATUM["CH1903+",
        ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Geodesy."],
        AREA["Liechtenstein; Switzerland."],
        BBOX[45.82,5.96,47.81,10.49]],
    ID["EPSG",4150]],
VERTCRS["Custom Vertical",
    BASEVERTCRS["LHN95 height",
        VDATUM["Landeshohennetz 1995"]],
    DERIVINGCONVERSION["vertical offset",
        METHOD["Vertical Offset",ID["EPSG",9616]],
        PARAMETER["Vertical Offset", 10,
            LENGTHUNIT["metre",1],
            ID["EPSG",8603]]],
    CS[vertical,1],
        AXIS["gravity-related height (H)",up,
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]]]]'

Problem description

If I apply a transformation like this:

echo 46.9524055555556 7.43958333333333 400 | cs2cs EPSG:4979 "$CRS" -d 12

I get 46.953728425886 7.440534370598 351.014348853429, but I was expecting 46.953728425886 7.440534370598 361.014348853429

projinfo show these 2 candidates transformations as the first ones:

Operation No. 1:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=ch_swisstopo_chgeo2004_ETRS89_LHN95.tif
        +multiplier=1
  +step +proj=push +v_3
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=674.374 +y=15.056 +z=405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=pop +v_3
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

Operation No. 2:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=ch_swisstopo_chgeo2004_ETRS89_LHN95.tif
        +multiplier=1
  +step +proj=geogoffset +dh=10
  +step +proj=push +v_3
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=674.374 +y=15.056 +z=405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=pop +v_3
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

Expected Output

I would expect the height to be higher according to the vertical offset and projinfo listing the 2nd operation as the first one.

Environment Information

  • PROJ version: 9.1.0
  • Operation System Information: Ubuntu 20.04

Installation method

Compiled from source.
@hernando hernando added the bug label Oct 21, 2022
@rouault rouault self-assigned this Oct 21, 2022
rouault added a commit to rouault/PROJ that referenced this issue Oct 21, 2022
@rouault
VerticalCRS::_isEquivalentTo(): do not consider VerticalCRS and Deriv…
92cc928 
…edVerticalCRS as equivalent

This cause incorrect result when transforming from/to a
DerivedVerticalCRS or a CompoundCRS made of a DerivedVerticalCRS

With that fix:
```
$ echo 46.9524055555556 7.43958333333333 400 | cs2cs -d 12 EPSG:4979 "$(cat test.wkt)"
46.953728423809	7.440534369109 361.014348853429
```

Fixes OSGeo#3407
@rouault rouault mentioned this issue Oct 21, 2022
Merged
@hernando
Copy link
Contributor Author
hernando commented Oct 22, 2022 •
edited
Loading

Thanks for the quick fix. Now I have a follow up, let me know if I should write a separate ticket for that.

The problem I'm facing now is that in fact, what I want to do is to declare a custom vertical for defining an offset on the ellipsoidal height, using as ellipsoid that one from the horizontal part of the compound.
If I try this:

TGT=$(cat <<EOF
COMPOUNDCRS["CH1903+ + custom height",
$(projinfo -o WKT2:2019 -q --single-line EPSG:4150),
VERTCRS["Custom Vertical",
    BASEVERTCRS["Ellipsoid",
        VDATUM["Ellipsoid"]],
    DERIVINGCONVERSION["vertical offset",
        METHOD["Vertical Offset",ID["EPSG",9616]],
        PARAMETER["Vertical Offset", 10,
            LENGTHUNIT["metre",1],
            ID["EPSG",8603]]],
    CS[vertical,1],
        AXIS["gravity-related height (H)",up,
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]]]]
EOF
)

echo 46.95 7.44 400 | cs2cs EPSG:4979 "$TGT" --3d -d 12

The result is very bad because it has ballpark transformations left and right and I get 46.95, 7.44, 400 instead of 46.953723 7.4405343 360.38.

Projinfo shows this pipeline as the first alternative:

+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=geogoffset +dh=10
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

But I would like to see something equivalent to this:

+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=cart +ellps=WGS84
  +step +inv +proj=helmert +x=674.374 +y=15.056 +z=405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=geogoffset +dh=10
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

Is it even possible express this in a WKT2 string?

Note: The "Ellipsoid" datum has been taken from the output of projinfo with --allow-ellipsoidal-height-as-vertical-crs for WKT1 format.

@rouault rouault mentioned this issue Oct 22, 2022
Merged
@rouault
Copy link
Member
rouault commented Oct 22, 2022

The problem I'm facing now is that in fact, what I want to do is to declare a custom vertical for defining an offset on the ellipsoidal height, using as ellipsoid that one from the horizontal part of the compound.

you need to create a DerivedGeographic 3D CRS. Cf #3411 for what you want and the fix.
Side note: it seems you have a lot of interesting advanced uses of PROJ. It would be really nice to contribute a documentation page to PROJ doc with different examples. Pretty sure it would help many people.

@hernando
Copy link
Contributor Author

The problem I'm facing now is that in fact, what I want to do is to declare a custom vertical for defining an offset on the ellipsoidal height, using as ellipsoid that one from the horizontal part of the compound.

you need to create a DerivedGeographic 3D CRS. Cf #3411 for what you want and the fix.

Great! That's exactly what I was looking for. In the beginning I tried something similar with a BOUNDCRS, but after reading the WKT2 standard, I realized that a vertical offset is not a valid transformation here.

Side note: it seems you have a lot of interesting advanced uses of PROJ. It would be really nice to contribute a documentation page to PROJ doc with different examples. Pretty sure it would help many people.

I'm happy to contribute something, but I can't commit to do it in a short time frame. The experiments I'm making are all around vertical transformations in situations where there's no geoid model grid available. I could describe use cases and include some C++ snippets with WKT2 strings.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@rouault rouault

Labels
bug
Projects
None yet
Milestone
No milestone
2 participants
@hernando @rouault
0