Bugfixes for transform_itrf
#45
Conversation
* Support plate-motion-model (PMM) coordinate propagation when there's no need to transform the ITRF. * Remove `+inv` from pmm step. I think this performs the reverse of the operation we want! * Update to the latest version of pyproj. This isn't a bugfix but I wanted to ensure I had the latest version of `pyproj` to test against.
| @@ -60,7 +62,7 @@ def transform_itrf( | |||
| transformed_chunks = [] | |||
| for source_itrf, chunk in data.groupby(by="ITRF"): | |||
| # If the source ITRF is the same as the target for this chunk, skip transformation. | |||
| if source_itrf == target_itrf: | |||
| if source_itrf == target_itrf and target_epoch is None: | |||
There was a problem hiding this comment.
We only skip transforming a chunk if both the source and target ITRF match, and no target epoch is given.
This allows us to do plate-motion coordinate propagation within a single ITRF (e.g., IS2 is in ITRF2014, and we want to compare to other data in ITRF2014, but we want to be able to change the epoch to account for plate motion over time).
| # provided plate motion model (PMM). An example is given in the | ||
| # message of this commit: | ||
| # https://github.com/OSGeo/PROJ/commit/403f930355926aced5caba5bfbcc230ad152cf86 | ||
| f"+step +init={target_itrf}:{plate} +t_epoch={target_epoch} " |
There was a problem hiding this comment.
I've removed the +inv step here. I think this is correct!
The +inv operator reverses the operation, which I think was placed here by mistake. We want to use the target ITRF's plate motion model to do a coordinate propagation to the target epoch, not the reverse.
See the linked commit's message for an example, which does not utilize the +inv operator: OSGeo/PROJ@403f930
There was a problem hiding this comment.
This is a complete guess, but maybe we put the +inv operator here because we have it in the ITRF transformation step:
+step +inv +init={target_itrf}:{source_itrf}
We do this for the ITRF transformation step because the +init= syntax looks up an init file with the corresponding entry (see https://proj.org/en/9.3/resource_files.html#init-files). There is not an init file for every ITRF, and for valkyrie we're mostly looking at doing ITRF conversions to ITRF2014.
So our pipeline is often constructed to look like this:
+step +inv +init=ITRF2014:ITRF2008
This looks up the ITRF2014 data file and finds the TIRF2008 transformation parameters (parameters to go from ITRF2014 to ITRF2008). By adding the +inv, we reverse the transformation so that our ITRF2008 data are converted to ITRF2014.
There was a problem hiding this comment.
Here's a proj mailing list thread about the correct syntax for plate motion models: https://lists.osgeo.org/pipermail/proj/2022-February/010527.html
In the example, a set of coordinates at epoch 2014.0 are propagated to 2022.0 using the ITRF2014:EURA model:
echo 3496746.5664 743256.4272 5264477.0556 2014.0 | cct +init=ITRF2014:EURA +t_epoch=2022.0
This would seem to confirm that the +inv should indeed be removed!
But later, in this message from Even: https://lists.osgeo.org/pipermail/proj/2022-February/010530.html
apply the ITRF2014:NOAM model using a target epoch of 2010 (hence the +inv, as naturally the transformation would use +t_epoch as the "source") for a point at observation epoch of 2022.1
echo -100 40 0 2022.1 | cct -d 8 +proj=pipeline \
+step +proj=cart +ellps+GRS80 \
+step +inv +init=ITRF2014:ITRF2000 \
+step +inv +init=ITRF2014:NOAM +t_epoch=2010 \
+step +proj=set +v_4=2010 \
+step +init=ITRF2014:ITRF2000 \
+step +inv +proj=cart +ellps+GRS80
| # https://github.com/OSGeo/PROJ/blob/master/data/ITRF2014). The | ||
| # `+inv` reverses the transformation. So `+init=ITRF2014:ITRF2008` | ||
| # performs a helmert transform from ITRF2008 to ITRF2014. | ||
| itrf_transformation_step = f"+step +inv +init={target_itrf}:{source_itrf} " |
There was a problem hiding this comment.
I want to think more about how we can support ITRF transformations more broadly. E.g., if someone passes in a target_itrf of ITRF2008 and a source_itrf of ITRF2014, this will fail because there is no ITRF2008 data file with a defined transformation to ITRF2014. But there is an ITRF2014 data file with transformation parameters for ITRF2008. We'd just drop the +inv here to support that.
I think the reason we did it this way in valkyrie was probably because we always expected to go from older -> newer, and e.g., IS2 was planned to use ITRF2014.
There was a problem hiding this comment.
We could have some code that tests each way for validity. E.g., if +init=ITRF2008:ITRF2014 fails, swap and add +inv: +inv +init=ITRF2014:ITRF2008.
This may not always work however. To make things a little more complicated, the data files that make the +init= approach work are limited to a selection of ITRFs and is dependent on the version of proj that gets installed.
The latest version of proj has a transformations file available for ITRF2020, for example, but the version I've been testing against does not.
A user could construct their own pipeline and pass in the helmert transformation parameters directly, or maybe we add that as an option here.
There was a problem hiding this comment.
WIP implementation that tests each way (with and without +inv) for validity: #46
| # Performs unit conversion from lon/lat degrees to radians. | ||
| # TODO: This step appears to be unnecessary. Removing it does not appear to | ||
| # affect the output. The following steps require that the | ||
| # coordinates be geodedic, which could be radians or degrees. | ||
| f"+step +proj=unitconvert +xy_in=deg +xy_out=rad " |
There was a problem hiding this comment.
I did some testing and this step does not appear to be necessary. Proj seems to assume lat/lon inputs by default and is "unit-aware" in that it will raise an error between steps if it detects an incompatibility.
I don't think it hurts to be here, but if not necessary I'm inclined to remove it to simplify the transformation.
$ echo "1 1 1 1" | cct -v -d 14 +proj=pipeline +ellps=WGS84 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=latlon +step +proj=cart +step +inv +init=ITRF2014:ITRF2008 +step +init=ITRF2014:ANTA +t_epoch=2019.7 +step +inv +proj=cart +step +proj=unitconvert +xy_in=rad +xy_out=deg
0.99961901119316 0.99981772606697 1.37719855178148 1.0000
is the same as:
$ echo "1 1 1 1" | cct -v -d 14 +proj=pipeline +ellps=WGS84 +step +proj=latlon +step +proj=cart +step +inv +init=ITRF2014:ITRF2008 +step +init=ITRF2014:ANTA +t_epoch=2019.7 +step +inv +proj=cart
0.99961901119316 0.99981772606697 1.37719855178148 1.0000
There was a problem hiding this comment.
I think we got the outline from our pipeline from this response from Even Rouault in a proj thread that I started: https://lists.osgeo.org/pipermail/proj/2019-May/008509.html
In a subsequent message Alan Snow notes:
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.
Which would seem to suggest that the unit conversion step could be skipped.
There was a problem hiding this comment.
The examples given in this PR to add the ITRF2020 params to proj shows the unit conversion steps. This may be a best-practice. See OSGeo/PROJ#4235.
| f"+step +proj=unitconvert +xy_in=deg +xy_out=rad " | ||
| # This step explicitly sets the projection as lat/lon. It won't | ||
| # change the coordinates, but they will be identified as geodetic, | ||
| # which is necessary for the next steps. | ||
| f"+step +proj=latlon " |
There was a problem hiding this comment.
This doesn't do anything aside from declare the coordinates to be geodetic, which is required by the subsequent steps. Not technically necessary but maybe better to be explicit in this case?
|
TODO: re-generate the example notebook |
|
+invfrom pmm step. I think this performs the reverse of the operation we want!pyprojto test against.📚 Documentation preview 📚: https://iceflow--45.org.readthedocs.build/en/45/