Merge pull request #79 from localdevices/issue_78

Issue 78

CHANGELOG.md

@@ -1,3 +1,16 @@
+## [0.2.3] - 2022-08-10
+### Added
+### Changed
+- pyorc.transect.get_q added method="log_interp" using a log-depth normalized velocity and linear interpolation 
+
+### Deprecated
+### Removed
+### Fixed
+- pyorc.transect.get_q improved method="log_fit" (former "log_profile"), so that it works if no dry bed points are found
+
+### Security
+
+
 ## [0.2.2] - 2022-08-01
 ### Added
 - pytest

docs/conf.py

@@ -53,7 +53,7 @@ def remove_dir_content(path: str) -> None:
 # The full version, including alpha/beta/rc tags
 # TODO: uncomment this as soon as we have a version number on the package within pypi
 # release = pkg_resources.get_distribution("ODMax").version
-release = '0.2.2'
+release = '0.2.3'
 
 # -- General configuration ---------------------------------------------------
 

pyorc/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.2.2"
+__version__ = "0.2.3"
 from .api.cameraconfig import CameraConfig, load_camera_config, get_camera_config
 from .api.video import Video
 from .api.frames import Frames

pyorc/api/cameraconfig.py

@@ -155,6 +155,45 @@ def get_depth(self, z, h_a=None):
         z_pressure = np.maximum(self.gcps["z_0"] - h_ref + h_a, z)
         return z_pressure - z
 
+
+    def get_dist_shore(self, x, y, z, h_a=None):
+        """Retrieve depth for measured bathymetry points using the camera configuration and an actual water level, measured
+        in local reference (e.g. staff gauge).
+
+        Parameters
+        ----------
+        z : list of floats
+            measured bathymetry point depths
+        h_a : float, optional
+            actual water level measured [m], if not set, assumption is that a single video
+            is processed and thus changes in water level are not relevant. (default: None)
+
+        Returns
+        -------
+        depths : list of floats
+
+        """
+        # retrieve depth
+        depth = self.get_depth(z, h_a=h_a)
+        if h_a is None:
+            assert(self.gcps["h_ref"] is None), "No actual water level is provided, but a reference water level is provided"
+            h_a = 0.
+            h_ref = 0.
+        else:
+            h_ref = self.gcps["h_ref"]
+        z_dry = depth <= 0
+        z_dry[[0, -1]] = True
+        # compute distance to nearest dry points with Pythagoras
+        dist_shore = np.array([(((x[z_dry] - _x) ** 2 + (y[z_dry] - _y) ** 2) ** 0.5).min() for _x, _y, in zip(x, y)])
+        return dist_shore
+
+
+    def get_dist_wall(self, x, y, z, h_a=None):
+        depth = self.get_depth(z, h_a=h_a)
+        dist_shore = self.get_dist_shore(x, y, z, h_a=h_a)
+        dist_wall = (dist_shore**2 + depth**2)**0.5
+        return dist_wall
+
     def z_to_h(self, z):
         """Convert z coordinates of bathymetry to height coordinates in local reference (e.g. staff gauge)
 

pyorc/api/transect.py

@@ -7,6 +7,7 @@
 from .plot import _Transect_PlotMethods
 from .orcbase import ORCBase
 
+
 @xr.register_dataset_accessor("transect")
 class Transect(ORCBase):
     """Transect functionalities that can be applied on ``xarray.Dataset``"""
@@ -57,14 +58,14 @@ def vector_to_scalar(self, v_x="v_x", v_y="v_y"):
         v_eff = np.cos(angle_diff) * v_scalar
         v_eff.attrs = {
             "standard_name": "velocity",
-            "long_name": "velocity in perpendicular direction of cross section, measured by angle in radians, measured from up-direction",
+            "long_name": "velocity in perpendicular direction of cross section, measured by angle in radians, "
+                         "measured from up-direction",
             "units": "m s-1",
         }
         # set name
         v_eff.name = "v_eff_nofill"  # there still may be gaps in this series
         self._obj["v_eff_nofill"] = v_eff
 
-
     def get_xyz_perspective(self, M=None, xs=None, ys=None, mask_outside=True):
         """Get camera-perspective column, row coordinates from cross-section locations.
         
@@ -154,7 +155,6 @@ def get_river_flow(self, q_name="q", Q_name="river_flow"):
         Q.name = "Q"
         self._obj[Q_name] = Q
 
-
     def get_q(self, v_corr=0.9, fill_method="zeros"):
         """Depth integrated velocity for quantiles of time series using a correction v_corr between surface velocity and
         depth-average velocity.
@@ -164,9 +164,10 @@ def get_q(self, v_corr=0.9, fill_method="zeros"):
         v_corr : float, optional
             correction factor (default: 0.9)
         fill_method : method to fill missing values. "zeros" fills NaNS with zeros, "interpolate" interpolates values
-            from nearest neighbour, "log_profile" fits a 4-parameter logarithmic profile with depth and with changing
-             velocities towards banks on known velocities, and fills missing with the fitted relationship.
-             (Default value = "zeros")
+            from nearest neighbour, "log_interp" interpolates values linearly with velocities scaled by the log of
+             depth over a roughness length, "log_fit" fits a 4-parameter logarithmic profile with depth and with
+             changing velocities towards banks on known velocities, and fills missing with the fitted relationship
+             (experimental) (Default value = "zeros").
 
         Returns
         -------
@@ -176,7 +177,8 @@ def get_q(self, v_corr=0.9, fill_method="zeros"):
 
         """
         # aggregate to a limited set of quantiles
-        assert(fill_method in ["zeros", "log_profile", "interpolate"]), f'fill_method must be "zeros", "log_profile", or "interpolate", instead "{fill_method}" given'
+        assert(fill_method in ["zeros", "log_fit", "log_interp", "interpolate"]),\
+            f'fill_method must be "zeros", "log_fit", "log_interp", or "interpolate", instead "{fill_method}" given'
         ds = self._obj
         x = ds["xcoords"].values
         y = ds["ycoords"].values
@@ -185,8 +187,12 @@ def get_q(self, v_corr=0.9, fill_method="zeros"):
         depth = self.camera_config.get_depth(z, self.h_a)
         if fill_method == "zeros":
             ds["v_eff"] = ds["v_eff_nofill"].fillna(0.)
-        elif fill_method == "log_profile":
-            ds["v_eff"] = helpers.velocity_fill(x, y, depth, ds["v_eff_nofill"], groupby="quantile")
+        elif fill_method == "log_fit":
+            dist_shore = self.camera_config.get_dist_shore(x, y, z, self.h_a)
+            ds["v_eff"] = helpers.velocity_log_fit(ds["v_eff_nofill"], depth, dist_shore, dim="quantile")
+        elif fill_method == "log_interp":
+            dist_wall = self.camera_config.get_dist_wall(x, y, z, self.h_a)
+            ds["v_eff"] = helpers.velocity_log_interp(ds["v_eff_nofill"], dist_wall, dim="quantile")
         elif fill_method == "interpolate":
             depth = ds.zcoords*0 + self.camera_config.get_depth(ds.zcoords, self.h_a)
             # interpolate gaps in between known values

pyorc/helpers.py

@@ -204,7 +204,7 @@ def log_profile(X, z0, k_max, s0=0., s1=0.):
     k_max: float
         maximum scale factor of log-profile function [-]
     s0: float, optional
-        distance from bank (defaukt: 0.) where k equals zero (and thus velocity is zero) [m]
+        distance from bank (default: 0.) where k equals zero (and thus velocity is zero) [m]
     s1: float, optional
         distance from bank (default: 0. meaning no difference over distance) where k=k_max (k cannot be larger than
         k_max) [m]
@@ -215,8 +215,8 @@ def log_profile(X, z0, k_max, s0=0., s1=0.):
         V values from log-profile, equal shape as arrays inside X [m s-1]
     """
     z, s = X
-    k = k_max * np.minimum(np.maximum((s-s0)/(s1-s0), 0), 1)
-    v = k*np.maximum(np.log(np.maximum(z, 1e-6)/z0), 0)
+    k = k_max * np.minimum(np.maximum((s - s0) / (s1 - s0), 0), 1)
+    v = k * np.maximum(np.log(np.maximum(z, 1e-6) / z0), 0)
     return v
 
 
@@ -265,11 +265,11 @@ def neighbour_stack(array, stride=1, missing=-9999.):
     array = copy.deepcopy(array)
     array[np.isnan(array)] = missing
     array_move = []
-    for vert in range(-stride, stride+1):
-        for horz in range(-stride, stride+1):
-            conv_arr = np.zeros((abs(vert)*2+1, abs(horz)*2+1))
-            _y = int(np.floor((abs(vert)*2+1)/2)) + vert
-            _x = int(np.floor((abs(horz)*2+1)/2)) + horz
+    for vert in range(-stride, stride + 1):
+        for horz in range(-stride, stride + 1):
+            conv_arr = np.zeros((abs(vert) * 2 + 1, abs(horz) * 2 + 1))
+            _y = int(np.floor((abs(vert) * 2 + 1) / 2)) + vert
+            _x = int(np.floor((abs(horz) * 2 + 1) / 2)) + horz
             conv_arr[_y, _x] = 1
             array_move.append(convolve2d(array, conv_arr, mode="same", fillvalue=np.nan))
     array_move = np.stack(array_move)
@@ -289,7 +289,6 @@ def optimize_log_profile(
         seed=0,
         **kwargs
 ):
-
     """optimize velocity log profile relation of v=k*max(z/z0) with k a function of distance to bank and k_max
     A differential evolution optimizer is used.
 
@@ -299,7 +298,7 @@ def optimize_log_profile(
         depths [m]
     v : list
         surface velocities [m s-1]
-    dist_bank : list
+    dist_bank : list, optional
         distances to bank [m]
     **kwargs : keyword arguments for scipy.optimize.differential_evolution
 
@@ -309,7 +308,7 @@ def optimize_log_profile(
         fitted parameters of log_profile {z_0, k_max, s0 and s1}
     """
     if dist_bank is None:
-        dist_bank = np.inf(len(v))
+        dist_bank = np.ones(len(v)) * np.inf
     v = np.array(v)
     z = np.array(z)
     X = (z, dist_bank)
@@ -361,20 +360,18 @@ def rotate_u_v(u, v, theta, deg=False):
     return u2, v2
 
 
-def velocity_fill(x, y, depth, v, groupby="quantile"):
+def velocity_log_fit(v, depth, dist_shore, dim="quantile"):
     """Fill missing surface velocities using a velocity depth profile with
 
     Parameters
     ----------
-    x : xr.DataArray (points)
-        x-coordinates of bathymetry depths (ref. CRS)
-    y : xr.DataArray (points)
-        y-coordinates of bathymetry depths (ref. CRS)
-    depth : xr.DataArray (points)
-        depths [m]
     v : xr.DataArray (time, points)
         effective velocity at surface [m s-1]
-    groupby: str, optional
+    depth : xr.DataArray (points)
+        bathymetry depths [m]
+    dist_shore : xr.DataArray (points)
+        shortest distance to a dry river bed point
+    dim: str, optional
         dimension over which data should be grouped, default: "quantile", dimension must exist in v, typically
         "quantile" or "time"
 
@@ -384,18 +381,65 @@ def velocity_fill(x, y, depth, v, groupby="quantile"):
     v_fill: xr.DataArray (quantile or time, points)
         filled surface velocities  [m s-1]
     """
-    def fit(_v):
-        pars = optimize_log_profile(depth[np.isfinite(_v).values], _v[np.isfinite(_v).values], dist_bank[np.isfinite(_v).values])
-        _v[np.isnan(_v).values] = log_profile((depth[np.isnan(_v).values], dist_bank[np.isnan(_v).values]), **pars)
+
+    def log_fit(_v):
+        pars = optimize_log_profile(
+            depth[np.isfinite(_v).values],
+            _v[np.isfinite(_v).values],
+            dist_shore[np.isfinite(_v).values]
+        )
+        _v[np.isnan(_v).values] = log_profile(
+            (
+                depth[np.isnan(_v).values],
+                dist_shore[np.isnan(_v).values]
+            ),
+            **pars
+        )
         # enforce that velocities are zero with zero depth
-        _v[depth<=0] = 0.
+        _v[depth <= 0] = 0.
         return np.maximum(_v, 0)
 
-    z_dry = depth <= 0
-    dist_bank = np.array([(((x[z_dry] - _x) ** 2 + (y[z_dry] - _y) ** 2) ** 0.5).min() for _x, _y, in zip(x, y)])
-    # per time slice or quantile, fill missings
-    v_group = copy.deepcopy(v).groupby(groupby)
-    return v_group.map(fit)
+    # fill per grouped dimension
+    v.load()
+    v_group = copy.deepcopy(v).groupby(dim)
+    return v_group.map(log_fit)
+
+
+def velocity_log_interp(v, dist_wall, d_0=0.1, dim="quantile"):
+    """
+
+    Parameters
+    ----------
+    v : xr.DataArray (time, points)
+        effective velocity at surface [m s-1]
+    dist_wall : xr.DataArray (points)
+        shortest distance to the river bed
+    d_0 : float, optional
+        roughness length (default: 0.1)
+    dim: str, optional
+        dimension over which data should be grouped, default: "quantile", dimension must exist in v, typically
+        "quantile" or "time"
+
+    Returns
+    -------
+
+    """
+
+    def log_interp(_v):
+        # scale with log depth
+        c = xr.DataArray(_v / np.log(np.maximum(dist_wall, d_0) / d_0))
+        # fill dry points with the nearest valid value for c
+        c[dist_wall == 0] = c.interpolate_na(dim="points", method="nearest", fill_value="extrapolate")[dist_wall == 0]
+        # interpolate with linear interpolation
+        c = c.interpolate_na(dim="points")
+        # use filled c to interpret missing v
+        _v[np.isnan(_v)] = (np.log(np.maximum(dist_wall, d_0) / d_0) * c)[np.isnan(_v)]
+        return _v
+
+    # fill per grouped dimension
+    v.load()
+    v_group = copy.deepcopy(v).groupby(dim)
+    return v_group.map(log_interp)
 
 
 def wrap_mse(pars_iter, *args):
@@ -475,6 +519,7 @@ def xy_angle(x, y):
     angles[-1] = np.arctan2(x[-1] - x[-2], y[-1] - y[-2])
     return angles
 
+
 def xy_to_perspective(x, y, resolution, M, reverse_y=None):
     """
     Back transform local meters-from-top-left coordinates from frame to original perspective of camera using M

setup.py

@@ -12,7 +12,7 @@
 setup(
     name="pyopenrivercam",
     description="pyopenrivercam (pyorc) is a front and backend to control river camera observation locations",
-    version="0.2.2",
+    version="0.2.3",
     long_description=readme + "\n\n",
     long_description_content_type="text/markdown",
     url="https://github.com/localdevices/pyorc",

tests/test_transect.py

@@ -14,7 +14,8 @@ def test_get_river_flow(piv_transect):
     "fill_method",
     [
         "zeros",
-        "log_profile",
+        "log_interp",
+        "log_fit",
         "interpolate",
     ]
 )