[PROJ] Transforming between ITRF realizations

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

[PROJ] Transforming between ITRF realizations

Trey Stafford
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
Reply | Threaded
Open this post in threaded view
|

Re: Transforming between ITRF realizations

Alan Snow
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?
Alan

On 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
Reply | Threaded
Open this post in threaded view
|

Re: Transforming between ITRF realizations

Even Rouault-2
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
Reply | Threaded
Open this post in threaded view
|

Re: 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
>
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> 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

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Transforming between ITRF realizations

Markus Metz-3
In reply to this post by Even Rouault-2


On Fri, May 3, 2019 at 8:58 PM Even Rouault <[hidden email]> wrote:
>
> On vendredi 3 mai 2019 13:01:01 CEST Alan Snow wrote:
[...]
>
> 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

that means we (developers) need to check with proj_angular_input() if such a unit conversion is needed? AFAIU, units are now degrees as in EPSG:4326, and in this case a conversion to radians is no longer required, as in proj-4.

Markus M


_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj