In many cases it is fine to first interpolate FESOM data to the regular grid and then plot them. Usually it will allow you to make plots faster. However sometimes you would like to work with original FESOM triangular mesh.
If you don't need interactivity on the plots, comment:
#%matplotlib notebook
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
from matplotlib import cm
from netCDF4 import Dataset
Load mesh
In [2]:
meshpath ='../../swift.dkrz.de/COREII'
mesh = pf.load_mesh(meshpath, usepickle=True)
Open the netCDF file that contains your data:
In [3]:
fl = Dataset('../../swift.dkrz.de/COREII_data/fesom.1951.oce.mean.nc')
There are several variables that we can plot:
In [4]:
fl.variables.keys()
Out[4]:
For every timestep there is a vector of values, each coresponding to the node of the model mesh.
In [5]:
fl.variables['temp'].shape
Out[5]:
In order to get values at some depth one should use get_data function from pyfesom. The input variables for this function are:
data : array
complete 3d data for one timestep
mesh : fesom_mesh object
mesh representation
depth : int
desired depth
As an output you will get:
level_data : array
2d array (actually vector) with data from the desired level.
elem_no_nan : array
array with triangles (defined as triplets of node indexes) with
not NaN elements.
The actual model level will be printed out.
Here is the result of requesting data from the level that is closest to 500 meters:
In [6]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][0,:], mesh, 500)
The information from level_data and elem_no_nan is enought to make the plot:
Define the map:
In [7]:
m = Basemap(projection='robin',lon_0=0, resolution='c')
x, y = m(mesh.x2, mesh.y2)
Plot the data from the last timestep of our dataset at 100m depth:
In [8]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][-1,:],mesh,100)
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
levels = np.arange(-3., 30., 1)
plt.tricontourf(x, y, elem_no_nan[::], level_data, levels = levels, cmap=cm.Spectral_r, extend='both')
cbar = plt.colorbar(orientation='horizontal', pad=0.03);
cbar.set_label("Temperature, $^{\circ}$C")
plt.title('Temperature at 100m depth')
plt.tight_layout()
If your color range covers the whole range of your data everything is fine. But if you try to narrow the range, the places with values outside your range become masked out:
In [9]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][-1,:],mesh,100)
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
levels = np.arange(-3., 20., 1)
plt.tricontourf(x, y, elem_no_nan[::], level_data, levels = levels, cmap=cm.Spectral_r, extend='both')
cbar = plt.colorbar(orientation='horizontal', pad=0.03);
cbar.set_label("Temperature, $^{\circ}$C")
plt.title('Temperature at 100m depth')
plt.tight_layout()
The way around it is to chnage your data :)
In [11]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][-1,:],mesh,100)
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
levels = np.arange(-3., 20., 1)
eps=(levels.max()-levels.min())/50.
level_data[level_data<=levels.min()]=levels.min()+eps
level_data[level_data>=levels.max()]=levels.max()-eps
plt.tricontourf(x, y, elem_no_nan[::], level_data, levels = levels, cmap=cm.Spectral_r, extend='both')
cbar = plt.colorbar(orientation='horizontal', pad=0.03);
cbar.set_label("Temperature, $^{\circ}$C")
plt.title('Temperature at 100m depth')
plt.tight_layout()
The things above were easy - since it was just doing contour plots. But sometimes it is nessesary to have a look at the data on original grid. Since properties in FESOM 1.4 are defined on the verticies of the triangles, the values for making the colors of the triangles will be calculated as a mean of the values in the verticies.
In [12]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][-1,:],mesh,100)
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
m.drawcoastlines()
plt.tripcolor(x, y, elem_no_nan, \
level_data, \
edgecolors='k',\
lw = 0.01,
cmap=cm.Spectral_r,
vmin = -3,
vmax = 30)
cbar = plt.colorbar(orientation='horizontal', pad=0.03);
cbar.set_label("Temperature, $^{\circ}$C")
plt.title('Temperature at 100m depth')
plt.tight_layout()
As you can see this takes considerably longer time. If you would like to have a look at the values exactly at the points were they are defined, you can try to use technique described in the fesom_with_datashader.ipynb tootorial.
In [13]:
m = Basemap(projection='mill',llcrnrlat=45,urcrnrlat=60,\
llcrnrlon=-10,urcrnrlon=15,resolution='c')
x2, y2 = m(mesh.x2, mesh.y2)
In [14]:
level_data, elem_no_nan = pf.get_data(fl.variables['temp'][-1,:],mesh,0)
plt.figure(figsize=(10,7))
m.drawmapboundary(fill_color='0.9')
#m.drawcoastlines()
plt.tripcolor(x2, y2, elem_no_nan, \
level_data, \
edgecolors='k',\
lw = 0.1,
cmap=cm.Spectral_r,
vmin = 5,
vmax = 15)
cbar = plt.colorbar(orientation='horizontal', pad=0.03);
cbar.set_label("Temperature, $^{\circ}$C")
plt.title('Temperature at 0m depth')
plt.tight_layout()
In [ ]: