# [PROJ] Transforming between ITRF realizations

6 messages
Open this post in threaded view
|

## [PROJ] Transforming between ITRF realizations

 Hello, I have a number of elevation datasets that each use a different ITRF realization. The coordinates in these datasets are referenced to the WGS84 ellipsoid. An example lon, lat, elev(meters above ellipsoid) tuple would be: -117.61748591770593, 59.496749040316935, 329.39612 According to the proj4 documentation on Helmert transformations (https://proj4.org/operations/transformations/helmert.html), I should be able to use something like this to convert between ITRF2000 and ITRF93: proj +proj=helmert +convention=position_vector +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195 +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001 +rx=-0.00039 +ry=0.00080 +rz=-0.00114 +drx=-0.00011 +dry=-0.00019 +drz=0.00007 +t_epoch=1988.0 However, entering my coordinates yield the following result: 0.05    -0.01 329.39612 I've used proj4 previously, but not to do transformations between different ITRF realizations. I suspect I am missing something, but I have not had revelations yet. Any suggestions or insight would be greatly appreciated! Thanks, Trey Stafford _______________________________________________ PROJ mailing list [hidden email] https://lists.osgeo.org/mailman/listinfo/proj
Open this post in threaded view
|

## Re: Transforming between ITRF realizations

 Hi Trey,I tried it with pyproj and I got one answer:>>> from pyproj import Transformer>>> trans = Transformer.from_pipeline("+init=ITRF2000:ITRF93")>>> trans.transform(-117.61748, 59.4967, 329.396)(-111.83893294174032, 59.90053033616967, 330.56817290107665)Then I used the `cct` tool and it seems to be slightly different:\$ cct --versioncct: Rel. 6.1.0, September 1st, 2019\$ cct +init=ITRF2000:ITRF93-117.61748, 59.4967, 329.396    -117.6048        59.5032      329.3751           infDo you have a result you are expecting?AlanOn Fri, May 3, 2019 at 12:33 PM Trey Stafford <[hidden email]> wrote:Hello, I have a number of elevation datasets that each use a different ITRF realization. The coordinates in these datasets are referenced to the WGS84 ellipsoid. An example lon, lat, elev(meters above ellipsoid) tuple would be: -117.61748591770593, 59.496749040316935, 329.39612 According to the proj4 documentation on Helmert transformations (https://proj4.org/operations/transformations/helmert.html), I should be able to use something like this to convert between ITRF2000 and ITRF93: proj +proj=helmert +convention=position_vector +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195 +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001 +rx=-0.00039 +ry=0.00080 +rz=-0.00114 +drx=-0.00011 +dry=-0.00019 +drz=0.00007 +t_epoch=1988.0 However, entering my coordinates yield the following result: 0.05    -0.01 329.39612 I've used proj4 previously, but not to do transformations between different ITRF realizations. I suspect I am missing something, but I have not had revelations yet. Any suggestions or insight would be greatly appreciated! Thanks, Trey Stafford _______________________________________________ PROJ mailing list [hidden email] https://lists.osgeo.org/mailman/listinfo/proj -- Alan Snow _______________________________________________ PROJ mailing list [hidden email] https://lists.osgeo.org/mailman/listinfo/proj
Open this post in threaded view
|

## Re: Transforming between ITRF realizations

 On vendredi 3 mai 2019 13:01:01 CEST Alan Snow wrote: > Hi Trey, > > I tried it with pyproj and I got one answer: > >>> from pyproj import Transformer > >>> trans = Transformer.from_pipeline("+init=ITRF2000:ITRF93") > >>> trans.transform(-117.61748, 59.4967, 329.396) > > (-111.83893294174032, 59.90053033616967, 330.56817290107665) > > Then I used the `cct` tool and it seems to be slightly different: > \$ cct --version > cct: Rel. 6.1.0, September 1st, 2019 > \$ cct +init=ITRF2000:ITRF93 > -117.61748, 59.4967, 329.396 >     -117.6048        59.5032      329.3751           inf > > Do you have a result you are expecting? None of the above make sense (not sure why you don't get the same result as pyproj as with cct. You must perform some additional conversion: unit ?). For transformations between ITRF realizations, the difference between the input and output coordinates should be very small (typically of the order of one metre or less, ie in the range of 1e-6 to 1e-5 degree) The +init=ITRF2000:ITRF93 definition uses a +proj=helmert operation which operates in the cartesian/geocentric coordinate space. So you need to convert from geographic to geocentric coordinates first. \$ echo "-117.61748, 59.4967, 329.396" | src/cct -d 7 +proj=pipeline \   +step +proj=unitconvert +xy_in=deg +xy_out=rad \   +step +proj=cart +ellps=GRS80 \   +step +init=ITRF2000:ITRF93 \   +step +inv +proj=cart +ellps=GRS80 \   +step +proj=unitconvert +xy_in=rad +xy_out=deg  -117.6174799     59.4967002   329.3845529           inf Note also that being a time-dependent operation, you should generally supply the observation time as a 4th coordinate. If you do not do so, +proj=helmert will assume that the observation time is the same as the +t_epoch parameter, here 1988. You can also get the same result with PROJ 6 using the EPSG dataset with: echo "59.4967 -117.61748 329.396 1988" | src/cs2cs -f "%.7f" ITRF2000 ITRF93 59.4967002 -117.6174799 329.3845529 1988 ITRF2000 will resolve to EPSG:7909 (ITRF2000 Geographic3D CRS) and ITRF93 to EPSG:7905 (ITRF93 Geographic3D CRS). Note the lat, long order due to using EPSG definitions. The explicit 1988 here is needed since otherwise it would default to 0 (bug I just fixed in PROJ master / 6.1dev) Even -- Spatialys - Geospatial professional services http://www.spatialys.com_______________________________________________ PROJ mailing list [hidden email] https://lists.osgeo.org/mailman/listinfo/proj
Open this post in threaded view
|

## Re: Transforming between ITRF realizations

 All, Thanks for the feedback, this is super helpful! I was definitely missing the conversion to geocentric coordinates. Note that the proj4 docs do actually have a discussion about using pipelines for exactly this purpose, which I neglected to see before submitting this question: https://proj4.org/usage/transformation.htmlThe only thing that I think needs to be changed in Even's example for my use case is the ellipsoid that is used. Instead of GRS80, my data are referenced to WGS84. I will also adjust my input to use the observation time as the fourth parameter. Each of my observations does have a datetime associated with it, so that should be easy enough to do. I'm going to do a bit more reading to ensure I understand each of these steps, but I think I have a better grasp now than I did! Thanks again, this is much appreciated! Trey On 5/3/19 12:58 PM, Even Rouault wrote: > On vendredi 3 mai 2019 13:01:01 CEST Alan Snow wrote: >> Hi Trey, >> >> I tried it with pyproj and I got one answer: >>>>> from pyproj import Transformer >>>>> trans = Transformer.from_pipeline("+init=ITRF2000:ITRF93") >>>>> trans.transform(-117.61748, 59.4967, 329.396) >> (-111.83893294174032, 59.90053033616967, 330.56817290107665) >> >> Then I used the `cct` tool and it seems to be slightly different: >> \$ cct --version >> cct: Rel. 6.1.0, September 1st, 2019 >> \$ cct +init=ITRF2000:ITRF93 >> -117.61748, 59.4967, 329.396 >>      -117.6048        59.5032      329.3751           inf >> >> Do you have a result you are expecting? > None of the above make sense (not sure why you don't get the same result as > pyproj as with cct. You must perform some additional conversion: unit ?). For > transformations between ITRF realizations, the difference between the input > and output coordinates should be very small (typically of the order of one > metre or less, ie in the range of 1e-6 to 1e-5 degree) > > The +init=ITRF2000:ITRF93 definition uses a +proj=helmert operation which > operates in the cartesian/geocentric coordinate space. So you need to convert > from geographic to geocentric coordinates first. > > \$ echo "-117.61748, 59.4967, 329.396" | src/cct -d 7 +proj=pipeline \ >    +step +proj=unitconvert +xy_in=deg +xy_out=rad \ >    +step +proj=cart +ellps=GRS80 \ >    +step +init=ITRF2000:ITRF93 \ >    +step +inv +proj=cart +ellps=GRS80 \ >    +step +proj=unitconvert +xy_in=rad +xy_out=deg > >   -117.6174799     59.4967002   329.3845529           inf > > Note also that being a time-dependent operation, you should generally supply > the observation time as a 4th coordinate. If you do not do so, +proj=helmert > will assume that the observation time is the same as the +t_epoch parameter, > here 1988. > > You can also get the same result with PROJ 6 using the EPSG dataset with: > > echo "59.4967 -117.61748 329.396 1988" | src/cs2cs -f "%.7f" ITRF2000 ITRF93 > 59.4967002 -117.6174799 329.3845529 1988 > > ITRF2000 will resolve to EPSG:7909 (ITRF2000 Geographic3D CRS) and ITRF93 to > EPSG:7905 (ITRF93 Geographic3D CRS). > Note the lat, long order due to using EPSG definitions. > The explicit 1988 here is needed since otherwise it would default to 0 (bug I > just fixed in PROJ master / 6.1dev) > > Even > _______________________________________________ PROJ mailing list [hidden email] https://lists.osgeo.org/mailman/listinfo/proj