search

mrJean1 / PyGeodesy

Pure Python geodesy tools
https://mrjean1.github.io/PyGeodesy/
297 stars 58 forks source link

Question about TRF conversions #80

Closed ggaessler closed 9 months ago

ggaessler commented 9 months ago

This is a great library, and a good documentation, however I struggled to find some examples using it. I want to convert coordinates from ITRF2014 at a specific epoch to ETRS89. Is this possible at all with this library down to a cm-level? If so, can this be achieved (using trf.trfXform) without the knowledge of additional transformation parameters? It seems like ETRS89 isn't implemented and each new transform requires a Helmert-7-tuple as input (which I don't have)?

Thanks a lot in advance.

mrJean1 commented 9 months ago

To convert between reference frames, get an ellipsoidal LatLon and use keyword arguments reframe and epoch with the source, "from" frame and the fractional calendar year, respectively. Then, use method toRefFrame to convert to convert that LatLon to the destination, "to" frame. Example 

from pygeodesy import ellipsoidalVincenty, RefFrames
p = ellipsoidalVincenty.LatLon(lat, lon, reframe=RefFrames.ITRF2014, epoch=2014.25)
q = p.toRefFrame(RefFrames.ETRS89)  # use ITRF**

For the highest precision use module ellipsoidalExact. See Karney's original C++ documentation for details about accuracy, better than 1 mm and as far as 1 nm, the roundoff.

**) There is no RefFrame named ETRS98, it is an ensemble of frames. However, there are frames like ITRF90, etc. till ITRF2014 in the RefFrames and conversion tables. The conversion coefficients for those were taken from here, for example and similar memos, see the comments in module trf.

Function trfXform is used to add a new entry to the internal conversion table in module trf. And that requires 2 reference frames plus the xform and rates coefficients, each in Helmert form.

ggaessler commented 9 months ago

Thanks a lot so far!

**) There is no RefFrame named ETRS98, it is an ensemble of frames. However, there are frames like ITRF90, etc. till ITRF2014 in the RefFrames and conversion tables. The conversion coefficients for those were taken from here, for example and similar memos, see the comments in module trf.

Regarding ETRS89, I usually have ETRS89@1989.0, which then should be ETRF1989 if I'm not mistaken. Would it be possible to convert to this frame without additional coefficients? I know it's not in the RefFrames list yet.

For testing purposes I have a dataset with exact coordinates in ETRS89@1989.0 and ITRF2014@2018.80 as ground truth, derived from the post-processing output report of a GNSS receiver vendor:

# ETRS89@1989.0
lat = dms_to_dd([48, 46, 36.89676])
lon = dms_to_dd([8, 55, 21.25713])
height = 476.049
print(f'ETRS89@1989.0\n{lat:.10f}°N, {lon:.10f}°E, +{height:.10f}m')

# ITRF2014@2018.80
lat = dms_to_dd([48, 46, 36.91314])
lon = dms_to_dd([8, 55, 21.28095])
height = 476.037
print(f'ITRF2014@2018.80\n{lat:.10f}°N, {lon:.10f}°E, +{height:.10f}m')

epoch = 2018.80

p = ellipsoidalExact.LatLon(lat, lon, height=height, reframe=RefFrames.ITRF2014, epoch=epoch)
q = p.toRefFrame(RefFrames.ETRF2000)
print(f'Conversion result\n{q.lat:.10f}°N, {q.lon:.10f}°E, +{q.height:.10f}m')

leads to

ETRS89@1989.0
48.7769157667°N, 8.9225714250°E, +476.0490000000m
ITRF2014@2018.80
48.7769203167°N, 8.9225780417°E, +476.0370000000m
Conversion result
48.7769154513°N, 8.9225706176°E, +476.0285140935m

I had to pick ETRF2000 here, this definitely points towards the correct direction for lat and lon, but of course isn't the exact frame I was looking for.

Converting to ETRS89@1989.0 involves converting to a different epoch, which I'm not sure is possible with toRefFrame since I can't specify a target epoch. This would would generally also be of interest, e.g. doing a conversion ITRF2014@2018.80 -> ITRF2008@2005.0.

mrJean1 commented 9 months ago

Without having the precise parameters for ETRF1989, it will be an approximation, unfortunately.

Your other point is very worthwhile, being able to convert from a frame@epoch to an other frame@epoch. That feature is missing in pygeodesy.trf. I'll be investigating this further in the next couple of days. Expect an update to this issue soon and keep it open.

ggaessler commented 9 months ago

Without having the precise parameters for ETRF1989, it will be an approximation, unfortunately.

I've kept on digging and found this document here: http://etrs89.ensg.ign.fr/pub/EUREF-TN-1.pdf Seems like it would provide the required parameters.

Your other point is very worthwhile, being able to convert from a frame@epoch to an other frame@epoch. That feature is missing in pygeodesy.trf. I'll be investigating this further in the next couple of days. Expect an update to this issue soon and keep it open.

That sounds great. Thanks a lot for looking into it.

mrJean1 commented 9 months ago

Great find!  All entries from those 4 tables not present in trf.py will be added./JeanOn Jan 23, 2024, at 3:32 AM, Gabriel Gaessler @.***> wrote:

Without having the precise parameters for ETRF1989, it will be an approximation, unfortunately. I've kept on digging and found this document here: http://etrs89.ensg.ign.fr/pub/EUREF-TN-1.pdf Seems like it would provide the required parameters.

Your other point is very worthwhile, being able to convert from a @. to an other @. That feature is missing in pygeodesy.trf. I'll be investigating this further in the next couple of days. Expect an update to this issue soon and keep it open. That sounds great. Thanks a lot for looking into it.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you commented.Message ID: @.***>

ggaessler commented 9 months ago

Great find!  All entries from those 4 tables not present in trf.py will be added.

Awesome - thanks a lot!

mrJean1 commented 9 months ago

One more question. The table of RefFrames needs to be extended as well, including the default epoch for each.

Perhaps, your exceptional search skills can find the proper epoch for all ETRF and ITRF 89 thru 2014 RefFrames?

mrJean1 commented 9 months ago

Below are some results from the enhanced trf.RefFrames and internal conversion table trf._trfXs.

In addition, there's a new keyword argument epoch2 for method toRefFrame (of ellipsoidal classes LatLon and Cartesian). Use epoch2 to specify a "to" epoch, different from the 'from" epoch given with keyword argument epoch**.

Using the example from your earlier post:

from pygeodesy.ellipsoidalExact import LatLon
from pygeodesy import RefFrames as RF

p = LatLon('48 46 36.89676N', '8 55 21.25713E', height=476.049, reframe=RF.ETRF89, epoch=1989)
q = p.toRefFrame(RF.ITRF2014, epoch2=2018.8)
q.toStr(form='D', prec=12)
'48.776915766667°N, 008.922571425°E, +476.05m'

p = LatLon('48 46 36.91314N', '8 55 21.28095E', height=476.049, reframe=RF.ITRF2014, epoch=2018.8)
q = p.toRefFrame(RF.ETRF89, epoch2=1989)
q.toStr(form='D', prec=12)
'48.776920316667°N, 008.922578041667°E, +476.05m'

If this works, I will push a release on GitHub. However, there are still loose ends, like the default epoch for the new RefFrames, more tests, etc.

**) The latter overrides the instance's epoch which in turn overrides the reference frame's epoch.

ggaessler commented 9 months ago

One more question. The table of RefFrames needs to be extended as well, including the default epoch for each. Perhaps, your exceptional search skills can find the proper epoch for all ETRF and ITRF 89 thru 2014 RefFrames?

I'd rather call it luck than exceptional search skills ;). To me it seems like while there are reference epochs for ITRF which are even documented in Wikipedia (https://en.wikipedia.org/wiki/International_Terrestrial_Reference_System_and_Frame) there might be no such thing for ETRF. In their documents they always talk about converting to central epochs instead and that seems to be a rather flexible approach.

While searching I also found a version of the above technical note that has been updated just a few weeks ago, also containing more recent data than the other one: http://etrs89.ensg.ign.fr/pub/EUREF-TN-24-01.pdf Comparing both documents one can see that even the euref guys don't use a fixed reference epoch for their ETRF2014 transformation parameters table - in one document it's 2010.0 and in the other 2015.0.

Below are some results from the enhanced trf.RefFrames and internal conversion table trf._trfXs.

In addition, there's a new keyword argument epoch2 for method toRefFrame (of ellipsoidal classes LatLon and Cartesian). Use epoch2 to specify a "to" epoch, different from the 'from" epoch given with keyword argument epoch**.

Using the example from your earlier post:

from pygeodesy.ellipsoidalExact import LatLon
from pygeodesy import RefFrames as RF

p = LatLon('48 46 36.89676N', '8 55 21.25713E', height=476.049, reframe=RF.ETRF89, epoch=1989)
q = p.toRefFrame(RF.ITRF2014, epoch2=2018.8)
q.toStr(form='D', prec=12)
'48.776915766667°N, 008.922571425°E, +476.05m'

p = LatLon('48 46 36.91314N', '8 55 21.28095E', height=476.049, reframe=RF.ITRF2014, epoch=2018.8)
q = p.toRefFrame(RF.ETRF89, epoch2=1989)
q.toStr(form='D', prec=12)
'48.776920316667°N, 008.922578041667°E, +476.05m'

If this works, I will push a release on GitHub. However, there are still loose ends, like the default epoch for the new RefFrames, more tests, etc.

**) The latter overrides the instance's epoch which in turn overrides the reference frame's epoch.

Sounds great, but however the results don't seem to fit.

'48.776915766667°N, 008.922571425°E, +476.05m'

is just the decimal degree representation of

p = LatLon('48 46 36.89676N', '8 55 21.25713E', height=476.049, reframe=RF.ETRF89, epoch=1989)

and no transformation has been applied if I'm not mistaken.

mrJean1 commented 9 months ago

Sorry, my mistake. Here are the results from the new test cases added to testTrf.py, on the lines marked # ****.

p = LatLon('48 46 36.89676N', '8 55 21.25713E', height=476.049, reframe=RefFrames.ETRF89, epoch=1989)
x = p.toRefFrame(RefFrames.ITRF2014, epoch2=2018.8)
self.test('Issue80', x.toStr(form='D', prec=10), '48.7769169572°N, 008.9225709519°E, +476.09m')  # ****
t = x.toRefFrame(RefFrames.ETRF89, epoch=1989)
self.test('Issue80', t.toStr(prec=6), '48°46′36.898876″N, 008°55′21.257278″E, +476.10m')

p = LatLon('48 46 36.91314N', '8 55 21.28095E', height=476.049, reframe=RefFrames.ITRF2014, epoch=2018.8)
x = p.toRefFrame(RefFrames.ETRF89, epoch2=1989)
self.test('Issue80', x.toStr(form='D', prec=10), '48.7769191262°N, 008.9225785148°E, +476.01m')  # ****
t = x.toRefFrame(RefFrames.ITRF2014, epoch=2018.8)
self.test('Issue80', t.toStr(prec=6), '48°46′36.929398″N, 008°55′21.308766″E, +476.05m')

Note, converting a LatLon to a different frame requires 2 additional conversions: toCartesian before and back toLatLon after the frame conversion.

mrJean1 commented 9 months ago

PyGeodesy 24.1.24 is still in the oven and not yet on PyPI. It includes the updates to RefFrames and the conversion table, however from the previous EUREF-TN-1 note, not yet the very latest EUREF-TN-24-1, forthcoming.

The toRefFrame method of all ellipsoidal LatLon and Cartesian classes has a new keyword argument epoch2 to specify a to epoch different from the from epoch.

Class RefFrame has a new method toRefFrame to convert an elipsoidal LatLon or Cartesian instance from one reframe @ epoch to an other reframe2 @ epoch2.

New function trfTransforms can be used to check whether conversion from reframe @ epoch to an other reframe2 @ epoch2 is available.

Lastly, there are still only a few test cases in test/testTrf.py. Those validate that converting from reframe @ epoch to an other reframe2 @ epoch2 and back to reframe @ epoch produces the expected result. However, there are still questions about the accuracy, especially for Cartesians.

HTH

mrJean1 commented 9 months ago

PyGeodesy 24.1.24 has been modified to handle additional test/testTrf.py test cases properly. Class Transform has a new property isunity. Cartesian and LatLon instances are converted to and from the ReFrame datum(s) involved in the RefFrame conversion.

mrJean1 commented 9 months ago

Release 24.2.4 on GitHub contains all the changes in module trf to improve accuracy and coverage. Following are the results for Example 1, Appendix B of EUREF-TN-1-31-Jan-2024.

Line 388 is the TransformXform instance (T) obtained from new function trfTransform0. Line 389 are the transformed X, Y and Z coordinates from method T.transform with the expected values (taken from EUREF-TN). Line 390 are the X, Y and Z differences and the absolute error, all below 0.1 mm in all cases.

Line 391 shows the velocities in mm/year computed by method T.velocities using equation (3) from EUREF-TN. The velocity errors in line 392 are too far off. It is unclear what the root cause of these discrepancies is.

388 Transform0: name='ITRF2020@2015xETRF2020@2010', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=8.75573508e-09, ry=5.28398431e-08, rz=-7.66635874e-08, s=0.0, sx=0.001806, sy=0.010899, sz=-0.015813
389 Transform0c: (4027893.9585, 307045.5550, 4919474.9620)  FAILED, KNOWN, expected (4027893.9585, 307045.5550, 4919474.9619)
390     Error0c: (-0.000016463, 0.0000336, 0.000055143) 6.664e-05
391 Transform0v: (0.08600, 0.51900, -0.75300)  FAILED, KNOWN, expected (0.00011, 0.00011, 0.00024)
392     Error0v: (0.08589, 0.51889, -0.75324) 0.9187

393 Transform0: name='-ITRF2014@2015xETRF2020@2010', tx=-0.0014, ty=-0.0004, tz=0.0004, s1=1.0, rx=-8.75573508e-09, ry=-5.28398431e-08, rz=7.66635874e-08, s=-0.00042, sx=-0.001806, sy=-0.010899, sz=0.015813
394 Transform0c: (4027893.6719, 307045.9063, 4919475.1704)  FAILED, KNOWN, expected (4027893.6719, 307045.9064, 4919475.1704)
395     Error0c: (0.000024786, -0.000062539, -0.000021305) 7.056e-05
396 Transform0v: (-0.08600, -0.51900, 0.75300)  FAILED, KNOWN, expected (-0.01361, 0.01676, 0.01044)
397     Error0v: (-0.072390003, -0.535759998, 0.742560001) 0.9185

398 Transform0: name='ITRF2014@2015xETRF2014@2010', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=8.65392421e-09, ry=5.40615736e-08, rz=-7.83943722e-08, s=0.0, sx=0.001785, sy=0.011151, sz=-0.01617
399 Transform0c: (4027893.9619, 307045.5481, 4919474.9553)  FAILED, KNOWN, expected (4027893.9620, 307045.5480, 4919474.9553)
400     Error0c: (-0.00007476, 0.000063039, 0.000002882) 9.783e-05
401 Transform0v: (0.08500, 0.53100, -0.77000)  FAILED, KNOWN, expected (0.00020, -0.00030, 0.00020)
402     Error0v: (0.0848, 0.5313, -0.7702) 0.9395

403 Transform0: name='-ITRF2000@2015xETRF2014@2010', tx=0.0007, ty=0.0012, tz=-0.0261, s1=1.0, rx=-8.65392421e-09, ry=-5.40615736e-08, rz=7.83943722e-08, s=0.00212, sx=-0.001785, sy=-0.011151, sz=0.01617
404 Transform0c: (4027893.6812, 307045.9082, 4919475.1547)
405     Error0c: (0.000013934, -0.000012081, 0.000026423) 3.222e-05
406 Transform0v: (-0.08500, -0.53100, 0.77000)  FAILED, KNOWN, expected (-0.01307, 0.01690, 0.00908)
407     Error0v: (-0.071929905, -0.547900009, 0.760920005) 0.9404

408 Transform0: name='ITRF2000@2015xETRF2000@2010', tx=0.054, ty=0.051, tz=-0.048, s1=1.0, rx=8.24668072e-09, ry=4.98873278e-08, rz=-8.06342114e-08, s=0.0, sx=0.001701, sy=0.01029, sz=-0.016632
409 Transform0c: (4027894.0054, 307045.5938, 4919474.9083)  FAILED, KNOWN, expected (4027894.0053, 307045.5939, 4919474.9083)
410     Error0c: (0.000077874, -0.000055372, -0.000008743) 9.595e-05
411 Transform0v: (0.08100, 0.49000, -0.79200)  FAILED, KNOWN, expected (-0.00020, -0.00050, -0.00036)
412     Error0v: (0.0812, 0.4905, -0.79164) 0.9348
mrJean1 commented 9 months ago

After adding an option to scale the velocities from mm/year (to meter/year by default) in the method TransformXform.velocities, the velocity errors are below 3 cm/year and less than 1 mm/year for most X, Y and Z velocities. For comparison ...

...
391 Transform0v: (0.00009, 0.00052, -0.00075)  FAILED, KNOWN, expected (0.00011, 0.00011, 0.00024)
392     Error0v: (-0.000024, 0.000409, -0.000993) 0.001074
...
396 Transform0v: (-0.00009, -0.00052, 0.00075)  FAILED, KNOWN, expected (-0.01361, 0.01676, 0.01044)
397     Error0v: (0.013524, -0.017279, -0.009687) 0.02399
...
401 Transform0v: (0.00009, 0.00053, -0.00077)  FAILED, KNOWN, expected (0.00020, -0.00030, 0.00020)
402     Error0v: (-0.000115, 0.000831, -0.00097) 0.001282
...
406 Transform0v: (-0.00008, -0.00053, 0.00077)  FAILED, KNOWN, expected (-0.01307, 0.01690, 0.00908)
407     Error0v: (0.012985, -0.017431, -0.00831) 0.02327
...
411 Transform0v: (0.00008, 0.00049, -0.00079)  FAILED, KNOWN, expected (-0.00020, -0.00050, -0.00036)
412     Error0v: (0.000281, 0.00099, -0.000432) 0.001116
mrJean1 commented 9 months ago

One usage example, in particular when velocities and several conversions are needed. 

>>> from pygeodesy import ellipsoidalExact as E, RefFrames as R, trfTransform0

>>> T = trfTransform0('ETRF89', R.ITRF2014, epoch2=2018.8)  # get the transform, once
TransformXform(name='-ITRF2014@2010xITRF89@1989+ITRF89@1989xETRF89@2018.800', tx=-0.0283, ty=-0.046, tz=0.0615, s1=1.0, rx=0.0, ry=0.0, rz=7.757e-10, s=-0.00567, sx=0.0, sy=0.0, sz=0.00016)

>>> T.toRefFrame(E.Cartesian(4160476.944, 653192.600, 4774604.455))  # convert a cartesian
Cartesian(4160476.892, 653192.554, 4774604.489)

>>> T.toRefFrame(E.LatLon('48 46 36.89676', '08 55 21.25713', height=467.049))  # convert a geodetic
LatLon(48°46′36.9″N, 008°55′21.26″E, +467.05m)

>>> print(T.velocities())  # get the velocities
(-0.00011, -0.00057, 0.00069)

PS) There is one issue: the Cartesian and LatLon above will have incorrect refname and epoch attributes. This has been fixed in release 24.2.6. Also, use release 24.2.6 or later for the improved method toRefFrame of class TransformXform.

ggaessler commented 9 months ago

Again, thanks a lot for all the work done here and all the communication off and on record here. While the accuracy is not where I initially hoped it would be this still helps a lot. I'll close the ticket for now and open a new one in case I have some new ideas regarding improving the accuracy at a later point in time.

mrJean1 commented 8 months ago

FYI, PyGeodesy 24.3.9 includes a new method transform2 for class trf.TransformXform to convert a position with velocities from one reference frame to an other.

The differences with the Examples 1 and 2 in Appendix B of EUREF-TN-1-31-Jan-2024 are all below 2 cm respectively 1 mm/year for smaller epoch ranges.