Downloading model fields with NCSS
1. What is NCSS?¶
In [1]:
# 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 [2]:
ncss.variables
Out[2]:
From here, we can build a query to ask for the data we want from the server.
In [3]:
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)
2. Making sense of netCDF¶
In [4]:
data
Out[4]:
We can use a library called XArray to make working with this a little simpler
In [5]:
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 [6]:
var = ds.metpy.parse_cf('Temperature_isobaric')
var
Out[6]:
XArray handles parsing things like dates, times, latitude, and longitude for us
In [7]:
latitude = var.metpy.y
longitude = var.metpy.x
Visualize the grid¶
In [8]:
%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.add_feature(cfeature.COASTLINE)
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 [9]:
# 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.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.add_feature(cfeature.BORDERS)
# add some lat/lon gridlines
ax.gridlines()
Out[9]:
SOLUTION
In [10]:
# %load solutions/map.py
# Cell content replaced by load magic replacement.
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.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.add_feature(cfeature.BORDERS)
# add some lat/lon gridlines
ax.gridlines()
Out[10]:
4. NCSS Point Request¶
We can also request data for a specfic lon/lat point, across vertical coordinates or times.
In [11]:
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 [12]:
from metpy.units import units
import numpy as np
# get netCDF variables
pressure = point_data.variables["isobaric"]
dname_temp = point_data.variables["Temperature_isobaric"].dimensions[2]
lev_temp = point_data.variables[dname_temp]
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()
# Due to a different number of vertical levels find where they are common
_, _, common_ind = np.intersect1d(pressure, lev_temp, return_indices=True)
T = T[common_ind]
We also need to calculate dewpoint from our relative humidity data:
In [13]:
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 [14]:
p.ito(units.millibar)
T.ito(units.degC)
Td.ito(units.degC)
u.ito(units.knot)
v.ito(units.knot)
In [15]:
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()
Out[15]:
