Merge pull request #125 from ASFHyP3/seperate-dems

Resample the Geocoding DEM

CHANGELOG.md

@@ -7,6 +7,10 @@ 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.6.2]
+### Changed
+* The geocoding DEM is now resampled to ~20m in the case of 5x1 looks.
+
 ## [0.6.1]
 ### Added
 * `apply_water_mask` optional argument to apply water mask in the wrapped and unwrapped phase geotiff files

src/hyp3_isce2/dem.py

@@ -56,11 +56,32 @@ def buffer_extent(extent: list, buffer: float) -> list:
     ]
 
 
+def distance_meters_to_degrees(distance_meters, latitude):
+    """Convert a distance from meters to degrees in longitude and latitude
+
+    Args:
+        distance_meters: Arc length in meters.
+        latitude: The line of latitude at which the calculation takes place.
+    Returns:
+        The length in degrees for longitude and latitude, respectively.
+    """
+    if np.abs(latitude) == 90:
+        # np.cos won't return exactly 0, so we must manually raise this exception.
+        raise ZeroDivisionError('A Latitude of 90 degrees results in dividing by zero.')
+    EARTHS_CIRCUMFERENCE = 40_030_173.59204114  # 2 * pi * 6,371,000.0 (Earth's average radius in meters)
+    latitude_circumference = EARTHS_CIRCUMFERENCE * np.cos(np.radians(latitude))
+    distance_degrees_lon = distance_meters / latitude_circumference * 360
+    distance_degrees_lat = distance_meters / EARTHS_CIRCUMFERENCE * 360
+    return (np.round(distance_degrees_lon, 15), np.round(distance_degrees_lat, 15))
+
+
 def download_dem_for_isce2(
         extent: list,
         dem_name: str = 'glo_30',
         dem_dir: Path = None,
-        buffer: float = .4) -> Path:
+        buffer: float = .4,
+        resample_20m: bool = False
+) -> Path:
     """Download the given DEM for the given extent.
 
     Args:
@@ -69,6 +90,7 @@ def download_dem_for_isce2(
         dem_dir: The output directory.
         buffer: The extent buffer in degrees, by default .4, which is about 44 km at the equator
                 (or about 2.5 bursts at the equator).
+        resample_20m: Whether or not the DEM should be resampled to 20 meters.
     Returns:
         The path to the downloaded DEM.
     """
@@ -77,13 +99,27 @@ def download_dem_for_isce2(
 
     extent_buffered = buffer_extent(extent, buffer)
 
-    dem_array, dem_profile = dem_stitcher.stitch_dem(
-        extent_buffered,
-        dem_name,
-        dst_ellipsoidal_height=True,
-        dst_area_or_point='Point',
-        n_threads_downloading=5,
-    )
+    if resample_20m:
+        res_degrees = distance_meters_to_degrees(20.0, extent_buffered[1])
+        dem_array, dem_profile = dem_stitcher.stitch_dem(
+            extent_buffered,
+            dem_name,
+            dst_ellipsoidal_height=True,
+            dst_area_or_point='Point',
+            n_threads_downloading=5,
+            dst_resolution=res_degrees
+        )
+        dem_path = dem_dir / 'full_res_geocode.dem.wgs84'
+    else:
+        dem_array, dem_profile = dem_stitcher.stitch_dem(
+            extent_buffered,
+            dem_name,
+            dst_ellipsoidal_height=True,
+            dst_area_or_point='Point',
+            n_threads_downloading=5,
+        )
+        dem_path = dem_dir / 'full_res.dem.wgs84'
+
     dem_array[np.isnan(dem_array)] = 0.
 
     dem_profile['nodata'] = None
@@ -93,7 +129,6 @@ def download_dem_for_isce2(
     for key in ['blockxsize', 'blockysize', 'compress', 'interleave', 'tiled']:
         del dem_profile[key]
 
-    dem_path = dem_dir / 'full_res.dem.wgs84'
     with rasterio.open(dem_path, 'w', **dem_profile) as ds:
         ds.write(dem_array, 1)
 

src/hyp3_isce2/insar_tops_burst.py

@@ -84,9 +84,26 @@ def insar_tops_burst(
     log.info(f'InSAR ROI: {insar_roi}')
     log.info(f'DEM ROI: {dem_roi}')
 
-    dem_path = download_dem_for_isce2(dem_roi, dem_name='glo_30', dem_dir=dem_dir, buffer=0)
+    dem_path = download_dem_for_isce2(
+        dem_roi,
+        dem_name='glo_30',
+        dem_dir=dem_dir,
+        buffer=0,
+        resample_20m=False
+    )
     download_aux_cal(aux_cal_dir)
 
+    if range_looks == 5:
+        geocode_dem_path = download_dem_for_isce2(
+            dem_roi,
+            dem_name='glo_30',
+            dem_dir=dem_dir,
+            buffer=0,
+            resample_20m=True
+        )
+    else:
+        geocode_dem_path = dem_path
+
     orbit_dir.mkdir(exist_ok=True, parents=True)
     for granule in (ref_params.granule, sec_params.granule):
         downloadSentinelOrbitFile(granule, str(orbit_dir))
@@ -98,6 +115,7 @@ def insar_tops_burst(
         aux_cal_directory=str(aux_cal_dir),
         roi=insar_roi,
         dem_filename=str(dem_path),
+        geocode_dem_filename=str(geocode_dem_path),
         swaths=swath_number,
         azimuth_looks=azimuth_looks,
         range_looks=range_looks,

src/hyp3_isce2/topsapp.py

@@ -51,6 +51,7 @@ def __init__(
         orbit_directory: str,
         aux_cal_directory: str,
         dem_filename: str,
+        geocode_dem_filename: str,
         roi: Iterable[float],
         swaths: int or Iterable[int] = [1, 2, 3],
         azimuth_looks: int = 4,
@@ -63,7 +64,7 @@ def __init__(
         self.aux_cal_directory = aux_cal_directory
         self.roi = [roi[1], roi[3], roi[0], roi[2]]
         self.dem_filename = dem_filename
-        self.geocode_dem_filename = dem_filename
+        self.geocode_dem_filename = geocode_dem_filename
         self.azimuth_looks = azimuth_looks
         self.range_looks = range_looks
         self.do_unwrap = do_unwrap

tests/test_dem.py

@@ -4,6 +4,7 @@
 import rasterio
 from affine import Affine
 from lxml import etree
+from pytest import raises
 from rasterio import CRS
 
 from hyp3_isce2 import dem
@@ -78,3 +79,16 @@ def test_buffer_extent():
     assert dem.buffer_extent(extent2, 0.1) == [-170, 53, -167, 55]
     assert dem.buffer_extent(extent2, 0.3) == [-170, 53, -167, 55]
     assert dem.buffer_extent(extent2, 0.4) == [-171, 52, -166, 56]
+
+
+def test_distance_meters_to_degrees():
+    assert dem.distance_meters_to_degrees(distance_meters=20, latitude=0) == (0.000179864321184, 0.000179864321184)
+    assert dem.distance_meters_to_degrees(distance_meters=20, latitude=45) == (0.000254366562405, 0.000179864321184)
+    assert dem.distance_meters_to_degrees(distance_meters=20, latitude=89.9) == (0.103054717208573, 0.000179864321184)
+    assert dem.distance_meters_to_degrees(distance_meters=20, latitude=-45) == (0.000254366562405, 0.000179864321184)
+    assert dem.distance_meters_to_degrees(distance_meters=20, latitude=-89.9) == (0.103054717208573, 0.000179864321184)
+    # This is since cos(90) = 0, leading to a divide by zero issue.
+    with raises(ZeroDivisionError):
+        dem.distance_meters_to_degrees(20, 90)
+    with raises(ZeroDivisionError):
+        dem.distance_meters_to_degrees(20, -90)

tests/test_topsapp.py

@@ -11,6 +11,7 @@ def test_topsapp_burst_config(tmp_path):
         aux_cal_directory='aux_cal',
         roi=[-118.0, 37.0, -117.0, 38.0],
         dem_filename='dem.tif',
+        geocode_dem_filename='dem_geocode.tif',
         swaths=1,
     )
 
@@ -25,6 +26,7 @@ def test_topsapp_burst_config(tmp_path):
         assert 'orbits' in template
         assert 'aux_cal' in template
         assert 'dem.tif' in template
+        assert 'dem_geocode.tif' in template
         assert '[37.0, 38.0, -118.0, -117.0]' in template
         assert '[1]' in template
 
@@ -55,6 +57,7 @@ def test_run_topsapp_burst(tmp_path, monkeypatch):
         aux_cal_directory='',
         roi=[0, 1, 2, 3],
         dem_filename='',
+        geocode_dem_filename='',
         swaths=1,
         azimuth_looks=1,
         range_looks=1,