Show FESOM mesh

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.

In [2]:
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.

In [3]:
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.

In [4]:
mesh
Out[4]:
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).

In [5]:
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

In [19]:
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 
In [20]:
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.

In [19]:
#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)
In [20]:
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.

North Pole stereographic projection

In [16]:
map = Basemap(projection='npstere',boundinglat=60,lon_0=0,resolution='l')
x, y = map(mesh.x2, mesh.y2)
In [18]:
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')

Mercator projection with etopo topography

In [29]:
map = Basemap(projection='merc',llcrnrlat=30,urcrnrlat=70,\
            llcrnrlon=-30,urcrnrlon=30,lat_ts=20,resolution='l')

x, y = map(mesh.x2, mesh.y2)
In [30]:
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()