

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


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

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


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 1e6 to 1e5 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 timedependent 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


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 1e6 to 1e5 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 timedependent 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


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
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 1e6 to 1e5 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 timedependent 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


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 proj4.
Markus M
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj

