Merge pull request #74 from ACCESS-Cloud-Based-InSAR/ionosphere

Adding Ionosphere Phase Estimation

.github/workflows/test.yml

@@ -26,12 +26,12 @@ jobs:
           username: ${{ secrets.EARTHDATA_USERNAME }}
           password: ${{ secrets.EARTHDATA_PASSWORD }}
 
-      - uses: conda-incubator/setup-miniconda@v2
+      - uses: mamba-org/provision-with-micromamba@main
         with:
-          mamba-version: "*"
-          python-version: ${{ matrix.python-version }}
-          activate-environment: topsapp_env
+          environment-name: topsapp_env
           environment-file: environment.yml
+          extra-specs: |
+            python=${{ matrix.python-version }}
 
       - name: Pytest in conda environment
         shell: bash -l {0}

.trufflehog.txt

@@ -2,4 +2,6 @@
 .*/templates/.*
 tests/sample_loc_metadata.json
 Dockerfile
+tests/test_data/sec_manifest.xml
+tests/test_data/ref_manifest.xml
 tests/test_data/

CHANGELOG.md

@@ -6,6 +6,14 @@ 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.2]
+
+### Added
+* Provide prototype (internal) for burst analysis thanks to Forrest Williams and Joseph Kennedy (see PR #73)
+* CLI (and API) can switch between burst and SLC ifg generation thanks to entry point magic (see PR #73 for details)
+* Exposes Ionosphere correction in CLI (and API)
+* Exposes ESD and ESD threshold in CLI (and API)
+
 ## [0.2.1]
 
 * Fixes write of start/stop sensing times due to changes in ASF Search v5.0.0 (see #79)

README.md

@@ -49,6 +49,21 @@ This is reflected in the [`sample_run.sh`](sample_run.sh).
 
 To be even more explicity, you can use [`tee`](https://en.wikipedia.org/wiki/Tee_(command)) to record output to both including `> >(tee -a topsapp_img.out) 2> >(tee -a topsapp_img.err >&2)`.
 
+### Customizations
+
+#### Estimating Ionospheric Phase Delay and ESD
+
+This example shows how to obtain a layer with ionsopheric phase delay. The SLCs are over the Arabian peninusula where the ionosphere can be seen:
+
+```
+isce2_topsapp --reference-scenes S1B_IW_SLC__1SDV_20171117T145926_20171117T145953_008323_00EBAB_AFB8 \
+              --secondary-scenes S1A_IW_SLC__1SDV_20171111T150004_20171111T150032_019219_0208AF_EE89 \
+              --estimate-ionosphere-delay True \
+              --do-esd True \
+              --esd-coherence-threshold .5 \
+              > topsapp_img.out 2> topsapp_img.err
+```
+
 # Running with Docker (locally or on a server)
 
 1. When running locally with root privileges (i.e. at your local workstation), build the docker image using:

environment.yml

@@ -34,6 +34,7 @@ dependencies:
  - pytest-cov
  - pytest-mock
  - rasterio
+ - rioxarray
  - setuptools
  - setuptools_scm
  - shapely

isce2_topsapp/__main__.py

@@ -88,8 +88,16 @@ def gunw_slc():
     parser.add_argument('--dry-run', action='store_true')
     parser.add_argument('--reference-scenes', type=str.split, nargs='+', required=True)
     parser.add_argument('--secondary-scenes', type=str.split, nargs='+', required=True)
+    parser.add_argument('--estimate-ionosphere-delay', type=bool, default=False)
+    parser.add_argument('--do-esd', type=bool, default=False)
+    parser.add_argument('--esd-coherence-threshold', type=float, default=-1)
     args = parser.parse_args()
 
+    do_esd_arg = (args.esd_coherence_threshold != -1) == args.do_esd
+    if not do_esd_arg:
+        raise ValueError('If ESD is turned on, specify esd_coherence_threshold between 0 and 1; '
+                         'Otherwise, do not or set the threshold to -1')
+
     ensure_earthdata_credentials(args.username, args.password)
 
     args.reference_scenes = [item for sublist in args.reference_scenes for item in sublist]
@@ -110,19 +118,28 @@ def gunw_slc():
                        secondary_slc_zips=loc_data['sec_paths'],
                        orbit_directory=loc_data['orbit_directory'],
                        extent=loc_data['extent'],
+                       estimate_ionosphere_delay=args.estimate_ionosphere_delay,
+                       do_esd=args.do_esd,
+                       esd_coherence_threshold=args.esd_coherence_threshold,
                        dem_for_proc=loc_data['full_res_dem_path'],
                        dem_for_geoc=loc_data['low_res_dem_path'],
-                       dry_run=args.dry_run
+                       dry_run=args.dry_run,
                        )
 
     ref_properties = loc_data['reference_properties']
     sec_properties = loc_data['secondary_properties']
     extent = loc_data['extent']
 
+    additional_2d_layers = []
+    if args.estimate_ionosphere_delay:
+        additional_2d_layers.append('ionosphere')
+
+    additional_2d_layers = additional_2d_layers or None
     nc_path = package_gunw_product(isce_data_directory=Path.cwd(),
                                    reference_properties=ref_properties,
                                    secondary_properties=sec_properties,
-                                   extent=extent
+                                   extent=extent,
+                                   additional_2d_layers=additional_2d_layers
                                    )
 
     # Move final product to current working directory
@@ -146,6 +163,7 @@ def gunw_burst():
     parser.add_argument('--burst-number', type=int, required=True)
     parser.add_argument('--azimuth-looks', type=int, default=2)
     parser.add_argument('--range-looks', type=int, default=10)
+    parser.add_argument('--estimate-ionosphere-delay', type=bool, default=False)
     args = parser.parse_args()
 
     ensure_earthdata_credentials(args.username, args.password)
@@ -184,6 +202,7 @@ def gunw_burst():
         extent=roi,
         dem_for_proc=dem['full_res_dem_path'],
         dem_for_geoc=dem['low_res_dem_path'],
+        estimate_ionosphere_delay=args.estimate_ionosphere_delay,
         azimuth_looks=args.azimuth_looks,
         range_looks=args.range_looks,
         swaths=[ref_burst.swath],

isce2_topsapp/packaging.py

@@ -5,15 +5,21 @@
 from pathlib import Path
 from typing import Union
 
+import h5py
 import numpy as np
 from dateparser import parse
 
 import isce2_topsapp
+from isce2_topsapp.packaging_utils.additional_layers import add_2d_layer
+from isce2_topsapp.packaging_utils.ionosphere import format_ionosphere_for_gunw
 from isce2_topsapp.templates import read_netcdf_packaging_template
 
 DATASET_VERSION = '2.0.6'
 
 
+PERMISSIBLE_2D_LAYERS = ['ionosphere']
+
+
 """Warning: the packaging scripts were written as command line scripts and
 are highly dependent on the current working directory and its structure.
 
@@ -211,24 +217,51 @@ def perform_netcdf_packaging(*,
     return out_nc_file
 
 
+def package_additional_layers_into_gunw(gunw_path: Path,
+                                        isce_data_directory: Path,
+                                        additional_2d_layers: list):
+    # Current workflow of additional layers
+    # 1. Do any additional processing/formatting outside of GUNW
+    # 2. Add layer into GUNW
+    # 3. Update Version
+    if not set(additional_2d_layers).issubset(set(PERMISSIBLE_2D_LAYERS)):
+        raise RuntimeError('Additional 2d layers must be subset of '
+                           f'{PERMISSIBLE_2D_LAYERS}')
+
+    if 'ionosphere' in additional_2d_layers:
+        # current working directory is ISCE directory
+        _ = format_ionosphere_for_gunw(isce_data_directory, gunw_path)
+
+    # Assumes ionosphere raster is written to specific path
+    [add_2d_layer(layer, gunw_path) for layer in additional_2d_layers]
+
+    # Update
+    with h5py.File(gunw_path, mode='a') as file:
+        file.attrs.modify('version', '1c')
+    return gunw_path
+
+
 def package_gunw_product(*,
                          isce_data_directory: Union[str, Path],
                          reference_properties: list,
                          secondary_properties: list,
-                         extent: list) -> Path:
+                         extent: list,
+                         additional_2d_layers: list = None) -> Path:
     """Creates a GUNW standard product netcdf from the ISCE outputs and some
     initial metadata.
 
     Parameters
     ----------
-    isce_data_directory
+    isce_data_directory: str | path
         Where the ISCE outputs are relative to current working directory
-    reference_properties
+    reference_properties: list
         Each item a dictionary per ASF API including ID, starttime, etc.
-    secondary_properties
+    secondary_properties: list
         Each item a dictionary per ASF API including ID, starttime, etc
-    extent
+    extent: list
         List of extents ([xmin, ymin, xmax, ymax])
+    additional_2d_layers: list
+        List of 2d layers to add. Currently, supported is ionosphere.
 
     Returns
     -------
@@ -244,4 +277,8 @@ def package_gunw_product(*,
     out_nc_file = perform_netcdf_packaging(isce_data_dir=isce_data_directory,
                                            gunw_id=gunw_id)
 
+    if additional_2d_layers is not None:
+        package_additional_layers_into_gunw(out_nc_file,
+                                            isce_data_directory,
+                                            additional_2d_layers)
     return out_nc_file

isce2_topsapp/packaging_utils/additional_layers.json

@@ -0,0 +1,11 @@
+{"ionosphere": {"dst_group": "/science/grids/corrections/derived/ionosphere",
+                "dst_variable": "ionosphere",
+                             "input_relative_path": "merged/ionosphere_for_gunw.geo",
+                             "attrs": {
+                                       "standard_name": "ionospherePhaseCorrection",
+                                       "long_name": "ionospherePhaseCorrection",
+                                       "description": "Estimated ionosphere phase correction"
+                                       }
+                             }
+
+}

isce2_topsapp/packaging_utils/additional_layers.py

@@ -0,0 +1,51 @@
+import json
+from pathlib import Path
+
+import h5py
+import xarray as xr
+
+LAYER_JSON = Path(__file__).parents[0] / 'additional_layers.json'
+ADDITIONAL_LAYERS = json.load(open(LAYER_JSON))
+
+
+def add_2d_layer(layer_name: str, gunw_netcdf_path: Path) -> Path:
+    """
+    Combines a lot of standard formatting of the netcdf via rioxarray and
+    deletes the previous placeholder (we assume it exists via the placeholder).
+
+    We also assume any additional processing specific to GUNW is done outside of
+    this function.
+    """
+
+    layer_data = ADDITIONAL_LAYERS[layer_name]
+    dst_group = layer_data['dst_group']
+    dst_variable = layer_data['dst_variable']
+
+    # 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]
+
+    ds = xr.open_dataset(layer_data['input_relative_path'],
+                         engine='rasterio')
+
+    # Renaming ensures correct geo-referencing with spatial_ref grid mapping
+    ds = ds.rename({'x': 'longitude',
+                    'y': 'latitude',
+                    'band_data': dst_variable})
+    ds['latitude'].attrs.update({'long_name': 'latitude',
+                                 'standard_name': 'latitude'})
+    ds['longitude'].attrs.update({'long_name': 'longitude',
+                                  'standard_name': 'longitude'})
+
+    # removes channel (aka band) dimension
+    ds = ds.squeeze(['band'], drop=True)
+    ds[layer_name].attrs.update(layer_data['attrs'])
+
+    ds.to_netcdf(gunw_netcdf_path,
+                 group=layer_data['dst_group'],
+                 mode='a')
+
+    return gunw_netcdf_path

isce2_topsapp/packaging_utils/ionosphere.py

@@ -0,0 +1,47 @@
+from pathlib import Path
+
+import numpy as np
+import rasterio
+from dem_stitcher.rio_tools import (reproject_arr_to_match_profile,
+                                    update_profile_resolution)
+
+
+def format_ionosphere_for_gunw(isce_directory: Path,
+                               gunw_netcdf_path: Path) -> Path:
+    with rasterio.open(isce_directory / "merged/topophase.ion.geo") as ds:
+        X_ion = ds.read(1)
+        p_ion = ds.profile
+
+    X_ion[X_ion == 0] = np.nan
+    p_ion['nodata'] = np.nan
+
+    # Get GUNW Mask
+    nc_path_str = (f'netcdf:{gunw_netcdf_path}:'
+                   '/science/grids/data/connectedComponents')
+    with rasterio.open(nc_path_str) as ds:
+        cc = ds.read(1)
+        p_cc = ds.profile
+    mask = (cc == -1).astype(np.int32)
+
+    # Lower resolution by factor of 11 to approximatly 990 m at equator
+    p_ion_low_res = update_profile_resolution(p_ion, 0.0091666666)
+    X_ion_low_res, _ = reproject_arr_to_match_profile(X_ion,
+                                                      p_ion,
+                                                      p_ion_low_res,
+                                                      resampling='bilinear')
+    p_cc['nodata'] = None
+    p_cc['dtype'] = np.int32
+    mask_low_res, _ = reproject_arr_to_match_profile(mask,
+                                                     p_cc,
+                                                     p_ion_low_res,
+                                                     resampling='nearest')
+
+    X_ion_low_res = X_ion_low_res[0, ...]
+    mask_low_res = mask_low_res[0, ...].astype(bool)
+    X_ion_low_res[mask_low_res] = np.nan
+
+    out_path = isce_directory / 'merged/ionosphere_for_gunw.geo'
+    with rasterio.open(out_path, 'w', **p_ion_low_res) as ds:
+        ds.write(X_ion_low_res, 1)
+
+    return out_path

isce2_topsapp/templates/topsapp_template.xml

@@ -22,9 +22,11 @@
     <property name="geocodeDemFilename">{{ geocodeDemFilename }}</property>
     <property name="do unwrap">True</property>
     <property name="unwrapper name">snaphu_mcf</property>
+    <property name="do ionosphere correction">{{ estimate_ionosphere_delay }}</property>
+    <property name="apply ionosphere correction">False</property>
     <property name="do ESD">{{ do_esd }}</property>
     <property name="ESD coherence threshold">{{ esd_coherence_threshold }}</property>
     <property name="use virtual files">{{ use_virtual_files }}</property>
-    <property name="geocode list">['merged/phsig.cor', 'merged/filt_topophase.unw', 'merged/los.rdr', 'merged/topophase.flat', 'merged/filt_topophase.flat', 'merged/filt_topophase_2stage.unw', 'merged/topophase.cor', 'merged/filt_topophase.unw.conncomp']</property>
+    <property name="geocode list">{{ geocode_list }}</property>
   </component>
 </topsApp>