Use rasterio to makde sardem pip-installable

README.md

@@ -24,6 +24,15 @@ pip install sardem
 ```
 This creates the command line executable `sardem`
 
+Using `uv`:
+```bash
+uv add sardem
+```
+You can also run the `sardem` command line tool using `uvx`:
+```bash
+uvx sardem --help
+```
+
 Alternatively, you can clone to build/install:
 
 ```bash
@@ -32,41 +41,23 @@ cd sardem
 # Install requirements using either pip or conda
 # conda install -c conda-forge --file environment.yml
 # pip install -r requirements.txt
-# the conda environment.yml is more complete, as GDAL is required for some of the functionality
+# the conda environment.yml is more complete, as rasterio is required for some of the functionality
 pip install -e .
 ```
 which will run `pip install --upgrade .` and create the command line script.
 
-
 ## Data sources
 The default data source, `--data-source NASA`, uses the SRTM 1 arcsecond data. You can also use the newer [Copernicus Digital Surface Model (DSM)](https://registry.opendata.aws/copernicus-dem/). 
 To see a comparison of the two, see the [srtm_copernicus_comparison](notebooks/srtm_copernicus_comparison.ipynb) notebook.
 
-**Note:** To convert the elevation values to heights about the WGS84 ellipsoid (which is the default), or to use the Copernicus data, **GDAL is required**. 
-For the Copernicus data, the minimum required GDAL version is 3.4.2; versions earlier than 3.4.0 seem to hang upon using `gdalwarp` on the global VRT, and <3.4.2 have an internal bug https://github.com/isce-framework/isce2/issues/556 .
-
+**Note:** To convert the elevation values to heights about the WGS84 ellipsoid (which is the default), or to use the Copernicus data, **rasterio is required**.
+For the Copernicus data, rasterio >= 1.0.0 is required for coordinate transformations and data warping.
 
 ## Bounding box convention
 
 `sardem` uses the gdal convention ("pixel is area") where `--bbox` points to the *edges* of the [left, bottom, right, top] pixels.
 I.e. (left, bottom) refers to the lower left corner of the lower left pixel.
 
-
-
-### Converting to WGS84 ellipsoidal heights from EGM96/EGM2008 geoid heights
-
-GDAL is required for the conversion, which is installed when using `conda install -c conda-forge sardem`.
-If you already are using an existing environment, make sure that the GDAL version is >=3.4.2.
-
-```bash 
-conda install -c conda-forge "gdal>=3.4.2"
-
-# or
-# conda install -c conda-forge mamba
-# mamba install -c conda-forge "gdal>=3.4.2"
-```
-
-
 ## Command Line Interface
 
 The full options for the command line tool in `sardem/cli.py` can be found using
@@ -83,7 +74,7 @@ Stiches SRTM .hgt files to make (upsampled) DEM
     the necessary SRTM1 tiles, stitch together, then upsample.
 
     The `--bbox` convention points to the *edges* of the [left, bottom, right, top]
-    pixels, following the "pixel is area" convention as used in gdal.
+    pixels, following the "pixel is area" convention as used in rasterio.
     I.e. (left, bottom) refers to the lower left corner of the lower left pixel.
 
     Usage Examples:

environment.yml

@@ -2,7 +2,7 @@ channels:
 - conda-forge
 dependencies:
 - python >= 3.6
-- gdal >= 3.4.2
+- rasterio >= 1.0.0
 - numpy >= 1.14.2
 - requests >= 2.20.0
 - pip

requirements.txt

@@ -1,2 +1,3 @@
 numpy>=1.14.2
 requests>=2.20.0
+rasterio>=1.0.0

sardem/cli.py

@@ -1,6 +1,7 @@
 """
 Command line interface for running sardem
 """
+
 import json
 from argparse import (
     ArgumentError,
@@ -29,7 +30,7 @@ def positive_small_int(argstring):
     the necessary SRTM1 tiles, stitch together, then upsample.
 
     The `--bbox` convention points to the *edges* of the [left, bottom, right, top]
-    pixels, following the "pixel is area" convention as used in gdal.
+    pixels, following the "pixel is area" convention as used in rasterio.
     I.e. (left, bottom) refers to the lower left corner of the lower left pixel.
 
     Usage Examples:
@@ -77,7 +78,7 @@ def get_cli_args():
         help="Bounding box of area of interest "
         " (e.g. --bbox -106.1 30.1 -103.1 33.1 ). \n"
         "--bbox points to the *edges* of the pixels, \n"
-        " following the 'pixel is area' convention as used in gdal. ",
+        " following the 'pixel is area' convention as used in rasterio. ",
     )
     parser.add_argument(
         "--geojson",

sardem/conversions.py

@@ -30,26 +30,41 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="eg
     # Target srs: WGS84 horizontal/vertical
     t_srs = '"epsg:4326"'
 
-    # Note: https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-tr
-    # If not specified, gdalwarp will generate an output raster with xsize=ysize
-    # We want it to match the input file
+    # Use rasterio for coordinate system transformation
+    import rasterio
+    from rasterio.warp import reproject, Resampling
+
     xsize, ysize = _get_size(filename)
-    cmd = (
-        'gdalwarp {overwrite} -s_srs {s_srs} -t_srs {t_srs}'
-        ' -of ROI_PAC -ts {xsize} {ysize} '
-        ' -multi -wo NUM_THREADS=4 -wm 4000 {inp} {out}'
-    )
-    cmd = cmd.format(
-        inp=filename,
-        out=output,
-        overwrite="-overwrite" if overwrite else "",
-        xsize=xsize,
-        ysize=ysize,
-        s_srs=s_srs,
-        t_srs=t_srs,
-    )
-    logger.info(cmd)
-    subprocess.run(cmd, check=True, shell=True)
+
+    with rasterio.open(filename) as src:
+        # Create output profile
+        profile = src.profile.copy()
+        profile.update(
+            {
+                "driver": "ENVI",  # ROI_PAC equivalent
+                "crs": t_srs.strip('"'),
+                "width": xsize,
+                "height": ysize,
+            }
+        )
+
+        if not overwrite and os.path.exists(output):
+            logger.warning("Output file {} exists and overwrite=False".format(output))
+            return output
+
+        logger.info("Converting {} from {} to {}".format(filename, s_srs, t_srs))
+        with rasterio.open(output, "w", **profile) as dst:
+            for i in range(1, src.count + 1):
+                reproject(
+                    source=rasterio.band(src, i),
+                    destination=rasterio.band(dst, i),
+                    src_transform=src.transform,
+                    src_crs=s_srs.strip('"'),
+                    dst_transform=dst.transform,
+                    dst_crs=t_srs.strip('"'),
+                    resampling=Resampling.nearest,
+                    num_threads=4,
+                )
 
     if copy_rsc:
         rsc_file = filename + ".rsc"
@@ -59,23 +74,19 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="eg
 
 
 def _get_size(filename):
-    """Retrieve the raster size from a gdal-readable file"""
-    from osgeo import gdal
+    """Retrieve the raster size from a rasterio-readable file"""
+    import rasterio
 
-    ds = gdal.Open(filename)
-    xsize, ysize = ds.RasterXSize, ds.RasterYSize
-    ds = None
+    with rasterio.open(filename) as ds:
+        xsize, ysize = ds.width, ds.height
     return xsize, ysize
 
 
 def convert_dem_to_wgs84(dem_filename, geoid="egm96"):
     """Convert the file `dem_filename` from EGM96 heights to WGS84 ellipsoidal heights
 
-    Overwrites file, requires GDAL to be installed
+    Overwrites file, uses rasterio for conversion
     """
-    if not utils._gdal_installed_correctly():
-        logger.error("GDAL required to convert DEM to WGS84")
-        return
 
     path_, fname = os.path.split(dem_filename)
     rsc_filename = os.path.join(path_, fname + ".rsc")
@@ -95,4 +106,4 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96"):
         logger.error("Failed to convert DEM:", exc_info=True)
         logger.error("Reverting back, using EGM dem as output")
         os.rename(output_egm, dem_filename)
-        os.rename(rsc_filename_egm, rsc_filename)
+        os.rename(rsc_filename_egm, rsc_filename)

sardem/cop_dem.py

@@ -21,7 +21,7 @@ def download_and_stitch(
     yrate=1,
     vrt_filename=None,
     output_format="ENVI",
-    output_type="int16"
+    output_type="int16",
 ):
     """Download the COP DEM from AWS.
 
@@ -32,9 +32,10 @@ def download_and_stitch(
         https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198
         https://copernicus-dem-30m.s3.amazonaws.com/readme.html
     """
-    from osgeo import gdal
+    import rasterio
+    import rasterio.warp
+    from rasterio.enums import Resampling
 
-    gdal.UseExceptions()
     # TODO: does downloading make it run any faster?
     # if download_vrt:
     #     cache_dir = utils.get_cache_dir()
@@ -56,74 +57,105 @@ def download_and_stitch(
     resamp = "bilinear" if (xrate > 1 or yrate > 1) else "nearest"
 
     # access_mode = "overwrite" if overwrite else None
-    option_dict = dict(
-        format=output_format,
-        outputBounds=bbox,
-        dstSRS=t_srs,
-        srcSRS=s_srs,
-        xRes=xres,
-        yRes=yres,
-        outputType=gdal.GetDataTypeByName(output_type.title()),
-        resampleAlg=resamp,
-        multithread=True,
-        warpMemoryLimit=5000,
-        warpOptions=["NUM_THREADS=4"],
-    )
-
-    # Used the __RETURN_OPTION_LIST__ to get the list of options for debugging
+    # Map resampling algorithm names
+    resamp_map = {
+        "bilinear": Resampling.bilinear,
+        "nearest": Resampling.nearest,
+        "cubic": Resampling.cubic,
+        "average": Resampling.average,
+    }
+    resample_alg = resamp_map.get(resamp, Resampling.bilinear)
+
     logger.info("Creating {}".format(output_name))
     logger.info("Fetching remote tiles...")
-    try:
-        cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_dict)
-        logger.info("Running GDAL command:")
-        logger.info(cmd)
-    except Exception:
-        # Can't form the cli version due to `deepcopy` Pickle error, just skip
-        logger.info("Running gdal.Warp with options:")
-        logger.info(option_dict)
-        pass
-    # Now convert to something GDAL can actually use
-    option_dict["callback"] = gdal.TermProgress
-    gdal.Warp(output_name, vrt_filename, options=gdal.WarpOptions(**option_dict))
-    return
-
 
-def _gdal_cmd_from_options(src, dst, option_dict):
-    from osgeo import gdal
-
-    opts = deepcopy(option_dict)
-    # To see what the list of cli options are (gdal >= 3.5.0)
-    opts["options"] = "__RETURN_OPTION_LIST__"
-    opt_list = gdal.WarpOptions(**opts)
-    out_opt_list = deepcopy(opt_list)
-    for idx, o in enumerate(opt_list):
-        # Wrap the srs option in quotes
-        if o.endswith("srs"):
-            out_opt_list[idx + 1] = '"{}"'.format(out_opt_list[idx + 1])
-    return "gdalwarp {} {} {}".format(src, dst, " ".join(out_opt_list))
+    # Use rasterio for warping
+    with rasterio.open(vrt_filename) as src:
+        # Calculate output transform and dimensions
+        dst_transform, dst_width, dst_height = (
+            rasterio.warp.calculate_default_transform(
+                src.crs if s_srs is None else s_srs,
+                t_srs if t_srs is not None else src.crs,
+                src.width,
+                src.height,
+                *src.bounds,
+                resolution=(xres, yres),
+            )
+        )
+
+        # If bbox is provided, override the transform calculation
+        if bbox:
+            left, bottom, right, top = bbox
+            dst_width = int(round((right - left) / xres))
+            dst_height = int(round((top - bottom) / abs(yres)))
+            dst_transform = rasterio.transform.from_bounds(
+                left, bottom, right, top, dst_width, dst_height
+            )
+
+        # Create output profile
+        profile = src.profile.copy()
+        profile.update(
+            {
+                "driver": "ENVI" if output_format == "ENVI" else output_format,
+                "height": dst_height,
+                "width": dst_width,
+                "transform": dst_transform,
+                "crs": t_srs if t_srs is not None else src.crs,
+                "dtype": output_type,
+            }
+        )
+
+        warp_opts = {"XSCALE": 1, "YSCALE": 1, "APPLY_VERTICAL_SHIFT": True}
+        with rasterio.open(output_name, "w", **profile) as dst:
+            for i in range(1, src.count + 1):
+                rasterio.warp.reproject(
+                    source=rasterio.band(src, i),
+                    destination=rasterio.band(dst, i),
+                    src_transform=src.transform,
+                    src_crs=src.crs if s_srs is None else s_srs,
+                    dst_transform=dst_transform,
+                    dst_crs=t_srs if t_srs is not None else src.crs,
+                    resampling=resample_alg,
+                    num_threads=4,
+                    **warp_opts,
+                )
+    return
 
 
 def make_cop_vrt(outname="copernicus_GLO_30_dem.vrt"):
     """Build a VRT from the Copernicus GLO 30m DEM COG dataset
 
     Note: this is a large VRT file, ~15MB, so it can many hours to build.
     """
-    from osgeo import gdal
-
-    gdal.UseExceptions()
+    import rasterio
+    from rasterio.vrt import WarpedVRT
 
     tile_list = get_tile_list()
     url_list = _make_url_list(tile_list)
-    vrt_options = gdal.BuildVRTOptions(
-        resampleAlg=gdal.GRIORA_NearestNeighbour,
-        outputBounds=[-180, -90, 180, 90],
-        resolution="highest",
-        outputSRS="EPSG:4326+3855",
-    )
+
     logger.info("Building VRT {}".format(outname))
-    vrt_file = gdal.BuildVRT(outname, url_list, options=vrt_options)
-    vrt_file.FlushCache()
-    vrt_file = None
+
+    # Create a simple VRT by writing XML directly (rasterio doesn't have BuildVRT equivalent)
+    vrt_xml = '<?xml version="1.0" encoding="UTF-8"?>\n'
+    vrt_xml += '<VRTDataset rasterXSize="1296000" rasterYSize="417600">\n'
+    vrt_xml += "  <SRS>EPSG:4326+3855</SRS>\n"
+    vrt_xml += "  <GeoTransform>-1.800000000000000e+02, 2.777777777777778e-04, 0.000000000000000e+00, 8.333333333333334e+01, 0.000000000000000e+00, -2.777777777777778e-04</GeoTransform>\n"
+    vrt_xml += '  <VRTRasterBand dataType="Int16" band="1">\n'
+    vrt_xml += "    <NoDataValue>-32768</NoDataValue>\n"
+
+    for url in url_list:
+        vrt_xml += f"    <SimpleSource>\n"
+        vrt_xml += f'      <SourceFilename relativeToVRT="0">{url}</SourceFilename>\n'
+        vrt_xml += f"      <SourceBand>1</SourceBand>\n"
+        vrt_xml += f"    </SimpleSource>\n"
+
+    vrt_xml += "  </VRTRasterBand>\n"
+    vrt_xml += "</VRTDataset>\n"
+
+    with open(outname, "w") as f:
+        f.write(vrt_xml)
+
+    logger.info(f"VRT written to {outname} with {len(url_list)} tiles")
 
 
 def get_tile_list():