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, it makes for an interesting case study. We will be meandering through,
To get our bearings let's see what there is inside our netCDF file. We will be using the xray library to dig inside our netCDF data. xray is similar to pandas, but for the Common Data Model. We could have just used the netcdf4-python library but xray has output that is more nicely formatted. Let's first import xray and open the dataset.
In [1]:
import xray
ds = xray.open_dataset('https://motherlode.ucar.edu/repository/opendap/41f2b38a-4e70-4135-8ff8-dbf3d1dcbfc1/entry.das',
decode_times=False)
print(ds)
As far as the dimensions and coordinates go, the most relevant and important coordinates variables are x
and y
. We can see the data variables, such as, temperature (t
), mixing ratio (mr
), and potential temperature (th
), are mostly on a 1901 x 1801 grid. There is also the mysterious nav
dimension and associated data variables which we will be examining later.
Let's set a goal of visualizing potential temperature with the Cartopy plotting package.
The first step is to get more information concerning the variables we are interested in. For example, let's look at potential temperature or th
.
In [2]:
print(ds['th'])
In [3]:
th = ds['th'].values[0][0]
print(th)
In order, to visualize the data that are contained within a two-dimensional array onto a map that represents a three-dimensional globe, we need to understand the projection of the data.
We can make an educated guess these are contained in the data variables with the nav
cooridinate variable.
Specifically,
grid_type
grid_type_code
x_dim
y_dim
Nx
Ny
La1
Lo1
LoV
Latin1
Latin2
Dx
Dy
But what are these??
In [4]:
print(ds['grid_type_code'])
A simple Google search of GRIB-1 GDS data representation type
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.
In [5]:
print(ds['grid_type_code'].values[0])
grid_type_code
of 5
?Let's look at Table 6 . A grid_type_code
of 5
corresponds to a projection of Polar Stereographic.
In [6]:
grid_type = ds['grid_type'].values
print('The grid type is ', grid_type[0])
In [7]:
nx, ny = ds['Nx'].values[0], ds['Ny'].values[0]
print(nx, ny)
In [8]:
la1, lo1 = ds['La1'].values[0], ds['Lo1'].values[0]
print(la1, lo1)
Latin1
and Latin2
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.
In [9]:
latin1, latin2 = ds['Latin1'].values[0], ds['Latin2'].values[0]
print(latin1, latin2)
In [10]:
lov = ds['LoV'].values[0]
print(lov)
In [11]:
print(ds['Dx'])
print(ds['Dy'])
In [12]:
dx,dy = ds['Dx'].values[0],ds['Dy'].values[0]
print(dx,dy)
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
).
In [13]:
%matplotlib inline
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib as mpl
In [14]:
proj = ccrs.LambertConformal(central_longitude=lov,standard_parallels=(latin1,latin2))
Remember, we have:
Nx
, Ny
)Dx
, Dy
)first latitude
, first longitude
).We have one of the corners in latitude and longitude, but we need to convert it LC coordinates and derive the other corner.
In [15]:
pc = ccrs.PlateCarree()
In [16]:
left,bottom = proj.transform_point(lo1,la1,pc)
print(left,bottom)
In [17]:
right,top = left + nx*dx,bottom + ny*dy
print(right,top)
In [18]:
#Define the figure
fig = plt.figure(figsize=(12, 12))
# Define the extents and add the data
ax = plt.axes(projection=proj)
extents = (left, right, bottom, top)
ax.contourf(th, origin='lower', extent=extents, transform=proj)
# Add bells and whistles
ax.coastlines(resolution='50m', color='black', linewidth=2)
ax.add_feature(ccrs.cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',facecolor='none'))
ax.add_feature(ccrs.cartopy.feature.BORDERS, linewidth='1', edgecolor='black')
ax.gridlines()
plt.show()
In [22]:
th.shape
Out[22]:
In [ ]: