Merge pull request #91 from ACCESS-Cloud-Based-InSAR/set

Solid Earth Tides

.gitignore

@@ -99,6 +99,7 @@ tests/out/*.ipynb
 !tests/localize-data.ipynb
 /isce2_topsapp.egg-info/
 !tests/prepare-for-delivery.ipynb
+!tests/solid_earth_tides.ipynb
 
 # Jsons in Test Directory
 tests/**/*.json

CHANGELOG.md

@@ -6,6 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.2.3]
+
+### Added
+* Added support to compute and embed solid earth tide correction layers into GUNW products (see PR #91)
+
 ## [0.2.2]
 
 ### Added

README.md

@@ -64,6 +64,17 @@ isce2_topsapp --reference-scenes S1B_IW_SLC__1SDV_20171117T145926_20171117T14595
 ```
 Not including `--esd-coherence-threshold` means no ESD correction will be applied. The ESD threshold refers to a coherence value and therefore must be in $[0, 1]$.
 
+### Apply Solid Earth Tide Correction
+
+```
+isce2_topsapp --reference-scenes S1B_IW_SLC__1SDV_20210723T014947_20210723T015014_027915_0354B4_B3A9 \
+              --secondary-scenes S1B_IW_SLC__1SDV_20210711T014922_20210711T014949_027740_034F80_859D \
+                                 S1B_IW_SLC__1SDV_20210711T014947_20210711T015013_027740_034F80_D404 \
+                                 S1B_IW_SLC__1SDV_20210711T015011_20210711T015038_027740_034F80_376C \
+              --compute-solid-earth-tide True \
+              > topsapp_img_set.out 2> topsapp_img_set.err
+```
+
 ### Using "fixed frames" (experimental)
 
 Sentinel-1 Frames are not constant over passes. We generate fixed frames [here](https://github.com/ACCESS-Cloud-Based-InSAR/s1-frame-generation) and enumerate interferograms using this [repo](https://github.com/ACCESS-Cloud-Based-InSAR/s1-frame-enumerator). This is highly experimental. We then ensure ISCE processes only over the frame. The key is overlap. We provide some examples of the additional options (you will need to run this in *two* separate directories because ISCE2 outputs are organized with respect to the working directory of the processing). For one frame over CA:

environment.yml

@@ -29,6 +29,7 @@ dependencies:
  - notebook
  - numpy<1.24
  - pandas
+ - pysolid
  - papermill
  - pytest
  - pytest-cov
@@ -39,6 +40,7 @@ dependencies:
  - scipy<1.10
  - setuptools
  - setuptools_scm
+ - scipy<1.10
  - shapely
  - tqdm
  - dem_stitcher>=2.3.1

isce2_topsapp/__main__.py

@@ -8,13 +8,15 @@
 from pathlib import Path
 from typing import Optional
 
+import h5py
 
 from isce2_topsapp import (BurstParams, aws, download_aux_cal, download_bursts,
                            download_dem_for_isce2, download_orbits,
                            download_slcs, get_asf_slc_objects, get_region_of_interest,
                            package_gunw_product, prepare_for_delivery,
                            topsapp_processing)
 from isce2_topsapp.json_encoder import MetadataEncoder
+from isce2_topsapp.solid_earth_tides import update_gunw_with_solid_earth_tide
 
 
 def localize_data(reference_scenes: list,
@@ -114,6 +116,7 @@ def gunw_slc():
     parser.add_argument('--secondary-scenes', type=str.split, nargs='+', required=True)
     parser.add_argument('--estimate-ionosphere-delay', type=true_false_string_argument, default=False)
     parser.add_argument('--frame-id', type=int, default=-1)
+    parser.add_argument('--compute-solid-earth-tide', type=true_false_string_argument, default=False)
     parser.add_argument('--esd-coherence-threshold', type=float, default=-1.)
     args = parser.parse_args()
 
@@ -162,9 +165,15 @@ def gunw_slc():
                                    reference_properties=ref_properties,
                                    secondary_properties=sec_properties,
                                    extent=extent,
-                                   additional_2d_layers=additional_2d_layers
+                                   additional_2d_layers=additional_2d_layers,
                                    )
 
+    if args.compute_solid_earth_tide:
+        nc_path = update_gunw_with_solid_earth_tide(nc_path)
+        # Update to 1c
+        with h5py.File(nc_path, mode='a') as file:
+            file.attrs.modify('version', '1c')
+
     # Move final product to current working directory
     final_directory = prepare_for_delivery(nc_path, loc_data)
 

isce2_topsapp/packaging_utils/additional_layers.py

@@ -24,9 +24,12 @@ def add_2d_layer(layer_name: str, gunw_netcdf_path: Path) -> Path:
     # The layers generally already exist within the file
     with h5py.File(gunw_netcdf_path, mode='a') as file:
         if dst_group in file:
-            for key in file[dst_group].keys():
-                del file[dst_group][key]
-            del file[dst_group]
+            # Delete the variable to be written to
+            if dst_variable in file[dst_group]:
+                del file[dst_group][dst_variable]
+            # Delete the group if there are no variables left
+            if len(file[dst_group].keys()) == 0:
+                del file[dst_group]
 
     ds = xr.open_dataset(layer_data['input_relative_path'],
                          engine='rasterio')

isce2_topsapp/solid_earth_tides.py

@@ -0,0 +1,167 @@
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# Authors: Simran Sangha, Charles Marshak, Brett Buzzanga, David Bekaert,
+#         Marin Govorcin, Zhang Yunjun
+# Copyright 2023, by the California Institute of Technology. ALL RIGHTS
+# RESERVED. United States Government Sponsorship acknowledged.
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+import datetime as dt
+from pathlib import Path
+from typing import Tuple
+
+import h5py
+import numpy as np
+import pysolid
+import xarray as xr
+
+
+def compute_solid_earth_tide_from_gunw(gunw_path: str) -> xr.Dataset:
+    """ Read GUNW and compute/export differential SET """
+
+    gunw_basename = Path(gunw_path).name
+    dates = gunw_basename.split('-')[6].split('_')
+    ref_time = gunw_basename.split('-')[7]
+
+    # get GUNW attributes
+    (wavelength,
+     z_meta,
+     inc_angle,
+     az_angle,
+     set_attrs) = get_gunw_attrs(gunw_path)
+
+    # compute differential SET ENU
+    # Use reference time twice (assuming exact same pass time)
+    tide_e_ref, tide_n_ref, tide_u_ref = compute_enu_se_tide(dates[0],
+                                                             ref_time,
+                                                             set_attrs)
+    tide_e_sec, tide_n_sec, tide_u_sec = compute_enu_se_tide(dates[1],
+                                                             ref_time,
+                                                             set_attrs)
+    tide_e = tide_e_sec - tide_e_ref
+    tide_n = tide_n_sec - tide_n_ref
+    tide_u = tide_u_sec - tide_u_ref
+
+    # compute SET LOS estimate for each height level
+    tide_los = np.zeros(inc_angle.shape)
+    for i in range(len(z_meta)):
+        tide_los[i] = compute_los_solid_earth_tide(tide_e,
+                                                   tide_n,
+                                                   tide_u,
+                                                   inc_angle[i],
+                                                   az_angle[i],
+                                                   wavelength)
+
+    solidtide_corr_ds = export_se_tides_to_dataset(gunw_path, tide_los)
+
+    return solidtide_corr_ds
+
+
+def get_gunw_attrs(gunw_path: str) -> Tuple[float, np.array, np.array,
+                                            np.array, xr.Dataset, dict]:
+    """ Access necessary GUNW attributes to compute SET """
+
+    group = 'science/radarMetaData'
+    with xr.open_dataset(gunw_path, group=group) as ds:
+        wavelength = ds['wavelength'].item()
+
+    group = 'science/grids/imagingGeometry'
+    with xr.open_dataset(gunw_path, group=group, engine='rasterio') as ds:
+        z_meta = ds.heightsMeta.data
+        lat_meta = ds.y.data
+        lon_meta = ds.x.data
+        # convert angles to rad
+        inc_angle = np.deg2rad(ds.incidenceAngle.data)
+        az_angle = np.deg2rad(ds.azimuthAngle.data-90)
+        solidtide_atr = {
+            'LENGTH': len(lat_meta),
+            'WIDTH': len(lon_meta),
+            'X_FIRST': lon_meta[0],
+            'Y_FIRST': lat_meta[0],
+            'X_STEP':  lon_meta[1] - lon_meta[0],
+            'Y_STEP': lat_meta[1] - lat_meta[0],
+        }
+
+    return wavelength, z_meta, inc_angle, az_angle, solidtide_atr
+
+
+def compute_enu_se_tide(acq_date: str,
+                        acq_time: str,
+                        solidtide_atr: dict) -> np.array:
+    """ Compute SET in ENU """
+    # source:
+    # https://github.com/insarlab/ and
+    # PySolid/blob/main/docs/plot_grid_SET.ipynb
+    dt_obj = dt.datetime.strptime(acq_date + '-' + acq_time,
+                                  "%Y%m%d-%H%M%S")
+    tide_e, tide_n, tide_u = pysolid.calc_solid_earth_tides_grid(dt_obj,
+                                                                 solidtide_atr,
+                                                                 display=False,
+                                                                 verbose=False)
+
+    return tide_e, tide_n, tide_u
+
+
+def compute_los_solid_earth_tide(tide_e: np.array,
+                                 tide_n: np.array,
+                                 tide_u: np.array,
+                                 inc_angle: np.array,
+                                 az_angle: np.array,
+                                 wavelength: float) -> np.array:
+    """ Compute SET in LOS """
+    # source:
+    # https://github.com/insarlab/
+    # PySolid/blob/main/docs/plot_grid_SET.ipynb
+
+    # project ENU to radar line-of-sight (LOS)
+    # with positive for motion towards satellite
+    tide_e_slant = tide_e * np.sin(inc_angle) * np.sin(az_angle) * -1
+    tide_n_slant = tide_n * np.sin(inc_angle) * np.cos(az_angle)
+    tide_u_slant = tide_u * np.cos(inc_angle)
+    tide_los = tide_e_slant + tide_n_slant + tide_u_slant
+
+    # Convert m to rad
+    tide_los = tide_los / (wavelength / (4*np.pi))
+
+    return tide_los
+
+
+def export_se_tides_to_dataset(gunw_path: str,
+                               tide_los: np.array,
+                               lyr_name='solidEarthTide') -> xr.Dataset:
+    """ Create SET array and populate with metadata leveraging
+    the same geodata from the imaging geometry in the gunw"""
+
+    # obtain affine transformation and coordinate metadata
+    group = 'science/grids/imagingGeometry'
+    solidtide_corr_ds = xr.open_dataset(gunw_path,
+                                        group=group)
+
+    solidtide_corr_ds = solidtide_corr_ds.rename({'longitudeMeta': 'longitude',
+                                                  'latitudeMeta': 'latitude',
+                                                  'heightsMeta': 'height'})
+    attrs = {'description': 'Solid Earth tide',
+             'units': 'radians',
+             'long_name': lyr_name,
+             'standard_name': lyr_name}
+
+    dim_order = ['height', 'latitude', 'longitude']
+    solidtide_corr_ds[lyr_name] = (dim_order, tide_los)
+    solidtide_corr_ds[lyr_name].attrs.update(attrs)
+    solidtide_corr_ds = solidtide_corr_ds[[lyr_name]]
+    solidtide_corr_ds = solidtide_corr_ds.astype(np.float32)
+    solidtide_corr_ds.rio.write_crs('epsg:4326', inplace=True)
+    solidtide_corr_ds['solidEarthTide'].rio.write_nodata(0, inplace=True)
+    return solidtide_corr_ds
+
+
+def update_gunw_with_solid_earth_tide(gunw_path: Path) -> Path:
+    se_tide_group = '/science/grids/corrections/external/tides'
+    with h5py.File(gunw_path, 'a') as file:
+        if se_tide_group in file:
+            del file[se_tide_group]
+    solid_earth_tide_ds = compute_solid_earth_tide_from_gunw(gunw_path)
+    solid_earth_tide_ds.to_netcdf(gunw_path,
+                                  mode='a',
+                                  group=se_tide_group)
+    return gunw_path