Add scripts to help build the hybrid elevation map

prepare_layers/README.md

@@ -0,0 +1,13 @@
+# Preparing input rasters for STAR
+
+## Elevation
+
+STAR uses [FABDEM v1.2](https://data.bris.ac.uk/data/dataset/s5hqmjcdj8yo2ibzi9b4ew3sn) for its elevation layer, however this raster needs to be processed before it can be used:
+
+* There are some errata tiles to correct mistakes in the area of Florida's Forgotten Coast. Unfortunately these are only available at time of writing via a [Google Drive link](https://drive.google.com/file/d/1DIAaheKT-thuWPhzXWpfVnTqRn13nqvi/view?usp=sharing).
+* FABDEM was created when certain tiles in the Copernicus GLO DEM were not available, leaving a gap around Azerbaijan and nearby countries.
+
+FABDEM itself is very large and slow to download, and so we leave that as an exercise for the user rather than automating it as part of the pipeline. Once that has downloaded and the google drive link has been followed and the tiles expanded, we provide two scripts to complete the job:
+
+* `fetch_cglo.py` - this will fetch the missing tiles from Copernicus GLO that are now available.
+* `make_hybrid_elevation_map.py` - this takes three inputs: a folder with the FABDEM v1.2 tiles, a folder with the errata files from Google Drive, and a folder with the additional CGLO tiles, and outputs the compiled hybrid elevation map. Note that this will be around 500 GB in size.

prepare_layers/fetch_cglo.py

@@ -0,0 +1,81 @@
+# STAR uses FABDEM, which is based on the Copernicus GLO 30 Digital Elevation Model (CGLO), but
+# was built at a time when CGLO was missing certain areas of the globe. This script downloads those
+# extra areas between 40.7371 to 50.9321 degrees longitude and 37.9296 to 45.7696 latitude.
+import argparse
+from pathlib import Path
+
+import boto3
+from botocore import UNSIGNED
+from botocore.config import Config
+
+def download_copernicus_dem_tiles(
+    min_lon: float,
+    max_lon: float,
+    min_lat: float,
+    max_lat: float,
+    output_dir: Path,
+) -> tuple[list[str],list[str]]:
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    s3 = boto3.client(
+        's3',
+        endpoint_url='https://opentopography.s3.sdsc.edu',
+        config=Config(signature_version=UNSIGNED)
+    )
+
+    bucket = 'raster'
+
+    lon_tiles = range(int(min_lon), int(max_lon) + 1)
+    lat_tiles = range(int(min_lat), int(max_lat) + 1)
+
+    downloaded = []
+    failed = []
+
+    for lat in lat_tiles:
+        for lon in lon_tiles:
+            lat_prefix = 'N' if lat >= 0 else 'S'
+            lon_prefix = 'E' if lon >= 0 else 'W'
+            lat_str = f"{abs(lat):02d}"
+            lon_str = f"{abs(lon):03d}"
+
+            tile_name = f"Copernicus_DSM_10_{lat_prefix}{lat_str}_00_{lon_prefix}{lon_str}_00_DEM"
+            s3_key = f"COP30/COP30_hh/{tile_name}.tif"
+            local_path = f"{output_dir}/{tile_name}.tif"
+
+            try:
+                print(f"Downloading {tile_name}.tif...")
+                s3.download_file(bucket, s3_key, local_path)
+                downloaded.append(tile_name)
+                print(f"  ✓ Saved to {local_path}")
+            except Exception as e: # pylint: disable=W0718
+                print(f"  ✗ Failed: {e}")
+                failed.append(tile_name)
+
+    print(f"\n{'='*60}")
+    print(f"Downloaded: {len(downloaded)} tiles")
+    if failed:
+        print(f"Failed: {len(failed)} tiles")
+        print("Failed tiles:", failed)
+
+    return downloaded, failed
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Convert IUCN crosswalk to minimal common format.")
+    parser.add_argument(
+        '--output',
+        type=Path,
+        help='Destination folder for tiles',
+        required=True,
+        dest='output_dir',
+    )
+    args = parser.parse_args()
+
+    min_lon = 40.7371
+    max_lon = 50.9321
+    min_lat = 37.9296
+    max_lat = 45.7696
+
+    download_copernicus_dem_tiles(min_lon, max_lon, min_lat, max_lat, args.output_dir)
+
+if __name__ == "__main__":
+    main()

prepare_layers/make_hybrid_elevation_map.py

@@ -0,0 +1,90 @@
+import argparse
+from pathlib import Path
+from alive_progress import alive_bar # type: ignore
+
+import yirgacheffe as yg
+from yirgacheffe.layers import RescaledRasterLayer
+
+def resize_cglo(
+    projection: yg.MapProjection,
+    cglo_path: Path,
+    output_path: Path,
+) -> None:
+    with yg.read_rasters(list(cglo_path.glob("*.tif"))) as cglo_30:
+        rescaled = RescaledRasterLayer(cglo_30, projection, nearest_neighbour=False)
+        with alive_bar(manual=True) as bar:
+            rescaled.to_geotiff(output_path, parallelism=True, callback=bar)
+
+def make_hybrid_elevation_map(
+    fabdem_path: Path,
+    fabdem_patch_path: Path,
+    cglo_path: Path,
+    output_path: Path,
+) -> None:
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+
+    tmpdir = output_path.parent
+
+    # The CGLO files are at a different resolution to FABDEM, so we first
+    # need to scale them.
+    resized_cglo_path = tmpdir / "cglo.tif"
+    if not resized_cglo_path.exists():
+        # Get the map projection and pixel scale for fabdem
+        fabdem_example_tile = list(fabdem_path.glob("*.tif"))[0]
+        with yg.read_raster(fabdem_example_tile) as example:
+            projection = example.map_projection
+
+        resize_cglo(projection, cglo_path, resized_cglo_path)
+
+    # Now we build up a large group layer, and rely on the fact that
+    # Yirgacheffe will render earlier layers over later layers
+    file_list = list(fabdem_patch_path.glob("*.tif")) + \
+        list(fabdem_path.glob("*.tif")) + \
+        [resized_cglo_path]
+
+    full_layer = yg.read_rasters(file_list)
+
+    with alive_bar(manual=True) as bar:
+        full_layer.to_geotiff(output_path, parallelism=256, callback=bar)
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Convert IUCN crosswalk to minimal common format.")
+    parser.add_argument(
+        '--fabdem_tiles',
+        type=Path,
+        help="Directory containing original FABDEM tiles",
+        required=True,
+        dest="fabdem_path",
+    )
+    parser.add_argument(
+        '--fabdem_patch_tiles',
+        type=Path,
+        help="Directory containing original FABDEM errata tiles",
+        required=True,
+        dest="fabdem_patch_path",
+    )
+    parser.add_argument(
+        '--cglo_tiles',
+        type=Path,
+        help="Directory containing missing_cglo tiles",
+        required=True,
+        dest="cglo_path",
+    )
+    parser.add_argument(
+        '--output',
+        type=Path,
+        help="Final output raster",
+        required=True,
+        dest='output_path',
+    )
+    args = parser.parse_args()
+
+    make_hybrid_elevation_map(
+        args.fabdem_path,
+        args.fabdem_patch_path,
+        args.cglo_path,
+        args.output_path,
+    )
+
+if __name__ == "__main__":
+    main()

requirements.txt

@@ -7,6 +7,7 @@ pymer4
 pyproj
 scikit-image
 redlistapi
+boto3
 yirgacheffe>=1.9
 aoh[validation]>=1.0
 

scripts/run.sh

@@ -83,11 +83,11 @@ if [[ ! -f "${DATADIR}"/elevation/elevation-max-1k.tif || ! -f "${DATADIR}"/elev
     fi
     if [ ! -f "${DATADIR}"/elevation/elevation-max-1k.tif ]; then
         echo "Generating elevation max layer..."
-        gdalwarp -t_srs ESRI:54009 -tr 1000 -1000 -r max -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation/elevation.tif "${DATADIR}"/elevation/elevation-max-1k.tif
+        gdalwarp -t_srs ESRI:54009 -tr 1000 -1000 -r max -tap -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation/elevation.tif "${DATADIR}"/elevation/elevation-max-1k.tif
     fi
     if [ ! -f "${DATADIR}"/elevation/elevation-min-1k.tif ]; then
         echo "Generating elevation min layer..."
-        gdalwarp -t_srs ESRI:54009 -tr 1000 -1000 -r min -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation/elevation.tif "${DATADIR}"/elevation/elevation-min-1k.tif
+        gdalwarp -t_srs ESRI:54009 -tr 1000 -1000 -r min -tap -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation/elevation.tif "${DATADIR}"/elevation/elevation-min-1k.tif
     fi
 fi