Fix matrix passing to C/C++ (#51)

* Fix matrix order

* Add test files

* Temporary CI branch rename

* Change order to 'F' in GridObject

* Add seed parameter to gen_random

* Add fillsinks tests

* Revise testing setup

* Add multiple python versions to test

* Fix python-version in matrix

* Add pytest import

* Change required python version

* Add first identifyflats test

* Change actions version to v3

* Add order='F' to read_tif

* Revert branch name to main

.github/workflows/ci.yaml

@@ -4,20 +4,29 @@ on:
   push:
     branches: ["main"]
 jobs:
-  build:
+  build-and-test:
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
+        python-version: ['3.10', '3.11']
+
     steps:
-      - name: Checkout package
-        uses: actions/checkout@v4
-      - name: Setup Python
-        uses: actions/setup-python@v5
+      - name: Checkout repository
+        uses: actions/checkout@v3
+
+      - name: Set up Python
+        uses: actions/setup-python@v3
         with:
-          python-version: '3.10'
-      - name: Install package
-        run: python -m pip install .[opensimplex]
-      - name: Test package
-        run: python -c "import topotoolbox as topo; dem = topo.gen_random(); assert (dem.fillsinks() >= dem).z.all()"
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install pytest
+          pip install .[opensimplex]
+
+      - name: Run tests
+        run: |
+          python -m pytest

.gitignore

@@ -14,4 +14,7 @@ docs/_build/*
 docs/_autosummary/*
 
 # ignore _temp
-docs/_temp
+docs/_temp
+
+# ignore pytest_cache
+.pytest_cache/

CMakeLists.txt

@@ -4,13 +4,13 @@ project(${SKBUILD_PROJECT_NAME}
   VERSION ${SKBUILD_PROJECT_VERSION}
   LANGUAGES CXX)
 
-# Find TopoToolbox somewhere
+# include libtopotoolbox
 include(FetchContent)
 FetchContent_Declare(
   topotoolbox
   GIT_REPOSITORY https://github.com/TopoToolbox/libtopotoolbox.git
-  GIT_TAG main # In the future, we should track a specific tag/commit
-	       # and bump it as we release versions
+  # GIT_TAG main
+  GIT_TAG 2024-W27
 )
 FetchContent_MakeAvailable(topotoolbox)
 

pyproject.toml

@@ -19,7 +19,7 @@ dependencies = [
     "matplotlib",
     "rasterio"
 ]
-requires-python = ">=3.7"
+requires-python = ">=3.10"
 keywords = ["TopoToolbox"]
 classifiers = []
 

src/topotoolbox/grid_object.py

@@ -22,11 +22,11 @@ def __init__(self) -> None:
         """
         # path to file
         self.path = ''
-        # name of dem
+        # name of DEM
         self.name = ''
 
         # raster metadata
-        self.z = np.empty(())
+        self.z = np.empty((), order='F')
         self.rows = 0
         self.columns = 0
         self.shape = self.z.shape
@@ -47,8 +47,7 @@ def fillsinks(self):
             The filled DEM.
         """
 
-        dem = self.z.astype(np.float32)
-
+        dem = self.z.astype(np.float32, order='F')
         output = np.zeros_like(dem)
 
         grid_fillsinks(output, dem, self.rows, self.columns)
@@ -67,7 +66,8 @@ def identifyflats(
         raw : bool, optional
             If True, returns the raw output grid as np.ndarray. 
             Defaults to False.
-        output : list of str, optional
+        output : list of str,
+                flat_neighbors = 0 optional
             List of strings indicating desired output types. Possible values 
             are 'sills', 'flats'. Defaults to ['sills', 'flats'].
 
@@ -87,8 +87,8 @@ def identifyflats(
         if output is None:
             output = ['sills', 'flats']
 
-        dem = self.z.astype(np.float32)
-        output_grid = np.zeros_like(dem).astype(np.int32)
+        dem = self.z.astype(np.float32, order='F')
+        output_grid = np.zeros_like(dem, dtype=np.int32)
 
         grid_identifyflats(output_grid, dem, self.rows, self.columns)
 
@@ -98,13 +98,13 @@ def identifyflats(
         result = []
         if 'flats' in output:
             flats = copy.copy(self)
-            flats.z = np.zeros_like(flats.z)
+            flats.z = np.zeros_like(flats.z, order='F')
             flats.z = np.where((output_grid & 1) == 1, 1, flats.z)
             result.append(flats)
 
         if 'sills' in output:
             sills = copy.copy(self)
-            sills.z = np.zeros_like(sills.z)
+            sills.z = np.zeros_like(sills.z, order='F')
             sills.z = np.where((output_grid & 2) == 2, 1, sills.z)
             result.append(sills)
 
@@ -122,11 +122,16 @@ def info(self):
         print(f"transform: {self.transform}")
         print(f"crs: {self.crs}")
 
-    def show(self):
+    def show(self, cmap='terrain'):
         """
         Display the GridObject instance as an image using Matplotlib.
+
+        Parameters
+        ----------
+        cmap : str, optional
+            Matplotlib colormap that will be used in the plot. 
         """
-        plt.imshow(self, cmap='terrain')
+        plt.imshow(self, cmap=cmap)
         plt.title(self.name)
         plt.colorbar()
         plt.tight_layout()

src/topotoolbox/utils.py

@@ -59,7 +59,7 @@ def write_tif(dem: GridObject, path: str) -> None:
         dataset.write(dem.z, 1)
 
 
-def show(*grid: GridObject, dpi: int = 100):
+def show(*grid: GridObject, dpi: int = 100, cmap: str = 'terrain'):
     """
     Display one or more GridObject instances using Matplotlib.
 
@@ -70,6 +70,8 @@ def show(*grid: GridObject, dpi: int = 100):
         should have an attribute `name` and be suitable for use with `imshow`.
     dpi : int, optional
         The resolution of the plots in dots per inch. Default is 100.
+    cmap : str, optional
+        Matplotlib colormap that will be used in the plot. 
 
     Notes
     -----
@@ -89,7 +91,7 @@ def show(*grid: GridObject, dpi: int = 100):
     fig, axes = plt.subplots(1, num_grids, figsize=(5*num_grids, 5), dpi=dpi)
     for i, dem in enumerate(grid):
         ax = axes[i] if num_grids > 1 else axes
-        im = ax.imshow(dem, cmap="terrain")
+        im = ax.imshow(dem, cmap=cmap)
         ax.set_title(dem.name)
         fig.colorbar(im, ax=ax, orientation='vertical')
 
@@ -125,7 +127,7 @@ def read_tif(path: str) -> GridObject:
         grid.path = path
         grid.name = os.path.splitext(os.path.basename(grid.path))[0]
 
-        grid.z = dataset.read(1).astype(np.float32)
+        grid.z = dataset.read(1).astype(np.float32, order='F')
         grid.rows = dataset.height
         grid.columns = dataset.width
         grid.shape = grid.z.shape
@@ -139,7 +141,7 @@ def read_tif(path: str) -> GridObject:
 
 
 def gen_random(hillsize: int = 24, rows: int = 128, columns: int = 128,
-               cellsize: float = 10.0) -> 'GridObject':
+               cellsize: float = 10.0, seed: int = 3) -> 'GridObject':
     """Generate a GridObject instance that is generated with OpenSimplex noise.
 
     Parameters
@@ -152,6 +154,8 @@ def gen_random(hillsize: int = 24, rows: int = 128, columns: int = 128,
         Number of columns. Defaults to 128.
     cellsize : float, optional
         Size of each cell in the grid. Defaults to 10.0.
+    seed : int, optional
+        Seed for the terrain generation. Defaults to 3
 
     Raises
     ------
@@ -171,7 +175,9 @@ def gen_random(hillsize: int = 24, rows: int = 128, columns: int = 128,
                "box[opensimplex]\" or \"pip install .[opensimplex]\"")
         raise ImportError(err) from None
 
-    noise_array = np.empty((rows, columns), dtype=np.float32)
+    noise_array = np.empty((rows, columns), dtype=np.float32, order='F')
+
+    simplex.seed(seed)
     for y in range(0, rows):
         for x in range(0, columns):
             value = simplex.noise4(x / hillsize, y / hillsize, 0.0, 0.0)
@@ -180,7 +186,6 @@ def gen_random(hillsize: int = 24, rows: int = 128, columns: int = 128,
 
     grid = GridObject()
 
-    grid.path = ''
     grid.z = noise_array
     grid.rows = rows
     grid.columns = columns

tests/test_grid_object.py

@@ -0,0 +1,85 @@
+import numpy as np
+import pytest
+
+import topotoolbox as topo
+
+
+@pytest.fixture
+def square_dem():
+    return topo.gen_random(rows=128, columns=128, seed=12)
+
+
+@pytest.fixture
+def wide_dem():
+    return topo.gen_random(rows=64, columns=128, seed=12)
+
+
+@pytest.fixture
+def tall_dem():
+    return topo.gen_random(rows=128, columns=64, seed=12)
+
+
+def test_fillsinks(square_dem, wide_dem, tall_dem):
+    for grid in [square_dem, wide_dem, tall_dem]:
+        # since grid is a fixture, it has to be assigned/called first
+        dem = grid
+        filled_dem = dem.fillsinks()
+
+        # Loop over all cells of the DEM
+        for i in range(dem.shape[0]):
+            for j in range(dem.shape[1]):
+
+                # Test: no filled cell lower than before calling fillsinks
+                assert dem[i, j] <= filled_dem[i, j]
+
+                # Test: cell isn't a sink
+                sink = 0
+                for i_offset, j_offset in [
+                        (-1, -1),
+                        (-1, 0),
+                        (-1, 1),
+                        (0, -1),
+                        (0, 1),
+                        (1, -1),
+                        (1, 0),
+                        (1, 1)]:
+
+                    i_neighbor = i + i_offset
+                    j_neighbor = j + j_offset
+
+                    if (i_neighbor < 0 or i_neighbor >= dem.z.shape[0]
+                            or j_neighbor < 0 or j_neighbor >= dem.z.shape[1]):
+                        continue
+
+                    if filled_dem[i_neighbor, j_neighbor] > filled_dem[i, j]:
+                        sink += 1
+
+                assert sink < 8
+
+
+def test_identifyflats(square_dem, wide_dem, tall_dem):
+    for dem in [square_dem, wide_dem, tall_dem]:
+        sills, flats = dem.identifyflats()
+
+        for i in range(dem.shape[0]):
+            for j in range(dem.shape[1]):
+
+                for i_offset, j_offset in [
+                        (-1, -1),
+                        (-1, 0),
+                        (-1, 1),
+                        (0, -1),
+                        (0, 1),
+                        (1, -1),
+                        (1, 0),
+                        (1, 1)]:
+
+                    i_neighbor = i + i_offset
+                    j_neighbor = j + j_offset
+
+                    if (i_neighbor < 0 or i_neighbor >= dem.z.shape[0]
+                            or j_neighbor < 0 or j_neighbor >= dem.z.shape[1]):
+                        continue
+
+                    if flats[i_neighbor, j_neighbor] < flats[i, j]:
+                        assert flats[i, j] == 1.0

tests/test_utils.py

@@ -0,0 +1,2 @@
+import pytest
+import topotoolbox.utils