[PROJ] Transforming between ITRF realizations

[Trey Stafford]

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.html

The 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
>