First we need some grids of values to work with. We can do this by dowloading information from the latest run of the GFS available on Unidata's THREDDS data server. First we access the catalog for the half-degree GFS output, and look for the dataset called the "Best GFS Half Degree Forecast Time Series". This dataset combines multiple sets of model runs to yield a time series of output with the shortest forecast offset.
In [ ]:
from siphon.catalog import TDSCatalog
cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/'
'NCEP/GFS/Global_0p5deg/catalog.xml')
best = cat.datasets['Best GFS Half Degree Forecast Time Series']
Next, we set up access to request subsets of data from the model. This uses the NetCDF Subset Service (NCSS) to make requests from the GRIB collection and get results in netCDF format.
In [ ]:
subset_access = best.subset()
query = subset_access.query()
Let's see what variables are available. Instead of just printing subset_access.variables, we can ask Python to only display variables that end with "isobaric", which is how the TDS denotes GRIB fields that are specified on isobaric levels.
In [ ]:
sorted(v for v in subset_access.variables if v.endswith('isobaric'))
Now we put together the "query"--the way we ask for data we want. We give ask for a wide box of data over the U.S. for the time step that's closest to now. We also request temperature, height, winds, and relative humidity. By asking for netCDF4 data, the result is compressed, so the download is smaller.
In [ ]:
from datetime import datetime
query.time(datetime.utcnow())
query.variables('Temperature_isobaric', 'Geopotential_height_isobaric',
'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric',
'Relative_humidity_isobaric')
query.lonlat_box(west=-130, east=-50, south=10, north=60)
query.accept('netcdf4')
Now all that's left is to actually make the request for data:
In [ ]:
nc = subset_access.get_data(query)
Open the returned netCDF data using XArray:
In [ ]:
from xarray.backends import NetCDF4DataStore
import xarray as xr
ds = xr.open_dataset(NetCDF4DataStore(nc))
Now let's take what we've downloaded, and use it to make an isentropic map. In this case, we're interpolating from one vertical coordinate, pressure, to another: potential temperature. MetPy has a function isentropic_interpolation that can do this for us. First, let's start with a few useful imports.
In [ ]:
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
Let's parse out the metadata for the isobaric temperature and get the projection information.
In [ ]:
dat = ds.metpy.parse_cf('Temperature_isobaric')
In [ ]:
data_proj = dat.metpy.cartopy_crs
Let's pull out the grids out into some shorter variable names. We also index with 0 to get the first, and only, time:
In [ ]:
lat = dat.metpy.y
lon = dat.metpy.x
press = ds['isobaric']
temperature = ds['Temperature_isobaric'][0]
# Need to adjust units on humidity because '%' causes problems
ds['Relative_humidity_isobaric'].attrs['units'] = 'percent'
rh = ds['Relative_humidity_isobaric'][0]
height = ds['Geopotential_height_isobaric'][0]
u = ds['u-component_of_wind_isobaric'][0]
v = ds['v-component_of_wind_isobaric'][0]
Next, we perform the isentropic interpolation. At a minimum, this must be given one or more isentropic levels, the 3-D temperature field, and the pressure levels of the original field; it then returns the 3D array of pressure values (2D slices for each isentropic level). You can also pass addition fields which will be interpolated to these levels as well. Below, we interpolate the winds (and pressure) to the 295K isentropic level:
In [ ]:
isen_level = np.array([320]) * units.kelvin
isen_press, isen_u, isen_v = mpcalc.isentropic_interpolation(isen_level, press,
temperature, u, v)
Let's plot the results and see what it looks like:
In [ ]:
# Need to squeeze() out the size-1 dimension for the isentropic level
isen_press = isen_press.squeeze()
isen_u = isen_u.squeeze()
isen_v = isen_v.squeeze()
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# Create a plot and basic map projection
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-100))
ax.coastlines()
# Contour the pressure values for the isentropic level. We keep the handle
# for the contour so that we can have matplotlib label the contours
levels = np.arange(300, 1000, 25)
cntr = ax.contour(lon, lat, isen_press, transform=data_proj,
colors='black', levels=levels)
ax.clabel(cntr, fmt='%d')
# Set up slices to subset the wind barbs--the slices below are the same as `::5`
# We put these here so that it's easy to change and keep all of the ones below matched
# up.
lon_slice = slice(None, None, 5)
lat_slice = slice(None, None, 3)
ax.barbs(lon[lon_slice], lat[lat_slice],
isen_u[lat_slice, lon_slice].to('knots').magnitude,
isen_v[lat_slice, lon_slice].to('knots').magnitude,
transform=ccrs.PlateCarree(), zorder=2)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=2)
ax.add_feature(cfeature.STATES, linestyle=':')
ax.set_extent((-120, -70, 25, 55), crs=ccrs.PlateCarree())
Let's add some moisture information to this plot. Feel free to choose a different isentropic level.
mpcalc)isentropic_interpolation with mixing ratio--you should copy the one from above and add mixing ratio to the call so that it interpolates everything.contour (in green) or contourf your moisture information on the map alongside pressureYou'll want to refer to the MetPy API documentation to see what calculation functions would help you.
In [ ]:
# Needed to make numpy broadcasting work between 1D pressure and other 3D arrays
# Use .metpy.unit_array to get numpy array with units rather than xarray DataArray
pressure_for_calc = press.metpy.unit_array[:, None, None]
#
# YOUR CODE: Calculate mixing ratio using something from mpcalc
#
# Take the return and convert manually to units of 'dimenionless'
#mixing.ito('dimensionless')
#
# YOUR CODE: Interpolate all the data
#
# Squeeze the returned arrays
#isen_press = isen_press.squeeze()
#isen_mixing = isen_mixing.squeeze()
#isen_u = isen_u.squeeze()
#isen_v = isen_v.squeeze()
# Create Plot -- same as before
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-100))
ax.coastlines()
levels = np.arange(300, 1000, 25)
cntr = ax.contour(lon, lat, isen_press, transform=ccrs.PlateCarree(),
colors='black', levels=levels)
ax.clabel(cntr, fmt='%d')
lon_slice = slice(None, None, 8)
lat_slice = slice(None, None, 8)
ax.barbs(lon[lon_slice], lat[lat_slice],
isen_u[lat_slice, lon_slice].to('knots').magnitude,
isen_v[lat_slice, lon_slice].to('knots').magnitude,
transform=ccrs.PlateCarree(), zorder=2)
#
# YOUR CODE: Contour/Contourf the mixing ratio values
#
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=2)
ax.add_feature(cfeature.STATES, linestyle=':')
ax.set_extent((-120, -70, 25, 55), crs=ccrs.PlateCarree())
In [ ]:
# %load solutions/mixing.py
In [ ]:
isen_press.units, isen_u.units, isen_v.units
In [ ]:
isen_press = mpcalc.smooth_gaussian(isen_press.squeeze(), 9)
isen_u = mpcalc.smooth_gaussian(isen_u.squeeze(), 9)
isen_v = mpcalc.smooth_gaussian(isen_v.squeeze(), 9)
In [ ]:
isen_press.units
Next, we need to take our grid point locations which are in degrees, and convert them to grid spacing in meters--this is what we need to pass to functions taking derivatives.
In [ ]:
# Use .values because we don't care about using DataArray
dx, dy = mpcalc.lat_lon_grid_deltas(lon.values, lat.values)
Now we can calculate the isentropic ascent. $\omega$ is given by:
$$\omega = \left(\frac{\partial P}{\partial t}\right)_\theta + \vec{V} \cdot \nabla P + \frac{\partial P}{\partial \theta}\frac{d\theta}{dt}$$Note, the second term of the above equation is just pressure advection (negated). Therefore, we can use MetPy to calculate this as:
In [ ]:
lift = -mpcalc.advection(isen_press, [isen_u, isen_v], [dx, dy], dim_order='yx')
In [ ]:
# YOUR CODE GOES HERE
In [ ]:
# %load solutions/lift.py