Merge pull request #74 from nsidc/41-support-other-ilvis2-lat-lon-elev-itrf-transforms

41 support other ilvis2 lat lon elev itrf transforms

Author: trey-stafford

CHANGELOG.md

@@ -1,3 +1,17 @@
+# v1.1.0
+
+- Support selecting alternative lat/lon/elev triplet as primary lat/lon/elev
+  fields for ILVIS2 data. By default, the `low_mode` coordinates are used, which
+  represent the center of the lowest detected mode within the waveform. This
+  replicates the behavior of the valkyrie service this code is based on. Users
+  may now choose between the following coordinate sets: `low_mode` (the
+  default), `high_mode`, `centroid` (ILVIS2 v1 only), and `highest_signal`
+  (ILVIS2 v2 only).
+- Update documentation around ILVIS2 datasets and their multiple tuplets of
+  lat/lon/elev fields.
+- Update data model for `IceflowDataFrame` to include `dataset`, which gives the
+  dataset short name and version as a string (e.g., ILVIS v2 is "ILVISv2").
+
 # v1.0.0
 
 - Update `transform_itrf` function to be more flexible. Both forward and reverse

docs/altimetry-data-overview.md

@@ -91,7 +91,7 @@ pulses are used to determine elevation data.
 ```{note}
 
 We recommend using the [_icepyx_](https://github.com/icesat2py/icepyx)
-library to access and interact with ICESat-2 data. Learn more about using `icepyx` with `icelfow` in the [Using iceflow with icepyx to Generate an Elevation Timeseries](notebooks/iceflow-with-icepyx) Jupyter notebook.
+library to access and interact with ICESat-2 data. Learn more about using `icepyx` with `iceflow` in the [Using iceflow with icepyx to Generate an Elevation Timeseries](notebooks/iceflow-with-icepyx) Jupyter notebook.
 
 ```
 
@@ -117,7 +117,8 @@ further mission information and documentation for each data set:
 | [BLATM L1B](https://nsidc.org/data/BLATM1B)               | South: N:-53, S: -90, E:180, W:-180 <br> North: N:90, S: 60, E:180, W:-180 | 23 Jun. 1993 - 30 Oct. 2008                    | Pre-IceBridge | ATM                                                  |
 | [ILATM L1B V1](https://nsidc.org/data/ILATM1B/versions/1) | South: N:-53, S: -90, E:180, W:-180 <br> North: N:90, S: 60, E:180, W:-180 | 31 Mar. 2009 - 8 Nov. 2012 <br> (updated 2013) | IceBridge     | ATM                                                  |
 | [ILATM L1B V2](https://nsidc.org/data/ILATM1B/versions/2) | South: N:-53, S: -90, E:180, W:-180 <br> North: N:90, S: 60, E:180, W:-180 | 20 Mar. 2013 - 16 May 2019 <br> (updated 2020) | IceBridge     | ATM                                                  |
-| [ILVIS2](https://nsidc.org/data/ILVIS2)                   | North: N:90, S: 60, E:180, W:-180                                          | 25 Aug. 2017 - 20 Sept. 2017                   | IceBridge     | ALTIMETERS, LASERS, LVIS                             |
+| [ILVIS2 v1](https://nsidc.org/data/ilvis2/versions/1)     | South: N:-53, S: -90, E:180, W:-180 <br> North: N:90, S: 60, E:180, W:-180 | 14 Apr. 2009 - 31 Oct. 2015                    | IceBridge     | ALTIMETERS, LASERS, LVIS                             |
+| [ILVIS2 v2](https://nsidc.org/data/ilvis2/versions/2)     | North: N:90, S: 60, E:180, W:-180                                          | 25 Aug. 2017 - 20 Sept. 2017                   | IceBridge     | ALTIMETERS, LASERS, LVIS                             |
 | [GLAH06](https://nsidc.org/data/GLAH06/)                  | Global: N:86, S: -86, E:180, W:-180                                        | 20 Feb. 2003 - 11 Oct. 2009                    | ICESat/GLAS   | ALTIMETERS, CD, GLAS, GPS, <br> GPS Receiver, LA, PC |
 
 ---
@@ -129,6 +130,36 @@ guides or contact NSIDC user services at [email protected]
 
 ```
 
+### ILVIS2 data
+
+ILVIS2 contain multiple sets of latitude/longitude/elevation values.
+
+- `GLAT`/`GLON`/`GZ` represent the center of the lowest mode in the waveform.
+- `HLAT`/`HLON`/`HZ` represent the center of the highest detected mode within
+  the waveform. Both of these sets of lat/lon/elev are available across v1 and
+  v2 ILIVS data.
+
+ILVIS V1 data:
+
+- `CLAT`/`CLON`/`ZC` represent the centroid of the corresponding LVIS Level-1B
+  waveform.
+
+ILVIS V2 data:
+
+- `TLAT`/`TLON`/`ZT`, which represent the highest detected signal.
+
+By default, `iceflow` will use `GLAT`/`GLON`/`GZ` as the primary
+latitude/longitude/elevation fields in `IceflowDataFrame`s. Use the
+`ilvis2_coordinate_set` kwarg on `read_iceflow_datafile(s)` or
+`make_iceflow_parquet` to select an different primary set of
+latitude/longitude/elevation fields. Alternatively, manually set the fields:
+
+```
+# 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"]]
+```
+
 ## Challenges
 
 The wealth of data from these missions presents an opportunity to study the
@@ -185,6 +216,9 @@ ICESat-2.
 - [OpenAltimetry](https://openaltimetry.earthdatacloud.nasa.gov/data/): Advanced
   discovery, processing, and visualization services for ICESat and ICESat-2
   altimeter data
+- [icepyx](https://icepyx.readthedocs.io/en/latest/): icepyx is both a software
+  library and a community composed of ICESat-2 data users, developers, and the
+  scientific community.
 - [ITS_LIVE](https://its-live.jpl.nasa.gov/): A NASA MEaSUREs project to provide
   automated, low latency, global glacier flow and elevation change data sets.
 

docs/getting-started.md

@@ -101,14 +101,15 @@ from nsidc.iceflow import make_iceflow_parquet
 parquet_path = make_iceflow_parquet(
     data_dir=Path("/path/to/data/dir/"),
     target_itrf="ITRF2014",
+    ilvis2_coordinate_set="low_mode",
 )
 df = dd.read_parquet(parquet_path)
 ```
 
 Note that `make_iceflow_parquet` creates a parquet datastore for the data in the
 provided `data_dir` with the data transformed into a common
-[ITRF](https://itrf.ign.fr/) to facilitate analysis. Only datetime, lat, lon,
-and elevation fields are preserved in the parquet datastore.
+[ITRF](https://itrf.ign.fr/) to facilitate analysis. Only datetime, latitude,
+longitude, elevation, and dataset fields are preserved in the parquet datastore.
 
 To access and analyze the full data record in the source files, use
 [`read_iceflow_datafiles`](nsidc.iceflow.read_iceflow_datafiles):
@@ -117,7 +118,10 @@ To access and analyze the full data record in the source files, use
 from nsidc.iceflow import read_iceflow_datafiles
 
 # Read all of the data in the source files - not just lat/lon/elev.
-df = read_iceflow_datafiles(downloaded_files)
+df = read_iceflow_datafiles(
+    downloaded_files,
+    ilvis2_coordinate_set="low_mode",
+)
 
 # Optional: transform lat/lon/elev to common ITRF:
 from nsidc.iceflow import transform_itrf
@@ -130,3 +134,12 @@ df = transform_itrf(
 Note that `read_iceflow_datafiles` reads all of the data from the given
 filepaths. This could be a large amount of data, and could cause your program to
 crash if physical memory limits are exceeded.
+
+#### Special considerations for ILVIS2 data
+
+Users of ILVIS2 data should be aware that ILVIS2 data contains multiple sets of
+lat/lon/elev that may be of interest. By default, the `low_mode` set is used as
+the primary set of latitude/longitude/elevation used by `iceflow`.
+
+See [ILVIS2 data](./altimetry-data-overview.md#ilvis2-data) for more
+information.

src/nsidc/iceflow/api.py

@@ -8,6 +8,7 @@
 import dask.dataframe as dd
 from loguru import logger
 
+from nsidc.iceflow.data.ilvis2 import ILVIS2_DEFAULT_COORDINATE_SET
 from nsidc.iceflow.data.read import read_iceflow_datafiles
 from nsidc.iceflow.data.supported_datasets import ALL_SUPPORTED_DATASETS
 from nsidc.iceflow.itrf.converter import transform_itrf
@@ -19,6 +20,7 @@ def make_iceflow_parquet(
     target_itrf: str,
     overwrite: bool = False,
     target_epoch: str | None = None,
+    ilvis2_coordinate_set=ILVIS2_DEFAULT_COORDINATE_SET,
 ) -> Path:
     """Create a parquet dataset containing the lat/lon/elev data in `data_dir`.
 
@@ -52,7 +54,10 @@ def make_iceflow_parquet(
     ]
     for subdir in all_subdirs:
         iceflow_filepaths = [path for path in subdir.iterdir() if path.is_file()]
-        iceflow_df = read_iceflow_datafiles(iceflow_filepaths)
+        iceflow_df = read_iceflow_datafiles(
+            iceflow_filepaths,
+            ilvis2_coordinate_set=ilvis2_coordinate_set,
+        )
 
         iceflow_df = transform_itrf(
             data=iceflow_df,
@@ -61,8 +66,6 @@ def make_iceflow_parquet(
         )
 
         # Add a string col w/ dataset name and version.
-        short_name, version = subdir.name.split("_")
-        iceflow_df["dataset"] = [f"{short_name}v{version}"] * len(iceflow_df.latitude)
         common_columns = ["latitude", "longitude", "elevation", "dataset"]
         common_dask_df = dd.from_pandas(iceflow_df[common_columns])
         if parquet_subdir.exists():

src/nsidc/iceflow/data/atm1b.py

@@ -440,15 +440,19 @@ def atm1b_data(filepath: Path) -> ATM1BDataFrame:
     if year >= 2013:
         file_date = _ilatm1b_date(filename)
         data = _ilatm1bv2_data(filepath, file_date)
+        dataset = "ILATM1Bv2"
     elif year >= 2009:
         file_date = _ilatm1b_date(filename)
         data = _atm1b_qfit_data(filepath, file_date)
+        dataset = "ILATM1Bv1"
     else:
         file_date = _blatm1bv1_date(filename)
         data = _atm1b_qfit_data(filepath, file_date)
+        dataset = "BLATM1Bv1"
 
     itrf = extract_itrf(filepath)
     data["ITRF"] = itrf
+    data["dataset"] = dataset
 
     data = data.set_index("utc_datetime")
 

src/nsidc/iceflow/data/glah06.py

@@ -253,6 +253,7 @@ def glah06_data(filepath: Path) -> GLAH06DataFrame:
     df["latitude"] = df["d_lat"]
     df["longitude"] = df["d_lon"]
     df["elevation"] = df["d_elev"]
+    df["dataset"] = "GLAH06v034"
 
     # We index the data by utc datetime.
     df = df.set_index("utc_datetime")

src/nsidc/iceflow/data/ilvis2.py

@@ -4,13 +4,47 @@
 import re
 from collections import namedtuple
 from pathlib import Path
+from typing import Literal
 
 import numpy as np
 import pandas as pd
 import pandera as pa
 
 from nsidc.iceflow.data.models import ILVIS2DataFrame
 
+ILVIS2_COORDINATE_SETS = Literal["low_mode", "high_mode", "centroid", "highest_signal"]
+ILVIS2_COORDINATE_SET_MAPPING: dict[ILVIS2_COORDINATE_SETS, dict[str, str]] = {
+    # The center of the lowest detected mode within the waveform. Available in
+    # both v1 and v2.
+    "low_mode": {
+        "latitude": "GLAT",
+        "longitude": "GLON",
+        "elevation": "ZG",
+    },
+    # The center of the highest detected mode within the waveform. Available in
+    # both v1 and v2.
+    "high_mode": {
+        "latitude": "HLAT",
+        "longitude": "HLON",
+        "elevation": "ZH",
+    },
+    # Centroid of the corresponding LVIS Level-1B waveform. v1 only.
+    "centroid": {
+        "latitude": "CLAT",
+        "longitude": "CLON",
+        "elevation": "ZC",
+    },
+    # Highest detected signal. v2 only.
+    "highest_signal": {
+        "latitude": "TLAT",
+        "longitude": "TLON",
+        "elevation": "ZT",
+    },
+}
+# In the valkyrie service that this code is based on, only the low_mode was
+# exposed to users. This default replicates that behavior.
+ILVIS2_DEFAULT_COORDINATE_SET: ILVIS2_COORDINATE_SETS = "low_mode"
+
 Field = namedtuple("Field", ["name", "type", "scale_factor"])
 
 """
@@ -160,7 +194,10 @@ def _ilvis2_data(filepath: Path, file_date: dt.date, fields) -> pd.DataFrame:
 
 
 @pa.check_types()
-def ilvis2_data(filepath: Path) -> ILVIS2DataFrame:
+def ilvis2_data(
+    filepath: Path,
+    coordinate_set: ILVIS2_COORDINATE_SETS = ILVIS2_DEFAULT_COORDINATE_SET,
+) -> ILVIS2DataFrame:
     """Return the ilvis2 data given a filepath.
 
     Parameters
@@ -187,34 +224,39 @@ def ilvis2_data(filepath: Path) -> ILVIS2DataFrame:
 
     if year < 2017:
         # This corresponds to ILVIS v1
+        dataset = "ILVIS2v1"
         the_fields = ILVIS2_V104_FIELDS
         # The user guide indicates ILVIS2 v1 data uses ITRF2000 as a reference frame:
         # https://nsidc.org/sites/default/files/documents/user-guide/ilvis2-v001-userguide.pdf
         itrf_str = "ITRF2000"
+        # Ensure that the highest_signal coordinate set has not been selected.
+        if coordinate_set == "highest_signal":
+            raise ValueError(
+                "ILVIS coordinate set 'highest_signal' is only available in v2 data."
+            )
     else:
         # This corresponds to ILVIS v2
+        dataset = "ILVIS2v2"
         the_fields = ILVIS2_V202b_FIELDS
         # The user guide indicates ILVIS2 v1 data uses ITRF2008 as a reference frame:
         # https://nsidc.org/sites/default/files/documents/user-guide/ilvis2-v002-userguide.pdf
         itrf_str = "ITRF2008"
+        # Ensure that the centroid coordinate set has not been selected.
+        if coordinate_set == "centroid":
+            raise ValueError(
+                "ILVIS coordinate set 'centroid' is only available in v1 data."
+            )
 
     file_date = _file_date(filename)
 
     data = _ilvis2_data(filepath, file_date, the_fields)
     data["ITRF"] = itrf_str
+    data["dataset"] = dataset
 
-    # TODO: this data does not have a single set of latitude, longitude, and
-    # elevation fields. Instead, it has e.g., "CLON" and "GLON" and "HLON". In
-    # the original `valkyrie` service code, it looks like "GLON", "GLAT", and
-    # "ZG" cols were used as for the points stored in the valkyrie database and
-    # transformed by the ITRF transformation service. Ideally, we support
-    # consistent transformation of the ITRF across all lat/lon/elev
-    # fields. E.g., a user may be more interested in looking at the "CLON",
-    # "CLAT", and "ZC" fields instead.
-    # For now, we will replicate the behavior of `valkyrie`:
-    data["latitude"] = data["GLAT"]
-    data["longitude"] = data["GLON"]
-    data["elevation"] = data["ZG"]
+    coordinate_fields = ILVIS2_COORDINATE_SET_MAPPING[coordinate_set]
+    data["latitude"] = data[coordinate_fields["latitude"]]
+    data["longitude"] = data[coordinate_fields["longitude"]]
+    data["elevation"] = data[coordinate_fields["elevation"]]
 
     data = data.set_index("utc_datetime")
 

src/nsidc/iceflow/data/models.py

@@ -17,6 +17,10 @@ class CommonDataColumnsSchema(pa.DataFrameModel):
     latitude: Series[float] = pa.Field(coerce=True)
     longitude: Series[float] = pa.Field(coerce=True)
     elevation: Series[float] = pa.Field(coerce=True)
+    # In practice, "dataset" is always available to provide provenance on where
+    # each point comes from, but we mark it as optional here because it is not
+    # strictly necessary for e.g., ITRF transformations.
+    dataset: str | None
 
 
 class ATM1BSchema(CommonDataColumnsSchema):
@@ -38,22 +42,41 @@ class ATM1BSchema(CommonDataColumnsSchema):
     pulse_width: Series[float] = pa.Field(nullable=True, coerce=True)
 
 
-# Note/TODO: the ILVIS2 data contain multiple sets of lat/lon/elev. The common
+# Note: the ILVIS2 data contain multiple sets of lat/lon/elev. The common
 # schema assumes one set of lat/lon/elev which is used for the ITRF
 # transformation code.
 class ILVIS2Schema(CommonDataColumnsSchema):
+    """ILVIS2 Data Schema.
+
+    Note that ILVIS2 data contain multiple sets of lat/lon/elev.
+
+    * GLAT/GLON/GZ represent the center of the lowest mode in the waveform.
+    * HLAT/HLON/HZ represent the center of the highest detected mode within the
+      waveform. Both of these sets of lat/lon/elev are available across v1 and
+      v2 ILVIS2 data.
+
+    ILVIS V1 data:
+    * CLAT/CLON/ZC represent the centroid of the corresponding LVIS Level-1B waveform.
+
+    ILVIS V2 data:
+    * TLAT/TLON/ZT, which represent the highest detected signal
+    """
+
     # Common columns
     LFID: Series[float] = pa.Field(nullable=True, coerce=True)
     SHOTNUMBER: Series[float] = pa.Field(nullable=True, coerce=True)
     TIME: Series[float] = pa.Field(nullable=True, coerce=True)
+    # ZG/GLAT/GLON: the center of the lowest detected mode within the waveform
     ZG: Series[float] = pa.Field(nullable=True, coerce=True)
     GLAT: Series[float] = pa.Field(nullable=True, coerce=True)
     GLON: Series[float] = pa.Field(nullable=True, coerce=True)
+    # HLAT/HLON/ZH: the center of the highest detected mode within the waveform
     HLAT: Series[float] = pa.Field(nullable=True, coerce=True)
     HLON: Series[float] = pa.Field(nullable=True, coerce=True)
     ZH: Series[float] = pa.Field(nullable=True, coerce=True)
 
     # V104-specific
+    # CLAT/CLON/ZC: Centroid of the corresponding LVIS Level-1B waveform
     CLAT: Series[float] = pa.Field(nullable=True, coerce=True)
     CLON: Series[float] = pa.Field(nullable=True, coerce=True)
     ZC: Series[float] = pa.Field(nullable=True, coerce=True)
@@ -66,6 +89,7 @@ class ILVIS2Schema(CommonDataColumnsSchema):
     COMPLEXITY: Series[float] = pa.Field(nullable=True, coerce=True)
     INCIDENT_ANGLE: Series[float] = pa.Field(nullable=True, coerce=True)
     RANGE: Series[float] = pa.Field(nullable=True, coerce=True)
+    # RH%%%: Height (relative to ZG) at which % of the waveform energy occurs
     RH10: Series[float] = pa.Field(nullable=True, coerce=True)
     RH15: Series[float] = pa.Field(nullable=True, coerce=True)
     RH20: Series[float] = pa.Field(nullable=True, coerce=True)
@@ -89,6 +113,7 @@ class ILVIS2Schema(CommonDataColumnsSchema):
     RH98: Series[float] = pa.Field(nullable=True, coerce=True)
     RH99: Series[float] = pa.Field(nullable=True, coerce=True)
     RH100: Series[float] = pa.Field(nullable=True, coerce=True)
+    # Highest detected signal
     TLAT: Series[float] = pa.Field(nullable=True, coerce=True)
     TLON: Series[float] = pa.Field(nullable=True, coerce=True)
     ZT: Series[float] = pa.Field(nullable=True, coerce=True)