Replace srtm4 dependency by SRTM DEM on GCS

[NBrown140]

Problem

srtm4 is a python library that downloads srtm files locally and then reads them to extract elevation data at specific latitudes and longitudes. There are two problems with srtm4:

1. The large srtm tiles are downloaded at every data-ingestion pipeline run. This takes time and is inefficient.
2. The srtm4 library requires special treatment to run locally on MacOS and complicates local development.

Solution

We already have a Copernicus DEM on GCS which dynamically builds up as new tiles are required. We could obtain elevation data by very quickly reading the elevation points from these cloud optimized geotifs.
I've downloaded the full SRTM v4.1 (90m) dataset at gs://overstory-dtms/srtm_v41_90m and built a vrt index. This can be used to very efficiently query for elevations at specific latitudes and longitudes.

To Do

• Download full or partial SRTM DEM dataset. Keep in a notebook/pipeline if we ever need to recreate or expand the dataset.
• Make a vrt file on GCS to link all the tiles together

Reproducibility

Here is a script demonstrating the difference between using SRTM ourselves instead of using srtm4 library:

from typing import List
import rasterio as rio
from rasterio.windows import Window
import srtm4

pyproj.network.set_network_enabled(active=True)
import os
os.environ["PROJ_UNSAFE_SSL"]="TRUE"


SRTM_VRT = "gs://overstory-dtms/srtm_v41_90m/index.vrt"
lons = [-72.2,  -102, -98, -115, -77]
lats = [45.4 , 50  , 35  , 33. , 37]

def get_srtm_elevations(lons: List[float], lats: List[float], convert_ellipsoidal: bool):
    with rio.open(SRTM_VRT) as src:
        assert src.crs == "EPSG:4326"
        elevations = [float(e[0]) for e in src.sample([(lon, lat) for lon, lat in zip(lons, lats)])]
        # Convert nodata
        elevations = [e if e != src.nodata else float("nan") for e in elevations]
        # Convert to ellipsoidal heights
        if convert_ellipsoidal:
            elevations = to_ellipsoid(lons, lats, elevations)
    return elevations

def to_ellipsoid(lons: List[float], lats: List[float], elevations: List[float]):
    """
    Convert geoidal heights to ellipsoidal heights.
    """
    # WGS84 with ellipsoid height as vertical axis
    ellipsoid = pyproj.CRS.from_epsg(4979)
    # WGS84 with Gravity-related height (EGM96)
    geoid = pyproj.CRS("EPSG:4326+5773")

    trf = pyproj.Transformer.from_crs(geoid, ellipsoid)
    new_elevations = trf.transform(lats, lons, elevations, errcheck=True)[-1]

    return [round(new_alt, 2) for new_alt in new_elevations]


old_elev = srtm4.srtm4(lons, lats)
print(old_elev)

new_elev = get_srtm_elevations(lons, lats, True)
print(new_elev)

Output is:

[266.777, 618.553, 342.51, 122.856, -5.45667]
[275.74, 618.55, 339.5, 123.86, -0.46]

The values are very similar, but not exactly the same. It's possible that srtm4 does a bit more than just converting to ellipsoidal height. Maybe it's related to the half-pixel offset due to the srtm geotifs having the wrong AREA_OR_POINT tiftag: https://github.com/centreborelli/srtm4/blob/master/srtm4/raster.py#L2-L6

[NBrown140]

I accidentally created a PR in the original rpcm repo: centreborelli#23
But I closed it immediately.

[NBrown140]

Remove data-ingestion dependencies on barely used libraries

[daviddemeij]

This is great, thanks 👍