testing calc_rmax

[oceandie]

I found some inconsistencies in the way calc_rmax works. I've developed a quick python code for:

1. read AMM7 bathymetry;
2. create an envelope bathymetry and smooth it: it implements exactly NEMO Fortran90 code (very inefficient but at least we know what it does)
3. compute the maximum rmax with three different codes: original numpy version of James; xarray version of James code; my old code (actually adapted from pyROMS, it very inefficient but I think it is clear what it does).

If we smooth the envelope with an rn_rmax threshold of 0.24, the following is the pring of the iterative smoothing:
python test_rmax_calc.py
Iter: 1 rmax: 1.0
Iter: 2 rmax: 0.9353992591835926
Iter: 3 rmax: 0.9187177073344386
Iter: 4 rmax: 0.8765126597150313
Iter: 5 rmax: 0.8185786808404448
Iter: 6 rmax: 0.7347370279088926
Iter: 7 rmax: 0.7347370279088925
Iter: 8 rmax: 0.6074021541852422
Iter: 9 rmax: 0.5717732932642535
Iter: 10 rmax: 0.4755012414806615
Iter: 11 rmax: 0.37086087298353765
Iter: 12 rmax: 0.24000000000000007

The rmax computed by the three codes are the following:

1. calc_rmax_xr: 0.34985075453036585 (xarray code)
2. calc_rmax_np: 0.34985075453036585 (numpy code)
3. SlopeParam: 0.24000000000000005 (ROMS code)

[jdha]

Probably should drop this note here too:

To reproduce your SlopeParam code @oceandie [which looks very much like my old Matlab ;) ] I did the following (see below). The large values seem to appear in the outer halo as they where included in the np.diff (your code only goes from the second index to the last but one) ... I also think I'd got the indices out of sync by 1 - but I'm not too sure ... I'll look again on Monday (and do an XR version of this).

If I diff the the results from the following code with SlopeParam they match.

def calc_rmax_np(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])
    Parameters
    ----------
    depth: float
            Bottom depth (units: m).
    Returns
    -------
    rmax: float
            Slope steepness value (units: None)
    """
    
    rmax_x, rmax_y = np.zeros_like(depth), np.zeros_like(depth)

    rmax_x[1:-1, 1:-1] = np.maximum (
        np.abs( np.diff(depth[1:-1, :-1], axis=1) / (depth[1:-1,  :-2] + depth[1:-1, 1:-1]) ),
        np.abs( np.diff(depth[1:-1,1:  ], axis=1) / (depth[1:-1, 1:-1] + depth[1:-1, 2:  ]) )
    )
    rmax_y[1:-1, 1:-1] = np.maximum (
        np.abs( np.diff(depth[ :-1, 1:-1], axis=0) / (depth[ :-2, 1:-1] + depth[1:-1, 1:-1]) ) ,
        np.abs( np.diff(depth[1:  , 1:-1], axis=0) / (depth[1:-1, 1:-1] + depth[2:,   1:-1]) )
    )

    return np.maximum(np.abs(rmax_x), np.abs(rmax_y))

[jdha]

OK - I had to return quickly to look at the xr - just to clear my head! This does the same as the np, but I'm not too sure why the outer halo -> rmax[0,:], rmax[-1,:], rmax[:,0], rmax[:,-1] - still has values in. Diffing rmax[1:-1,1:-1] with the np and SlopeParam gives the same results

@malmans2 perhaps it has something to do with
rmax[{dim: [0, -1]}] = 0
inside the dim loop and that it should really be applied to:
np.maximum(*both_rmax)
? Sorry - still learning the finer details of xr

def calc_rmax_xr(depth):
    """
    Calculate rmax: measure of steepness
    This function returns the maximum slope paramater

    rmax = abs(Hb - Ha) / (Ha + Hb)

    where Ha and Hb are the depths of adjacent grid cells (Mellor et al 1998).

    Reference:
    *) Mellor, Oey & Ezer, J Atm. Oce. Tech. 15(5):1122-1131, 1998.

    Parameters
    ----------
    depth: DataArray
        Bottom depth (units: m).

    Returns
    -------
    DataArray
        2D maximum slope parameter (units: None)

    """
    both_rmax = []

    for dim in depth.dims:

        # |Ha - Hb| / (Ha + Hb)
        rolled = depth.rolling({dim: 2}).construct("window_dim")
        # Ha - Hb: diff is computed at U/V points
        diff = rolled.diff("window_dim").squeeze("window_dim")
        # rmax is computed at U/V points
        rmax = np.abs(diff) / rolled.sum("window_dim")
        # find max either side of t-point
        rolled = rmax.rolling({dim: 2}).construct("window_dim")
        rmax = rolled.max("window_dim", skipna=True)
        
        # 1. Place on the correct index (shift -1)
        # 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)

    return np.maximum(*both_rmax)

[jdha]

@oceandie if you're happy with these solutions, I'll push back to the branch

[malmans2]

Looks like pyroms is also using the "conservative approach" we discussed yesterday. That's already implemented so you just have to use the latest version of #40.

There are two other minor sources of error at the boundaries and near coast, although I'm not convinced pyroms method is better (see notebook below)... Anyways, as @jdha mentioned we just have to mask land and boundaries.

Here is the code to compare and reproduce pyroms implementation: https://nbviewer.jupyter.org/gist/malmans2/f6f1ee945fd0c5125b87cef275e86e2b

[oceandie]

Hi @jdha , thanks for looking at it!

Glad the SlopeParam is similar to your matlab code, but as I specified in the doc of the function and in the description of the issue that is NOT my original function, I adapted it (very little) from pyROMS (if you are curious, see here ).

So if I undersatnd correctly, there are two changes that are needed to replicate SlopeParam with calc_rmax_np:

1. substitute the average with the maximum (conservative approach)
2. exclude the outer halo in the computation.

[oceandie]

Looks like pyroms is also using the "conservative approach" we discussed yesterday. That's already implemented so you just have to use the latest version of #40.

There are two other minor sources of error at the boundaries and near coast, although I'm not convinced pyroms method is better (see notebook below)... Anyways, as @jdha mentioned we just have to mask land and boundaries.

Here is the code to compare and reproduce pyroms implementation: https://nbviewer.jupyter.org/gist/malmans2/f6f1ee945fd0c5125b87cef275e86e2b

I think that more importantly, pyROMS code and our "conservative approach" are consistent with the way NEMO implement the smoothing algorithm ...

[jdha]

Hi @jdha , thanks for looking at it!

Glad the SlopeParam is similar to your matlab code, but as I specified in the doc of the function and in the description of the issue that is NOT my original function, I adapted it (very little) from pyROMS (if you are curious, see here ).

So if I undersatnd correctly, there are two changes that are needed to replicate SlopeParam with calc_rmax_np:

1. substitute the average with the maximum (conservative approach)
2. exclude the outer halo in the computation.

for calc_rmax_np, I re-wrote it above - and satisfies 1 & 2

[oceandie]

Hi @jdha , thanks for looking at it!
Glad the SlopeParam is similar to your matlab code, but as I specified in the doc of the function and in the description of the issue that is NOT my original function, I adapted it (very little) from pyROMS (if you are curious, see here ).
So if I undersatnd correctly, there are two changes that are needed to replicate SlopeParam with calc_rmax_np:

1. substitute the average with the maximum (conservative approach)
2. exclude the outer halo in the computation.

for calc_rmax_np, I re-wrote it above - and satisfies 1 & 2

Yes, yes I saw it ... was just a recap of the diffs ;)

[oceandie]

Looks like pyroms is also using the "conservative approach" we discussed yesterday. That's already implemented so you just have to use the latest version of #40.

There are two other minor sources of error at the boundaries and near coast, although I'm not convinced pyroms method is better (see notebook below)... Anyways, as @jdha mentioned we just have to mask land and boundaries.

Here is the code to compare and reproduce pyroms implementation: https://nbviewer.jupyter.org/gist/malmans2/f6f1ee945fd0c5125b87cef275e86e2b

Regarding NaNs on the land:

I think pyROMS and NEMO approaches are more correct - we don't want to include land points in the computation since

1. as @jdha pointed out they are never used in computing hpg
2. the smoothing algorithm is an iterative process, if we include also land points they will sligthly change also the points close to the coast. For the same reason, NEMO code, before smoothing, change from land to wet all the land points adjacent to coastal wet cells, see this very slow loop ...