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()

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() |

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

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.

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.

How do you install the Basemap in windows 10

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

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

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