MAINT: Switch from gdal to rasterio

.github/workflows/code-checks.yml

@@ -33,8 +33,8 @@ jobs:
 
         - name: Install dependencies
           run: |
-            sudo apt-get update -y
-            sudo apt-get install build-essential make openjdk-11-jre-headless
+            sudo apt-get update
+            sudo apt-get install build-essential openjdk-11-jre-headless
             python -m venv test
             . test/bin/activate
             python -m pip install --upgrade pip

.github/workflows/doc.yml

@@ -33,14 +33,11 @@ jobs:
 
         - name: Install dependencies
           run: |
-            sudo apt-add-repository ppa:ubuntugis/ubuntugis-unstable
             sudo apt-get update
-            sudo apt-get install build-essential make gcc gfortran gdal-bin libgdal-dev
+            sudo apt-get install build-essential gfortran
             python -m venv test
             . test/bin/activate
             python -m pip install --upgrade pip
-            pip install wheel numpy
-            pip install gdal==$(gdal-config --version)
             pip install -r requirements-dev.txt
 
         - name: Build

.github/workflows/unit-tests.yml

@@ -33,14 +33,11 @@ jobs:
 
         - name: Install dependencies
           run: |
-            sudo apt-add-repository ppa:ubuntugis/ubuntugis-unstable
             sudo apt-get update
-            sudo apt-get install build-essential make gcc gfortran gdal-bin libgdal-dev
+            sudo apt-get install build-essential gfortran
             python -m venv test
             . test/bin/activate
             python -m pip install --upgrade pip
-            pip install wheel numpy
-            pip install gdal==$(gdal-config --version)
             pip install -r requirements-dev.txt
 
         - name: Build

doc/source/release/1.0.0-notes.rst

@@ -59,8 +59,11 @@ This release was made possible thanks to the contributions of:
 Compatibilities
 ---------------
 
-- The `wheel <https://wheel.readthedocs.io/en/stable/>`__ package has been added to the list of dependencies
-  to build GDAL.
+- the `gdal <https://gdal.org/api/python_bindings.html>`__ package has been replaced by
+  `rasterio <https://rasterio.readthedocs.io/en/stable/>`__. This change has been made purely to improve
+  the distribution of `smash`. Despite the better performance in terms of raster reading time, gdal does
+  not provide a binary wheel directly, unlike rasterio
+  (see `discussion <https://github.com/OSGeo/gdal/issues/3060>`__).
 
 - The `pytest-cov <https://pytest-cov.readthedocs.io/en/latest/>`__ plugin has been added to the list of
   dependencies to perform test coverage.

environment-dev.yml

@@ -11,10 +11,9 @@ dependencies:
     - gfortran
 
     # required dependencies
-    - wheel
     - numpy>=1.13
     - f90wrap
-    - gdal
+    - rasterio
     - pandas
     - h5py
     - tqdm

environment.yml

@@ -11,10 +11,9 @@ dependencies:
     - gfortran
 
     # required dependencies
-    - wheel
     - numpy>=1.13
     - f90wrap
-    - gdal
+    - rasterio
     - pandas
     - h5py
     - tqdm

requirements.txt

@@ -1,8 +1,7 @@
 # required dependencies
-wheel
 numpy >= 1.13
 f90wrap
-gdal
+rasterio
 pandas
 h5py
 tqdm

smash/core/model/_build_model.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import pandas as pd
+import rasterio
 
 from f90wrap.runtime import FortranDerivedType
 from smash._constant import (
@@ -84,24 +85,25 @@ def _build_input_data(setup: SetupDT, mesh: MeshDT, input_data: Input_DataDT):
     if setup.read_qobs:
         _read_qobs(setup, mesh, input_data)
 
-    if setup.read_prcp:
-        _read_prcp(setup, mesh, input_data)
+    with rasterio.Env():
+        if setup.read_prcp:
+            _read_prcp(setup, mesh, input_data)
 
-    if setup.read_pet:
-        _read_pet(setup, mesh, input_data)
+        if setup.read_pet:
+            _read_pet(setup, mesh, input_data)
 
-    if setup.read_snow:
-        _read_snow(setup, mesh, input_data)
+        if setup.read_snow:
+            _read_snow(setup, mesh, input_data)
 
-    if setup.read_temp:
-        _read_temp(setup, mesh, input_data)
+        if setup.read_temp:
+            _read_temp(setup, mesh, input_data)
+
+        if setup.read_descriptor:
+            _read_descriptor(setup, mesh, input_data)
 
     if setup.prcp_partitioning:
         wrap_compute_prcp_partitioning(setup, mesh, input_data)  # % Fortran subroutine
 
-    if setup.read_descriptor:
-        _read_descriptor(setup, mesh, input_data)
-
     if setup.compute_mean_atmos:
         print("</> Computing mean atmospheric data")
         wrap_compute_mean_atmos(setup, mesh, input_data)  # % Fortran subroutine

smash/core/model/_read_input_data.py

@@ -1,13 +1,14 @@
 from __future__ import annotations
 
 import glob
+import os
 import re
 import warnings
 from typing import TYPE_CHECKING
 
 import numpy as np
 import pandas as pd
-from osgeo import gdal
+import rasterio
 from tqdm import tqdm
 
 from smash._constant import RATIO_PET_HOURLY
@@ -47,7 +48,7 @@ def _find_index_files_containing_date(
     ind = -1
     regex_pattern = _get_date_regex_pattern(dt, daily_interannual)
     for i, f in enumerate(files):
-        re_match = re.search(regex_pattern, f)
+        re_match = re.search(regex_pattern, os.path.basename(f))
         if daily_interannual:
             fdate = pd.Timestamp(f"{date.year}{re_match.group()}")
         else:
@@ -85,7 +86,7 @@ def _get_atmos_files(
         # % Adjust list by removing files that are ealier than start_time
         regex_pattern = _get_date_regex_pattern(date_range.freq.n, daily_interannual)
         for i, f in enumerate(files):
-            re_match = re.search(regex_pattern, f)
+            re_match = re.search(regex_pattern, os.path.basename(f))
             fdate = pd.Timestamp(re_match.group())
             if fdate >= date_range[0]:
                 return files[i:]
@@ -96,87 +97,88 @@ def _get_atmos_files(
 
 
 def _read_windowed_raster(path: FilePath, mesh: MeshDT) -> tuple[np.ndarray, dict[str, int]]:
-    gdal.UseExceptions()
-
     warning = {"res": 0, "overlap": 0, "outofbound": 0}
 
     # % Get raster information
-    ds = gdal.Open(path)
-    band = ds.GetRasterBand(1)
-    transform = ds.GetGeoTransform()
-    xmin = transform[0]
-    ymax = transform[3]
-    xres = transform[1]
-    yres = -transform[5]
-
-    # % Manage absolute tolerance wrt to xres or yres value
-    atol = 1e-2
-    xatol = atol * 10 ** min(0, np.floor(np.log10(np.abs(xres))))
-    yatol = atol * 10 ** min(0, np.floor(np.log10(np.abs(yres))))
-
-    # % Resolution missmatch
-    if not np.isclose(mesh.xres, xres, atol=xatol) or not np.isclose(mesh.yres, yres, atol=yatol):
-        warning["res"] = 1
-
-    # % Overlap missmatch
-    dxmin = mesh.xmin - xmin
-    dymax = ymax - mesh.ymax
-    xol_match = np.abs(dxmin / xres - np.round(dxmin / xres))
-    yol_match = np.abs(dymax / yres - np.round(dymax / yres))
-    if not np.isclose(xol_match, 0, atol=xatol) or not np.isclose(yol_match, 0, atol=yatol):
-        warning["overlap"] = 1
-
-    # % Allocate buffer
-    arr = np.zeros(shape=(mesh.nrow, mesh.ncol), dtype=np.float32)
-    arr.fill(np.float32(-99))
-
-    # % Pad offset to the nearest integer
-    xoff = np.rint((mesh.xmin - xmin) / xres)
-    yoff = np.rint((ymax - mesh.ymax) / yres)
-
-    # % Resolution ratio
-    xres_ratio = mesh.xres / xres
-    yres_ratio = mesh.yres / yres
-
-    # Reading window
-    win_xsize = np.rint(mesh.ncol * xres_ratio)
-    win_ysize = np.rint(mesh.nrow * yres_ratio)
-
-    # % Totally out of bound
-    # % Return the buffer with no data
-    if xoff >= band.XSize or yoff >= band.YSize or xoff + win_xsize <= 0 or yoff + win_ysize <= 0:
-        warning["outofbound"] = 2
-        return arr, warning
-
-    # % Partially out of bound
-    if xoff < 0 or yoff < 0 or xoff + win_xsize > band.XSize or yoff + win_ysize > band.YSize:
-        warning["outofbound"] = 1
-
-    # % Readjust offset
-    pxoff = max(0, xoff)
-    pyoff = max(0, yoff)
-
-    # % Readjust reading window
-    pwin_xsize = min(win_xsize + min(0, xoff), band.XSize - pxoff)
-    pwin_ysize = min(win_ysize + min(0, yoff), band.YSize - pyoff)
-
-    # % Buffer slices
-    xs = pxoff - xoff
-    ys = pyoff - yoff
-    xe = xs + pwin_xsize
-    ye = ys + pwin_ysize
-    arr_xslice = slice(int(xs * 1 / xres_ratio), int(xe * 1 / xres_ratio))
-    arr_yslice = slice(int(ys * 1 / yres_ratio), int(ye * 1 / yres_ratio))
-
-    # % Reading and writting into buffer
-    band.ReadAsArray(pxoff, pyoff, pwin_xsize, pwin_ysize, buf_obj=arr[arr_yslice, arr_xslice])
+    with rasterio.open(path) as ds:
+        transform = ds.get_transform()
+        xmin = transform[0]
+        ymax = transform[3]
+        xres = transform[1]
+        yres = -transform[5]
+
+        # % Manage absolute tolerance wrt to xres or yres value
+        atol = 1e-2
+        xatol = atol * 10 ** min(0, np.floor(np.log10(np.abs(xres))))
+        yatol = atol * 10 ** min(0, np.floor(np.log10(np.abs(yres))))
+
+        # % Resolution missmatch
+        if not np.isclose(mesh.xres, xres, atol=xatol) or not np.isclose(mesh.yres, yres, atol=yatol):
+            warning["res"] = 1
+
+        # % Overlap missmatch
+        dxmin = mesh.xmin - xmin
+        dymax = ymax - mesh.ymax
+        xol_match = np.abs(dxmin / xres - np.round(dxmin / xres))
+        yol_match = np.abs(dymax / yres - np.round(dymax / yres))
+        if not np.isclose(xol_match, 0, atol=xatol) or not np.isclose(yol_match, 0, atol=yatol):
+            warning["overlap"] = 1
+
+        # # % Allocate buffer
+        arr = np.zeros(shape=(mesh.nrow, mesh.ncol), dtype=np.float32)
+        arr.fill(np.float32(-99))
+
+        # % Pad offset to the nearest integer
+        xoff = np.rint((mesh.xmin - xmin) / xres)
+        yoff = np.rint((ymax - mesh.ymax) / yres)
+
+        # % Resolution ratio
+        xres_ratio = mesh.xres / xres
+        yres_ratio = mesh.yres / yres
+
+        # Reading window
+        win_xsize = np.rint(mesh.ncol * xres_ratio)
+        win_ysize = np.rint(mesh.nrow * yres_ratio)
+
+        # % Totally out of bound
+        # % Return the buffer with no data
+        if xoff >= ds.width or yoff >= ds.height or xoff + win_xsize <= 0 or yoff + win_ysize <= 0:
+            warning["outofbound"] = 2
+            return arr, warning
+
+        # % Partially out of bound
+        if xoff < 0 or yoff < 0 or xoff + win_xsize > ds.width or yoff + win_ysize > ds.height:
+            warning["outofbound"] = 1
+
+        # % Readjust offset
+        pxoff = max(0, xoff)
+        pyoff = max(0, yoff)
+
+        # % Readjust reading window
+        pwin_xsize = min(win_xsize + min(0, xoff), ds.width - pxoff)
+        pwin_ysize = min(win_ysize + min(0, yoff), ds.height - pyoff)
+
+        # % Buffer slices
+        xs = pxoff - xoff
+        ys = pyoff - yoff
+        xe = xs + pwin_xsize
+        ye = ys + pwin_ysize
+        arr_xslice = slice(int(xs * 1 / xres_ratio), int(xe * 1 / xres_ratio))
+        arr_yslice = slice(int(ys * 1 / yres_ratio), int(ye * 1 / yres_ratio))
+
+        # % Reading and writting into buffer
+        ds.read(
+            1,
+            window=rasterio.windows.Window(pxoff, pyoff, pwin_xsize, pwin_ysize),
+            out=arr[arr_yslice, arr_xslice],
+        )
 
-    # % Manage NoData
-    nodata = band.GetNoDataValue()
-    if nodata is not None:
-        arr = np.where(arr == nodata, -99.0, arr)
+        # % Manage NoData
+        nodata = ds.nodata
+        if nodata is not None:
+            arr = np.where(arr == nodata, -99.0, arr)
 
-    return arr, warning
+        return arr, warning
 
 
 def _get_reading_warning_message(reading_warning: dict[str, int | list[str]]) -> str:

smash/factory/mesh/_standardize.py

@@ -6,7 +6,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from osgeo import gdal
+import rasterio
 
 from smash.factory.mesh._tools import _get_transform
 
@@ -28,7 +28,7 @@ def _standardize_generate_mesh_flwdir_path(flwdir_path: FilePath) -> str:
     return flwdir_path
 
 
-def _standardize_generate_mesh_bbox(flwdir_dataset: gdal.Dataset, bbox: ListLike) -> np.ndarray:
+def _standardize_generate_mesh_bbox(flwdir_dataset: rasterio.DatasetReader, bbox: ListLike) -> np.ndarray:
     # % Bounding Box (xmin, xmax, ymin, ymax)
 
     if not isinstance(bbox, (list, tuple, np.ndarray)):
@@ -91,7 +91,7 @@ def _standardize_generate_mesh_bbox(flwdir_dataset: gdal.Dataset, bbox: ListLike
 
 
 def _standardize_generate_mesh_x_y_area(
-    flwdir_dataset: gdal.Dataset,
+    flwdir_dataset: rasterio.DatasetReader,
     x: Numeric | ListLike,
     y: Numeric | ListLike,
     area: Numeric | ListLike,
@@ -183,11 +183,9 @@ def _standardize_generate_mesh_args(
     max_depth: Numeric,
     epsg: AlphaNumeric | None,
 ) -> AnyTuple:
-    gdal.UseExceptions()
-
     flwdir_path = _standardize_generate_mesh_flwdir_path(flwdir_path)
 
-    flwdir_dataset = gdal.Open(flwdir_path)
+    flwdir_dataset = rasterio.open(flwdir_path)
 
     if x is None and bbox is None:
         raise ValueError("bbox argument or (x, y, area) arguments must be defined")