Add provisional data and update isce3

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.3.1]
+
+### Changed
+- Updated ISCE3 version and included provisional collection for GUNWs.
+
 ## [0.3.0]
 
 > [!IMPORTANT]

pyproject.toml

@@ -128,7 +128,7 @@ numpy = "<2.4"
 opencv = "*"
 pandas = ">=1.4,<3.0"
 pyaps3="*"
-pydap="*"
+pydap="<3.5.11"
 pygrib="*"
 pyproj = ">=3.3"
 python-dateutil = "*"
@@ -170,7 +170,7 @@ pyre = "*"
 pyre = "*"
 
 [tool.pixi.pypi-dependencies]
-snaphu = "*"
+whirlwind-insar = ">=0.10.0"
 raider = {git = "https://github.com/dbekaert/RAiDER.git", rev = "v0.5.3"}
 
 [tool.pixi.feature.pycharm.dependencies]
@@ -192,15 +192,19 @@ py313 = { features = ["develop", "py313"]}
 py314 = { features = ["develop", "py314"], solve-group = "default"}
 
 
+[tool.pixi.tasks.install-snaphu]
+cmd = "pip install snaphu --no-cache-dir"
+env = { SKBUILD_STRICT_CONFIG = "0", CMAKE_ARGS = "-DCMAKE_C_FLAGS=-Wno-error=incompatible-pointer-types" }
+
 [tool.pixi.tasks.install-isce3]
 cmd = [
     "python",
     "-m",
     "pip",
     "install",
-    "git+https://github.com/isce-framework/isce3.git@bd58fcd3dad94f7f06871f927b902eca20680439"
+    "git+https://github.com/isce-framework/isce3.git@4d2e5f43e6e873972e59a21123f917d1d5dc794a"
 ]
-depends-on = []
+depends-on = ["install-snaphu"]
 
 [tool.pixi.tasks.install-editable]
 cmd = [

src/hyp3_isce3/__main__.py

@@ -73,7 +73,7 @@ def main() -> None:
 
         if len(args.subset) == 0:
             args.subset = None
-        elif not len(args.subset) == 4:
+        elif len(args.subset) != 4:
             raise ValueError('The number of coordinates is not four')
 
     username = os.getenv('EARTHDATA_USERNAME')

src/hyp3_isce3/process.py

@@ -3,17 +3,16 @@
 import argparse
 import logging
 import zipfile
-from datetime import datetime, timedelta
+from datetime import UTC, datetime, timedelta
 from pathlib import Path
 
 import asf_search as asf
 import earthaccess
 import utm
 import yaml
-from hyp3lib.dem import prepare_dem_geotiff
 from nisar.workflows import h5_prep, insar, stage_dem
 from nisar.workflows.insar_runconfig import InsarRunConfig
-from osgeo import ogr, osr
+from osgeo import gdal
 
 import hyp3_isce3
 from hyp3_isce3.crop_rslc import crop_streamed, geocode_subset_box, stream_skeleton
@@ -32,6 +31,7 @@ def get_config(
     secondary_orbit: str,
     reference_tropo: str,
     secondary_tropo: str,
+    dem_path: str,
     tec_path: str,
     watermask: str,
     template_yaml: Path,
@@ -47,6 +47,7 @@ def get_config(
         secondary_orbit: Path of the secondary orbit.
         reference_tropo: Path of the ECMWF file for the reference scene.
         secondary_tropo: Path of the ECMWF file for the secondary scene.
+        dem_path: Path to DEM file.
         tec_path: Path of the TEC file for the reference scene.
         watermask: Path of the water mask file.
         template_yaml: Path of the downloaded JPL runconfig (from :func:`download_yaml`) used as the tail/template.
@@ -86,6 +87,8 @@ def get_config(
                 newstring += line.replace('ref_tropo', reference_tropo)
             elif 'secondary_tropo' in line:
                 newstring += line.replace('sec_tropo', secondary_tropo)
+            elif 'dem_file' in line:
+                newstring += line.replace('dem_path', dem_path)
             elif 'tec_file' in line:
                 newstring += line.replace('tec_path', tec_path)
             elif 'watermask' in line:
@@ -143,9 +146,15 @@ def download_yaml(reference_path: str) -> Path:
     Returns:
         tmp_path: Path of the yaml file.
     """
-    short_name = 'NISAR_L2_GUNW_BETA_V1'
-    keyword = '_'.join(reference_path.split('_')[4:8])
-    results = earthaccess.search_data(short_name=short_name, granule_name=f'*{keyword}*')
+    keyword = '_'.join(reference_path.split('_')[5:8])
+    # Prefer the PROVISIONAL template (current production settings); fall back to BETA.
+    short_names = ['NISAR_L2_GUNW_PROVISIONAL_V1', 'NISAR_L2_GUNW_BETA_V1']
+    for short_name in short_names:
+        results = earthaccess.search_data(short_name=short_name, granule_name=f'*{keyword}*')
+        if results:
+            break
+    else:
+        raise ValueError(f'No GUNW granule found for {keyword} in {short_names}')
     gunw = results[0].data_links()[0].split('/')[-2]
     res = asf.granule_search(gunw)
     yaml_url = res.find_urls(pattern=r'.yaml')[0]
@@ -187,8 +196,8 @@ def get_orbit(scene_name: str) -> str:
         orbit_path: Path of the orbit file.
     """
     short_name = 'NISAR_OE'
-    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S')
-    end_date = datetime.strptime(scene_name.split('_')[12], '%Y%m%dT%H%M%S')
+    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S').replace(tzinfo=UTC)
+    end_date = datetime.strptime(scene_name.split('_')[12], '%Y%m%dT%H%M%S').replace(tzinfo=UTC)
     temporal = (start_date.strftime('%Y-%m-%d %H:%M:%S'), end_date.strftime('%Y-%m-%d %H:%M:%S'))
     results = earthaccess.search_data(short_name=short_name, granule_name='*POE*', temporal=temporal)
     if len(results) == 0:
@@ -214,8 +223,8 @@ def get_tropo(scene_name: str) -> str:
         tropo_path: Path of the file.
     """
     short_name = 'ASF_ECMWF_TROP'
-    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S')
-    day = datetime(start_date.year, start_date.month, start_date.day)
+    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S').replace(tzinfo=UTC)
+    day = datetime(start_date.year, start_date.month, start_date.day, tzinfo=UTC)
     if start_date.hour % 6 < 3:
         tropo_date = day + timedelta(hours=int(start_date.hour / 6) * 6)
     else:
@@ -239,8 +248,8 @@ def get_tec(scene_name: str) -> str:
         tropo_path: Path of the file.
     """
     short_name = 'NISAR_TEC'
-    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S')
-    end_date = datetime.strptime(scene_name.split('_')[12], '%Y%m%dT%H%M%S')
+    start_date = datetime.strptime(scene_name.split('_')[11], '%Y%m%dT%H%M%S').replace(tzinfo=UTC)
+    end_date = datetime.strptime(scene_name.split('_')[12], '%Y%m%dT%H%M%S').replace(tzinfo=UTC)
     temporal = (start_date.strftime('%Y-%m-%d %H:%M:%S'), end_date.strftime('%Y-%m-%d %H:%M:%S'))
     results = earthaccess.search_data(short_name=short_name, temporal=temporal)
     files = sorted(earthaccess.download(results))
@@ -261,37 +270,65 @@ def get_watermask(reference_path: str, subset: list[float] | None = None) -> str
         tropo_path: Path of the file.
     """
     short_name = 'NISAR_WATERMASK'
-    if subset is None:
+    if subset:
+        bbox = tuple(subset)
+    else:
         poly, _ = stage_dem.determine_polygon(reference_path, bbox=None, bbox_epsg='4326')
         bbox = poly.bounds
-    else:
-        bbox = subset
     bbox = (bbox[0] - 1, bbox[1] - 1, bbox[2] + 1, bbox[3] + 1)
     results = earthaccess.search_data(short_name=short_name, bounding_box=bbox)
     files = sorted(earthaccess.download(results))
-
-    return str(files[-1])
+    folder = files[0].parent
+    files = [f.resolve() for f in files if '.tif' in f.name]
+
+    output_raster = folder / 'watermask.tif'
+    vrt_dataset = gdal.BuildVRT(str(folder / 'mosaic.vrt'), files)
+    vrt_dataset = None
+    vrt_dataset = folder / 'mosaic.vrt'
+    # gdal.Warp(str(output_raster), vrt_dataset, dstSRS="EPSG:4326")
+    gdal.Warp(str(output_raster), vrt_dataset)
+
+    if output_raster.exists():
+        return str(output_raster)
+    else:
+        raise RuntimeError('watermask could not be downloaded')
 
 
-def get_dem(scene_poly: ogr.Geometry, epsg_code: int, dem_path: str = 'dem.tif') -> str:
-    """Download DEM for a given polygon.
+def get_dem(reference_path: str, subset: list[float] | None = None) -> str:
+    """Download the NISAR DEM tiles covering the scene and mosaic them into one GeoTIFF.
 
     Args:
-        scene_poly: Scene polygon.
-        epsg_code: EPSG code for the output projection.
-        dem_path: Output path for the DEM.
+        reference_path: Path of the reference scene.
+        subset: Optional AOI [lon_min, lat_min, lon_max, lat_max]; when set, the DEM
+            is fetched over the AOI instead of the whole frame (the 1-degree buffer
+            below covers the crop margin).
 
     Returns:
-        dem_path: Path of the DEM file.
+        dem_path: Path of the mosaicked DEM file.
     """
-    return str(
-        prepare_dem_geotiff(
-            output_name=dem_path,
-            geometry=scene_poly,
-            epsg_code=4326,
-            pixel_size=0.001,
-        )
-    )
+    short_name = 'NISAR_DEM'
+    if subset:
+        bbox = tuple(subset)
+    else:
+        poly, _ = stage_dem.determine_polygon(reference_path, bbox=None, bbox_epsg='4326')
+        bbox = poly.bounds
+    bbox = (bbox[0] - 1, bbox[1] - 1, bbox[2] + 1, bbox[3] + 1)
+    results = earthaccess.search_data(short_name=short_name, bounding_box=bbox)
+    files = sorted(earthaccess.download(results))
+    folder = files[0].parent
+    files = [f.resolve() for f in files if '.tif' in f.name]
+
+    output_raster = folder / 'dem.tif'
+    vrt_dataset = gdal.BuildVRT(str(folder / 'mosaic.vrt'), files)
+    vrt_dataset = None
+    vrt_dataset = folder / 'mosaic.vrt'
+    # gdal.Warp(str(output_raster), vrt_dataset, dstSRS="EPSG:4326")
+    gdal.Warp(str(output_raster), vrt_dataset)
+
+    if output_raster.exists():
+        return str(output_raster)
+    else:
+        raise RuntimeError('DEM could not be downloaded')
 
 
 def get_epsg(lat: float, lon: float) -> int:
@@ -312,36 +349,20 @@ def get_epsg(lat: float, lon: float) -> int:
     return epsg_base + zone_number
 
 
-def get_scene_polygon(reference_path: str, subset: list[float] | None = None) -> ogr.Geometry:
-    """Get Polygon for reference scene.
+def get_scene_epsg(reference_path: str) -> int:
+    """Get the UTM EPSG code for the reference scene.
+
+    Taken from the full scene's footprint centroid, so a subset run uses the same
+    projection as a full-frame run.
 
     Args:
-        reference_path: Path of the downloaded h5 file.
-        subset: Optional AOI [lon_min, lat_min, lon_max, lat_max]; when set, the DEM is
-            staged over the AOI (plus a buffer for the crop margin and radar-processing
-            edges) instead of the whole frame. EPSG is still taken from the full scene.
+        reference_path: Path of the reference scene (full product or skeleton).
 
     Returns:
-        geom: Polygon of the reference scene.
+        epsg_code: UTM EPSG code of the scene centroid.
     """
     poly, _ = stage_dem.determine_polygon(reference_path, bbox=None, bbox_epsg='4326')
-    epsg_code = get_epsg(poly.centroid.y, poly.centroid.x)
-    if subset is None:
-        poly, _ = stage_dem.determine_polygon(reference_path, bbox=None, bbox_epsg=str(epsg_code))
-    else:
-        # Buffer the AOI past the 512-px crop margin's ground extent (~5-6 km); the
-        # extra apply_margin_to_geographic_box 5 km below then adds further headroom.
-        buf = 0.1  # degrees (~11 km)
-        bbox = [subset[0] - buf, subset[1] - buf, subset[2] + buf, subset[3] + buf]
-        poly, _ = stage_dem.determine_polygon(reference_path, bbox=bbox, bbox_epsg='4326')
-    poly = stage_dem.apply_margin_to_geographic_box(poly)
-    geom = ogr.CreateGeometryFromWkt(str(poly))
-
-    srs = osr.SpatialReference()
-    srs.ImportFromEPSG(epsg_code)
-    geom.AssignSpatialReference(srs)
-
-    return geom, epsg_code
+    return get_epsg(poly.centroid.y, poly.centroid.x)
 
 
 def get_product_id(reference_scene: str, secondary_scene: str) -> str:
@@ -404,8 +425,8 @@ def process_isce3(reference_scene: str, secondary_scene: str, subset: list[float
 
     tec_path = get_tec(reference_scene)
 
-    scene_polygon, epsg_code = get_scene_polygon(reference_path, subset)
-    dem_path = get_dem(scene_polygon, epsg_code)
+    epsg_code = get_scene_epsg(reference_path)
+    dem_path = get_dem(reference_path, subset)
 
     # The JPL runconfig is both our config template (its tail) and the source of the
     # crossmul looks the crop aligns to; download once and reuse for both.
@@ -436,6 +457,7 @@ def process_isce3(reference_scene: str, secondary_scene: str, subset: list[float
         secondary_orbit,
         reference_tropo,
         secondary_tropo,
+        dem_path,
         tec_path,
         watermask,
         template_yaml,

src/hyp3_isce3/schemas/insar.yaml

@@ -11,7 +11,7 @@ runconfig:
             qa_runw_input_file: output/RUNW_product.h5
             qa_gunw_input_file: output/GUNW_product.h5
         dynamic_ancillary_file_group:
-            dem_file: dem.tif
+            dem_file: dem_path
             dem_file_description: Digital Elevation Model (DEM) for the 
                 NASA-ISRO NISAR mission, version 1.2, based on the Copernicus 
                 DEM 30-m (2023_1), referenced to the WGS84 ellipsoid. This 

tests/test_process.py

@@ -14,6 +14,8 @@ def test_get_config(monkeypatch, tmp_path):
     reference_tropo = 'REFERENCE_TROPO.nc'
     secondary_tropo = 'SECONDARY_TROPO.nc'
 
+    dem_path = 'DEM.tif'
+
     tec_path = 'TEC.json'
     watermask = 'WATERMASK.vrt'
 
@@ -27,6 +29,7 @@ def test_get_config(monkeypatch, tmp_path):
         secondary_orbit,
         reference_tropo,
         secondary_tropo,
+        dem_path,
         tec_path,
         watermask,
         temp_yaml,
@@ -112,6 +115,7 @@ def test_get_config_subset(monkeypatch, tmp_path):
         'SECORB.xml',
         'REFTROP.nc',
         'SECTROP.nc',
+        'DEM.tif',
         'TEC.json',
         'WMASK.vrt',
         Path('temp.yaml'),