Using GeoCAT to plot MPAS on native mesh

[RCarpenter]

Has anyone used GeoCAT to plot MPAS on its native mesh? I know that capability exists in NCL, but I don't know if it has been ported to GeoCAT.

 

[mcurry]

Hi RCarpenter, great question! We have not yet here, but this is a good question to prompt us to do so. I'll investigate over the next few days and see if I can generate some examples for you!

 

[mcurry]

Hi RCarpenter, sorry about the very slow reply, but I finally have some examples using some of the Geocat-Viz package, Cartopy and NetCDF4-Python

The first, is a general scrip that takes a variable name, vertical level, and an MPAS output file (that contains the requested variable and LatCell and LonCell), and produces a black and white contour plot, similar to the one produced in this GeoCat example plot. Instead of contour, it uses tricontour.

Here is that plotting script, and an example output (they are both also attached below):

import argparse
import sys

import numpy as np
from netCDF4 import Dataset
import cartopy.crs as ccrs

import matplotlib.pyplot as plt
import matplotlib.tri as tri

from scipy.spatial import Delaunay
import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter

import geocat.viz as viz
import geocat.viz.util as gvutil


def normalize_coords(lat, lon):
    lat = np.rad2deg(lat)
    lon = np.rad2deg(lon)
    lon[lon > 180.0] -= 360
    return lat, lon

parser = argparse.ArgumentParser()
parser.add_argument('var', type=str, help="Desired variable to plot")
parser.add_argument('vlevel', type=int, help="Desired vertical level")
parser.add_argument('file',
                    type=str,
                    help='''File you want to plot from''')
args = parser.parse_args()

mesh = Dataset(args.file)
varName = args.var
vlevel = args.vlevel

var = mesh.variables[varName]
field = var[:]
field = field.squeeze()

latCell = mesh.variables['latCell'][:]
lonCell = mesh.variables['lonCell'][:]
zgrid = mesh.variables['zgrid'][:]
latCell, lonCell = normalize_coords(latCell, lonCell)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
# set_global() is currently needed for global plots. It seems Cartopy (or Matplotlib) currently has
# issues with setting the extent. See https://github.com/SciTools/cartopy/issues/1345
ax.set_global()
ax.add_feature(cfeature.LAND, color='silver')

gvutil.set_axes_limits_and_ticks(ax,
                                 xlim=(-180, 180),
                                 ylim=(-90, 90),
                                 xticks=np.linspace(-180,180,13),
                                 yticks=np.linspace(-90,90,7))

gvutil.add_major_minor_ticks(ax)
gvutil.add_lat_lon_ticklabels(ax)

ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))
ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))

gvutil.set_titles_and_labels(ax,
                             maintitle="{0} at {1:.0f} m".format(var.name, zgrid[0,vlevel]),
                             lefttitle=var.long_name,
                             lefttitlefontsize=13,
                             righttitle=var.units,
                             righttitlefontsize=13,
                             xlabel="",
                             ylabel="")

contour = ax.tricontour(lonCell, latCell, field[:,vlevel], 
                        colors='black',
                        levels=7,
                        linewidths=0.5,
                        linestyles='-',
                        transform=ccrs.PlateCarree())

ax.clabel(contour, contour.levels, fontsize=12, fmt="%.0f")

plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.tight_layout()
plt.savefig('mpas_{0}_l{1}.png'.format(var.name, vlevel))

To me, this plot looks a little funny, and that might be because I'm plotting uZonal winds at a height level, rather than a pressure level, but, its a start, and you can adapt the script for your needs.

You can easily adapt the above script to use tricontourf, rather than tricontour. But here is another example of plotting the terrain field in case you are curious:

import sys

import numpy as np
from netCDF4 import Dataset
import cartopy.crs as ccrs

import matplotlib.pyplot as plt
import matplotlib.tri as tri

from scipy.spatial import Delaunay
import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter

import geocat.viz as viz
import geocat.viz.util as gvutil


def normalize_coords(lat, lon):
    lat = np.rad2deg(lat)
    lon = np.rad2deg(lon)
    lon[lon > 180.0] -= 360
    return lat, lon

print('Argv[1]: ',sys.argv[1])
mesh = Dataset(sys.argv[1])

ter = mesh.variables['ter']
latCell = mesh.variables['latCell'][:]
lonCell = mesh.variables['lonCell'][:]
latCell, lonCell = normalize_coords(latCell, lonCell)

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()

gvutil.set_axes_limits_and_ticks(ax,
                                 xlim=(-180, 180),
                                 ylim=(-90, 90),
                                 xticks=np.linspace(-180,180,13),
                                 yticks=np.linspace(-90,90,7))

gvutil.add_major_minor_ticks(ax, labelsize=10)
gvutil.add_lat_lon_ticklabels(ax)

gvutil.set_titles_and_labels(ax,
                             maintitle="Terrain Height",
                             lefttitle=ter.long_name,
                             lefttitlefontsize=13,
                             righttitle=ter.units,
                             righttitlefontsize=13,
                             xlabel="",
                             ylabel="")

contour = ax.tricontourf(lonCell, latCell, ter[:], 
                         levels=100,
                         cmap='terrain', 
                         transform=ccrs.PlateCarree())

fig.colorbar(contour, orientation="horizontal")

plt.savefig('mpas_terrain.png')

If you notice there are some missing data on the upper corners of the map. I belive GeoCat-Viz has a utility xr_add_cyclic_longitudes that you can use to fix this. Instead of using NetCDF4, you'll need to use XArray to open your NetCDF file.


Lastly, if you would like to use colored filled plots (and not contour filled plots), you can use the mesh_to_triangles.py function found in the MPAS-Tools repository to allow you to use tripcolor. You can see some examples of using mesh_to_triangles.py in PR #319. The examples in that PR use the Ocean core of the MPAS-Model, but it is very easy to substitute SST for an MPAS-Atmosphere variable.