In two previous posts I showed how to create a NetCDF file from a mess of ASCII data, and how to reduce its size by storing only the lat-lon grid cells for which data exists. This post will show a brief example of reading and plotting data from a NetCDF file using Matplotlib with the Basemap toolkit.
1. Read in NetCDF variables
This is an easy one—just open the file and assign the values you want to plot.
root_grp = Dataset('climatology-nc/vic_LHS_climatology_cells.nc')
vic_runoff = root_grp.variables['vic_runoff']
lat = root_grp.variables['latitude'][:]
lon = root_grp.variables['longitude'][:]
2. Create the Basemap object
Next we create the Basemap object, m. I do this by calling a custom mapformat() function, which returns the object. Here is the function in its entirety:
from __future__ import division import numpy as np from mpl_toolkits.basemap import Basemap def mapformat(): m = Basemap(projection='robin', lon_0=0,resolution='c') # resolution c, l, i, h, f in that order m.drawmapboundary(fill_color='white', zorder=-1) m.fillcontinents(color='0.8', lake_color='white', zorder=0) m.drawcoastlines(color='0.6', linewidth=0.5) m.drawcountries(color='0.6', linewidth=0.5) m.drawparallels(np.arange(-90.,91.,30.), labels=[1,0,0,1], dashes=[1,1], linewidth=0.25, color='0.5') m.drawmeridians(np.arange(0., 360., 60.), labels=[1,0,0,1], dashes=[1,1], linewidth=0.25, color='0.5') return m
I did this because I wanted to modularize the projection and color scheme that I found to look “nice” for plotting global data. This gives light gray continents with white oceans, so that the emphasis is on the plotted data rather than the map background; your preferences may vary.
3. Populate a gridded array for plotting
Next, we’ll take our NetCDF variable of interest and put the relevant values into a single 2D array to be plotted. Recall in the previous post we arranged this file into a list format, so the NetCDF data is not gridded by default—we need to do a little bit of work to stick it into a 2D grid.
array = np.empty((180,360))
array[:] = np.NAN;
x = np.arange(0.5, 360.5, 1.0)
y = np.arange(-90, 91, 1.0)
for i in xrange(0, lat.size):
ilat = int(lat[i] + 90.0)
ilon = int(lon[i] - 0.5)
if ilon < 179 or ilon > 181: # HACK: to avoid date-line wraparound problem
array[ilat,ilon] = vic_runoff[i,ensemble,month]
array, x = shiftgrid(180, array, x)
array_mask = np.ma.masked_invalid(array)
x,y = np.meshgrid(x,y)
x,y = m(x,y)
Here, x and y are the coordinates we’re going to use for plotting. We loop through our list of grid cells and assign the values of vic_runoff to the array variable (which is just a 2D Numpy array). Note the creation of array_mask, which removes the NaN values from array.
4. Create and customize plot
The goal here is to call a function of the form plot(x,y,array). Specifically I’ll use the m.pcolormesh function, which is part of the Basemap library. Creating the plot looks like this:
rbb = np.loadtxt('cmaps/runoff-brownblue.txt')/255;
rbb = mpl.colors.ListedColormap(rbb)
m.pcolormesh(x,y,array_mask,vmin=0.0,vmax=100.0, cmap=rbb, rasterized=False, edgecolor='0.6', linewidth=0)
cbar = m.colorbar()
cbar.solids.set_edgecolor("face")
cbar.set_ticks([0,100])
plt.title("VIC Runoff (mm/month)")
plt.show()
And the result:
Shameless plug: I created this beige-blue colormap using Colormap. I think it does a nice job of conveying low-to-high runoff values—intuitively, the transition between a “desert” and a “lake”. Also it’s monotonically decreasing in grayscale, so it will print correctly in black and white. None of this is required for the task at hand, but it’s something to think about.
So there you have it, NetCDF makes it easy to plot data because you don’t need to spend time redefining all of the columns of data that you would normally read in from ASCII files. Thanks for reading, here is the full gist if you’re interested:
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import matplotlib.cm as cm | |
| from mpl_toolkits.basemap import Basemap,shiftgrid | |
| import matplotlib as mpl | |
| from mapformat import mapformat | |
| from netCDF4 import Dataset | |
| root_grp = Dataset('climatology-nc/vic_LHS_climatology_cells.nc') | |
| vic_runoff = root_grp.variables['vic_runoff'] | |
| lat = root_grp.variables['latitude'][:] | |
| lon = root_grp.variables['longitude'][:] | |
| ensemble = 5000 | |
| month = 7 | |
| m = mapformat() | |
| array = np.empty((180,360)) | |
| array[:] = np.NAN; | |
| x = np.arange(0.5, 360.5, 1.0) | |
| y = np.arange(-90, 91, 1.0) | |
| for i in xrange(0, lat.size): | |
| ilat = int(lat[i] + 90.0) | |
| ilon = int(lon[i] - 0.5) | |
| if ilon < 179 or ilon > 181: # HACK: to avoid date-line wraparound problem | |
| array[ilat,ilon] = vic_runoff[i,ensemble,month] | |
| array, x = shiftgrid(180, array, x) | |
| array_mask = np.ma.masked_invalid(array) | |
| x,y = np.meshgrid(x,y) | |
| x,y = m(x,y) | |
| rbb = np.loadtxt('cmaps/runoff-brownblue.txt')/255; | |
| rbb = mpl.colors.ListedColormap(rbb) | |
| m.pcolormesh(x,y,array_mask,vmin=0.0,vmax=100.0, cmap=rbb, rasterized=False, edgecolor='0.6', linewidth=0) | |
| cbar = m.colorbar() | |
| cbar.solids.set_edgecolor("face") | |
| cbar.set_ticks([0,100]) | |
| plt.title("VIC Runoff (mm/month)") | |
| plt.show() | |
Water Programming: A Collaborative Research Blog 