Plotting data on a map with python

There are two main python packages that may be used to plot data on a map using Python: Basemap, Iris/Cartopy, and PyNGL (the python interface to NCL). 

This is a brief tutorial on how to use the Basemap package.

  • First thing, import the module!

from mpl_toolkits.basemap import Basemap

  • make sure you have also imported numpy
  • Now read in some of your WRF netcdf data (use whichever package you want: Nio, netCDF4, Make sure to read in your latitude and longitude values, as you need these in order to have georeferenced data. After you've done this, the next step is to set up your map. 
Which projection are you going to use?
 Since we set up our original WRF domain as a Lambert Conic Conformal (lcc) projection, it makes sense to use that projection. To set up the projection, we'll have to specify many of the same things we did in our namelist.wps (see API documentation for more details:
  • lat_1: first standard parallel
  • lat_2: second standard parallel
  • lat_0, lon_0: latitude and longitude of central point
  • projection name: here will be 'lcc'
  • area_thresh: may be specified if you want to ignore plotting of coastline features less than a certain area
  • resolution: specifies how coarse or detailed you want the plotted boundaries. May be specified as
    • 'c': coarse
    • 'l': low
    • 'i': intermediate
    • 'h': high
    • 'f': full
  • llcrnrlat: latitude of lower-left corner of your map
  • urcrnrlat: latitude of upper-right corner of your map
  • llcrnrlon: longitude of lower-left corner of your map
  • urcrnrlon: latitude of upper-right corner of your map

#these next four parameters are specified the same as in our namelist.wps
lon_0 = -95.  
lat_0 = 40.
lat_1 = 45.
lat_2 = 28.

#set lat/lons for lower-left and upper-right corners of our map
#this assumes that you have read in your 2-D (no time dimension) lat/lon values 
#(here, names of arrays are 'lat' and 'lon')

imax,jmax = lat.shape
minlat = lat[0,0]
maxlat = lat[imax-1,jmax-1]
minlon = lon[0,0]
maxlon = lon[imax-1,jmax-1]

#set up the map
m = Basemap(projection='lcc',lon_0=lon_0,lat_0=lat_0,\
              llcrnrlat=minlat,urcrnrlat=maxlat, llcrnrlon=minlon,urcrnrlon=maxlon,\

        You can of course change the lat/lon values of your corners if you want to zoom in on an area, say Oklahoma. I've found that:
minlat = 33. 
maxlat = 38. 
minlon = -103. 
maxlon = -93.

        works well for plotting  Oklahoma.
  • In order to be able to plot data on our map, we have to convert our georeferenced coordinates to map projection coordinates. This is made simple in basemap with a simple function call:
x, y = m(lon,lat)

        Our longitude values will be our x coordinate, and latitude the y coordinate.