Do not drop indexes when computing rmax #40

malmans2 · 2021-06-21T09:17:15Z

Closes calc_rmax: Should we compute the rolling average using absolute values? #39
Closes testing calc_rmax #48
Passes pre-commit run --all-files
Project, label, and assignee tabs are populated

codecov · 2021-06-21T09:18:45Z

Codecov Report

Merging #40 (6e8ed5e) into main (f68612d) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##             main      #40   +/-   ##
=======================================
  Coverage   94.85%   94.85%           
=======================================
  Files           5        5           
  Lines         214      214           
=======================================
  Hits          203      203           
  Misses         11       11

Flag	Coverage Δ
unittests	`94.85% <ø> (ø)`

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
pydomcfg/domzgr/zgr.py	`100.00% <ø> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f68612d...6e8ed5e. Read the comment docs.

pydomcfg/tests/bathymetry.py

malmans2 · 2021-06-21T11:27:09Z

Ready for review!

malmans2 · 2021-06-21T11:40:41Z

@oceandie @jdha moving the discussion here.
What is you preference for land points (0/0):

Return inf
Return NaN
Return 0

Also what about the rolling average. I.e., what should be the output of coastal points. Say val0=0.1, val1=0/0, the rolling average should return:

0.1
0.1 / 2
Others? (0, or NaN, or inf, ....)?

jdha · 2021-06-21T11:54:26Z

Sorry, I'm probably a little slow on the uptake but re land points: a DS is passed (and therefore a mask) so the land isn't used in the calc, is that correct?

Re: coastal point. If we're adopting @malmans2 solution, then the output at coastal points is just the maximum rmax of the surrounding u/v points, no? (u/v being zero adjacent to land)

oceandie · 2021-06-21T13:13:34Z

Sorry, I'm probably a little slow on the uptake but re land points: a DS is passed (and therefore a mask) so the land isn't used in the calc, is that correct?

Re: coastal point. If we're adopting @malmans2 solution, then the output at coastal points is just the maximum rmax of the surrounding u/v points, no? (u/v being zero adjacent to land)

Right now the input is a DataArray (DA) and we don't specify anywhere that has to be masked. Also, creating the mask inside the function, doesn't really work, since the input DA not necesserely correspond to the original model bathymetry (i.e. original land-sea mask) : for example NEMO set first land point adjacent to a wet cell to min_dep before the smoothing as this needs to be included in smoothing (see here)

I think that even implementing #39 we still have division by zero on the land, which will give nans and then will create a problem since np.maximum propagates nans, isn't it?

malmans2 · 2021-06-21T13:22:55Z

I think @jdha is right, but currently we are not expecting users to pass a mask along with Bathymetry (I wasn't sure if we should expect mask from users, so the Bathymetry adds a mask but it is never used).

So I think we just have to add depth = depth.where(depth > 0, 0) at the beginning, which infers sea points and replaces all land points with zeros.

These are the steps for each dimensions of the current implementation:

Move to U/V points: rmax_vel = abs(H[1] - H[0]) / (H[0] + H[1]) -> 0/0 gives NaN
Move back to T points: rmax = mean(rmax_vel[0], rmax_vel[1]) -> By default, skipna=True, so mean([1, NaN]) == 1
replace NaN with 0 (NaN are either land points or first/last row/column)

Finally, rmax = max(rmax_x, rmax_y)

My question was mainly about how to handle 2 in this comment, when we move back to T points. I think there are cases where it would make a difference if we replace NaN with 0. So mean(1, NaN) != mean(1, 0), which could affect the final result near land. It makes more sense to me to use NaN, but I'm not 100% sure. Does it make sense?

malmans2 · 2021-06-21T13:33:40Z

These are the steps for each dimensions of the current implementation:

With current implementation I mean this PR... I think I see now the difference in this implementation that fixes issues with land:

pyDOMCFG/pydomcfg/tests/bathymetry.py

Line 173 in 9d928d9

rmax = rmax.shift({dim: -1}).fillna(0)

We are now filling all nan with zeros (i.e., land and first/last row/column), before we were just padding the boundaries with zeros.

jdha · 2021-06-21T13:34:23Z

I think @jdha is right, but currently we are not expecting users to pass a mask along with Bathymetry (I wasn't sure if we should expect mask from users, so the Bathymetry adds a mask but it is never used).

So I think we just have to add depth = depth.where(depth > 0, 0) at the beginning, which infers sea points and replaces all land points with zeros.

These are the steps for each dimensions of the current implementation:

Move to U/V points: rmax_vel = abs(H[1] - H[0]) / (H[0] + H[1]) -> 0/0 gives NaN

Move back to T points: rmax = mean(rmax_vel[0], rmax_vel[1]) -> By default, skipna=True, so mean([1, NaN]) == 1

replace NaN with 0 (NaN are either land points or first/last row/column)

Finally, rmax = max(rmax_x, rmax_y)

My question was mainly about how to handle 2 in this comment, when we move back to T points. I think there are cases where it would make a difference if we replace NaN with 0. So mean(1, NaN) != mean(1, 0), which could affect the final result near land. It makes more sense to me to use NaN, but I'm not 100% sure. Does it make sense?

@malmans2 I would definitely use NaN, as that value/point is never used to calculate pressure gradients in NEMO, so is of no interest.

oceandie · 2021-06-21T13:40:17Z

Since we are not dropping nans anymore an given the default behaviour of mean (which solves the problem with maximum I was mentioning before), yes I would agree to leave nan to avoind weird results with zeros is the way to go.

oceandie · 2021-06-21T13:44:35Z

However, I still think the input should be a DA, not a DS., what do you think?

malmans2 · 2021-06-21T13:53:03Z

However, I still think the input should be a DA, not a DS., what do you think?

Yes, the bathymetry is the only DataArray needed to compute rmax. The mask is inferred from positive values, and anything that is land (not positive) is replaced with zeros.

I also added skipna=True to make explicit what's going on.

malmans2 · 2021-06-21T14:02:42Z

pre-commit.ci autofix

for more information, see https://pre-commit.ci

oceandie · 2021-06-21T14:02:57Z

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)


Do we really need this here? I think when bathy or envelope are passed to this function they should have already zeros for land

What happen if the user bathymetry has negative values or NaN on land?
Should we remove it and clarify in the documentation that Bathymetry must be >= 0?

It's not a bad idea to get rid of this because it is potentially a very expensive but also useless operation...

I would say yes, since I think we should write a function that put all land values to zero when interpolating the input bathymetry onto the model grid. So the model bathymetry will be positive defined by definition, with zeros for land (also I think NEMO will fail if this conditions are not met, but not sure about this last point)

Sorry, I've been dipping in and out of this today (so no continuity in my thought process). @oceandie re NEMO will fail, NEMO as far as I understand doesn't use bathymetry field, just scale factors and wet cells [i.e. if top_level is equal to zero it's land]. So in that sense we can define land values in bathymetry as we see fit. I'll look over the python again tomorrow (and you may have already included some of this), but the key user defined vars are rn_sbot_min (for ln_sco) and rn_hmin (for ln_zco, ln_zps; which take different meanings depending on whether they are -ve or +ve: min number of levels or min depth). So when calculating rmax this user choices will affect the outcome (although, not in the case of tests).

from the Fortran (where the bathy is updated):

IF ( .not. ln_sco ) THEN !== set a minimum depth ==! IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth ENDIF zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans END WHERE IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik ENDIF

and

bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) DO jj = 1, jpj DO ji = 1, jpi IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) END DO END DO

note there is also a rn_sbot_max

Apologies if I'm going over old ground...

HI @jdha, yes sorry, I was still thinking about NEMO3.6 or DOMCFG, that we are actually trying to replicate :)
Regarding all the params and bits of code you mentioned, yes I started to consider them - soon I will push where I am now in sco_dev , then any feedback is more than welcome ;)

OK, this is helpful! (maybe open an issue about this as this comment will get kind of lost when we merge this PR?)

This PR should be just the starting point, and _calc_rmax will become the general function to compute rmax. Once we merge this PR, I think @oceandie will merge main into #33 and will move the whole _calc_rmax function into utils. I guess vcoord-specific parameters such as rn_sbot_min will be handled internally by our domzgr classes (but all classes will make use of the same _calc_rmax function).

Hi @malmans2 and @jdha , how can I see the original code for calc_rmax by @jdha (the one using numpy) ?

This PR started from the xarray version. You'd have to look in James' PR (already merged). Should be this one: #17

I'd go there and look at the appropriate commit (e.g., click on commits and select the first commit).

pydomcfg/tests/bathymetry.py

malmans2

I think it's good to go now... I changed stuff back and forward quite a bit, so it's probably better if you compare the changes from all commits (button at the top of this page if you didn't use it before).

pydomcfg/tests/bathymetry.py

malmans2 · 2021-06-26T10:38:56Z

This is now using the same implementation in pyroms

Do not drop indexes when computing rmax

34f32e4

malmans2 added the enhancement New feature or request label Jun 21, 2021

malmans2 requested a review from jdha June 21, 2021 09:17

malmans2 added this to In progress in domzgr.F90 -> pyDOMCFG via automation Jun 21, 2021

malmans2 commented Jun 21, 2021

View reviewed changes

pydomcfg/tests/bathymetry.py Outdated Show resolved Hide resolved

implement pyNEMO#39

9d928d9

oceandie mentioned this pull request Jun 21, 2021

Developing sco class #33

Draft

6 tasks

replace land with 0s

9a382a8

[pre-commit.ci] auto fixes from pre-commit.com hooks

56f2193

for more information, see https://pre-commit.ci

oceandie reviewed Jun 21, 2021

View reviewed changes

malmans2 commented Jun 21, 2021

View reviewed changes

pydomcfg/tests/bathymetry.py Outdated Show resolved Hide resolved

malmans2 and others added 3 commits June 21, 2021 19:19

handle land properly

d93b66b

Merge branch 'pyNEMO:main' into rmax_indexes

f0ae17f

use max rather than mean and use xarray max instead of np maximum

944226c

malmans2 commented Jun 25, 2021

View reviewed changes

pydomcfg/tests/bathymetry.py Show resolved Hide resolved

pydomcfg/tests/bathymetry.py Show resolved Hide resolved

malmans2 mentioned this pull request Jun 26, 2021

testing calc_rmax #48

Closed

same as pyroms

6e8ed5e

domzgr.F90 -> pyDOMCFG automation moved this from In progress to Reviewer approved Jun 27, 2021

oceandie approved these changes Jun 27, 2021

View reviewed changes

oceandie merged commit 9002bb0 into pyNEMO:main Jun 27, 2021

domzgr.F90 -> pyDOMCFG automation moved this from Reviewer approved to Done Jun 27, 2021

malmans2 deleted the rmax_indexes branch June 30, 2021 07:12

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Do not drop indexes when computing rmax #40

Do not drop indexes when computing rmax #40

malmans2 commented Jun 21, 2021 •

edited

codecov bot commented Jun 21, 2021 •

edited

malmans2 commented Jun 21, 2021

malmans2 commented Jun 21, 2021 •

edited

jdha commented Jun 21, 2021

oceandie commented Jun 21, 2021

malmans2 commented Jun 21, 2021

malmans2 commented Jun 21, 2021 •

edited

jdha commented Jun 21, 2021

oceandie commented Jun 21, 2021

oceandie commented Jun 21, 2021

malmans2 commented Jun 21, 2021 •

edited

malmans2 commented Jun 21, 2021

oceandie Jun 21, 2021

malmans2 Jun 21, 2021 •

edited

oceandie Jun 21, 2021 •

edited

jdha Jun 21, 2021 •

edited by malmans2

oceandie Jun 21, 2021

malmans2 Jun 21, 2021

oceandie Jun 25, 2021

malmans2 Jun 25, 2021

oceandie Jun 25, 2021

malmans2 left a comment

malmans2 commented Jun 26, 2021 •

edited

Do not drop indexes when computing rmax #40

Do not drop indexes when computing rmax #40

Conversation

malmans2 commented Jun 21, 2021 • edited

codecov bot commented Jun 21, 2021 • edited

Codecov Report

malmans2 commented Jun 21, 2021

malmans2 commented Jun 21, 2021 • edited

jdha commented Jun 21, 2021

oceandie commented Jun 21, 2021

malmans2 commented Jun 21, 2021

malmans2 commented Jun 21, 2021 • edited

jdha commented Jun 21, 2021

oceandie commented Jun 21, 2021

oceandie commented Jun 21, 2021

malmans2 commented Jun 21, 2021 • edited

malmans2 commented Jun 21, 2021

oceandie Jun 21, 2021

Choose a reason for hiding this comment

malmans2 Jun 21, 2021 • edited

Choose a reason for hiding this comment

oceandie Jun 21, 2021 • edited

Choose a reason for hiding this comment

jdha Jun 21, 2021 • edited by malmans2

Choose a reason for hiding this comment

oceandie Jun 21, 2021

Choose a reason for hiding this comment

malmans2 Jun 21, 2021

Choose a reason for hiding this comment

oceandie Jun 25, 2021

Choose a reason for hiding this comment

malmans2 Jun 25, 2021

Choose a reason for hiding this comment

oceandie Jun 25, 2021

Choose a reason for hiding this comment

malmans2 left a comment

Choose a reason for hiding this comment

malmans2 commented Jun 26, 2021 • edited

malmans2 commented Jun 21, 2021 •

edited

codecov bot commented Jun 21, 2021 •

edited

malmans2 commented Jun 21, 2021 •

edited

malmans2 commented Jun 21, 2021 •

edited

malmans2 commented Jun 21, 2021 •

edited

malmans2 Jun 21, 2021 •

edited

oceandie Jun 21, 2021 •

edited

jdha Jun 21, 2021 •

edited by malmans2

malmans2 commented Jun 26, 2021 •

edited