[PROJ] Transforming between ITRF realizations

[Even Rouault]

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