Downloading model fields using netCDF Subset Service (NCSS)

Unidata Python Workshop


Overview:

  • Teaching: 15 minutes
  • Exercises: 15 minutes

Questions

  1. What is the netCDF Subset Service (NCSS)?
  2. How can I use Siphon to make an NCSS request?
  3. How do I plot gridded fields using CartoPy?

Objectives

  1. Use siphon to make a request using NCSS
  2. Making sense of netCDF
  3. Plot on a map
  4. Requesting for a single point

1. What is NCSS?


In [ ]:
# Resolve the latest GFS dataset
import metpy
from siphon.catalog import TDSCatalog

# Set up access via NCSS
gfs_catalog = ('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
               'Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')
cat = TDSCatalog(gfs_catalog)
ncss = cat.datasets[0].subset()

We can see what variables are available from ncss as well:


In [ ]:
ncss.variables

From here, we can build a query to ask for the data we want from the server.


In [ ]:
from datetime import datetime, timedelta

# Create a new NCSS query
query = ncss.query()

# Request data in netCDF format
query.accept('netcdf')

# Ask for our variable
query.variables('Temperature_isobaric')

# Ask for the 500 hPa surface
query.vertical_level(50000)

# Set the time range of data we want
now = datetime.utcnow()
query.time_range(now, now + timedelta(days=1))

# Set the spatial limits
query.lonlat_box(west=-110, east=-45, north=50, south=10)

# get the data!
data = ncss.get_data(query)

Top


2. Making sense of netCDF


In [ ]:
data

We can use a library called XArray to make working with this a little simpler


In [ ]:
from xarray.backends import NetCDF4DataStore
import xarray as xr

# We need the datastore so that we can open the existing netcdf dataset we downloaded
ds = xr.open_dataset(NetCDF4DataStore(data))

In [ ]:
var = ds.metpy.parse_cf('Temperature_isobaric')
var

XArray handles parsing things like date and times for us


In [ ]:
latitude = var.metpy.y
longitude = var.metpy.x

Top


Visualize the grid


In [ ]:
%matplotlib inline
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

# GFS uses lon/lat grid
data_projection = ccrs.PlateCarree()

# Make it easy to change what time step we look at
t_step = 0

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())
mesh = ax.pcolormesh(longitude, latitude, var[t_step].squeeze(),
                     transform=data_projection, zorder=0)

# add some common geographic features
ax.coastlines(resolution='10m', color='black')
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.add_feature(cfeature.BORDERS)

# add some lat/lon gridlines
ax.gridlines()

# add a colorbar
cax = fig.colorbar(mesh)
cax.set_label(var.attrs['units'])
EXERCISE:
  • Extend the plot above by plotting contours of 500 hPa geopotential heights
  • Add a title to the plot with the correct time

In [ ]:
# Set up an NCSS query from thredds using siphon
query = ncss.query()

#
# SET UP QUERY
#

# Download data using NCSS
#ncss.get_data(query)

data_projection = ccrs.PlateCarree()

# Plot using CartoPy and Matplotlib
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

#
# YOUR PLOT HERE
#

# add some common geographic features
ax.coastlines(resolution='10m', color='black')
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.add_feature(cfeature.BORDERS)

# add some lat/lon gridlines
ax.gridlines()

import numpy as np

\# Set up an NCSS query from thredds using siphon
query = ncss.query()
query.accept('netcdf4')
query.variables('Temperature_isobaric', 'Geopotential_height_isobaric')
query.vertical_level(50000)
now = datetime.utcnow()
query.time_range(now, now + timedelta(days=1))
query.lonlat_box(west=-110, east=-45, north=50, south=10)

\# Download data using NCSS
data = ncss.get_data(query)
ds = xr.open_dataset(NetCDF4DataStore(data))

temp_var = ds.metpy.parse_cf('Temperature_isobaric')
height_var = ds.metpy.parse_cf('Geopotential_height_isobaric')
longitude = temp_var.metpy.x
latitude = temp_var.metpy.y
time_index = 0

\# Plot using CartoPy and Matplotlib
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

contours = np.arange(5000, 6000, 80)
ax.pcolormesh(longitude, latitude, temp_var[time_index].squeeze(),
              transform=data_projection, zorder=0)
ax.contour(longitude, latitude, height_var[time_index].squeeze(), contours, colors='k',
           transform=data_projection, linewidths=2, zorder=1)
ax.set_title(temp_var.metpy.time[time_index].values)

\# add some common geographic features
ax.coastlines(resolution='10m', color='black')
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.add_feature(cfeature.BORDERS)

\# add some lat/lon gridlines
ax.gridlines()

Top


4. NCSS Point Request

We can also request data for a specfic lon/lat point, across vertical coordinates or times.


In [ ]:
cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
                 'Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')
ncss = cat.datasets[0].subset()

point_query = ncss.query()
point_query.time(datetime.utcnow())
point_query.accept('netcdf4')
point_query.variables('Temperature_isobaric', 'Relative_humidity_isobaric')
point_query.variables('u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
point_query.lonlat_point(-101.877, 33.583)

# get the data! Unfortunately, xarray does not quite like what comes out of thredds
point_data = ncss.get_data(point_query)

Skew-T diagrams typical use specific units. First, let's assign units to the variables we requested:


In [ ]:
from metpy.units import units

# get netCDF variables
pressure = point_data.variables["isobaric"]
temp = point_data.variables["Temperature_isobaric"]
u_cmp = point_data.variables["u-component_of_wind_isobaric"]
v_cmp = point_data.variables["v-component_of_wind_isobaric"]
relh = point_data.variables["Relative_humidity_isobaric"]

# download data and assign the units based on the variables metadata
# Need to put units on the left to assure things work properly with masked arrays
p = units(pressure.units) * pressure[:].squeeze()
T = units(temp.units) * temp[:].squeeze()
u = units(u_cmp.units) * u_cmp[:].squeeze()
v = units(v_cmp.units) * v_cmp[:].squeeze()
relh = units('percent') * relh[:].squeeze()

We also need to calculate dewpoint from our relative humidity data:


In [ ]:
import metpy.calc as mpcalc

Td = mpcalc.dewpoint_rh(T, relh)

Now, let's change those units to what we typically see used in Skew-T diagrams. We use ito to do this in-place rather than manually reassigning to the same variable.


In [ ]:
p.ito(units.millibar)
T.ito(units.degC)
Td.ito(units.degC)
u.ito(units.knot)
v.ito(units.knot)

In [ ]:
from metpy.plots import SkewT

# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=45)

# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'tab:red')
skew.plot(p, Td, 'blue')
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)

# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()