Break area request into tiles that want hit opentopography request co…

bris-adapt/bris_adapt/orography/area.py

@@ -0,0 +1,38 @@
+import math
+
+EARTH_MEAN_RADIUS_KM = 6371.0  # mean Earth radius in km
+
+
+def norm_lon(lon: float) -> float:
+    """Map lon to [-180, 180]."""
+    if math.fabs(lon) <= 180.0:
+        return math.copysign(((lon + 180.0) % 360.0) - 180.0, lon)
+    else:
+        return ((lon + 180.0) % 360.0) - 180.0
+
+
+def bbox_area_km2(north: float, west: float, south: float, east: float) -> float:
+    '''Return area of bbox in km² using spherical Earth approximation.
+    Handles bounding boxes that cross the antimeridian.
+    Args:
+      north (float): northern latitude in degrees
+      south (float): southern latitude in degrees
+      west (float): western longitude in degrees
+      east (float): eastern longitude in degrees
+    Returns:
+      Area in km² (float)
+    '''
+    west = norm_lon(west)
+    east = norm_lon(east)
+
+    if west > east:  # crosses ±180 -> split
+        return (bbox_area_km2(north,  west, south, 180.0) +
+                bbox_area_km2(north, -180, south, east))
+
+    lat_south, lat_north = map(math.radians, (south, north))
+    lon_west, lon_east = map(math.radians, (west,  east))
+    area = (EARTH_MEAN_RADIUS_KM**2) * abs(math.sin(lat_north) -
+                                           math.sin(lat_south)) * abs(lon_east - lon_west)
+    # print(
+    #     f"bbox_area_km2 computed: {north}, {west}, {south}, {east} -> {area} km²")
+    return area

bris-adapt/bris_adapt/orography/download.py

@@ -1,38 +1,146 @@
+import shutil
 import requests
 from typing import BinaryIO
+from .tiles import Tiler
+from .area import bbox_area_km2
+from .merge import merge_tiles
+from .ot import get_dataset
+import tempfile
+import os
+
 
 def download(
-  area_latlon: tuple[float|str, float|str, float|str, float|str],
-  dest_stream: BinaryIO,
-  api_key: str,
-  dem_type: str = "SRTMGL3"
+    area_latlon: tuple[float | str, float | str, float | str, float | str],
+    dest_stream: BinaryIO,
+    api_key: str,
+    high_res: bool = False
+) -> None:
+    """
+    Download topography data from OpenTopography Global DEM API.
+
+    Args:
+      area_latlon (tuple): Bounding box coordinates as (north, west, south, east).
+      dest_stream (BinaryIO): A writable binary file-like object to save the downloaded data.
+      high_res (bool): If True, use high-resolution DEMs (30m), otherwise use low-resolution (90m).
+    """
+    with tempfile.NamedTemporaryFile(prefix="orography_", suffix=".tif", delete=False) as temp_file:
+        filename = temp_file.name
+        download_to_file(
+            area_latlon=area_latlon,
+            filename=filename,
+            api_key=api_key,
+            high_res=high_res
+        )
+        os.close(temp_file.fileno())
+
+        with open(filename, "rb") as src_file:
+            shutil.copyfileobj(src_file, dest_stream)
+
+
+def download_to_file(
+    area_latlon: tuple[float | str, float | str, float | str, float | str],
+    filename: str,
+    api_key: str,
+    high_res: bool = False,
+    delete_temp: bool = True
+) -> str | None:
+    '''
+      Download topography data from OpenTopography Global DEM API.
+
+      Args:
+        area_latlon (tuple): Bounding box coordinates as (north, west, south, east).
+        dest_stream (BinaryIO or str): A writable binary file-like object to save the downloaded data.
+        high_res (bool): If True, use high-resolution DEMs (30m), otherwise use low-resolution (90m).
+        delete_temp (bool): If True, delete temporary files after merging. If false the temporary directory path is returned.
+
+     Returns:
+        The path to the temporary directory containing downloaded tiles if delete_temp is False, else None
+    '''
+    print(
+        f"Creating tiles for area ({area_latlon[0]}, {area_latlon[2]}, {area_latlon[1]}, {area_latlon[3]})")
+
+    dataset_info = get_dataset(
+        north=float(area_latlon[0]),
+        west=float(area_latlon[1]),
+        south=float(area_latlon[2]),
+        east=float(area_latlon[3]),
+        high_res=high_res
+    )
+    if dataset_info is None:
+        raise ValueError("No suitable dataset found.")
+
+    tiler = Tiler(
+        north=float(area_latlon[0]),
+        south=float(area_latlon[2]),
+        west=float(area_latlon[1]),
+        east=float(area_latlon[3]),
+        max_km2=dataset_info[1],
+        dlat=None,
+        dlon=None
+    )
+    tiles = tiler.create_area_limited_tiles()
+
+    temp_dir = tempfile.mkdtemp(prefix="dem_download_")
+    print(f"Created temporary directory: {temp_dir}")
+
+    for tile in tiles:
+        tile_filename = _create_temp_file(tile, temp_dir)
+        with open(tile_filename, "wb") as dest_stream:
+            _download(tile, dest_stream, api_key, dataset_info)
+
+    print(f"All tiles downloaded to temporary directory: {temp_dir}")
+    merge_tiles(tilesdir=temp_dir, pattern="*.tif", output=filename)
+    print(f"Merged DEM saved to: {filename}")
+    if delete_temp:
+        print(f"Deleting temporary directory: {temp_dir}")
+        if os.path.exists(temp_dir):
+            shutil.rmtree(temp_dir)
+        return None
+    return temp_dir
+
+
+def _create_temp_file(tile: tuple[float, float, float, float], dir_path: str, prefix: str = "dem_tile_", suffix: str = ".tif") -> str:
+    tile_filename = os.path.join(
+        dir_path,
+        f"{prefix}{int(tile[0])}_{int(tile[1])}_{int(tile[2])}_{int(tile[3])}{suffix}"
+    )
+    return tile_filename
+
+
+def _download(
+    area_latlon: tuple[float | str, float | str, float | str, float | str],
+    dest_stream: BinaryIO,
+    api_key: str,
+    dem_type: str = "SRTMGL3"
 ) -> None:
-  """
-  Download topography data from OpenTopography Global DEM API.
-
-  Args:
-    area_latlon (tuple): Bounding box coordinates as (north, west, south, east).
-    dest_stream (BinaryIO): A writable binary file-like object to save the downloaded data.
-    dem_type (str): DEM type (default: 'SRTMGL3').
-  """
-  url = "https://portal.opentopography.org/API/globaldem"
-  params = {
-    "demtype": dem_type,
-    "south": area_latlon[2],
-    "north": area_latlon[0],
-    "west": area_latlon[1],
-    "east": area_latlon[3],
-    "outputFormat": "GTiff",
-    "API_Key": api_key
-  }
-  response = requests.get(url, params=params, stream=True)
-  response.raise_for_status()
-  print(f"Downloading DEM of type {dem_type} for area {area_latlon}...")
-  print("Downloading...", end="", flush=True)
-  count = 0
-  for chunk in response.iter_content(chunk_size=8192):
-    count += 1
-    if count % 100 == 0:
-      print(".", end="", flush=True)
-    dest_stream.write(chunk)
-  print(f"\nDownloaded DEM to stream")
+    """
+    Download topography data from OpenTopography Global DEM API.
+
+    Args:
+      area_latlon (tuple): Bounding box coordinates as (north, south, west, east).
+      dest_stream (BinaryIO): A writable binary file-like object to save the downloaded data.
+      dem_type (str): DEM type (default: 'SRTMGL3').
+    """
+    url = "https://portal.opentopography.org/API/globaldem"
+    print(
+        f"Preparing download for area {area_latlon}...({area_latlon[0]}, {area_latlon[1]}, {area_latlon[2]}, {area_latlon[3]}) Area: {bbox_area_km2(*area_latlon):_.0f} km²")
+    params = {
+        "demtype": dem_type,
+        "south": area_latlon[2],
+        "north": area_latlon[0],
+        "west": area_latlon[1],
+        "east": area_latlon[3],
+        "outputFormat": "GTiff",
+        "API_Key": api_key
+    }
+    response = requests.get(url, params=params, stream=True)
+    response.raise_for_status()
+    print(f"Downloading DEM of type {dem_type} for area {area_latlon}...")
+    print("Downloading...", end="", flush=True)
+    count = 0
+    for chunk in response.iter_content(chunk_size=8192):
+        count += 1
+        if count % 100 == 0:
+            print(".", end="", flush=True)
+        dest_stream.write(chunk)
+    print(f"\nDownloaded DEM to stream")

bris-adapt/bris_adapt/orography/merge.py

@@ -0,0 +1,53 @@
+import rasterio
+from rasterio.merge import merge
+from glob import glob
+
+
+def merge_tiles(tilesdir: str, pattern: str = "*.tif", output: str = "merged_dem.tif", nodata: int | None = None) -> None:
+    '''
+    Merge all GeoTIFF tiles in 'tilesdir/pattern' into a single DEM 'merged_dem.tif'.
+    In case of overlapping areas, the value from the last file (alphabetically)
+    is used.
+
+    Args:
+        tilesdir (str): Directory containing the tiles.
+        pattern (str): Glob pattern to match tile files. Default is "*.tif".
+        output (str): Output merged DEM file path. Default is "merged_dem.tif".
+        nodata (int | None): NoData value for the DEM. If None, defaults to -32768.
+    '''
+
+    # order matters: later = higher priority
+    paths = sorted(glob(f"{tilesdir}/{pattern}"))
+    srcs = [rasterio.open(p) for p in paths]
+
+    if len(srcs) == 0:
+        print(f"No tiles found in '{tilesdir}' matching pattern '{pattern}'.")
+        return
+
+    if nodata is None:
+        nodata = srcs[0].nodata
+
+    mosaic, transform = merge(
+        srcs,
+        nodata=nodata,            # your DEM NoData
+        method="last",            # <<< last file wins in overlaps
+        resampling=rasterio.enums.Resampling.bilinear  # or .nearest
+    )
+
+    profile = srcs[0].profile
+    profile.update(
+        driver="GTiff",
+        height=mosaic.shape[1],
+        width=mosaic.shape[2],
+        transform=transform,
+        compress="lzw",
+        tiled=True
+    )
+
+    with rasterio.open(output, "w", **profile) as dst:
+        dst.write(mosaic)
+
+    for s in srcs:
+        s.close()
+
+    print(f"Merged {len(paths)} tiles into '{output}'")

bris-adapt/bris_adapt/orography/ot.py

@@ -0,0 +1,61 @@
+"""Orography datasets information and selection logic for opentopography.org.
+There are many other datasets available, but these are most interesting for our purposes.
+In fact we use only the datasets SRTMGL3 (90m), SRTMGL1 (30m), and COP30 (30m) and COP90 (90m).
+"""
+
+MAX_AREA_KM2_DEFAULTS = [
+    (('SRTMGL3', 'COP90'), 4_050_000),
+    (('GEDI_L3', 'SRTM15Plus', 'GEBCOIceTopo', 'GEBCOSubIceTopo', ''), 125_000_000)
+]
+
+MAX_AREA_KM2 = {
+    'SRTMGL3': 4_050_000,
+    'COP90': 4_050_000,
+    'GEDI_L3': 125_000_000,
+    'SRTM15Plus': 125_000_000,
+    'GEBCOIceTopo': 125_000_000,
+    'GEBCOSubIceTopo': 125_000_000
+}
+
+GLOBAL_DEMS_BONDING_BOX = {
+    'SRTMGL3': (60, -180, -56, 180),  # 90m
+    'SRTMGL1': (60, -180, -56, 180),  # 30m
+    'COP90': (84, -180, -85, 180),  # 90m
+    'COP30': (84, -180, -85, 180),  # 30m
+    'GEDI_L3': (52, -180, -52, 180),
+    'SRTM15Plus': (85, -180, -85, 180),
+    'GEBCOIceTopo': (85, -180, -85, 180),
+    'GEBCOSubIceTopo': (85, -180, -85, 180)
+}
+
+MAX_AREA_KM2_DEFAULT = 450_000
+
+
+def get_dataset(north: float, west: float, south: float, east: float, high_res: bool = False) -> tuple[str, int] | None:
+    """Get the suitable dataset for the given bounding box.
+
+    Args:
+        north (float): northern latitude in degrees
+        south (float): southern latitude in degrees
+        west (float): western longitude in degrees
+        east (float): eastern longitude in degrees
+        high_res (bool): if True, prefer high-resolution datasets (30m) over lower-resolution ones (90m)
+
+        Returns:
+            A tuple of (dataset_name, max_area_km2) if a suitable dataset is found, else None.
+    """
+    only_dems = []
+    if high_res:
+        only_dems = ['SRTMGL1', 'COP30']
+    else:
+        only_dems = ['SRTMGL3', 'COP90']
+
+    for dataset_name in only_dems:
+        bbox = GLOBAL_DEMS_BONDING_BOX.get(dataset_name)
+
+        bbox_north, bbox_west, bbox_south, bbox_east = bbox  # type: ignore
+        # Check if incoming bbox is fully contained in dataset bbox
+        if (north <= bbox_north and south >= bbox_south and
+                west >= bbox_west and east <= bbox_east):
+            return (dataset_name, MAX_AREA_KM2.get(dataset_name, MAX_AREA_KM2_DEFAULT))
+    return None