import pyfesom2 as pf
import xarray as xr
import numpy as np
import matplotlib.cm as cm

mesh_path = '/path/to/meshes/SO3/'
mesh = pf.load_mesh(mesh_path)
diag = xr.open_dataset(f'{mesh_path}/fesom.mesh.diag.nc')
elem_area = diag.elem_area
resolution = np.sqrt((elem_area*2))/1.09/1e3
resolution_nodes = pf.tonodes(resolution, mesh)
pf.plot(mesh, resolution_nodes, levels=(2.5, 25, 19), cmap=cm.magma_r, units='km')