Skip to content

Commit a881e3e

Browse files
update based on review feedback
1 parent 31a7a7c commit a881e3e

File tree

2 files changed

+22
-21
lines changed

2 files changed

+22
-21
lines changed

src/hyp3_opera_rtc/dem.py

Lines changed: 20 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
1-
import os
21
from collections.abc import Callable
32
from pathlib import Path
3+
from tempfile import NamedTemporaryFile
44

5-
import backoff
65
import numpy as np
76
import shapely
87
import shapely.ops
98
import shapely.wkt
9+
from hyp3lib.util import GDALConfigManager
1010
from osgeo import gdal, osr
1111
from shapely.geometry import LinearRing, Polygon, box
1212

@@ -49,7 +49,7 @@ def polygon_from_bounds(bounds: tuple[float, float, float, float]) -> Polygon:
4949
return poly
5050

5151

52-
def check_antimeridean(poly: Polygon) -> list[Polygon]:
52+
def split_antimeridean(poly: Polygon) -> list[Polygon]:
5353
"""Check if the provided polygon crosses the antimeridian and split it if it does."""
5454
x_min, _, x_max, _ = poly.bounds
5555

@@ -92,14 +92,13 @@ def snap_coord(val: float, snap: float, offset: float, round_func: Callable) ->
9292
return round_func(float(val - offset) / snap) * snap + offset
9393

9494

95-
@backoff.on_exception(backoff.expo, Exception, max_time=600, max_value=32)
9695
def translate_dem(vrt_filename: str, output_path: str, bounds: tuple[float, float, float, float]) -> None:
97-
"""Translate the OPERA DEM from S3 to a region matching the provided boundaries.
96+
"""Write a local subset of the OPERA DEM for a region matching the provided bounds.
9897
9998
Params:
10099
vrt_filename: Path to the input VRT file
101100
output_path: Path to the translated output GTiff file
102-
bounds: tuple of (x_min, x_max, y_min, y_max)
101+
bounds: Bounding box in the form of (lon_min, lat_min, lon_max, lat_max)
103102
"""
104103
ds = gdal.Open(vrt_filename, gdal.GA_ReadOnly)
105104

@@ -159,21 +158,23 @@ def download_opera_dem_for_footprint(outfile: Path, bounds: tuple[float, float,
159158
"""Download a DEM from the specified S3 bucket.
160159
161160
Params:
162-
polys: List of shapely polygons.
163161
outfile: Path to the where the output DEM file is to be staged.
162+
bounds: Bounding box in the form of (lon_min, lat_min, lon_max, lat_max).
164163
"""
165164
poly = polygon_from_bounds(bounds)
166-
polys = check_antimeridean(poly)
165+
polys = split_antimeridean(poly)
167166
dem_list = []
168167

169-
os.environ['GDAL_HTTP_COOKIEJAR'] = 'cookies.txt'
170-
os.environ['GDAL_HTTP_COOKIEFILE'] = 'cookies.txt'
171-
os.environ['GDAL_DISABLE_READDIR_ON_OPEN'] = 'EMPTY_DIR'
172-
173-
vrt_filename = '/vsicurl/https://nisar.asf.earthdatacloud.nasa.gov/STATIC/DEM/v1.2/EPSG4326/EPSG4326.vrt'
174-
for idx, poly in enumerate(polys):
175-
output_path = f'{outfile.stem}_{idx}.tif'
176-
dem_list.append(output_path)
177-
translate_dem(vrt_filename, output_path, bounds)
178-
179-
gdal.BuildVRT(str(outfile), dem_list)
168+
with NamedTemporaryFile(suffix='.txt') as cookie_file:
169+
with GDALConfigManager(
170+
GDAL_HTTP_COOKIEJAR=cookie_file.name,
171+
GDAL_HTTP_COOKIEFILE=cookie_file.name,
172+
GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
173+
):
174+
vrt_filename = '/vsicurl/https://nisar.asf.earthdatacloud.nasa.gov/STATIC/DEM/v1.2/EPSG4326/EPSG4326.vrt'
175+
for idx, poly in enumerate(polys):
176+
output_path = f'{outfile.stem}_{idx}.tif'
177+
dem_list.append(output_path)
178+
translate_dem(vrt_filename, output_path, bounds)
179+
180+
gdal.BuildVRT(str(outfile), dem_list)

tests/test_dem.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,12 +29,12 @@ def test_polygon_from_bounds():
2929

3030
def test_check_antimeridean():
3131
no_cross = box(-1, -1, 0, 0)
32-
polys = dem.check_antimeridean(no_cross)
32+
polys = dem.split_antimeridean(no_cross)
3333
assert len(polys) == 1
3434
assert polys[0].equals(no_cross)
3535

3636
cross = box(179, -1, 181, 0)
37-
polys = dem.check_antimeridean(cross)
37+
polys = dem.split_antimeridean(cross)
3838
negative_side = box(-180, -1, -179, 0)
3939
positive_side = box(179, -1, 180, 0)
4040
assert len(polys) == 2

0 commit comments

Comments
 (0)