-
Notifications
You must be signed in to change notification settings - Fork 11
Description
- 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")

