make continent masking optional for non-modern land-sea distribution?

[christian-stepanek]

• pyfesom2 version: 0.1.0 (as per init.py)
• Python version: 3.10
• Operating System: Linux 4.18.0-372.13.1.el8_6.x86_64 GNU/Linux

Description

When plotting FESOM2 output for meshes where the land sea mask differs significantly from the modern one, masking of continents that is applied per standard based on Natural Earth (defined in ut.py) hides data where model's and Natural Earth's land sea masks differ (ocean in model vs. land in Natural Earth).

Masking of the land sea mask is at the moment hard-coded in functions get_vector_forplot and plot of file pyfesom2/plotting.py as far as I see, potentially in other functions as well.

    302     m2 = mask_ne(lonreg2, latreg2)
    303     u_int = u_int[0]
    304     v_int = v_int[0]
    305 
    306     u_int = np.ma.masked_where(m2, u_int)
    307     u_int = np.ma.masked_equal(u_int, 0)
    308 
    309     v_int = np.ma.masked_where(m2, v_int)
    310     v_int = np.ma.masked_equal(v_int, 0)
    311 
    312     return u_int, v_int, lonreg2, latreg2

This results in the attached figure midPliocene_surface_velocity.png, where data is hidden in regions that are land in Natural Earth but ocean in the mesh; for reference please also see attached figure midPliocene_surface_velocity_not_masked.png that results from commenting code that applies the masking and provides the full information from the simulation. Note the difference in plotting parameters, only the difference in data coverage in the figure is relevant to this issue.

First question: Is there necessity to do the masking at all? (For example towards deriving correct interpolation results at the coast line?). Or is the masking "just" towards plot quality?
Second question: Doesn't the mesh itself contain all information relevant to do the masking? If so, could the masking be easily done based on the mesh itself? This would make the plotting code more versatile and more robust in case of plotting data for modified meshes.

What I Did

Standard use of pyfesom2 data reading, interpolation, and plotting routines:

import pyfesom2 as pf
import matplotlib.pylab as plt
import xarray as xr
meshpath='Pliocene_3Ma/fesom2/midpli2/'
alpha, beta, gamma=[0, 0, 0]
mesh=pf.load_mesh(meshpath, abg=[alpha, beta, gamma], usepickle = False, usejoblib=True)
datapath = 'midPliocene/outdata/fesom/'
u = pf.get_data(datapath, 'u', [2249], mesh, depth = 0)
v = pf.get_data(datapath, 'v', [2249], mesh, depth = 0)
u_int, v_int, lonreg2, latreg2 = pf.get_vector_forplot(u,v, mesh)
box=[-180, 180, -90, -60]
u_int, v_int, lonreg2, latreg2 = pf.get_vector_forplot(u,v, mesh, box=box, res = [360, 20*10])
pf.plot_vector(u_int,  v_int, lonreg2, latreg2, box, mapproj='sp', sstep=1, scale=3, coastlines=False)
plt.savefig("midPliocene_surface_velocity.png")

[koldunovn]

Hi @christian-stepanek There is a possibility to turn off NE masking in the latest version of the code with no_pi_mask argument:

    m2 = mask_ne(lonreg2, latreg2)
    u_int = u_int[0]
    v_int = v_int[0]

    if not no_pi_mask:
        u_int = np.ma.masked_where(m2, u_int)
    u_int = np.ma.masked_equal(u_int, 0)

    if not no_pi_mask:
        v_int = np.ma.masked_where(m2, v_int)
    v_int = np.ma.masked_equal(v_int, 0)

    return u_int, v_int, lonreg2, latreg2

same is possible for plot function.

It's not necessary to do the masking - basically it is there just to make figures look nice. As most of the time we are doing nearest neighbour, and "radius of influence" can be quite substantial, resulting in ugly coastlines we have this feature, but one can turn it off.

Mesh itself contain the infrormation, but one have to do much slower interpolation to do it right (e.g. here https://github.com/FESOM/fint/blob/6ab4ddb47d7fcefade3fefaef7cfa3369ec593a5/src/fint/fint.py#L222). We plan to extend the pyfesom2 with methods from fint, but the timeline for it is not clear :(

[christian-stepanek]

Dear @koldunovn, thanks a lot, your suggestions are highly appreciated. Thanks a lot for developing pyfesom2.