In [1]:
import sys
sys.path.append("../")
import pyfesom as pf
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
import numpy as np
#%matplotlib notebook
%matplotlib inline
from matplotlib import cm
from netCDF4 import Dataset, MFDataset
from matplotlib.colors import Normalize
First, as usuall load the mesh:
In [2]:
meshpath ='../../swift.dkrz.de/COREII'
mesh = pf.load_mesh(meshpath, usepickle=True)
Load data for one year:
In [3]:
fl = Dataset('../../swift.dkrz.de/COREII_data/fesom.1951.oce.mean.nc')
In [4]:
fl.variables['u'].shape
Out[4]:
Vector data in FESOM are internally treated on the rotated mesh (it does not matter if the mesh you work with is "rotated" or not). So you have to convert vector data from rotated coordinates to geographical coordinates. We use function vec_rotate_r2g for this. The angles 50, 15, -90 in practice do not change, so it is save to use them.
In [5]:
level_data_u, elem_no_nan = pf.get_data(fl.variables['u'][0,:],mesh,100)
level_data_v, elem_no_nan = pf.get_data(fl.variables['v'][0,:],mesh,100)
u, v = pf.vec_rotate_r2g(50, 15, -90, mesh.x2, mesh.y2, level_data_u, level_data_v, flag=1)
So now, the u and v are the vectors fields in geographical coordinates that we have to work with.
For details see Basemap documentation. Note, that we rotate vectors to the map projection with rotate_vector. Parameter to play with to make the size of arrows better is scale. More information on different options in matplotlib qiver documentation.
In [6]:
m = Basemap(projection='merc',llcrnrlat=20,urcrnrlat=60,\
llcrnrlon=-100,urcrnrlon=0,lat_ts=20,resolution='c')
x, y = m(mesh.x2, mesh.y2)
In [7]:
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
u10_rot, v10_rot, x, y = m.rotate_vector(u, v, mesh.x2, mesh.y2, returnxy=True)
plt.quiver(x, y, u10_rot, v10_rot, scale=10, width=0.002)
plt.title('Speed at 100m depth')
plt.tight_layout()
Same thing but for the Arctic view
In [8]:
m = Basemap(projection='npstere',lon_0=0, boundinglat=70, resolution='c')
x, y = m(mesh.x2, mesh.y2)
In [9]:
plt.figure(figsize=(10,10))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
u10_rot, v10_rot, x, y = m.rotate_vector(u, v, mesh.x2, mesh.y2, returnxy=True)
m.quiver(x, y, u10_rot, v10_rot, scale=3, width=0.002,)
plt.title('Speed at 100m depth')
plt.tight_layout()
Sometimes it make sence to plot not all point, but some subset:
In [10]:
plt.figure(figsize=(10,10))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
u10_rot, v10_rot, x, y = m.rotate_vector(u, v, mesh.x2, mesh.y2, returnxy=True)
m.quiver(x[::2], y[::2], u10_rot[::2], v10_rot[::2], scale=3, width=0.002)
plt.title('Speed at 100m depth')
plt.tight_layout()
Please have a look at streamplot documentation in Basemap and in matplotlib.
Streamplot only works with regular data and preferably that it is regular for the map projection you choose. So first we create a regular grid for the progection:
In [11]:
m = Basemap(projection='npstere',lon_0=0, boundinglat=70, resolution='c')
x, y = m(mesh.x2, mesh.y2)
In [12]:
lons, lats, x, y = m.makegrid(300, 300, returnxy=True)
Then we interpolate our u and v values to this regular grid:
In [13]:
u_int = pf.fesom2regular(u, mesh, lons, lats, radius_of_influence=80000)
v_int = pf.fesom2regular(v, mesh, lons, lats, radius_of_influence=80000)
Vectors have to be rotated in order to have proper orientation for the map projection:
In [14]:
u10_rot, v10_rot, x, y = m.rotate_vector(u_int, v_int, lons, lats, returnxy=True)
And we calculate the speed to colorise our streamlines:
In [15]:
speed = np.hypot(u10_rot, v10_rot)
In [16]:
plt.figure(figsize=(10,10))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
m.streamplot(x, y, u10_rot, v10_rot, color=speed, density=5, norm = Normalize(vmin=0.01, vmax=0.05),
cmap=cm.Spectral_r)
cb = plt.colorbar(orientation = 'horizontal', pad=0.03)
cb.set_label('m/s')
In [ ]:
In [17]:
m = Basemap(projection='merc',llcrnrlat=20,urcrnrlat=60,\
llcrnrlon=-100,urcrnrlon=0,lat_ts=20,resolution='c')
x, y = m(mesh.x2, mesh.y2)
lons, lats, x, y = m.makegrid(300, 300, returnxy=True)
u_int = pf.fesom2regular(u, mesh, lons, lats, radius_of_influence=80000)
v_int = pf.fesom2regular(v, mesh, lons, lats, radius_of_influence=80000)
u10_rot, v10_rot, x, y = m.rotate_vector(u_int, v_int, lons, lats, returnxy=True)
speed = np.hypot(u10_rot, v10_rot)
plt.figure(figsize=(10,10))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
m.streamplot(x, y, u10_rot, v10_rot, color=speed, density=5, norm = Normalize(vmin=0.01, vmax=0.1),
cmap=cm.Spectral_r)
cb = plt.colorbar(orientation = 'horizontal', pad=0.03)
cb.set_label('m/s')
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: