(0.5.0) Better experiment.bathymetry method

[navidcy]

use kwargs

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[navidcy]

@ashjbarnes, what is the purpose of make_topography = True?
all experiment.bathymetry seems to be indented into an if make_topography: which beats the purpose?

[navidcy]

why multiplying by -1 ensures that? there is an assumption that it wasn't already positive down which implies that the user was careful enough, right?

[ashjbarnes]

No the user can provide their bathymetry without having to pre-multiply it. Passing the positive_down argument just specifies what the raw data is and the code adjusts by multiplying by -1 or not. gebco for example has positive values for land and neg for ocean whereas MOM6 wants depth to be positive

[navidcy]

My point was that this implies that the user knows what “positive down” is. We can alternatively have the code here test to figure out whether the provided bathymetry needs to be multiplied w -1 or not and do that and print an info message to the user?

[navidcy]

@ashjbarnes question: Is this method here an old name? make_bathymetry?

regional-mom6/regional_mom6/regional_mom6.py

Line 1287 in 4ceaf39

print("No topography file! Need to run make_bathymetry first")

[ashjbarnes]

@ashjbarnes, what is the purpose of make_topography = True? all experiment.bathymetry seems to be indented into an if make_topography: which beats the purpose?

yeah that was an old thing! I'm fixing it now

[ashjbarnes]

what does the * do here? I'm not familiar with this notation

    def setup_bathymetry(
        self,
        *,
        bathymetry_path,
        coordinate_names,
        fill_channels=False,
        minimum_layers=3,
        positive_down=False,
        chunks="auto",
    ):

[ashjbarnes]

If @navidcy is happy with the changes I've made then this is good to go. I've switched to calling everything 'bathymetry', added more docstrings, removed redundant make_topog argument and switched the coordinate names dict to kwargs instead

[navidcy]

what does the * do here? I'm not familiar with this notation

    def setup_bathymetry(
        self,
        *,
        bathymetry_path,
        coordinate_names,
        fill_channels=False,
        minimum_layers=3,
        positive_down=False,
        chunks="auto",
    ):

it splits the args from kwargs; self needs to be an argument (similarly it __init__)

[navidcy]

@ashjbarnes, I'd like to automate this multiplication and remove the positive_down kwarg. I believe (but I'm not sure that it's true) that there should be a way so that the user can be completely naïve about this positive up/down situation and whatever forcing the user provide the package deals with it and converts the vertical coordinate to how MOM6 expects it to be. Could we do that?

Could you give me an example of a vertical coordinate that it is positive_down and one that it is not?

[ashjbarnes]

@ashjbarnes, I'd like to automate this multiplication and remove the positive_down kwarg. I believe (but I'm not sure that it's true) that there should be a way so that the user can be completely naïve about this positive up/down situation and whatever forcing the user provide the package deals with it and converts the vertical coordinate to how MOM6 expects it to be. Could we do that?

Could you give me an example of a vertical coordinate that it is positive_down and one that it is not?

That's what I'd hoped originally, but what ended up with me accidentally inverting the bathymetry! If the user gives the bathymetry for an existing MOM6 run, then it will have depths of 0 -> 6000. If you give it GEBCO it'll have depths of -6000-> 0 and positive values correspond to land. What do you think is a robust way for the model to assume correctly the orientation of the vertical coordinate?

[navidcy]

The vertical coordinate could be one of:

• [0, 10, ..., 6000] (eg MOM6)
• [0, -10, ..., -6000] (eg GEBCO?)
• [-6000, ..., -10, 0] (eg Oceananigans)
• [6000, ..., 10, 0] (I have no example in mind)

Whatever the vertical coordinate convention is we need to convert it to MOM6 convention, right?

[navidcy]

@ashjbarnes what do you think of the following?

import numpy as np

def convert_to_MOM6_convention(vcoord):
    '''
    Converts the vertical coordinate (``vcoord``) to have the convention expected by
    MOM6, that is positive values that depict the distance below sea level and with
    ascending order. For example ``vcoord = [0, 10, 50, 100]``.
    '''

    # for single vertical layer config don't do anything
    if len(vcoord)==1:
        return vcoord

    # ensure vcoord is numpy.array
    vcoord = np.array(vcoord)
    
    # if in descending order then flip them
    if vcoord[1] > vcoord[0]:
        ascending = True
    elif vcoord[1] < vcoord[0]:
        return convert_to_MOM6_convention(np.flip(vcoord))
    else:
        raise ValueError("vcoord seems constant...")

    # if values are negative, convert to positive
    if vcoord[1] > 0 and vcoord[0] >= 0:
        is_positive = True
    elif vcoord[1] < 0 and vcoord[0] <= 0:
        return convert_to_MOM6_convention(-1 * vcoord)
    else:
        raise ValueError("vcoord has both positive and negative values...")

    # if vcoord follows MOM6 convention then all is well
    if is_positive and ascending:
        return vcoord

When I try it out:

# test that convert_to_MOM6 does what is supposed to do
vcoords = ([0, 1, 2, 10],
           [0, -1, -2, -10],
           [10, 2, 1, 0],
           [-10, -2, -1, 0],
           np.array([0, 1, 2, 10]),
           np.array([0, -1, -2, -10]),
           np.array([10, 2, 1, 0]),
           np.array([-10, -2, -1, 0]),
           [0.2, 1, 2, 10],
           [-0.2, -1, -2, -10],
           [10, 2, 1, 0.2],
           [-10, -2, -1, -0.2],
)

for vcoord in vcoords:
    print(convert_to_MOM6_convention(vcoord))

gives

[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]

demonstrating that it does return things how MOM6 expect them?

[navidcy]

In MOM6 what values does land have?

[navidcy]

I decided to merge this and potentially deal with the positive_down kwarg and its potential automation in a different PR.

[ashjbarnes]

@ashjbarnes what do you think of the following?

import numpy as np

def convert_to_MOM6_convention(vcoord):
    '''
    Converts the vertical coordinate (``vcoord``) to have the convention expected by
    MOM6, that is positive values that depict the distance below sea level and with
    ascending order. For example ``vcoord = [0, 10, 50, 100]``.
    '''

    # for single vertical layer config don't do anything
    if len(vcoord)==1:
        return vcoord

    # ensure vcoord is numpy.array
    vcoord = np.array(vcoord)
    
    # if in descending order then flip them
    if vcoord[1] > vcoord[0]:
        ascending = True
    elif vcoord[1] < vcoord[0]:
        return convert_to_MOM6_convention(np.flip(vcoord))
    else:
        raise ValueError("vcoord seems constant...")

    # if values are negative, convert to positive
    if vcoord[1] > 0 and vcoord[0] >= 0:
        is_positive = True
    elif vcoord[1] < 0 and vcoord[0] <= 0:
        return convert_to_MOM6_convention(-1 * vcoord)
    else:
        raise ValueError("vcoord has both positive and negative values...")

    # if vcoord follows MOM6 convention then all is well
    if is_positive and ascending:
        return vcoord

When I try it out:

# test that convert_to_MOM6 does what is supposed to do
vcoords = ([0, 1, 2, 10],
           [0, -1, -2, -10],
           [10, 2, 1, 0],
           [-10, -2, -1, 0],
           np.array([0, 1, 2, 10]),
           np.array([0, -1, -2, -10]),
           np.array([10, 2, 1, 0]),
           np.array([-10, -2, -1, 0]),
           [0.2, 1, 2, 10],
           [-0.2, -1, -2, -10],
           [10, 2, 1, 0.2],
           [-10, -2, -1, -0.2],
)

for vcoord in vcoords:
    print(convert_to_MOM6_convention(vcoord))

gives

[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0  1  2 10]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]
[ 0.2  1.   2.  10. ]

demonstrating that it does return things how MOM6 expect them?

Gebco is more like [-5000...0...5000] How do we automatically tell which bits are land and which are ocean? That's the issue. I also don't think that whether it's 0,1,2,3 or 3,2,1,0 matters since it gets interpolated onto the base grid anyway

[navidcy]

Let's continue the discussion about potential elimination of positive_down kwarg in #150