Reference corrections to point, convert to meters (#377)

* pass through the `ReferencePoint` found from timeseries

* convert `troposphere` to output delay in meters

* convert ionosphere to meters

* convert troposphere to meters

* add subtraction of reference_point to `ionosphere` and `troposphere`

* add functions for units and description

* add `get_` functions, export to `__all__`

* add units to `write_arr` for iono/tropo

* remove chatting wgs84 warning, add tropo logging

* use `timeseries_paths` for corrections, not ifg paths

* fix new test

src/dolphin/atmosphere/_netcdf.py

@@ -75,8 +75,8 @@ def delay_from_netcdf(
         try:
             wm_proj = CRS.from_wkt(ds["proj"].attrs["crs_wkt"])
         except KeyError:
-            logger.warning(
-                "WARNING: I can't find a CRS in the weather model file, "
+            logger.debug(
+                "I can't find a CRS in the weather model file, "
                 "so I will assume you are using WGS84"
             )
             wm_proj = CRS.from_epsg(4326)
@@ -265,7 +265,6 @@ def _build_cube(xpts, ypts, zpts, model_crs, pts_crs, interpolators):
 
     # Loop over heights and compute delays
     for ii, ht in enumerate(zpts):
-
         # pts is in weather model system;
         if model_crs != pts_crs:
             # lat / lon / height for hrrr

src/dolphin/atmosphere/ionosphere.py

@@ -16,6 +16,7 @@
 
 from dolphin import io
 from dolphin._types import Bbox, Filename
+from dolphin.timeseries import ReferencePoint
 from dolphin.utils import format_date_pair
 
 logger = logging.getLogger(__name__)
@@ -35,11 +36,12 @@ def estimate_ionospheric_delay(
     slc_files: Mapping[tuple[datetime.datetime], Sequence[Filename]],
     tec_files: Mapping[tuple[datetime.datetime], Sequence[Filename]],
     geom_files: dict[str, Path],
+    reference_point: ReferencePoint | None,
     output_dir: Path,
     epsg: int,
     bounds: Bbox,
 ) -> list[Path]:
-    """Estimate the range delay caused by ionosphere for each interferogram.
+    """Estimate the range delay (in meters) caused by ionosphere for each interferogram.
 
     Parameters
     ----------
@@ -52,6 +54,9 @@ def estimate_ionospheric_delay(
     geom_files : list[Path]
         Dictionary of geometry files with height and incidence angle, or
         los_east and los_north.
+    reference_point : tuple[int, int], optional
+        Reference point (row, col) used during time series inversion.
+        If not provided, no spatial referencing is performed.
     output_dir : Path
         Output directory.
     epsg : int
@@ -166,13 +171,20 @@ def estimate_ionospheric_delay(
             secondary_vtec, iono_inc_angle, freq, obs_type="phase"
         )
 
-        ifg_iono_range_delay = range_delay_reference - range_delay_secondary
+        ifg_iono_range_delay_radians = range_delay_reference - range_delay_secondary
+        # Convert to meters, where positive corresponds to motion toward the satellite
+        ifg_iono_range_delay = -wavelength / (4 * np.pi) * ifg_iono_range_delay_radians
+
+        if reference_point is not None:
+            ref_row, ref_col = reference_point
+            ifg_iono_range_delay -= ifg_iono_range_delay[ref_row, ref_col]
 
         # Write 2D tropospheric correction layer to disc
         io.write_arr(
             arr=ifg_iono_range_delay,
             output_name=iono_delay_product_name,
             like_filename=ifg,
+            units="meters",
         )
 
     return output_paths

src/dolphin/atmosphere/troposphere.py

@@ -17,6 +17,7 @@
 
 from dolphin import io
 from dolphin._types import Bbox, Filename, TropoModel, TropoType
+from dolphin.timeseries import ReferencePoint
 from dolphin.utils import format_date_pair
 
 from ._netcdf import delay_from_netcdf, group_netcdf_by_date
@@ -53,9 +54,6 @@ class DelayParams:
     geotransform: list[float]
     """Sequence of geotransformation parameters."""
 
-    wavelength: float
-    """Radar wavelength."""
-
     tropo_model: TropoModel
     """Model used for tropospheric correction."""
 
@@ -83,14 +81,15 @@ def estimate_tropospheric_delay(
     slc_files: Mapping[tuple[datetime.datetime], Sequence[Filename]],
     troposphere_files: Sequence[Filename],
     geom_files: dict[str, Path],
+    reference_point: ReferencePoint | None,
     output_dir: Path,
     file_date_fmt: str,
     tropo_model: TropoModel,
     tropo_delay_type: TropoType,
     epsg: int,
     bounds: Bbox,
 ) -> list[Path]:
-    """Estimate the tropospheric delay corrections for each interferogram.
+    """Estimate the tropospheric delay corrections (in meters) for each interferogram.
 
     Parameters
     ----------
@@ -103,6 +102,9 @@ def estimate_tropospheric_delay(
     geom_files : dict[str, Path]
         Dictionary of geometry files with height and incidence angle, or
         los_east and los_north.
+    reference_point : tuple[int, int], optional
+        Reference point (row, col) used during time series inversion.
+        If not provided, no spatial referencing is performed.
     output_dir : Path
         Output directory.
     file_date_fmt : str
@@ -121,6 +123,7 @@ def estimate_tropospheric_delay(
     -------
     list[Path]
         List of newly created tropospheric phase delay geotiffs.
+        Units are in meters.
 
     """
     # Read geogrid data
@@ -137,13 +140,6 @@ def estimate_tropospheric_delay(
 
     tropo_height_levels = np.concatenate(([-100], np.arange(0, 9000, 500)))
 
-    for key in slc_files:
-        if "compressed" not in str(slc_files[key][0]).lower():
-            one_of_slcs = slc_files[key][0]
-            break
-
-    wavelength = oput.get_radar_wavelength(one_of_slcs)
-
     delay_type = tropo_delay_type.value
 
     if str(troposphere_files[0]).endswith(".nc"):
@@ -216,7 +212,6 @@ def estimate_tropospheric_delay(
             epsg=epsg,
             tropo_model=tropo_model,
             delay_type=delay_type,
-            wavelength=wavelength,
             shape=(ysize, xsize),
             geotransform=gt,
             reference_file=tropo_files[closest_date_ref][0],
@@ -226,15 +221,21 @@ def estimate_tropospheric_delay(
             interferogram=format_date_pair(ref_date, sec_date),
         )
 
+        logger.info(f"Running troposphere computation for {ref_date}, {sec_date}")
         delay_datacube = tropo_run(delay_parameters)
 
         tropo_delay_2d = compute_2d_delay(delay_parameters, delay_datacube, geom_files)
 
+        if reference_point is not None:
+            ref_row, ref_col = reference_point
+            tropo_delay_2d -= tropo_delay_2d[ref_row, ref_col]
+
         # Write 2D tropospheric correction layer to disc
         io.write_arr(
             arr=tropo_delay_2d,
             output_name=tropo_delay_product_path,
             like_filename=ifg,
+            units="meters",
         )
 
     return output_paths
@@ -304,9 +305,8 @@ def compute_pyaps(delay_parameters: DelayParams) -> np.ndarray:
         phs_second = second_aps_estimator.getdelay()
 
         # Convert the delay in meters to radians
-        tropo_delay_datacube_list.append(
-            -(phs_ref - phs_second) * 4.0 * np.pi / delay_parameters.wavelength
-        )
+        relative_delay = phs_second - phs_ref
+        tropo_delay_datacube_list.append(relative_delay)
 
     # Tropo delay datacube
     tropo_delay_datacube = np.stack(tropo_delay_datacube_list)
@@ -315,7 +315,7 @@ def compute_pyaps(delay_parameters: DelayParams) -> np.ndarray:
 
 
 def compute_tropo_delay_from_netcdf(delay_parameters: DelayParams) -> np.ndarray:
-    """Compute tropospheric delay datacube from netcdf tropo file.
+    """Compute tropospheric delay (in meters) datacube from netcdf tropo file.
 
     Parameters
     ----------
@@ -364,8 +364,8 @@ def compute_tropo_delay_from_netcdf(delay_parameters: DelayParams) -> np.ndarray
             - tropo_delay_secondary[delay_parameters.delay_type]
         )
 
-    # Convert it to radians units
-    tropo_delay_datacube = -tropo_delay * 4.0 * np.pi / delay_parameters.wavelength
+    # Convert it to convention where positive means toward the satellite
+    tropo_delay_datacube = -1 * tropo_delay
 
     # Create a masked datacube that excludes the NaN values
     tropo_delay_datacube_masked = np.ma.masked_invalid(tropo_delay_datacube)
@@ -413,7 +413,7 @@ def compute_2d_delay(
     Returns
     -------
     np.ndarray
-        Computed 2D delay.
+        Computed 2D delay in meters.
 
     """
     dem_file = Path(geo_files["height"])

src/dolphin/io/_core.py

@@ -38,6 +38,10 @@
     "get_raster_bounds",
     "get_raster_dtype",
     "get_raster_metadata",
+    "get_raster_units",
+    "set_raster_units",
+    "get_raster_description",
+    "set_raster_description",
     "get_raster_gt",
     "get_raster_driver",
     "get_raster_chunk_size",
@@ -406,6 +410,81 @@ def set_raster_metadata(
     ds = None
 
 
+def get_raster_description(filename: Filename, band: int = 1):
+    """Get description of a raster band.
+
+    Parameters
+    ----------
+    filename : Filename
+        Path to the file to load.
+    band : int, optional
+        Band to get description for. Default is 1.
+
+    """
+    ds = gdal.Open(fspath(filename), gdal.GA_Update)
+    bnd = ds.GetRasterBand(band)
+    return bnd.GetDescription()
+
+
+def set_raster_description(filename: Filename, description: str, band: int = 1):
+    """Set description on a raster band.
+
+    Parameters
+    ----------
+    filename : Filename
+        Path to the file to load.
+    description : str
+        Description to set.
+    band : int, optional
+        Band to set description for. Default is 1.
+
+    """
+    ds = gdal.Open(fspath(filename), gdal.GA_Update)
+    bnd = ds.GetRasterBand(band)
+    bnd.SetDescription(description)
+    bnd.FlushCache()
+    ds = None
+
+
+def get_raster_units(filename: Filename, band: int = 1) -> str | None:
+    """Get units of a raster band.
+
+    Parameters
+    ----------
+    filename : Filename
+        Path to the file to load.
+    band : int
+        Band to get units for.
+        Default is 1.
+
+    """
+    ds = gdal.Open(fspath(filename))
+    bnd = ds.GetRasterBand(band)
+    return bnd.GetUnitType() or None
+
+
+def set_raster_units(filename: Filename, units: str, band: int | None = None) -> None:
+    """Set units on a raster band.
+
+    Parameters
+    ----------
+    filename : Filename
+        Path to the file to load.
+    units : str
+        Units to set.
+    band : int, optional
+        Band to set units for. Default is None, which sets for all bands.
+
+    """
+    ds = gdal.Open(fspath(filename), gdal.GA_Update)
+    if band is None:
+        bands = range(1, ds.RasterCount + 1)
+    for i in bands:
+        bnd = ds.GetRasterBand(i)
+        bnd.SetUnitType(units)
+        bnd.FlushCache()
+
+
 def rowcol_to_xy(
     row: int,
     col: int,
@@ -467,6 +546,8 @@ def write_arr(
     strides: Optional[dict[str, int]] = None,
     projection: Optional[Any] = None,
     nodata: Optional[float] = None,
+    units: Optional[str] = None,
+    description: Optional[str] = None,
 ):
     """Save an array to `output_name`.
 
@@ -507,6 +588,11 @@ def write_arr(
     nodata : float, optional
         Nodata value to save.
         Default is the nodata of band 1 of `like_filename` (if provided), or None.
+    units : str, optional
+        Units of the data. Default is None.
+        Value is stored in the metadata as "units".
+    description : str, optional
+        Description of the raster bands stored in the metadata.
 
     """
     fi = FileInfo.from_user_inputs(
@@ -552,12 +638,17 @@ def write_arr(
             bnd = ds_out.GetRasterBand(i + 1)
             bnd.WriteArray(arr[i])
 
-    # Set the nodata value for each band
-    if fi.nodata is not None:
-        for i in range(fi.nbands):
-            logger.debug(f"Setting nodata for band {i+1}/{fi.nbands}")
-            bnd = ds_out.GetRasterBand(i + 1)
+    # Set the nodata/units/description for each band
+    for i in range(fi.nbands):
+        logger.debug(f"Setting nodata for band {i+1}/{fi.nbands}")
+        bnd = ds_out.GetRasterBand(i + 1)
+        # Note: right now we're assuming the nodata/units/description
+        if fi.nodata is not None:
             bnd.SetNoDataValue(fi.nodata)
+        if units is not None:
+            bnd.SetUnitType(units)
+        if description is not None:
+            bnd.SetDescription(description)
 
     ds_out.FlushCache()
     ds_out = None