From 34f32e49fdc0f9b31f6dbebf594935568dbf0444 Mon Sep 17 00:00:00 2001 From: malmans2 Date: Mon, 21 Jun 2021 10:14:40 +0100 Subject: [PATCH 1/7] Do not drop indexes when computing rmax --- pydomcfg/tests/bathymetry.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index a135a6a..ca6b7fc 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -147,28 +147,30 @@ def _calc_rmax(depth): Parameters ---------- depth: float - Bottom depth (units: m). + Bottom depth (units: m). Returns ------- rmax: float - Slope steepness value (units: None) + Slope steepness value (units: None) """ - depth = depth.reset_index(list(depth.dims)) both_rmax = [] for dim in depth.dims: - # (H[0] - H[1]) / (H[0] + H[1]) - depth_diff = depth.diff(dim) - depth_rolling_sum = depth.rolling({dim: 2}).sum().dropna(dim) - rmax = depth_diff / depth_rolling_sum + rolled = depth.rolling({dim: 2}).construct("tmp_dim") + + # |(H[0] - H[1])| / (H[0] + H[1]) + # First value is NaN + diff = rolled.diff("tmp_dim").squeeze("tmp_dim") + rmax = diff / rolled.sum("tmp_dim") - # (R[0] + R[1]) / 2 - rmax = rmax.rolling({dim: 2}).mean().dropna(dim) + # (rmax[0] + rmax[1]) / 2 + # First two values are NaN + rmax = rmax.rolling({dim: 2}).mean() - # Fill first row and column - rmax = rmax.pad({dim: (1, 1)}, constant_values=0) + # First and last values are zero + rmax = rmax.shift({dim: -1}).fillna(0) both_rmax.append(np.abs(rmax)) From 9d928d9757fba02d6923513c662dd67d62ed59f8 Mon Sep 17 00:00:00 2001 From: malmans2 Date: Mon, 21 Jun 2021 12:18:15 +0100 Subject: [PATCH 2/7] implement #39 --- pydomcfg/tests/bathymetry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index ca6b7fc..d276135 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -163,7 +163,7 @@ def _calc_rmax(depth): # |(H[0] - H[1])| / (H[0] + H[1]) # First value is NaN diff = rolled.diff("tmp_dim").squeeze("tmp_dim") - rmax = diff / rolled.sum("tmp_dim") + rmax = np.abs(diff) / rolled.sum("tmp_dim") # (rmax[0] + rmax[1]) / 2 # First two values are NaN @@ -172,6 +172,6 @@ def _calc_rmax(depth): # First and last values are zero rmax = rmax.shift({dim: -1}).fillna(0) - both_rmax.append(np.abs(rmax)) + both_rmax.append(rmax) return np.maximum(*both_rmax) From 9a382a84c2b1a65b7c53ea2336455cd23b02572f Mon Sep 17 00:00:00 2001 From: malmans2 Date: Mon, 21 Jun 2021 14:49:56 +0100 Subject: [PATCH 3/7] replace land with 0s --- pydomcfg/tests/bathymetry.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index d276135..237d4ba 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -155,6 +155,9 @@ def _calc_rmax(depth): Slope steepness value (units: None) """ + # Replace land with zeros + depth = depth.where(depth > 0, 0) + both_rmax = [] for dim in depth.dims: @@ -167,7 +170,7 @@ def _calc_rmax(depth): # (rmax[0] + rmax[1]) / 2 # First two values are NaN - rmax = rmax.rolling({dim: 2}).mean() + rmax = rmax.rolling({dim: 2}).mean(skipna=True) # First and last values are zero rmax = rmax.shift({dim: -1}).fillna(0) From 56f2193d08cb1d4a46ce87a8c871138e0581ca30 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 21 Jun 2021 14:03:22 +0000 Subject: [PATCH 4/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- pydomcfg/tests/bathymetry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index 237d4ba..52d0f6e 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -157,7 +157,7 @@ def _calc_rmax(depth): # Replace land with zeros depth = depth.where(depth > 0, 0) - + both_rmax = [] for dim in depth.dims: From d93b66b9c690940f2046ca9564e5b988b530f0d4 Mon Sep 17 00:00:00 2001 From: malmans2 Date: Mon, 21 Jun 2021 19:19:09 +0100 Subject: [PATCH 5/7] handle land properly --- pydomcfg/domzgr/zgr.py | 3 ++- pydomcfg/tests/bathymetry.py | 27 +++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/pydomcfg/domzgr/zgr.py b/pydomcfg/domzgr/zgr.py index 5910107..d24b2a9 100644 --- a/pydomcfg/domzgr/zgr.py +++ b/pydomcfg/domzgr/zgr.py @@ -20,7 +20,8 @@ def __init__(self, ds_bathy: Dataset, jpk: int): Parameters ---------- ds_bathy: Dataset - xarray dataset including grid coordinates and bathymetry + xarray dataset including grid coordinates and bathymetry. + All bathymetry values MUST be >= 0, where 0 is land. jpk: int number of computational levels """ diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index 52d0f6e..7ce31f2 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -155,25 +155,24 @@ def _calc_rmax(depth): Slope steepness value (units: None) """ - # Replace land with zeros - depth = depth.where(depth > 0, 0) - both_rmax = [] for dim in depth.dims: - rolled = depth.rolling({dim: 2}).construct("tmp_dim") - - # |(H[0] - H[1])| / (H[0] + H[1]) - # First value is NaN - diff = rolled.diff("tmp_dim").squeeze("tmp_dim") - rmax = np.abs(diff) / rolled.sum("tmp_dim") + # |H[0] - H[1]| / (H[0] + H[1]) + rolled = depth.rolling({dim: 2}).construct("window_dim") + diff = rolled.diff("window_dim").squeeze("window_dim") + rmax = np.abs(diff) / rolled.sum("window_dim") # (rmax[0] + rmax[1]) / 2 - # First two values are NaN - rmax = rmax.rolling({dim: 2}).mean(skipna=True) - - # First and last values are zero - rmax = rmax.shift({dim: -1}).fillna(0) + rolled = rmax.rolling({dim: 2}).construct("window_dim") + rmax = rolled.mean("window_dim", skipna=True) + + # 1. Place on the correct index (shift -1 as we rolled twice) + # 2. Force first/last values = 0 + # 3. Replace land values with 0 + rmax = rmax.shift({dim: -1}) + rmax[{dim: [0, -1]}] = 0 + rmax = rmax.fillna(0) both_rmax.append(rmax) From 944226c3579b12c1c5d9d10386cf787835cf79c4 Mon Sep 17 00:00:00 2001 From: malmans2 Date: Fri, 25 Jun 2021 16:45:16 +0100 Subject: [PATCH 6/7] use max rather than mean and use xarray max instead of np maximum --- pydomcfg/tests/bathymetry.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index 7ce31f2..a1f2530 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -142,7 +142,7 @@ def _calc_rmax(depth): Calculate rmax: measure of steepness This function returns the slope steepness criteria rmax, which is simply - (H[0] - H[1]) / (H[0] + H[1]) + |H[0] - H[1]| / (H[0] + H[1]) Parameters ---------- @@ -153,27 +153,31 @@ def _calc_rmax(depth): ------- rmax: float Slope steepness value (units: None) + + Notes + ----- + This function uses a "conservative approach" and rmax is overestimated. + rmax at T points is the maximum rmax estimated at any adjacent U/V point. """ + # Loop over x and y both_rmax = [] for dim in depth.dims: - # |H[0] - H[1]| / (H[0] + H[1]) + # Compute rmax rolled = depth.rolling({dim: 2}).construct("window_dim") diff = rolled.diff("window_dim").squeeze("window_dim") rmax = np.abs(diff) / rolled.sum("window_dim") - # (rmax[0] + rmax[1]) / 2 - rolled = rmax.rolling({dim: 2}).construct("window_dim") - rmax = rolled.mean("window_dim", skipna=True) - - # 1. Place on the correct index (shift -1 as we rolled twice) - # 2. Force first/last values = 0 - # 3. Replace land values with 0 + # Construct dimension with velocity points adjacent to any T point + # We need to shift as we rolled twice and to force boundaries = NaN + rmax = rmax.rolling({dim: 2}).construct("vel_points") rmax = rmax.shift({dim: -1}) - rmax[{dim: [0, -1]}] = 0 - rmax = rmax.fillna(0) + rmax[{dim: [0, -1]}] = None both_rmax.append(rmax) - return np.maximum(*both_rmax) + rmax = xr.concat(both_rmax, "vel_points") + rmax = rmax.max("vel_points", skipna=True) + + return rmax.fillna(0) From 6e8ed5ec724e4fc9fc0e1f2084ec8f989426fc45 Mon Sep 17 00:00:00 2001 From: malmans2 Date: Sat, 26 Jun 2021 11:33:59 +0100 Subject: [PATCH 7/7] same as pyroms --- pydomcfg/tests/bathymetry.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/pydomcfg/tests/bathymetry.py b/pydomcfg/tests/bathymetry.py index a1f2530..ea2bcd0 100644 --- a/pydomcfg/tests/bathymetry.py +++ b/pydomcfg/tests/bathymetry.py @@ -160,6 +160,9 @@ def _calc_rmax(depth): rmax at T points is the maximum rmax estimated at any adjacent U/V point. """ + # Mask land + depth = depth.where(depth > 0) + # Loop over x and y both_rmax = [] for dim in depth.dims: @@ -170,14 +173,18 @@ def _calc_rmax(depth): rmax = np.abs(diff) / rolled.sum("window_dim") # Construct dimension with velocity points adjacent to any T point - # We need to shift as we rolled twice and to force boundaries = NaN + # We need to shift as we rolled twice rmax = rmax.rolling({dim: 2}).construct("vel_points") rmax = rmax.shift({dim: -1}) - rmax[{dim: [0, -1]}] = None both_rmax.append(rmax) + # Find maximum rmax at adjacent U/V points rmax = xr.concat(both_rmax, "vel_points") rmax = rmax.max("vel_points", skipna=True) + # Mask halo points + for dim in rmax.dims: + rmax[{dim: [0, -1]}] = 0 + return rmax.fillna(0)