Skip to content
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.

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

createOperations(): fix issues when transforming between Geog3D and DerivedGeog3D CRS with Geographic3D offsets method #3411

Merged
merged 1 commit into from
Oct 23, 2022

Conversation

rouault
Copy link
Member

@rouault rouault commented Oct 22, 2022

Given test.wkt with:

GEOGCRS["CH1903+ with 10m offset on ellipsoidal height",
    BASEGEOGCRS["CH1903+",
        DATUM["CH1903+",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6150]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    DERIVINGCONVERSION["Offset on ellipsoidal height",
        METHOD["Geographic3D offsets",
            ID["EPSG",9660]],
        PARAMETER["Latitude offset",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8601]],
        PARAMETER["Longitude offset",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8602]],
        PARAMETER["Vertical Offset",10,
            LENGTHUNIT["metre",1],
            ID["EPSG",8603]]],
    CS[ellipsoidal,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
$ projinfo -s EPSG:4979 -t @test.wkt --spatial-test intersects -o PROJ --hide-ballpark
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of CH1903+ to WGS 84 (1) + Offset on ellipsoidal height, 1 m, Liechtenstein; Switzerland.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=geogoffset +dlat=0 +dlon=0 +dh=10
  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
  +step +proj=axisswap +order=2,1

Fixes #3407 (comment)

…erivedGeog3D CRS with Geographic3D offsets method

Given test.wkt with:
```
GEOGCRS["CH1903+ with 10m offset on ellipsoidal height",
    BASEGEOGCRS["CH1903+",
        DATUM["CH1903+",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6150]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    DERIVINGCONVERSION["Offset on ellipsoidal height",
        METHOD["Geographic3D offsets",
            ID["EPSG",9660]],
        PARAMETER["Latitude offset",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8601]],
        PARAMETER["Longitude offset",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8602]],
        PARAMETER["Vertical Offset",10,
            LENGTHUNIT["metre",1],
            ID["EPSG",8603]]],
    CS[ellipsoidal,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
```

```
$ projinfo -s EPSG:4979 -t @test.wkt --spatial-test intersects -o PROJ --hide-ballpark
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of CH1903+ to WGS 84 (1) + Offset on ellipsoidal height, 1 m, Liechtenstein; Switzerland.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=geogoffset +dlat=0 +dlon=0 +dh=10
  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
  +step +proj=axisswap +order=2,1
```
@rouault rouault merged commit cb759f7 into OSGeo:master Oct 23, 2022
github-actions bot pushed a commit that referenced this pull request Oct 23, 2022
createOperations(): fix issues when transforming between Geog3D and DerivedGeog3D CRS with Geographic3D offsets method
@hernando
Copy link
Contributor

hernando commented Oct 23, 2022

There's still an issue, it seems the CRS above can't be created programmatically using the C++ API. The factory function DerivedGeographicCRS::create takes a ConversionNNPtr but the coordinate operation "Geographic3D offsets" is a Transformation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Bad transformation using a derived vertical CRS with base vertical that has a geoid
2 participants