I was asked by a colleague to visualize data contained within this netCDF file (OPeNDAP link) with Python. What follows is an exploration of how I achieved that objective. Because this exercise touches upon many technologies related to Unidata, I want to share this notebook with the community. We will be meandering through,
To get our bearings let's first open our netCDF file with ToolsUI. ToolsUI is a Unidata application to analyze a variety of geoscience data file formats including netCDF. It is entirely possible to glean all the information using only Python data exploration, but ToolsUI is an effective, mature, advanced tool for metadata exploration. Find a screenshot of ToolsUI below.
This step reveals a number of data variables such as temperature (t
), mixing ratio (mr
), and potential temperature (th
). Let's set a goal of visualizing potential temperature in Python with Basemap.
# For ToolsUI screenshot
from IPython.display import Image
Image(filename='tu1.png')
Because the record
and z
dimensions are 1 as shown by ToolsUI, we are basically looking at a two-dimensional 1801 x 1901 grid for the data variables of interest. The next challenge to visualize these data will require us to understand the map projection for the grid that contain these data. The information to derive the map projection is available in ToolsUI and is outlined in red below.
Image(filename='tu2.png')
This part of the exercise is rather disorienting if you do not have a familiarity with the map projection metadata this file is making use of. Fortunately, we are given a hint by the GRIB-1 GDS data representation type
description of the grid_type_code
variable.
A simple Google search of that description takes us to A GUIDE TO THE CODE FORM FM 92-IX Ext. GRIB Edition 1 from 1994 document. Therein one can find an explanation of the variables needed to understand the map projection. Let's review these variables. We will need to import netcdf4-python first and open the file with that API.
import netCDF4
f = netCDF4.Dataset('https://motherlode.ucar.edu/repository/opendap/41f2b38a-4e70-4135-8ff8-dbf3d1dcbfc1/entry.das')
grid_type_code = f.variables['grid_type_code'][0]
print 'grid type code is ' + str(grid_type_code)
grid type code is 5
A grid_type_code
of 5
corresponds to a projection of Polar Stereographic
. (See Table 6 .)
grid_type = f.variables['grid_type']
print 'The grid type is ' + ''.join(grid_type[0,:])
The grid type is secant lambert conformal
Note that this newest piece of information relating to a Lambert Conformal projection disagrees with the earlier projection information about a Polar Stereographic projection. There is a bug in the metadata description of the projection, but let's move on anyway. Next let's grab the number of grid points along the x and y axes as given by Nx
and Ny
.
nx, ny = f.variables['Nx'][0], f.variables['Ny'][0]
print nx, ny
1901 1801
Next let's get first latitude
, and first longitude
which is probably the latitude and longitude for one of the corners of the grid. ToolsUI shows all the latitude and longitude measurements that follow are in degrees.
la1, lo1 = f.variables['La1'][0], f.variables['Lo1'][0]
print la1, lo1
10.1608 78.9861
Next up are the rather mysteriously named Latin1
and Latin2
variables. When I first saw these identifiers, I thought they referred to a Unicode block, but in fact they relate to the secants of the projection cone. I do not know why they are called "Latin" and this name is confusing. At any rate, we can feel comfortable that we are dealing with Lambert Conformal rather than Polar Stereographic.
latin1, latin2 = f.variables['Latin1'][0], f.variables['Latin2'][0]
print latin1, latin2
30.0 60.0
If we are defining a Lambert Conformal projection, we will require the central meridian that the GRIB documentation refers to as LoV
.
lov = f.variables['LoV'][0]
print lov
102.0
Finally, let's look at the grid increments. The units are in meters
as shown by ToolsUI.
dx, dy = f.variables['Dx'][0], f.variables['Dy'][0]
print dx, dy
3000.0 3000.0
We now have all the information we need to understand the Lambert projection:
Latin1
, Latin2
)LoV
)Moreover, we have additional information that shows how the data grid relates to the projection:
Nx
, Ny
)Dx
, Dy
)first latitude
, first longitude
).We will assume a spherical Earth with a 6,371,200 meter radius.
In theory, we are ready to try to plot with Basemap
, but as we will see in a moment, we are not quite out of the woods yet. The Basemap
constructor for the Lambert projection will look something like this:
Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,
urcrnrlat=urcrnrlat, projection='lcc',lat_1=latin1,lat_2=latin2,lon_0=lov,
resolution ='l')
Here are the parameters in greater detail:
projection
the abbreviation of the projection we will be usingllcrnrlon
longitude of lower left hand cornerllcrnrlat
latitude of lower left hand cornerurcrnrlon
longitude of upper right hand cornerurcrnrlat
latitude of upper right hand cornerlat_1
the first standard parallel for lambert conformal projectionlat_2
the second standard parallel for lambert conformal projectionlon_0
the central meridianThe main difficulty is we do not have one of the corners of the grid though we do have llcrnrlat
and llcrnrlon
.
PROJ.4 is a library for handling cartographic projections. pyproj is a Python API that wraps the PROJ.4 library.
We will be using pyproj
to calculate the remaining corner of the grid. A bit of research (i.e., a lot of Google searches) reveals that we should be constructing the Lambert Conformal projection with pyproj
in this manner:
from pyproj import Proj
prj = Proj(proj='lcc', R=6371200, lat_1=latin1, lat_2=latin2, lat_0=la1 , lon_0=lov)
The parameter list is similar to the parameter list for the Basemap
constructor, but with a couple of differences.
R
is the Earth radius in meterslat_0
is the latitude of the originAt this point, we are ready to take our Proj
object instance for a test drive. Let's find out where the first latitude
and first longitude
are in x y space in meters relative to the origin of the projection. These x y values will turn out to be important as we will see in a minute.
prj(lo1, la1) # note, longitude is first!
(-2850539.9857709, 412501.2590335927)
Remember what are objectives are; to get the corners of the grids to make Basemap
happy. We have one corner; la1
and la0
, and actually we have the other too, but in x
y
distance relative to la1
and la0
. Let's redefine our Proj
object relative to la1
and la0
by adding so-called "false eastings" and "false northings"
fe, fn = prj(lo1, la1) # note, longitude is first!
prj = Proj(proj='lcc', R=6371200, lat_1=latin1, lat_2=latin2, lat_0=la1 , lon_0=lov, x_0=-fe, y_0=-fn)
Let's sanity check our our new Proj
object by making sure la1
and la0
are now at (or very close to) the origin of the projection.
prj(lo1, la1) # note, longitude is first!
(-8.996576070785522e-07, -4.0727900341153145e-07)
Good, it appears the results of our calculations are now going to be relative to la1
and la0
. The last part is now easy. We know the grid increment and number of grid points so we know the location in x y of the opposite corner. Note the inverse=True
parameter to back calculate the latitude and longitude from x y.
urcrnrlon, urcrnrlat = prj(nx*dx,ny*dy, inverse=True)
# And for parameter convenience
llcrnrlon, llcrnrlat = lo1,la1
print llcrnrlat, llcrnrlon, urcrnrlat, urcrnrlon
10.1608 78.9861 54.0011247454 149.41758739
We are now finally ready to plot our data, but first let's simply plot the projection we have worked so hard to obtain using Basemap
. The secants and central meridian are in red.
from mpl_toolkits.basemap import Basemap
fig = plt.figure(figsize=(10,10))
bm = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,
urcrnrlat=urcrnrlat,
projection='lcc',lat_1=latin1,lat_2=latin2,lon_0=lov,
resolution ='l')
bm.drawcoastlines()
bm.drawcountries()
bm.drawmapboundary(fill_color='#99ffff')
bm.fillcontinents(color='#cc9966',lake_color='#99ffff')
parallels = np.arange(0.,81,10.)
bm.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
bm.drawmeridians(meridians,labels=[True,False,False,True])
secants = [60,30]
bm.drawparallels(secants,labels=[False,False,False,False],
dashes=[1000, 1],linewidth=3,linestyle='solid',color = 'red')
central_meridian = [lov]
bm.drawmeridians(central_meridian,labels=[False,False,False,False],
dashes=[1000, 1],linewidth=3,linestyle='solid',color = 'red')
plt.show()
Now let's try to plot a field such as temperature.
# Get the temperature from the netCDF file
t = f.variables['th']
fig = plt.figure(figsize=(10,10))
bm = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,
urcrnrlat=urcrnrlat,
projection='lcc',lat_1=latin1,lat_2=latin2,lon_0=lov,
resolution ='l')
lons, lats = bm.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.
x, y = bm(lons, lats) # compute map proj coordinates.
cs = bm.contourf(x,y,t[0][0], zorder=1)
bm.drawcoastlines()
bm.drawcountries()
bm.drawmapboundary(fill_color='#99ffff')
bm.fillcontinents(color='#cc9966',lake_color='#99ffff', zorder=0)
parallels = np.arange(0.,81,10.)
bm.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
bm.drawmeridians(meridians,labels=[True,False,False,True])
cbar = bm.colorbar(cs,location='bottom',pad="5%")
cbar.set_label('kelvin')
plt.show()
That's it. The hard part was decrypting some obscure GRIB-1 metadata, and finding the corners of the grid with Proj
. After that, everything fell into place with Basemap
.
Thanks to Rich Signell, Ryan May, and Yuan Ho