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

Add support for transformations involving coordinate epoch changes (PointMotionOperation), mostly for Canadian NAD83(CSRS) case #3884

Merged
merged 29 commits into from
Sep 15, 2023

Conversation

rouault
Copy link
Member

@rouault rouault commented Sep 5, 2023

Summary

This pull request contains API, utility and database enhancements to better support dynamic CRS, PointMotionOperation and CoordinateMetadata concepts. All those changes address mostly changes of epoch within the same dynamic CRS, typically NAD83(CSRS)v7 using the NAD83v70VG grid.

"CRS A @ epoch X" to "CRS B @ epoch Y" is implemented, for geodetic CRS, currently by doing "CRS A @ epoch X" -> "CRS A @ epoch Y" -> "CRS B @ epoch Y".

Note: time-specific Helmert, which might be needed in some situations, is not implemented

Detailed contents

  • PointMotionOperation: implement constructor, import/export WKT and PROJJSON, and instanciation from database
  • Add AuthorityFactory::createGeodeticCRSFromDatum() utility class, adapted/moved from coordinateoperationfactory logic
  • Add AuthorityFactory::getPointMotionOperationsFor() and allow a CoordinateMetadata to be built for CRS that are not officially dynamic but have a PointMotionOperation
  • proj_crs_promote_to_3D(): make it work on CoordinateMetadata as well
  • createFromUserInput(): allow 'EPSG:XXXX @ YYYY'
  • createOperations(): use point motion operations for change of epoch between CoordinateMetadata
  • proj_create_operations(): allow CoordinateMetadata <--> CoordinateMetadata transformations
  • projinfo: allow CoordinateMetadata <--> CoordinateMetadata transformations
  • Database: add 'CGVD28 height to CGVD2013a(2010) height' with NAD83(CSRS)v7 as interpolation CRS
  • projinfo and cs2cs: add --s_epoch and --t_epoch
  • Database: add a 'coordinate_metadata' table and a AuthorityFactory::createCoordinateMetadata() method (Bump database layout version to 1.3)
  • WKT/PROJJSON importer: allow to build CoordinateMetadata with CRS like NAD83(CSRS)v7
  • createFromUserInput(): make it recognize urn:ogc:def:coordinateMetadata:AUTH_NAME::CODE syntax
  • Database: add a NRCAN authority with CoordinateMetadata records
  • createOperations(): compoundToCompound: tune behaviour with point motion operation, and records under NRCAN authority
  • Implement 'Geographic3D Offset by velocity grid (NRCan byn)' transformation method (e.g for NAD83(CSRS)v7 to NAD83(CSRS)v3)
  • createOperations(): compoundToCompound: tune behaviour when one of the vertical CRS has a geoid model

Improvements done that can be useful in other contexts:

  • createOperations(): take into account GEOIDMODEL[] referenced by ID
  • CoordinateOperationFactory: make sure to preserve interpolation CRS when substituting source/target CRS
  • CoordinateOperationFactory::Private::createOperationsWithDatumPivot(): use candidateSrcGeod and candidateDstGeod of same type in priority

Tests

Note: all the needed grids are already in PROJ-data (NAD83v70VG, HT2_1997/2002/2010 geoid grids for CGVD28 height, CGG2013an83 geoid grid for CGVD2013a(1997/2002/2010) height, and HT2_xxxx_CGG2013a grids for CGVD28 to CGVD2013a heights

1. Change of epoch on geographic coordinates on NAD83(CSRS)v7 using velocity grid NAD83v70VG

$ bin/projinfo -s "NAD83(CSRS)v7" --s_epoch 2010 -t "NAD83(CSRS)v7" --t_epoch 2002
+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=GRS80
+step +proj=set +v_4=2010 +omit_fwd
+step +proj=deformation +dt=-8 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2002 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1

$ echo 60 -100 0 | PROJ_NETWORK=ON bin/cct -d 8 +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=GRS80
+step +proj=set +v_4=2010 +omit_fwd
+step +proj=deformation +dt=-8 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2002 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
60.00000011 -100.00000045 -0.06764700 2002.0000

$ echo 60 -100 0 | PROJ_NETWORK=ON bin/cs2cs -d 8 "NAD83(CSRS)v7 @ 2010.0" "NAD83(CSRS)v7 @ 2002.0"
60.00000011 -100.00000045 -0.06764700

Matches results of the official TRX tool (https://webapp.csrs-scrs.nrcan-rncan.gc.ca/geod/tools-outils/trx.php?locale=fr):

screenshot

2. Another instance

$ echo 45.42936526 -75.70165558 39.524 | PROJ_NETWORK=ON PROJ_DATA=data bin/cs2cs -d 8 "NAD83(CSRS)v7 @ 1997" "NAD83(CSRS)v7 @ 2010"
45.42936503 -75.70165526 39.55114382

Matches results of the official TRX tool (https://webapp.csrs-scrs.nrcan-rncan.gc.ca/geod/tools-outils/trx.php?locale=fr):

screenshot2

3. CompoundCRS to CompoundCRS: MTM projection + CGVD28 height @ 1997 to UTM projection + CGVD2013a(2010) height @ 2010

==> Chains HT2_1997, NAD83v70VG and CGG2013an83

$ bin/projinfo -s urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_1997_MTM7_HT2_1997 -t urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM19_CGVD2013_2010 --single-line

+proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0
+ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_HT2_1997.tif +multiplier=1
+step +proj=cart +ellps=GRS80
+step +proj=set +v_4=1997 +omit_fwd
+step +proj=deformation +dt=13 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2010 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=utm +zone=19 +ellps=GRS80

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cct -d 3 +proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0
+ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_HT2_1997.tif +multiplier=1
+step +proj=cart +ellps=GRS80
+step +proj=set +v_4=1997 +omit_fwd
+step +proj=deformation +dt=13 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2010 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=utm +zone=19 +ellps=GRS80
356670.104 5540546.584 -0.306 2010.0000

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cs2cs -d 3 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_1997_MTM7_HT2_1997 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM19_CGVD2013_2010
356670.104 5540546.584 -0.306

4. CompoundCRS to CompoundCRS: UTM projection + CGVD2013a(2002) height @ 2002 to UTM projection + CGVD2013a(2010) height @ 2010

==> Chains CGG2013an83, NAD83v70VG and CGG2013an83

$ bin/projinfo -s urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2002_UTM17_CGVD2013_2002 -t urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM17_CGVD2013_2010

+proj=pipeline
+step +inv +proj=utm +zone=17 +ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=cart +ellps=GRS80
+step +proj=set +v_4=2002 +omit_fwd
+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2010 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=utm +zone=17 +ellps=GRS80

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cct -d 3 +proj=pipeline
+step +inv +proj=utm +zone=17 +ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=cart +ellps=GRS80
+step +proj=set +v_4=2002 +omit_fwd
+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2010 +omit_inv
+step +inv +proj=cart +ellps=GRS80
+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif +multiplier=1
+step +proj=utm +zone=17 +ellps=GRS80
268955.015 5540411.986 0.047 2010.0000

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cs2cs -d 3 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2002_UTM17_CGVD2013_2002 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM17_CGVD2013_2010
268955.015 5540411.986 0.047

5. CompoundCRS to CompoundCRS: UTM projection + CGVD28 height @ 2010 to UTM projection + CGVD2013a(2010) height @ 2010

==> Uses HT2_2010v70_CGG2013a grid

$ bin/projinfo -s urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_MTM7_HT2_2010 -t urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM19_CGVD2013_2010

+proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0
+ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_HT2_2010v70_CGG2013a.tif +multiplier=-1
+step +proj=utm +zone=19 +ellps=GRS80

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cct -d 3
+proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0
+ellps=GRS80
+step +proj=vgridshift +grids=ca_nrc_HT2_2010v70_CGG2013a.tif +multiplier=-1
+step +proj=utm +zone=19 +ellps=GRS80
356670.074 5540546.618 -0.306 inf

$ echo 268955 5540412 0 | PROJ_NETWORK=ON bin/cs2cs -d 3 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_MTM7_HT2_2010 urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_2010_UTM19_CGVD2013_2010
356670.074 5540546.618 -0.306

6. Change of epoch on geographic coordinates between NAD83(CSRS)v7 and ITRF2014 using velocity grid NAD83v70VG and time-dependent Helmert transformation

$ bin/projinfo -s "NAD83(CSRS)v7 @ 2005.0" -t "ITRF2014 @ 2020.0"

+proj=pipeline
+step +proj=set +v_4=2005
+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=GRS80
+step +proj=set +v_4=2005 +omit_fwd
+step +proj=deformation +dt=15+grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80
+step +proj=set +v_4=2020 +omit_inv
+step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=-0.02678138
+ry=0.00042027 +rz=-0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006
+dz=-0.00144 +drx=-6.667e-05 +dry=0.00075744 +drz=5.133e-05
+ds=-7.201e-05 +t_epoch=2010 +convention=position_vector
+step +inv +proj=cart +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
+step +proj=set +v_4=2020

$ echo 60 -100 0 | PROJ_NETWORK=ON bin/cs2cs -d 8 "NAD83(CSRS)v7 @ 2005.0" "ITRF2014 @ 2020.0"
60.00000772 -100.00002159 -0.24708346

Matches TRX:

screenshot3


Funded by Natural Resources Canada

@rouault rouault added this to the 9.4.0 milestone Sep 5, 2023
rouault added a commit to rouault/PDAL that referenced this pull request Sep 6, 2023
…nce built from SetFromUserInput()

This is useful to get the epoch for a SpatialReference built from a
string like 'urn:ogc:def:coordinateMetadata:AUTHORTIY::CODE code' (since
OSGeo/gdal#8340 and OSGeo/PROJ#3884), which
resolves to a CoordinateMetadata object with a CRS and potentially a
coordinate epoch
@rouault rouault force-pushed the pointmotionoperation branch 2 times, most recently from 69a769b to 47ad6a3 Compare September 6, 2023 12:56
@kbevers
Copy link
Member

kbevers commented Sep 8, 2023

@rouault impressive work! I haven't had time to give it more than a shallow look at this point but I look forward to go over the PR properly when I have a bit more time on my hands.

rouault added a commit to rouault/PDAL that referenced this pull request Sep 11, 2023
…nce built from SetFromUserInput()

This is useful to get the epoch for a SpatialReference built from a
string like 'urn:ogc:def:coordinateMetadata:AUTHORTIY::CODE code' (since
OSGeo/gdal#8340 and OSGeo/PROJ#3884), which
resolves to a CoordinateMetadata object with a CRS and potentially a
coordinate epoch
Copy link
Member

@kbevers kbevers left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As far as I can tell this is solid. I did stop for a second when I noticed that the coordinate epochs (4th component) were deliberately overwritten by some of the pipelines but as far as I understand that is how point motion operations are supposed to work.

When I did the time-varying Helmert implementations I decided on not changing the coordinate epochs since it's difficult to transform your coordinates back to where they came from, but having both options available is good. Depending on the use case on can be preferable over the other.

@rouault
Copy link
Member Author

rouault commented Sep 14, 2023

As far as I can tell this is solid. I did stop for a second when I noticed that the coordinate epochs (4th component) were deliberately overwritten by some of the pipelines but as far as I understand that is how point motion operations are supposed to work.

For a transformation like "NAD83(CSRS)v7 @ 2005.0" to "ITRF2014 @ 2020.0", you need to do first a point motion operation within NAD83(CSRS)v7 to go from 2005 to 2020, and then apply the time-dependent Helmert transformation from NAD83(CSRS)v7 to ITRF2014 at epoch 2020. So the point motion operation has to set the epoch to 2020 after it has completed.

@kbevers
Copy link
Member

kbevers commented Sep 14, 2023

So the point motion operation has to set the epoch to 2020 after it has completed.

I guess that is a matter of taste. The NKG-transformations are quite similar to the examples in this PR, as far as I can tell. In those transformations we keep the original epoch of the coordinate. A difference might be that the pipelines are hand-written which may make things simpler. I can see how a more generic approach benefits from the coordinate epochs changing.

In the end what matters is how the Canadians have defined their transformations - if they have specified that epochs should change then that is how it should be implemented :-)

…mation method (e.g for NAD83(CSRS)v7 to NAD83(CSRS)v3)
…: use candidateSrcGeod and candidateDstGeod of same type in priority
…CSRS_1997_yyyy_CGVD2013_1997' or 'NAD83_CSRS_2002_xxxx_HT2_2002 to NAD83_CSRS_2002_yyyy_CGVD2013_2002'
…s, i.e. when the point motion operation is on the end, and also make sure that the time-dependent Helmert transformation receives coordinates in the right epoch
@rouault rouault merged commit d4ebd92 into OSGeo:master Sep 15, 2023
18 checks passed
rouault added a commit to rouault/PDAL that referenced this pull request Sep 15, 2023
…nce built from SetFromUserInput()

This is useful to get the epoch for a SpatialReference built from a
string like 'urn:ogc:def:coordinateMetadata:AUTHORTIY::CODE code' (since
OSGeo/gdal#8340 and OSGeo/PROJ#3884), which
resolves to a CoordinateMetadata object with a CRS and potentially a
coordinate epoch
hobu pushed a commit to PDAL/PDAL that referenced this pull request Sep 18, 2023
* SpatialReference::set(): store coordinate epoch from OGRSpatialReference built from SetFromUserInput()

This is useful to get the epoch for a SpatialReference built from a
string like 'urn:ogc:def:coordinateMetadata:AUTHORTIY::CODE code' (since
OSGeo/gdal#8340 and OSGeo/PROJ#3884), which
resolves to a CoordinateMetadata object with a CRS and potentially a
coordinate epoch

* SrsTransform: make sure to propagate coordinate epoch when rebuilding OGRSpatialReference objects

* SpatialReference: add a m_wkt2 member filled in set() and returned by getWKT2()

* SrsTransform: use SpatialReference::getWKT2()

This limits the loss of information, particularly for CRS like
'urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_1997_MTM7_HT2_1997'
which include a GEOIDMODEL[] node, that would be lost otherwise

* SpatialReference::isWKT2(): add missing keywords, and fix typo VERITCALCRS -> VERTICALCRS
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.

2 participants