From fc1f7008effa22674edd3f3ea7331467ada8c585 Mon Sep 17 00:00:00 2001 From: Trey Stafford Date: Mon, 30 Jun 2025 11:53:22 -0600 Subject: [PATCH 1/4] Extract PMM transform step into function --- src/nsidc/iceflow/itrf/converter.py | 71 +++++++++++++++++------------ 1 file changed, 43 insertions(+), 28 deletions(-) diff --git a/src/nsidc/iceflow/itrf/converter.py b/src/nsidc/iceflow/itrf/converter.py index 4220b23..bfdc229 100644 --- a/src/nsidc/iceflow/itrf/converter.py +++ b/src/nsidc/iceflow/itrf/converter.py @@ -95,6 +95,42 @@ def _itrf_transformation_step(source_itrf: str, target_itrf: str) -> str: raise RuntimeError(err_msg) +def _get_pmm_transform_step( + data: IceflowDataFrame, + target_itrf: str, + target_epoch: str | None, + plate: str | None, +): + """Perform coordinate propagation to the target epoch using a proj-provided plate + motion model (PMM). + + This step uses the target_itrf's init file to lookup the associated plate's + PMM parameters. For example, ITRF2014 parameters are defined here: + https://github.com/OSGeo/PROJ/blob/8b65d5b14e2a8fbb8198335019488a2b2968df5c/data/ITRF2014. + The step is inverted because proj defined `t_epoch` as the "central epoch" - + not the "target epoch". The transformation uses a delta defined by + `t_observed - t_epoch` that are applied to the PMM's rate of change to + propagate the point into the past/future. See + https://proj.org/en/9.5/operations/transformations/helmert.html#mathematical-description + for more information. For example, if a point's observation date + (`t_observed`) is 1993.0 (1993-01-01T00:00:00) and the t_epoch is 2011.0 + (2011-01-01T00:00:00), then the delta is 1993 - 2011: -18. We need to invert + the step so that the point is propagated forward in time, from 1993 to 2011. + """ + plate_model_step = "" + if target_epoch: + if not plate: + plate = plate_name(Point(data.longitude.mean(), data.latitude.mean())) + plate_model_step = ( + f"+step +inv +init={target_itrf}:{plate} +t_epoch={target_epoch}" + ) + if not _check_valid_proj_step(plate_model_step): + err_msg = f"Failed to find pre-defined plate-model parameters for {target_itrf}:{plate}" + raise RuntimeError(err_msg) + + return plate_model_step + + @pa.check_types() def transform_itrf( data: IceflowDataFrame, @@ -152,34 +188,13 @@ def transform_itrf( transformed_chunks.append(chunk) continue - plate_model_step = "" - if target_epoch: - if not plate: - plate = plate_name(Point(chunk.longitude.mean(), chunk.latitude.mean())) - plate_model_step = ( - # Perform coordinate propagation to the target epoch using the - # provided plate motion model (PMM). - # This step uses the target_itrf's init file to lookup the - # associated plate's PMM parameters. For example, ITRF2014 - # parameters are defined here: - # https://github.com/OSGeo/PROJ/blob/8b65d5b14e2a8fbb8198335019488a2b2968df5c/data/ITRF2014. - # The step is inverted because proj defined `t_epoch` as the - # "central epoch" - not the "target epoch. The transformation - # uses a delta defined by `t_observed - t_epoch` that are - # applied to the PMM's rate of change to propagate the point - # into the past/future. See - # https://proj.org/en/9.5/operations/transformations/helmert.html#mathematical-description - # for more information. - # For example, if a point's observation date - # (`t_observed`) is 1993.0 (1993-01-01T00:00:00) and the t_epoch - # is 2011.0 (2011-01-01T00:00:00), then the delta is 1993 - - # 2011: -18. We need to invert the step so that the point is - # propagated forward in time, from 1993 to 2011. - f"+step +inv +init={target_itrf}:{plate} +t_epoch={target_epoch}" - ) - if not _check_valid_proj_step(plate_model_step): - err_msg = f"Failed to find pre-defined plate-model parameters for {target_itrf}:{plate}" - raise RuntimeError(err_msg) + chunk = cast(IceflowDataFrame, chunk) + plate_model_step = _get_pmm_transform_step( + chunk, + target_itrf, + target_epoch, + plate, + ) itrf_transformation_step = _itrf_transformation_step(source_itrf, target_itrf) From 87bf5ef128525a3a90eb4bfa6ee9c8f9869d72f6 Mon Sep 17 00:00:00 2001 From: Trey Stafford Date: Mon, 30 Jun 2025 17:52:59 -0600 Subject: [PATCH 2/4] WIP support using ORB params when available --- src/nsidc/iceflow/itrf/converter.py | 187 ++++++++++++++++++++-------- 1 file changed, 132 insertions(+), 55 deletions(-) diff --git a/src/nsidc/iceflow/itrf/converter.py b/src/nsidc/iceflow/itrf/converter.py index bfdc229..fb0b0c1 100644 --- a/src/nsidc/iceflow/itrf/converter.py +++ b/src/nsidc/iceflow/itrf/converter.py @@ -100,6 +100,7 @@ def _get_pmm_transform_step( target_itrf: str, target_epoch: str | None, plate: str | None, + pmm_use_orb: bool, ): """Perform coordinate propagation to the target epoch using a proj-provided plate motion model (PMM). @@ -116,21 +117,100 @@ def _get_pmm_transform_step( (`t_observed`) is 1993.0 (1993-01-01T00:00:00) and the t_epoch is 2011.0 (2011-01-01T00:00:00), then the delta is 1993 - 2011: -18. We need to invert the step so that the point is propagated forward in time, from 1993 to 2011. + + If `pmm_use_orb` is `True`, then this function will use Origin Rate Bias (ORB) + parameters for the selected PMM, if available (e.g., ITRF2020 and ITRF2008 + init files have ORB parameters defined; ITRF2014 has no ORB parameter sets). """ plate_model_step = "" - if target_epoch: - if not plate: - plate = plate_name(Point(data.longitude.mean(), data.latitude.mean())) + if not target_epoch: + return plate_model_step + + if not plate: + plate = plate_name(Point(data.longitude.mean(), data.latitude.mean())) + + if pmm_use_orb: + plate_model_name = f"{plate}_T" plate_model_step = ( - f"+step +inv +init={target_itrf}:{plate} +t_epoch={target_epoch}" + f"+step +inv +init={target_itrf}:{plate_model_name} +t_epoch={target_epoch}" ) - if not _check_valid_proj_step(plate_model_step): - err_msg = f"Failed to find pre-defined plate-model parameters for {target_itrf}:{plate}" - raise RuntimeError(err_msg) + + # If the step is valid, then return the step as-is. Otherwise, fallback + # on the non-ORB step. + if _check_valid_proj_step(plate_model_step): + return plate_model_step + + plate_model_name = f"{plate}" + plate_model_step = ( + f"+step +inv +init={target_itrf}:{plate_model_name} +t_epoch={target_epoch}" + ) + + if not _check_valid_proj_step(plate_model_step): + err_msg = f"Failed to find pre-defined plate-model parameters for {target_itrf}:{plate_model_name}" + raise RuntimeError(err_msg) return plate_model_step +def _run_itrf_pipeline( + data: IceflowDataFrame, + itrf_transformation_step: str, + plate_model_step: str, + target_itrf: str, +): + 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 " + # 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 " + # 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 " + # 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" + ) + + transformer = pyproj.Transformer.from_pipeline(pipeline) + + decimalyears = ( + data.reset_index().utc_datetime.apply(_datetime_to_decimal_year).to_numpy() + ) + # TODO: Should we create a new decimalyears when doing an epoch + # propagation since PROJ doesn't do this? + + lons, lats, elevs, _ = transformer.transform( + data.longitude, + data.latitude, + data.elevation, + decimalyears, + ) + + transformed_data = data.copy() + transformed_data["latitude"] = lats + transformed_data["longitude"] = lons + transformed_data["elevation"] = elevs + transformed_data["ITRF"] = target_itrf + + return transformed_data + + @pa.check_types() def transform_itrf( data: IceflowDataFrame, @@ -140,6 +220,10 @@ def transform_itrf( # is given but the plate name is not, each source ITRF is grouped together # and the mean of that chunk is used to determine the plate name. plate: str | None = None, + # Optional toggle for using ORB parameters, when available. + # Defaults to `True` because this is recommended per the PMM authors. See + # https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2023GL106373. + pmm_use_orb: bool = True, ) -> IceflowDataFrame: """Transform the data's latitude/longitude/elevation variables from the source ITRF to the target ITRF. @@ -154,6 +238,13 @@ def transform_itrf( defined for the target_epoch, so it is likely to be most accurate for points observed near the ITRF's defined epoch. + If `pmm_use_orb` is `True`, then this function will use Origin Rate Bias (ORB) + parameters for the selected PMM, if available (e.g., ITRF2020 and ITRF2008 + init files have ORB parameters defined; ITRF2014 has no ORB parameter sets), + for the horizontal component of each point (`longitude` and `latitude`). The + vertical component (`elevation`) will be transformed without ORB parameters + applied. + All ITRF and PMM transformations are dependent on the user's `proj` installation's ITRF init files (see https://proj.org/en/9.3/resource_files.html#init-files). For example, @@ -171,6 +262,7 @@ def transform_itrf( # TLAT/TLON/TZ are only available in ILVIS2v2 data: sel_ilvis2v2 = data.dataset == "ILVIS2v2" data.loc[sel_ilvis2v2, ["latitude", "longitude", "elevation"]] = data.loc[sel_ilvis2v2, ["TLAT", "TLON", "ZT"]] + ``` """ if not check_itrf(target_itrf): @@ -188,65 +280,50 @@ def transform_itrf( transformed_chunks.append(chunk) continue + itrf_transformation_step = _itrf_transformation_step(source_itrf, target_itrf) + chunk = cast(IceflowDataFrame, chunk) + plate_model_step = _get_pmm_transform_step( chunk, target_itrf, target_epoch, plate, + pmm_use_orb=False, ) - itrf_transformation_step = _itrf_transformation_step(source_itrf, target_itrf) - - 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 " - # 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 " - # 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 " - # 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" + transformed_chunk = _run_itrf_pipeline( + chunk, + itrf_transformation_step, + plate_model_step, + target_itrf, ) - transformer = pyproj.Transformer.from_pipeline(pipeline) - - decimalyears = ( - chunk.reset_index().utc_datetime.apply(_datetime_to_decimal_year).to_numpy() - ) - # TODO: Should we create a new decimalyears when doing an epoch - # propagation since PROJ doesn't do this? - - lons, lats, elevs, _ = transformer.transform( - chunk.longitude, - chunk.latitude, - chunk.elevation, - decimalyears, - ) + if pmm_use_orb: + plate_model_step_with_orb = _get_pmm_transform_step( + chunk, + target_itrf, + target_epoch, + plate, + pmm_use_orb=True, + ) + + # The ORB params are available if the resulting step is different. + if plate_model_step_with_orb != plate_model_step: + orb_transformed_chunk = _run_itrf_pipeline( + chunk, + itrf_transformation_step, + plate_model_step_with_orb, + target_itrf, + ) + # Transform only the lat/lon elements, per + # Z. Altamimi et. al. 2023 (https://doi.org/10.1029/2023GL106373): + # > it is recommended to ignore the artifactual vertical component of the + # > predicted velocities resulting from the addition of the ORB, and to only + # > consider the horizontal components of the predicted velocities. + transformed_chunk["latitude"] = orb_transformed_chunk.latitude + transformed_chunk["longitude"] = orb_transformed_chunk.longitude - transformed_chunk = chunk.copy() - transformed_chunk["latitude"] = lats - transformed_chunk["longitude"] = lons - transformed_chunk["elevation"] = elevs - transformed_chunk["ITRF"] = target_itrf transformed_chunks.append(transformed_chunk) transformed_df = pd.concat(transformed_chunks) From e576b53af173e15545837e48825ca05508ec076a Mon Sep 17 00:00:00 2001 From: Trey Stafford Date: Tue, 1 Jul 2025 16:54:44 -0600 Subject: [PATCH 3/4] CHANGELOG for v1.2.0 --- CHANGELOG.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 39f58b4..3bc381e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,12 @@ +# v1.2.0 + +- Add support to `transform_itrf` for PMM models using Origin Rate Bias (ORB) + parameters. By default, `transform_itrf` will use ORB parameters for the + selected PMM, if available (e.g., ITRF2020 and ITRF2008 init files have ORB + parameters defined; ITRF2014 has no ORB parameter sets), for the horizontal + component of each point (`longitude` and `latitude`). The vertical component + (`elevation`) will be transformed without ORB parameters applied. + # v1.1.0 - Support selecting alternative lat/lon/elev triplet as primary lat/lon/elev From 124645075c0ce9aa45e7d7dcc2c90abf5268d3cc Mon Sep 17 00:00:00 2001 From: Trey Stafford Date: Tue, 1 Jul 2025 17:02:28 -0600 Subject: [PATCH 4/4] Bumpversion v1.1.0 -> v1.2.0 --- pyproject.toml | 4 ++-- src/nsidc/iceflow/__init__.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 3723baf..a9c6ce7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "nsidc-iceflow" -version = "v1.1.0" +version = "v1.2.0" authors = [ { name = "NSIDC", email = "nsidc@nsidc.org" }, ] @@ -191,7 +191,7 @@ messages_control.disable = [ ] [tool.bumpversion] -current_version = "1.1.0" +current_version = "1.2.0" commit = false tag = false diff --git a/src/nsidc/iceflow/__init__.py b/src/nsidc/iceflow/__init__.py index 09414bd..a5e070c 100644 --- a/src/nsidc/iceflow/__init__.py +++ b/src/nsidc/iceflow/__init__.py @@ -33,7 +33,7 @@ ) from nsidc.iceflow.itrf.converter import transform_itrf -__version__ = "v1.1.0" +__version__ = "v1.2.0" __all__ = [