This notebook shows how to visualize FESOM grid. We are going to use %matplotlib notebook option in Jupyter notebooks in order to make the plots interactive.
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
First we have to load the FESOM mesh. We only need information about surface nodes and elements (triangles), so we put get3d to False.
meshpath ='../../../../FESOM/mesh/'
mesh = pf.fesom_mesh(meshpath, get3d=False)
Load 2d part of the grid in 1 second(s)
You can have a look at the basic information about the mesh. Note that the information about 3d nods are not loaded.
mesh
FESOM mesh: path = /home/magik/TOOLING/FESOM/mesh alpha, beta, gamma = 50, 15, -90 number of 2d nodes = 126859 number of 2d elements = 244660 number of 3d nodes = 0
Here we remove cyclic elements from the variables that contain information about elements and volumes (2d volumes, really areas).
elem2=mesh.elem[mesh.no_cyclic_elem,:]
voltri = mesh.voltri[mesh.no_cyclic_elem]
Setup a map and convert positions of the elements to the map projection
map = Basemap(projection='robin',lon_0=0, resolution='l')
x, y = map(mesh.x2, mesh.y2)
Plot triangular mesh:
x - lon (x) coordinates
y - lat (y) coordinates
elem2 - 3d array, where each raw have trhee indexes of the x and y. So, each raw describes triangle
plt.figure(figsize=(9,5))
map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.triplot(x, y, elem2, lw=0.2,color='k');
The resulting plot is interactive. You can zoom to specific areas and do some other basic manupulations by pressing icons on the lower left.
IMPORTANT - when the figure is in the interactive mode, it considered to be the active one, so the output of the following ploting functions will apear on this figure as well. To avoid this when you finish to work with the figure, you have to press the "switch off" button at the upper right of the figure.
One can plot color codded volume of the triangles as well. This is is much more computationally expensive.
#this cell can be executed only once if you don't plan to change the Basemap settings
map = Basemap(projection='robin',lon_0=0, resolution='l')
x, y = map(mesh.x2, mesh.y2)
map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.tripcolor(x, y, elem2, \
facecolors=mesh.voltri[mesh.no_cyclic_elem], \
edgecolors='k',\
cmap=cm.Spectral_r)
plt.colorbar(orientation='horizontal', pad=0.03);
The Basemap can use many different projections and one of course can define the region to be plotted. Please refer to the Basemap documentation for details.
map = Basemap(projection='npstere',boundinglat=60,lon_0=0,resolution='l')
x, y = map(mesh.x2, mesh.y2)
map.drawmapboundary(fill_color='0.9')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
#note that here we use *mesh.elem* directly since the projection can adequatelly process it.
im=plt.triplot(x, y, mesh.elem, lw=0.2,color='k')
map = Basemap(projection='merc',llcrnrlat=30,urcrnrlat=70,\
llcrnrlon=-30,urcrnrlon=30,lat_ts=20,resolution='l')
x, y = map(mesh.x2, mesh.y2)
plt.figure(figsize=(8,9))
map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.triplot(x, y, elem2, lw=0.2,color='k');
map.etopo()
<matplotlib.image.AxesImage at 0x7f37a0412e50>