Plotting a map of NetCDF data with Matplotlib/Basemap

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:
vic_runoff_example

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:

Advertisements

6 thoughts on “Plotting a map of NetCDF data with Matplotlib/Basemap

  1. Pingback: Basic NetCDF processing in Matlab | Water Programming: A Collaborative Research Blog

  2. Thank you very much for the program. I am facing error in adopting the code for a different netcdf file especially on the step 3. It will be much useful to have example file (vic_LHS_climatology_cells.nc) to learn further. Kindly provide a sample file.

  3. Thanks for reading. I won’t be able to provide the file because it’s fairly large (around 7 GB). This example may not generalize to netCDF files in other formats … this particular file was created manually as described here:
    https://waterprogramming.wordpress.com/2014/04/30/reducing-size-of-netcdf-spatial-data-with-list-representation/

    where the data values are just a 2-D array that you need to index manually using the lat/lon coordinates. If your netCDF file is already georeferenced, you may not need this step, you just need to get your data into a 2-D array somehow to plot it.

    Sorry I can’t be more help but you’ll need to look into how to read in your particular file.

  4. Pingback: Making Movies of Time-Evolving Global Maps with Python – Water Programming: A Collaborative Research Blog

  5. Pingback: Making Watershed Maps in Python – Water Programming: A Collaborative Research Blog

  6. Pingback: Water Programming Blog Guide (3) – Water Programming: A Collaborative Research Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s