From 4af3a28a8a554c746c39fc418614c83f1adaf692 Mon Sep 17 00:00:00 2001 From: Andrew Date: Thu, 26 Jan 2023 11:54:38 -0600 Subject: [PATCH 1/3] draft of two metrics on topset. both use equally-spaced radial sections to sample the data. --- deltametrics/plan.py | 146 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) diff --git a/deltametrics/plan.py b/deltametrics/plan.py index 43f18eca..5e457e33 100644 --- a/deltametrics/plan.py +++ b/deltametrics/plan.py @@ -1660,6 +1660,152 @@ def compute_shoreline_distance(shore_mask, origin=(0, 0), return_distances=False return np.nanmean(_dists), np.nanstd(_dists) +def _determine_equally_spaced_azimuths(*args, **kwargs): + """Helper to determine equally spaced azimuths for other computations. + + Parameters + ---------- + num : :obj:`int` or :obj:`float` + Number of azimuths to return. Coerced to `int`. + + start : :obj:`float` + Starting azimuth. + + end : :obj:`float` + Ending azimuth. + + buffered + Whether `start` and `end` represent the bounds for the calculation and + should be buffered before the first and last azimuth, or if `start` + and `end` are the actual end points of the azimuth array. This is + helpful when `start` and `end` represent model or experimental domain + edges. Default is `True`, include buffer. + + .. doctest:: + + >>> _determine_equally_spaced_azimuths(3, 0, 180, buffered=False) + [ 0. 90. 180.] + + >>> _determine_equally_spaced_azimuths(3, 0, 180, buffered=True) + [ 45. 90. 135.] + + >>> _determine_equally_spaced_azimuths( + ... num=5, start=22.5, end=157.5, buffered=True) + [ 45. 67.5 90. 112.5 135. ] + + >>> _determine_equally_spaced_azimuths() + [ 15. 30. 45. 60. 75. 90. 105. 120. 135. 150. 165.] + + """ + # process the input arguments + if len(args) > 0: + # must specify all four as args + if len(args) != 4: + raise ValueError( + "Must specify all four arguments as positional, " + "if specifying any as positional." + ) + # unpack args + num, start, end, buffered = args[0], args[1], args[2], args[3] + else: + # arguments have been specified as keywords or use defaults + num = kwargs.pop("num", 11) + start = kwargs.pop("start", 0) + end = kwargs.pop("end", 180) + buffered = kwargs.pop("buffered", True) + + # create the azimuths + if buffered: + azimuths = np.linspace(start, end, num + 2, dtype=float) + azimuths = azimuths[1:-1] + else: + azimuths = np.linspace(start, end, num, dtype=float) + return azimuths + + +def compute_shoreline_radius(shore_mask, origin=[0, 0], return_radii=False, **kwargs): + """Radially-averaged shoreline radius. + + .. note:: + + In implementation, this metric uses the last (farthest) intersection + of the radial section at `azimuth` and the shoreline mask. + + See also: + :obj:`compute_shoreline_distance`. + """ + azimuths = _determine_equally_spaced_azimuths(**kwargs) + radii = np.zeros((len(azimuths),)) + # note: + for a, azimuth in enumerate(azimuths): + a_section = dm.section.RadialSection( + cube, azimuth=azimuth, origin_idx=[0, eta_diff.shape[1] // 2] + ) + + # find where interescts shoreline mask + shoreline_mask_alongsection = deposit_percentile_line[ + a_section.trace_idx[:, 0], a_section.trace_idx[:, 1] + ] + where_intersects = np.nonzero(shoreline_mask_alongsection)[0] + if where_intersects.size > 0: + radii[a] = a_section.s[where_intersects[-1]] + else: + radii[a] = np.nan + + if return_radii: + return np.nanmean(slopes), np.nanstd(slopes), slopes + else: + return np.nanmean(slopes), np.nanstd(slopes) + + +def compute_topset_slope( + shore_mask, + elevation_data, + origin=[0, 0], + elevation_threshold=0, + point_threshold=5, + return_slopes=False, +): + """ + + Parameters + ---------- + + point_threshold + How many points above sea level must be identified in order to fit a + line. + """ + azimuths = _determine_equally_spaced_azimuths(**kwargs) + slopes = np.zeros((len(azimuths),)) + for a, azimuth in enumerate(azimuths): + a_section = dm.section.RadialSection( + cube, azimuth=azimuth, origin_idx=[0, eta_diff.shape[1] // 2] + ) + + _HSL = cube.meta["H_SL"][t_idx].data + yvalues = np.array(a_section["eta"][t_idx, :]) + xvalues = np.array(a_section.s) + if np.isnan(_HSL): + # is no sea level, take locations where deposited + _keep_bool = yvalues > a_section["eta"][0, :] + else: + # if sea level, take locations where above sl + _keep_bool = yvalues > (_HSL - 0.25) + yvalues_above = yvalues[_keep_bool] + xvalues_above = xvalues[_keep_bool] + + if xvalues_above.size > point_threshold: + m, b = np.polyfit(xvalues_above, yvalues_above, 1) + else: + m = np.nan + slopes[a] = m + + if return_slopes: + return np.nanmean(slopes), np.nanstd(slopes), slopes + else: + return np.nanmean(slopes), np.nanstd(slopes) + + @njit(parallel=True) def _compute_angles_between(test_set_points, query_set_points, numviews): """Private helper for shaw_opening_angle_method. From 3c2f0cb5e6b96537f95515849858057674538aef Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 3 Jun 2022 15:20:56 -0500 Subject: [PATCH 2/3] add draft of shoreline rugosity. --- deltametrics/plan.py | 83 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 78 insertions(+), 5 deletions(-) diff --git a/deltametrics/plan.py b/deltametrics/plan.py index 5e457e33..5b0c9eed 100644 --- a/deltametrics/plan.py +++ b/deltametrics/plan.py @@ -1323,6 +1323,73 @@ def compute_shoreline_roughness(shore_mask, land_mask, **kwargs): return rough +def compute_shoreline_rugosity(shore_mask, **kwargs): + """Compute shoreline rugosity. + + Computes the shoreline rugosity metric: + + .. math:: + + R_j = \\sqrt{1/N \\sum_{i=1}^N \\left( \\frac{r_{i,j}-\\bar{r}}{\\bar{r}}\\right)^2} + + where ... + + Parameters + ---------- + shore_mask : :obj:`~deltametrics.mask.ShorelineMask`, :obj:`ndarray` + Shoreline mask. Can be a :obj:`~deltametrics.mask.ShorelineMask` object, + or a binarized array. + + **kwargs + Keyword argument are passed to :obj:`compute_shoreline_length` + internally. + + Returns + ------- + rugosity : :obj:`float` + Shoreline rugosity, computed as described above. + + Examples + -------- + Compare the rugosity of the shoreline early in the model simulation with + the rugosity later. Here, we use the `elevation_offset` parameter (passed + to :obj:`~deltametrics.mask.ElevationMask`) to better capture the + topography of the `pyDeltaRCM` model results. + + after Barefoot + """ + # extract data from masks + if isinstance(shore_mask, mask.ShorelineMask): + shore_mask = shore_mask.mask + _sm = shore_mask.values + _dx = float(shore_mask[shore_mask.dims[0]][1] - + shore_mask[shore_mask.dims[0]][0]) + elif isinstance(shore_mask, xr.core.dataarray.DataArray): + _sm = shore_mask.values + _dx = float(shore_mask[shore_mask.dims[0]][1] - + shore_mask[shore_mask.dims[0]][0]) + elif isinstance(shore_mask, np.ndarray): + _sm = shore_mask + _dx = 1 + else: + raise TypeError('Invalid type {0}'.format(type(shore_mask))) + + # _ = kwargs.pop('return_line', None) # trash this variable if passed + # shorelength = compute_shoreline_length( + # shore_mask, return_line=False, **kwargs) + # find where the mask is True (all x-y pairs along shore) + _y, _x = np.argwhere(_sm).T + + N = np.sum(_sm) + if N > 0: + # compute rugosity + rugosity = np.sqrt(1/N * np.sum()) + else: + raise ValueError('No pixels in land mask.') + + return rugosity + + def compute_shoreline_length(shore_mask, origin=(0, 0), return_line=False): """Compute the length of a shoreline from a mask of the shoreline. @@ -2591,6 +2658,11 @@ def compute_surface_deposit_age(data, surface_idx, **kwargs): function :obj:`_compute_surface_deposit_time_from_etas`. Hint: you may want to control the `stasis_tol` parameter. + Returns + ------- + sfc_age : :obj:`np.ndarray` + Age of surface (in indices). + Examples -------- @@ -2606,12 +2678,12 @@ def compute_surface_deposit_age(data, surface_idx, **kwargs): >>> fig, ax = plt.subplots() >>> _ = ax.imshow(sfc_time, cmap='YlGn_r') """ - sfc_date = compute_surface_deposit_time(data, surface_idx, **kwargs) + sfc_time = compute_surface_deposit_time(data, surface_idx, **kwargs) # handle indices less than 0 if surface_idx < 0: # wrap to find the index surface_idx = surface_idx % data.shape[0] - return surface_idx - sfc_date + return surface_idx - sfc_time def _compute_surface_deposit_time_from_etas(etas, stasis_tol=0.01): @@ -2636,7 +2708,8 @@ def _compute_surface_deposit_time_from_etas(etas, stasis_tol=0.01): Returns ------- - sfc_date : obj:`ndarray` + sfc_time : :obj:`ndarray` + Time idx of surface. """ if not (stasis_tol > 0): raise ValueError( @@ -2646,6 +2719,6 @@ def _compute_surface_deposit_time_from_etas(etas, stasis_tol=0.01): etaf = np.array(etas[-1, :, :]) whr = np.abs(etaf - etas) < stasis_tol - sfc_date = np.argmax(whr, axis=0) + sfc_time = np.argmax(whr, axis=0) - return sfc_date + return sfc_time From 042fada0fd936be11bda12a5ea8e13d2dab015fc Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 3 Jun 2022 15:21:33 -0500 Subject: [PATCH 3/3] add some documentation for surface time. --- docs/source/reference/plan/index.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/reference/plan/index.rst b/docs/source/reference/plan/index.rst index bc0de150..dc42d5cd 100644 --- a/docs/source/reference/plan/index.rst +++ b/docs/source/reference/plan/index.rst @@ -49,3 +49,4 @@ Functions compute_surface_deposit_age shaw_opening_angle_method morphological_closing_method + _compute_surface_deposit_time_from_etas