Analytical overhaul

Makefile

@@ -9,7 +9,7 @@ qa:
 	pre-commit run --all-files
 
 unit-tests:
-	python -m pytest -vv --cov=seaduck --cov-report=$(COV_REPORT) --doctest-glob="*.md" --doctest-glob="*.rst"
+	python -m pytest -vv --cov=seaduck --cov-report=$(COV_REPORT) --doctest-glob="*.md" --doctest-glob="*.rst" -W ignore::RuntimeWarning
 
 type-check:
 	python -m mypy .

pyproject.toml

@@ -18,7 +18,7 @@ dependencies = [
   "numpy",
   "pandas",
   "scipy",
-  "dask",
+  "dask[array]",
   "xarray"
 ]
 description = "A python package that interpolates data from ocean dataset from both Eulerian and Lagrangian perspective. "
@@ -39,7 +39,8 @@ module = [
   "matplotlib.*",
   "numba",
   "pooch",
-  "scipy"
+  "scipy",
+  "scipy.interpolate"
 ]
 
 [tool.ruff]

seaduck/eulerian.py

@@ -1,4 +1,5 @@
 import copy
+import logging
 
 import numpy as np
 
@@ -112,8 +113,10 @@ class Position:
     To actually do interpolation, use from_latlon method to tell the ducks where they are.
     """
 
-    def __init__(self):
-        self.rel = RelCoord()
+    def __new__(cls, *arg, **kwarg):
+        new_position = object.__new__(cls)
+        new_position.rel = RelCoord()
+        return new_position
 
     def __getattr__(self, attr):
         if attr == "rel":
@@ -161,30 +164,30 @@ def from_latlon(self, x=None, y=None, z=None, t=None, data=None):
         if any(i != self.N for i in length if i > 1):
             raise ValueError("Shapes of input coordinates are not compatible")
 
-        if isinstance(x, float):
-            x = np.array([1.0]) * x
-        if isinstance(y, float):
-            y = np.array([1.0]) * y
-        if isinstance(z, float):
-            z = np.array([1.0]) * z
-        if isinstance(t, float):
-            t = np.array([1.0]) * t
+        if isinstance(x, (int, float, np.floating)):
+            x = np.ones(self.N, float) * x
+        if isinstance(y, (int, float, np.floating)):
+            y = np.ones(self.N, float) * y
+        if isinstance(z, (int, float, np.floating)):
+            z = np.ones(self.N, float) * z
+        if isinstance(t, (int, float, np.floating)):
+            t = np.ones(self.N, float) * t
 
         for thing in [x, y, z, t]:
             if thing is None:
                 continue
             if len(thing.shape) > 1:
                 raise ValueError("Input need to be 1D numpy arrays")
         if (x is not None) and (y is not None):
-            self.lon = x
-            self.lat = y
-            self.rel.update(self.ocedata.find_rel_h(x, y))
+            self.lon = copy.deepcopy(x)
+            self.lat = copy.deepcopy(y)
+            self.rel.update(self.ocedata.find_rel_h(self.lon, self.lat))
         else:
             self.rel.update(HRel._make([None for i in range(11)]))
             self.lon = None
             self.lat = None
         if z is not None:
-            self.dep = z
+            self.dep = copy.deepcopy(z)
             if self.ocedata.readiness["Z"]:
                 self.rel.update(self.ocedata.find_rel_v(z))
             else:
@@ -199,7 +202,7 @@ def from_latlon(self, x=None, y=None, z=None, t=None, data=None):
             self.dep = None
 
         if t is not None:
-            self.t = t
+            self.t = copy.deepcopy(t)
             if self.ocedata.readiness["time"]:
                 self.rel.update(self.ocedata.find_rel_t(t))
             else:
@@ -225,23 +228,47 @@ def subset(self, which):
         the_subset: Position object
             The selected Positions.
         """
-        p = Position()
+        p = object.__new__(type(self))
         vardict = vars(self)
         keys = vardict.keys()
-        for i in keys:
-            item = vardict[i]
+        for key in keys:
+            item = vardict[key]
             if isinstance(item, np.ndarray):
                 if len(item.shape) == 1:
-                    setattr(p, i, item[which])
-                    p.N = len(getattr(p, i))
+                    setattr(p, key, item[which])
+                    p.N = len(getattr(p, key))
+                elif key in ["px", "py"]:
+                    setattr(p, key, item[:, which])
                 else:
-                    setattr(p, i, item)
+                    setattr(p, key, item)
             elif isinstance(item, RelCoord):
-                setattr(p, i, item.subset(which))
+                setattr(p, key, item.subset(which))
             else:
-                setattr(p, i, item)
+                setattr(p, key, item)
         return p
 
+    def update_from_subset(self, sub, which):
+        vardict = vars(sub)
+        keys = vardict.keys()
+        for key in keys:
+            item = vardict[key]
+            if not hasattr(self, key):
+                logging.warning(
+                    f"A new attribute {key} defined" "after updating from subset"
+                )
+                setattr(self, key, item)
+            if getattr(self, key) is None:
+                continue
+            if isinstance(item, np.ndarray):
+                if len(item.shape) == 1:
+                    getattr(self, key)[which] = item
+                elif key in ["px", "py"]:
+                    getattr(self, key)[:, which] = item
+            elif isinstance(item, RelCoord):
+                self.rel.update_from_subset(item, which)
+            elif isinstance(item, list):
+                setattr(self, key, item)
+
     def fatten_h(self, knw, ind_moves_kwarg={}):
         """Fatten horizontal indices.
 
@@ -262,18 +289,18 @@ def fatten_h(self, knw, ind_moves_kwarg={}):
             Read Topology.ind_moves for more detail.
         """
         #         self.ind_h_dict
-        kernel = knw.kernel
+        kernel = knw.kernel.astype(int)
         kernel_tends = [_translate_to_tendency(k) for k in kernel]
         m = len(kernel_tends)
         n = len(self.iy)
         tp = self.ocedata.tp
 
         # the arrays we are going to return
         if self.face is not None:
-            n_faces = np.zeros((n, m))
+            n_faces = np.zeros((n, m), int)
             n_faces.T[:] = self.face
-        n_iys = np.zeros((n, m))
-        n_ixs = np.zeros((n, m))
+        n_iys = np.zeros((n, m), int)
+        n_ixs = np.zeros((n, m), int)
 
         # first try to fatten it naively(fast and vectorized)
         for i, node in enumerate(kernel):
@@ -302,9 +329,9 @@ def fatten_h(self, knw, ind_moves_kwarg={}):
             else:
                 n_iys[j, i], n_ixs[j, i] = n_ind
         if self.face is not None:
-            return n_faces.astype("int"), n_iys.astype("int"), n_ixs.astype("int")
+            return n_faces, n_iys, n_ixs
         else:
-            return None, n_iys.astype("int"), n_ixs.astype("int")
+            return None, n_iys, n_ixs
 
     def fatten_v(self, knw):
         """Fatten in vertical center coord.
@@ -319,13 +346,13 @@ def fatten_v(self, knw):
         if self.iz is None:
             return None
         if knw.vkernel == "nearest":
-            return copy.deepcopy(self.iz.astype(int))
+            return copy.deepcopy(self.iz)
         elif knw.vkernel in ["dz", "linear"]:
             try:
                 self.iz_lin
             except AttributeError:
                 self.rel.update(self.ocedata.find_rel_v_lin(self.dep))
-            return np.vstack([self.iz_lin.astype(int), self.iz_lin.astype(int) - 1]).T
+            return np.vstack([self.iz_lin, self.iz_lin - 1]).T
         else:
             raise ValueError("vkernel not supported")
 
@@ -342,13 +369,13 @@ def fatten_vl(self, knw):
         if self.izl is None:
             return None
         if knw.vkernel == "nearest":
-            return copy.deepcopy(self.izl.astype(int))
+            return copy.deepcopy(self.izl)
         elif knw.vkernel in ["dz", "linear"]:
             try:
                 self.izl_lin
             except AttributeError:
                 self.rel.update(self.ocedata.find_rel_vl_lin(self.dep))
-            return np.vstack([self.izl_lin.astype(int), self.izl_lin.astype(int) - 1]).T
+            return np.vstack([self.izl_lin, self.izl_lin - 1]).T
         else:
             raise ValueError("vkernel not supported")
 
@@ -365,13 +392,13 @@ def fatten_t(self, knw):
         if self.it is None:
             return None
         if knw.tkernel == "nearest":
-            return copy.deepcopy(self.it.astype(int))
+            return copy.deepcopy(self.it)
         elif knw.tkernel in ["dt", "linear"]:
             try:
                 self.izl_lin
             except AttributeError:
                 self.rel.update(self.ocedata.find_rel_t_lin(self.t))
-            return np.vstack([self.it_lin.astype(int), self.it_lin.astype(int) + 1]).T
+            return np.vstack([self.it_lin, self.it_lin + 1]).T
         else:
             raise ValueError("tkernel not supported")
 

seaduck/get_masks.py

@@ -107,7 +107,7 @@ def mask_w_node(maskC, tp=None):
     return maskW
 
 
-def get_masks(od, tp):
+def get_mask_arrays(od):
     """Mask all staggered valocity points.
 
     A wrapper around mask_u_node, mask_v_node, mask_w_node.
@@ -125,6 +125,7 @@ def get_masks(od, tp):
     maskC,maskU,maskV,maskW: numpy.ndarray
         masks at center points, U-walls, V-walls, W-walls respectively.
     """
+    tp = od.tp
     keys = od._ds.keys()
     if "maskC" not in keys:
         warnings.warn("no maskC in the dataset, assuming nothing is masked.")
@@ -217,10 +218,13 @@ def get_masked(od, ind, cuvwg="C"):
 
 def which_not_stuck(p):
     """Investigate which points are in land mask."""
+    ind = []
+    if p.izl_lin is not None:
+        ind.append(p.izl_lin - 1)
     if p.face is not None:
-        ind = (p.izl_lin - 1, p.face, p.iy, p.ix)
-    else:
-        ind = (p.izl_lin - 1, p.iy, p.ix)
+        ind.append(p.face)
+    ind += [p.iy, p.ix]
+    ind = tuple(ind)
     return get_masked(p.ocedata, ind).astype(bool)
 
 

seaduck/kernel_weight.py

@@ -121,8 +121,8 @@ def kernel_weight_x(kernel, ktype="interp", order=0):
         function to calculate the hotizontal interpolation/derivative
         weight
     """
-    xs = np.array(list(set(kernel.T[0]))).astype(float)
-    ys = np.array(list(set(kernel.T[1]))).astype(float)
+    xs = np.array(list(set(kernel.T[0])), dtype=float)
+    ys = np.array(list(set(kernel.T[1])), dtype=float)
 
     # if you the kernel is a line rather than a cross
     if len(xs) == 1:
@@ -149,8 +149,8 @@ def kernel_weight_x(kernel, ktype="interp", order=0):
             y_poly.append(
                 list(combinations([i for i in ys if i != ay], len(ys) - 1 - order))
             )
-    x_poly = np.array(x_poly).astype(float)
-    y_poly = np.array(y_poly).astype(float)
+    x_poly = np.array(x_poly, dtype=float)
+    y_poly = np.array(y_poly, dtype=float)
 
     @compileable
     def the_interp_func(rx, ry):
@@ -335,8 +335,8 @@ def kernel_weight_s(kernel, xorder=0, yorder=0):
         function to calculate the hotizontal interpolation/derivative
         weight
     """
-    xs = np.array(list(set(kernel.T[0]))).astype(float)
-    ys = np.array(list(set(kernel.T[1]))).astype(float)
+    xs = np.array(list(set(kernel.T[0])), dtype=float)
+    ys = np.array(list(set(kernel.T[1])), dtype=float)
     xmaxorder = False
     ymaxorder = False
     if xorder < len(xs) - 1:
@@ -363,8 +363,8 @@ def kernel_weight_s(kernel, xorder=0, yorder=0):
         y_poly.append(
             list(combinations([i for i in ys if i != ay], len(ys) - 1 - yorder))
         )
-    x_poly = np.array(x_poly).astype(float)
-    y_poly = np.array(y_poly).astype(float)
+    x_poly = np.array(x_poly, dtype=float)
+    y_poly = np.array(y_poly, dtype=float)
 
     @compileable
     def the_square_func(rx, ry):
@@ -722,8 +722,8 @@ def same_hsize(self, other):
         if not type_same:
             raise TypeError("the argument is not a KnW object")
         try:
-            return (self.kernel == other.kernel).all()
-        except AttributeError:
+            return np.allclose(self.kernel, other.kernel)
+        except (ValueError, AttributeError):
             return False
 
     def same_size(self, other):
@@ -738,9 +738,7 @@ def __eq__(self, other):
         type_same = isinstance(other, type(self))
         if not type_same:
             return False
-        shpe_same = (
-            self.kernel == other.kernel
-        ).all() and self.inheritance == other.inheritance
+        shpe_same = self.same_hsize(other) and self.inheritance == other.inheritance
         diff_same = (
             (self.hkernel == other.hkernel)
             and (self.vkernel == other.vkernel)