fix overflow errors numpy v2

pyflwdir/basins.py

@@ -89,11 +89,11 @@ def subbasins_streamorder(idxs_ds, seq, strord, mask=None, min_sto=-2):
         map with unique IDs for stream_order>=min_sto subbasins
     """
     if min_sto < 0:
-        min_sto = strord.max() + min_sto
+        min_sto = int(strord.max()) + min_sto
     subbas = np.full(idxs_ds.shape, 0, dtype=np.int32)
     idxs = []
     for idx0 in seq[::-1]:  # up- to downstream
-        if (mask is not None and mask[idx0] == False) or strord[idx0] < min_sto:
+        if (mask is not None and mask[idx0] is False) or strord[idx0] < min_sto:
             continue
         idx_ds = idxs_ds[idx0]
         if strord[idx0] != strord[idx_ds] or idx_ds == idx0:

pyflwdir/core.py

@@ -277,8 +277,8 @@ def _d8_idx(idx0, shape):
     """Returns linear indices of eight neighboring cells"""
     nrow, ncol = shape
     # assume c-style row-major
-    r = idx0 // ncol
-    c = idx0 % ncol
+    r = int(idx0 // ncol)
+    c = int(idx0 % ncol)
     idxs_lst = list()
     for dr in range(-1, 2):
         for dc in range(-1, 2):

pyflwdir/core_d8.py

@@ -2,8 +2,10 @@
 """Description of D8 flow direction type and methods to convert to/from general
 nextidx."""
 
-from numba import njit, vectorize
+from typing import Tuple
 import numpy as np
+from numba import njit
+
 from . import core
 
 __all__ = []
@@ -50,8 +52,8 @@ def from_array(flwdir, _mv=_mv, dtype=np.intp):
         if flwdir_flat[idx0] == _mv:
             continue
         dr, dc = drdc(flwdir_flat[idx0])
-        r_ds = int(idx0 // ncol + dr)
-        c_ds = int(idx0 % ncol + dc)
+        r_ds = int(idx0 // ncol) + int(dr)
+        c_ds = int(idx0 % ncol) + int(dc)
         pit = dr == 0 and dc == 0
         outside = r_ds >= nrow or c_ds >= ncol or r_ds < 0 or c_ds < 0
         idx_ds = c_ds + r_ds * ncol
@@ -82,25 +84,25 @@ def _downstream_idx(idx0, flwdir_flat, shape, mv=core._mv):
 
 # general
 @njit
-def to_array(idxs_ds, shape, mv=core._mv):
+def to_array(idxs_ds: np.ndarray[np.uint64], shape: Tuple[int, int], mv=core._mv):
     """convert downstream linear indices to dense D8 raster"""
     ncol = shape[1]
     flwdir = np.full(idxs_ds.size, _mv, dtype=np.uint8)
     for idx0 in range(idxs_ds.size):
         idx_ds = idxs_ds[idx0]
         if idx_ds == mv:
             continue
-        dr = (idx_ds // ncol) - (idx0 // ncol)
-        dc = (idx_ds % ncol) - (idx0 % ncol)
+        dr: np.int32 = np.int32(idx_ds // ncol) - np.int32(idx0 // ncol)
+        dc: np.int32 = np.int32(idx_ds % ncol) - np.int32(idx0 % ncol)
         if dr >= -1 and dr <= 1 and dc >= -1 and dc <= 1:
-            dd = _ds[dr + 1, dc + 1]
+            dd: np.uint8 = _ds[dr + 1, dc + 1]
         else:
-            raise ValueError("Invalid data downstream index outside 8 neighbours.")
+            raise ValueError("Invalid data downstream index outside 8 neighbors.")
         flwdir[idx0] = dd
     return flwdir.reshape(shape)
 
 
-def isvalid(flwdir, _all=_all):
+def isvalid(flwdir: np.uint8, _all: np.ndarray[np.uint8] = _all) -> bool:
     """True if 2D D8 raster is valid"""
     return (
         isinstance(flwdir, np.ndarray)

pyflwdir/core_ldd.py

@@ -17,20 +17,24 @@
 _all = np.array([7, 8, 9, 4, 5, 6, 1, 2, 3, 255], dtype=np.uint8)
 
 
+from numba import njit
+import numpy as np
+
+
 @njit("Tuple((int8, int8))(uint8)")
 def drdc(dd):
     """convert ldd value to delta row/col"""
     dr, dc = np.int8(0), np.int8(0)
     if dd >= np.uint8(4):  # W / PIT / E / NW / N / NE
-        if dd >= np.uint(7):  # NW / N / NE
+        if dd >= np.uint8(7):  # NW / N / NE
             dr = np.int8(-1)
-            dc = np.int8(dd - 8)
+            dc = np.int8(dd) - np.int8(8)
         else:  # W / PIT / E
             dr = np.int8(0)
-            dc = np.int8(dd - 5)
+            dc = np.int8(dd) - np.int8(5)
     else:  # SW / S / SE
         dr = np.int8(1)
-        dc = np.int8(dd - 2)
+        dc = np.int8(dd) - np.int8(2)
     return dr, dc
 
 
@@ -87,12 +91,12 @@ def to_array(idxs_ds, shape, mv=core._mv):
         idx_ds = idxs_ds[idx0]
         if idx_ds == mv:
             continue
-        dr = (idx_ds // ncol) - (idx0 // ncol)
-        dc = (idx_ds % ncol) - (idx0 % ncol)
+        dr: np.int32 = np.int32(idx_ds // ncol) - np.int32(idx0 // ncol)
+        dc: np.int32 = np.int32(idx_ds % ncol) - np.int32(idx0 % ncol)
         if dr >= -1 and dr <= 1 and dc >= -1 and dc <= 1:
             dd = _ds[dr + 1, dc + 1]
         else:
-            raise ValueError("Invalid data downstream index outside 8 neighbours.")
+            raise ValueError("Invalid data downstream index outside 8 neighbors.")
         flwdir[idx0] = dd
     return flwdir.reshape(shape)
 

pyflwdir/flwdir.py

@@ -22,20 +22,44 @@
 logger = logging.getLogger(__name__)
 
 
-@njit
-def get_lin_indices(idxs, idxs_ds):
-    idxs_ds0 = np.arange(idxs.size, dtype=idxs.dtype)
-    for j, idx_ds in enumerate(idxs_ds):
-        for i, idx in enumerate(idxs):
-            if idx == idx_ds:
-                idxs_ds0[j] = i
-    return idxs_ds0
+# @njit(cache=True)
+def get_loc_idx(idxs: np.ndarray, idxs_ds: np.ndarray) -> np.ndarray:
+    """Get location indices of downstream cells for all downstream cells"""
+    max_val = idxs.max()
+    idx_map = np.empty(max_val + 1, dtype=np.int_)
+    idx_map.fill(-1)  # Fill with -1 to indicate missing values
+    idx_map[idxs] = np.arange(idxs.size)
+    return idx_map[idxs_ds]
+
+
+def from_dataframe(df: "pandas.DataFrame", ds_col="idx_ds") -> "Flwdir":
+    """Create a Flwdir object from a dataframe with flow direction data.
+
+    Parameters
+    ----------
+    df : pandas.DataFrame
+        dataframe with flow directions
+    ds_col : str, optional
+        name of column with downstream indices, by default "idx_ds"
+
+    Returns
+    -------
+    Flwdir
+        flow direction object
+    """
 
-
-def from_dataframe(df, ds_col="idx_ds"):
     idxs_ds = df[ds_col].values
     idxs = df.index.values
-    return Flwdir(idxs_ds=get_lin_indices(idxs=idxs, idxs_ds=idxs_ds))
+    return Flwdir(idxs_ds=get_loc_idx(idxs=idxs, idxs_ds=idxs_ds))
+
+
+# def _get_idxs_ds_upstream(idxs: np.ndarray, idxs_up: np.ndarray) -> np.ndarray:
+#     idxs_ds0 = np.arange(idxs.size, dtype=idxs.dtype)
+#     for j, idx_up in enumerate(idxs_up):
+#         for i, idx in enumerate(idxs):
+#             if idx == idx_up:
+#                 idxs_up0[j] = i
+#     return idxs_ds0
 
 
 class Flwdir(object):
@@ -96,6 +120,8 @@ def __init__(
         if self.idxs_pit.size == 0:
             raise ValueError("Invalid FlwdirRaster: no pits found")
 
+    ### REPRESENTATION ###
+
     def __str__(self):
         return pprint.pformat(self._dict)
 

pyflwdir/gis_utils.py

@@ -99,7 +99,7 @@ def spread2d(obs, msk=None, nodata=0, frc=None, latlon=False, transform=IDENTITY
             for dc in range(-1, 2):
                 if dr == 0 and dc == 0:
                     continue
-                r1, c1 = r + dr, c + dc
+                r1, c1 = int(r) + dr, int(c) + dc
                 outside = r1 < 0 or r1 >= nrow or c1 < 0 or c1 >= ncol
                 if outside or (msk is not None and ~msk[r1, c1]):
                     continue
@@ -472,10 +472,10 @@ def distance(idx0, idx1, ncol, latlon=False, transform=IDENTITY):
     """
     xres, yres, north = transform[0], transform[4], transform[5]
     # compute delta row, col
-    r0 = idx0 // ncol
-    r1 = idx1 // ncol
+    r0 = int(idx0 // ncol)
+    r1 = int(idx1 // ncol)
     dr = abs(r1 - r0)
-    dc = abs((idx1 % ncol) - (idx0 % ncol))
+    dc = abs(int(idx1 % ncol) - int(idx0 % ncol))
     if latlon:  # calculate cell size in metres
         lat = north + (r0 + r1) / 2.0 * yres
         dy = 0.0 if dr == 0 else degree_metres_y(lat) * yres

pyflwdir/streams.py

@@ -302,7 +302,7 @@ def stream_distance(
     for idx0 in seq:  # down- to upstream
         idx_ds = idxs_ds[idx0]
         # sum distances; skip if at pit or mask is True
-        if idx0 == idx_ds or (mask is not None and mask[idx0] == True):
+        if idx0 == idx_ds or (mask is not None and mask[idx0] is True):
             continue
         if real_length:
             d = gis_utils.distance(idx0, idx_ds, ncol, latlon, transform)

pyflwdir/subgrid.py

@@ -465,7 +465,7 @@ def segment_slope(
                 idx1 == mv
                 or idx1 == idx
                 or outlets[idx1]
-                or (mask is not None and mask[idx1] == False)
+                or (mask is not None and mask[idx1] is False)
             ):
                 break
             idxs.append(idx1)
@@ -535,15 +535,15 @@ def fixed_length_slope(
         x1 = distnc[idx0] + length / 2
         while distnc[idx0] > x0:
             idx_ds = idxs_ds[idx0]
-            if idx_ds == idx0 or (mask is not None and mask[idx0] == False):
+            if idx_ds == idx0 or (mask is not None and mask[idx0] is False):
                 break
             idx0 = idx_ds
         # move upstream and collect x & z
         xs = [distnc[idx0]]
         zs = [elevtn[idx0]]
         while distnc[idx0] < x1:
             idx_us = idxs_us_main[idx0]
-            if idx_us == mv or (mask is not None and mask[idx_us] == False):
+            if idx_us == mv or (mask is not None and mask[idx_us] is False):
                 break
             xs.append(distnc[idx_us])
             zs.append(elevtn[idx_us])

pyflwdir/upscale.py

@@ -22,16 +22,16 @@
 @njit
 def subidx_2_idx(subidx, subncol, cellsize, ncol):
     """Returns the lowres index <idx> of highres cell index <subidx>."""
-    r = (subidx // subncol) // cellsize
-    c = (subidx % subncol) // cellsize
+    r = int(subidx // subncol) // cellsize
+    c = int(subidx % subncol) // cellsize
     return r * ncol + c
 
 
 @njit
 def in_d8(idx0, idx_ds, ncol):
     """Returns True if inside 3x3 (current and 8 neighboring) cells."""
-    cond1 = abs(idx_ds % ncol - idx0 % ncol) <= 1  # west - east
-    cond2 = abs(idx_ds // ncol - idx0 // ncol) <= 1  # south - north
+    cond1 = abs(int(idx_ds % ncol) - int(idx0 % ncol)) <= 1  # west - east
+    cond2 = abs(int(idx_ds // ncol) - int(idx0 // ncol)) <= 1  # south - north
     return cond1 and cond2
 
 
@@ -1068,8 +1068,8 @@ def ihu_minimize_error(
         check_pit = pit_out_of_cell > 0 and subidx_ds == subidx
         if check_pit:
             idx1 = subidx_2_idx(subidx_ds, subncol, cellsize, ncol)
-            dr = (idx1 % ncol) - (idx0 % ncol)
-            dc = (idx1 // ncol) - (idx0 // ncol)
+            dr = int(idx1 % ncol) - int(idx0 % ncol)
+            dc = int(idx1 // ncol) - int(idx0 // ncol)
             check_pit = abs(dr) <= pit_out_of_cell and abs(dc) <= pit_out_of_cell
         # not outlet cells and at pit -> reset outlet to pit
         if check_pit and (subidx_ds == subidx0 or len(idxs) == 0):
@@ -1359,7 +1359,7 @@ def upscale_error(subidxs_out, idxs_ds, subidxs_ds, mv=_mv):
                 # next iter
                 subidx = subidx1
         else:
-            connect_map[idx0] = np.uint8(-1)
+            connect_map[idx0] = np.uint8(255)  # -1
     return connect_map, np.array(idxs_fix_lst, dtype=idxs_ds.dtype)
 
 

tests/test_core.py

@@ -3,42 +3,16 @@
 
 import pytest
 import numpy as np
-import os
 
-from pyflwdir import core, core_d8, streams, from_dem
+from pyflwdir import core, streams
 
-# _mv = core._mv
 
-# # test data
-# # slice from real data
-# testdir = os.path.dirname(__file__)
-# flwdir0 = np.loadtxt(os.path.join(testdir, "flwdir.asc"), dtype=np.uint8)
-# idxs_ds0, idxs_pit0, _ = core_d8.from_array(flwdir0, dtype=np.uint32)
-# rank0, n0 = core.rank(idxs_ds0, mv=np.uint32(_mv))
-# seq0 = np.argsort(rank0)[-n0:]
-# parsed0 = (idxs_ds0, idxs_pit0, seq0, rank0, np.uint32(_mv))
-# # random data
-# # np.random.seed(2345)
-# flwdir1 = from_dem(np.random.rand(15, 10)).to_array("d8")
-# idxs_ds1, idxs_pit1, _ = core_d8.from_array(flwdir1, dtype=np.uint32)
-# rank1, n1 = core.rank(idxs_ds1, mv=np.uint32(_mv))
-# seq1 = np.argsort(rank1)[-n1:]
-# parsed1 = (idxs_ds1, idxs_pit1, seq1, rank1, np.uint32(_mv))
-# # combined
-# test_data = [(parsed0, flwdir0), (parsed1, flwdir1)]
-
-
-def test_downstream0(test_data0, flwdir0):
-    _test_downstream(test_data0, flwdir0)
-
-
-def test_downstream1(test_data1, flwdir1):
-    _test_downstream(test_data1, flwdir1)
-
-
-# @pytest.mark.parametrize("parsed, flwdir", test_data)
-# def test_downstream(test_data0, flwdir0):
-def _test_downstream(test_data, flwdir):
+@pytest.mark.parametrize(
+    "test_data, flwdir", [("test_data0", "flwdir0"), ("test_data0", "flwdir0")]
+)
+def test_downstream(test_data, flwdir, request):
+    test_data = request.getfixturevalue(test_data)
+    flwdir = request.getfixturevalue(flwdir)
     idxs_ds, idxs_pit, seq, rank, mv = [p.copy() for p in test_data]
     n, ncol = np.sum(idxs_ds != mv), flwdir.shape[1]
     # rank
@@ -86,17 +60,13 @@ def _test_downstream(test_data, flwdir):
         assert np.all(rank3[idxs1] == n_up)
 
 
-def test_upstream0(test_data0, flwdir0):
-    _test_upstream(test_data0, flwdir0)
-
-
-def test_upstream1(test_data1, flwdir1):
-    _test_upstream(test_data1, flwdir1)
-
-
-# @pytest.mark.parametrize("parsed, flwdir", test_data)
-def _test_upstream(parsed, flwdir):
-    idxs_ds, idxs_pit, seq, rank, mv = [p.copy() for p in parsed]
+@pytest.mark.parametrize(
+    "test_data, flwdir", [("test_data0", "flwdir0"), ("test_data0", "flwdir0")]
+)
+def test_upstream(test_data, flwdir, request):
+    test_data = request.getfixturevalue(test_data)
+    flwdir = request.getfixturevalue(flwdir)
+    idxs_ds, idxs_pit, seq, rank, mv = [p.copy() for p in test_data]
     idxs_ds[rank == -1] = mv
     n, ncol = np.sum(idxs_ds != mv), flwdir.shape[1]
     upa = streams.upstream_area(idxs_ds, seq, ncol, dtype=np.int32)

tests/test_core_xx.py

@@ -11,10 +11,9 @@
 @pytest.mark.parametrize("fd", [core_nextxy, core_d8, core_ldd])
 def test_core(fd):
     """test core_x.py submodules based on _us definitions"""
-    _name = fd._ftype
     # test isvalid
     assert fd.isvalid(fd._us)
-    assert not fd.isvalid(fd._us * -1)
+    assert not fd.isvalid(fd._us * 3)
     # test ispit
     assert np.all(fd.ispit(fd._pv))
     # test isnodata