Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfixes for transform_itrf #45

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Bugfixes for transform_itrf #45

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0d3366f
Bugfixes for `transform_itrf`
trey-stafford Oct 18, 2024
26a0cb6
Update ITRF-related docstring and comments for pipeline steps
trey-stafford Oct 21, 2024
da12291
Update comment w/ link
trey-stafford Oct 21, 2024
27d164c
Rerender iceflow example notebook
trey-stafford Oct 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ dependencies = [
"pandas >=2.2",
"h5py >=3.11",
"gps-timemachine >=1.1.4",
"pyproj >=3.6.1",
"pyproj >=3.7.0",
"shapely >=2.0.5",
"pandera[mypy] >= 0.20.3",
"pydantic >=2.8.2",
Expand Down
45 changes: 40 additions & 5 deletions src/nsidc/iceflow/itrf/converter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,10 @@ def transform_itrf(
# and the mean of that chunk is used to determine the plate name.
plate: str | None = None,
) -> IceflowDataFrame:
"""Pipeline string for proj to transform from the source to the target
ITRF frame and, optionally, epoch.
"""Transform the data's lon/lat/elev from the source ITRF to the target ITRF.

If a `target_epoch` is given, coordinate propagation is performed via a
plate motion model.
"""
if not check_itrf(target_itrf):
err_msg = (
Expand All @@ -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:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

transformed_chunks.append(chunk)
continue

Expand All @@ -69,17 +71,50 @@ def transform_itrf(
if not plate:
plate = plate_name(Point(chunk.longitude.mean(), chunk.latitude.mean()))
plate_model_step = (
f"+step +inv +init={target_itrf}:{plate} +t_epoch={target_epoch} "
# Perform coordinate propagation to the given epoch using the
# 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} "
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

)

itrf_transformation_step = ""
if source_itrf != target_itrf:
# This performs a helmert transform (see
# https://proj.org/en/9.4/operations/transformations/helmert.html). `+init=ITRF2014:ITRF2008`
# looks up the ITRF2008 helmert transformation step in the ITRF2014
# data file (see e.g.,
# 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} "
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

@trey-stafford trey-stafford Oct 22, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

WIP implementation that tests each way (with and without +inv) for validity: #46


pipeline = (
# This initializes the pipeline and declares the use of the WGS84
# ellipsoid for all of the following steps. See
# https://proj.org/en/9.5/operations/pipeline.html.
f"+proj=pipeline +ellps=WGS84 "
# 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 "
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

# 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 "
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

# Convert from lat/lon/elev geodetic coordinates to cartesian
# coordinates, which are required for the following steps.
# See: https://proj.org/en/9.5/operations/conversions/cart.html
f"+step +proj=cart "
f"+step +inv +init={target_itrf}:{source_itrf} "
# ITRF transformation. See above for definition.
f"{itrf_transformation_step}"
# See above for definition.
f"{plate_model_step}"
# Convert back from cartesian to lat/lon coordinates
f"+step +inv +proj=cart "
# Convert lon/lat from radians back to degrees.
# TODO: remove this if the initial conversion to radians above is not needed
f"+step +proj=unitconvert +xy_in=rad +xy_out=deg"
)

Expand Down