# Geoscience visualization with matplotlib using netCDF and THREDDS

• Sophisticated 2D and 3D plotting library.
• We'll give you a taste of the possibilities, but covering < 1%.
• Mature and can produce publication quality graphics.
• Static images, nothing interactive yet.
• We'll focus on geoscience, but matplotlib can likely meet most of your scientific plotting requirements.
• The best is to find an example of what you want in a gallery and emulate that.

## Basic example

$$f(x) = x^{3}$$

In [1]:
# Importing libraries we will need.
import numpy as np
import matplotlib.pyplot as plt

# Calculating X and Y coordinates
# linspace returns evenly spaced numbers over a specified interval.
x = np.linspace(-10,10,25)
# Raising x to third power, element-wise
y = x**3

# Make and show the plot
plt.plot(x,y)
plt.show()


## Plotting netCDF data

### Basic Plot

• Plotting temperature as a function of depth as predicted by the RTOFS model
In [2]:
# Importing libraries we will need.
import netCDF4
import matplotlib.pyplot as plt

# Open the netCDF file
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')

# Getting the n-dimensional data
tempv = f.variables['temperature']
depth = f.variables['Depth']

print("The temperature variable...\n")
# Note the temperature variable has time, depth, y and x dimensions
print(tempv)
print("The dimensions...\n")
print(tempv.dimensions)

The temperature variable...

<type 'netCDF4.Variable'>
float32 temperature(MT, Depth, Y, X)
coordinates: Longitude Latitude Date
standard_name: sea_water_potential_temperature
units: degC
_FillValue: 1.26765e+30
valid_range: [ -5.07860279  11.14989948]
long_name:   temp [90.9H]
unlimited dimensions: MT
current shape = (1, 10, 850, 712)

The dimensions...

(u'MT', u'Depth', u'Y', u'X')

In [3]:
# Continued from previous cell..

# Our goal is temperature as a function of depth so slicing along the depth axis
# at a specific time and place
temp = tempv[0,:,123,486]

print("The masked array containing the temperature data...")
print(repr(temp))

# trick for filtering out good values

print("Plotting...")
# plot and show data
plt.plot(y,x)
plt.show()

# close netCDF
f.close()

The masked array containing the temperature data...
masked_array(data = [6.485864639282227 4.6258392333984375 4.010849952697754 3.8229074478149414
3.4448373317718506 2.8652758598327637 1.785945177078247 1.333146333694458
-- --],
mask = [False False False False False False False False  True  True],
fill_value = 1.26765e+30)

Plotting...


### Let's build upon the previous plot into something that is ready for publication.

• Title
• Axis Legends
• Markers

Peruse matplotlib gallery and see and emulate what you like.

In [4]:
# Importing libraries we will need.
import netCDF4
import matplotlib.pyplot as plt

# Get figure hook to manipulate our plot
fig = plt.figure()
desc ='Figure 1. Temperature as a function of ocean depth as\n'\
'predicted by the RTOFS model'
plt.figtext(.5,.15,desc,fontsize=12,ha='center')
fig.set_size_inches(7,7)

# Improve axes
# Get axis hook to manipulate our plot
ax.set_xlabel('Depth (m)', fontweight='bold')
ax.set_ylabel('Temperature ($^\circ$C)', fontweight='bold')
# Don't show top and right axis
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
# Define ticks
ax.tick_params(axis='both', direction='out')
ax.get_xaxis().tick_bottom() # remove unneeded ticks
ax.get_yaxis().tick_left()

# Getting the data as we did before
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = f.variables['temperature']
depth = f.variables['Depth']
temp = tempv[0,:,123,486]
x = temp[~temp.mask] #trick for getting data

# Plotting line with triangle markers, and red line color.
plt.plot(y,x, marker=r'^', color='r', markersize=10, clip_on=False,linewidth=2.0)
plt.show()
f.close()


## Exercise

• Create a new notebook cell here.
• Based on cell above, plot salinity as a function of depth.
• Try using a different marker.
• Try using a different line color.

## Matplotlib Basemap

• High level API for dealing with maps
• Basemap allows you to plot data on a 2D map.
• Support 25 different map projections
• Support for shapefiles from the GIS world
• Cartopy may be a better choice eventually.

### Let's plot some data from RTOFS model on a map

• Data are coming from the local file system
In [5]:
#importing libraries
from mpl_toolkits.basemap import Basemap
import netCDF4
import matplotlib.pyplot as plt
import numpy as np

# Open and read netCDF variables
nc = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = nc.variables['temperature']
lat = nc.variables['Latitude'][:]
lon = nc.variables['Longitude'][:]
data = tempv[0,0,:,:]

# Construct our basemap object with lat/lon bounds, resolution, projection,
m = Basemap(llcrnrlon=-179,llcrnrlat=20,urcrnrlon=-100,urcrnrlat=85,
resolution='l',projection='stere',
lat_0=60,lon_0=-120.)

# Setting the plot size and text
fig = plt.figure(figsize=(10,8))
plt.figtext(.5,.15,'Figure 1. Sea surface temperature as predicted by the RTOFS model',fontsize=12,ha='center')

# compute map proj coordinates.
x, y = m(lon,lat)

# define color map
cmap = plt.cm.hsv

#Coloring the data

# Nice high-level, human-readable abstractions for dealing with maps.
m.drawcoastlines()
m.fillcontinents(color='#989865',lake_color='#336699')
m.drawmapboundary(fill_color='#336699')
m.drawparallels(np.arange(-80.,81.,20.),labels=[0,1,1,0])
m.drawmeridians(np.arange(-180.,180.,20.),labels=[1,0,0,1])

# Color bar
cbar.set_label('$^{o}C$')

nc.close()


## Exercise

• Create a new notebook cell here.
• Based on cell above, create a plot showing salinity on map.
• Try a different projection.
• Try a different color map (examples).

## Accessing and plotting data from the TDS

### Global 1/2 degree GFS example

• Plot highs and lows on a sea-level pressure map
• Note the data are coming from the TDS
In [8]:
#importing libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import minimum_filter, maximum_filter
from netCDF4 import Dataset
from pyudl.tds import get_latest_dods_url

# Function to determine local extrema based on a window size
def extrema(mat,mode='wrap',window=10):
"""find the indices of local extrema (min and max)
in the input array."""
mn = minimum_filter(mat, size=window, mode=mode)
mx = maximum_filter(mat, size=window, mode=mode)
# (mat == mx) true if pixel is equal to the local max
# (mat == mn) true if pixel is equal to the local in
# Return the indices of the maxima, minima
return np.nonzero(mat == mn), np.nonzero(mat == mx)

# Getting data from the TDS with our get_latest_dods_url helper function
gfs_data_url = \
"http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p5deg/catalog.xml"
latest = get_latest_dods_url(gfs_data_url)
data = Dataset(latest)

# Not the best way to get time. Assumes wanting to get the first time.
# And making fragile assumptions on the appearance of the unit string.
# Is there a better way of getting the time stamp?
time = data.variables['time'].units.split()[2]

lats = data.variables['lat'][:]
lons1 = data.variables['lon'][:]
nlats = len(lats)
nlons = len(lons1)

# read prmsl, convert to hPa (mb).
prmsl = 0.01*data.variables['Pressure_reduced_to_MSL_msl'][0]

# the window parameter controls the number of highs and lows detected.
# (higher value, fewer highs and lows)
local_min, local_max = extrema(prmsl, mode='wrap', window=50)

# create Basemap instance for the Miller Cylindrical projection
m = \
Basemap(llcrnrlon=0,llcrnrlat=-80,urcrnrlon=360,urcrnrlat=80,projection='mill')

# add wrap-around point in longitude.

# contour levels
clevs = np.arange(900,1100.,5.)

# find x,y of map projection grid.
# Meshgird returns coordinate matrices from two coordinate vectors.
lons, lats = np.meshgrid(lons, lats)
x, y = m(lons, lats)

# create figure.
fig=plt.figure(figsize=(8,4.5))
cs = m.contour(x,y,prmsl,clevs,colors='k',linewidths=1.)
m.drawcoastlines(linewidth=1.25)
m.fillcontinents(color='0.8')
m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
xlows = x[local_min]; xhighs = x[local_max]
ylows = y[local_min]; yhighs = y[local_max]
lowvals = prmsl[local_min]; highvals = prmsl[local_max]

# plot lows as blue L's, with min pressure value underneath.
xyplotted = []

# don't plot if there is already a L or H within dmin meters.
yoffset = 0.022*(m.ymax-m.ymin)
dmin = yoffset
for x,y,p in zip(xlows, ylows, lowvals):
if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin:
dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted]
if not dist or min(dist) > dmin:
plt.text(x,y,'L',fontsize=14,fontweight='bold',
ha='center',va='center',color='b')
plt.text(x,y-yoffset,repr(int(p)),fontsize=9,
ha='center',va='top',color='b',
bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5)))
xyplotted.append((x,y))

# plot highs as red H's, with max pressure value underneath.
xyplotted = []
for x,y,p in zip(xhighs, yhighs, highvals):
if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin:
dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted]
if not dist or min(dist) > dmin:
plt.text(x,y,'H',fontsize=14,fontweight='bold',
ha='center',va='center',color='r')
plt.text(x,y-yoffset,repr(int(p)),fontsize=9,
ha='center',va='top',color='r',
bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5)))
xyplotted.append((x,y))

plt.title('Mean Sea-Level Pressure (with Highs and Lows) ' + time)
plt.show()

In []: