Bugfixes for transform_itrf

notebooks/iceflow-example.ipynb

@@ -26,7 +26,7 @@
     "import datetime as dt\n",
     "\n",
     "from nsidc.iceflow.api import fetch_iceflow_df\n",
-    "from nsidc.iceflow.data.models import DatasetSearchParameters, BoundingBox, ATM1BDataset\n",
+    "from nsidc.iceflow.data.models import DatasetSearchParameters, BoundingBox, ILATM1BDataset\n",
     "\n",
     "# Downloaded data will go here.\n",
     "data_path = Path(\"./downloaded-data/\")\n",
@@ -35,8 +35,8 @@
     "# ILATM1B data can be very large, so for the purposes of this example we will focus on just a small area with a manageable amount of data.\n",
     "BBOX = BoundingBox(lower_left_lon=-103.125559, lower_left_lat=-75.180563, upper_right_lon=-102.677327, upper_right_lat=-74.798063)\n",
     "\n",
-    "# Define the dataset that we want to search for. ATM1B version 1 is the only version available at the time of writing.\n",
-    "atm1b_v1_dataset = ATM1BDataset(version=\"1\")\n",
+    "# Define the dataset that we want to search for.\n",
+    "atm1b_v1_dataset = ILATM1BDataset(version=\"1\")\n",
     "\n",
     "# We will define a short date range in 2009 to search for data. Again, we choose this primarily to keep the amount of data manageable.\n",
     "date_range = (dt.date(2009, 11, 1), dt.date(2009, 12, 31))\n",

pyproject.toml

@@ -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",

src/nsidc/iceflow/itrf/converter.py

@@ -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 = (
@@ -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:
             transformed_chunks.append(chunk)
             continue
 
@@ -69,17 +71,51 @@ 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} "
             )
 
+        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
+            # https://proj.org/en/9.3/resource_files.html#init-files and 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} "
+
         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 "
-            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"
         )