[PROJ] Transforming between ITRF realizations

[Alan Snow]

Thanks for the clarification Even, that helped make this much more clear. I
have much to learn!

It seems everything works properly in pyproj with the year added and
upgrading to the newest version of PROJ.
Just an update.
>>> from pyproj import Transformer

Testing the correct solution for this thread:
>>> proj_str = (
...     "+proj=pipeline "
...     "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
...     "+step +proj=cart +ellps=WGS84 "
...     "+step +init=ITRF2000:ITRF93 "
...     "+step +inv +proj=cart +ellps=WGS84 "
...     "+step +proj=unitconvert +xy_in=rad +xy_out=deg"
... )
>>> trans = Transformer.from_pipeline(proj_str)
>>> trans.transform(-117.61748591770593, 59.496749040316935, 329.39612,
1988)
(-117.61748584159756, 59.49674923319888, 329.3846729239449, 1988.0)

Just testing the cs2cs functionality (needed to add support for ITRF
initialization in pyproj - PR for that submitted):
>>> crs_trans = Transformer.from_crs("ITRF2000", "ITRF93")
>>> crs_trans.transform(59.496749040316935, -117.61748591770593, 329.39612,
1988)
(59.49674923319888, -117.61748584159756, 329.3846729239449, 1988.0)

The updated version of the idea from earlier with strange results:
>>> trans = Transformer.from_pipeline("+init=ITRF2000:ITRF93")
>>> trans.transform(-117.61748591770593, 59.496749040316935, 329.39612, 0)
(-111.83893885919258, 59.900579376362764, 330.5682929011372, 0.0)
>>> trans.transform(-117.61748591770593, 59.496749040316935, 329.39612,
1988)
(-117.60478454066292, 59.50325042920626, 329.37522098600806, 1988.0)

Very informative thread.

Best,
Alan

On Fri, May 3, 2019 at 3:10 PM Trey Stafford <trey.stafford at nsidc.org>
wrote:

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


-- 
Alan Snow
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190504/9f9ea5cb/attachment.html>