.. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_PV_baroclinic_isobaric.py: Baroclinic Potential Vorticity Analysis, Divergence, and Winds ============================================================== Classic baroclinic potential vorticity plot at 250 hPa using GFS analysis file. This example uses example data from the GFS analysis for 12 UTC 31 October 2016 and uses xarray as the main read source with using MetPy to calculate the baroclinic potential vorticity, divergence and wind speed with geographic plotting using Cartopy for a CONUS view of the 250-hPa surface with divergence and wind barbs. By: Kevin Goebbert .. code-block:: default import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt import metpy.calc as mpcalc from metpy.units import units import numpy as np import xarray as xr The following code reads the example data using the xarray open_dataset function and prints the coordinate values that are associated with the various variables contained within the file. This code uses the metpy accessor to parse file to make it easy to pull data using common coordinate names (e.g., vertical) and attach units. .. code-block:: default ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/' 'casestudies/python-gallery/GFS_20101026_1200.nc').metpy.parse_cf() Data Retrieval -------------- This code retrieves the necessary data from the file and completes some smoothing of the temperature, and wind fields using the MetPy function smooth_n_point. A nicely formated valid time (vtime) variable is also created. .. code-block:: default # Set subset slice for the geographic extent of data to limit download lon_slice = slice(200, 350) lat_slice = slice(85, 10) # Grab lat/lon values (GFS will be 1D) lats = ds.lat.sel(lat=lat_slice).values lons = ds.lon.sel(lon=lon_slice).values # Grab the pressure levels and select the data to be imported # Need all pressure levels for Temperatures, U and V Wind, and Rel. Humidity # Smooth with the gaussian filter from scipy pres = ds['isobaric3'].values[:] * units('Pa') tmpk_var = ds['Temperature_isobaric'].metpy.sel(lat=lat_slice, lon=lon_slice).squeeze() tmpk = mpcalc.smooth_n_point(tmpk_var, 9, 2) thta = mpcalc.potential_temperature(pres[:, None, None], tmpk) uwnd_var = ds['u-component_of_wind_isobaric'].metpy.sel(lat=lat_slice, lon=lon_slice).squeeze() vwnd_var = ds['v-component_of_wind_isobaric'].metpy.sel(lat=lat_slice, lon=lon_slice).squeeze() uwnd = mpcalc.smooth_n_point(uwnd_var, 9, 2) vwnd = mpcalc.smooth_n_point(vwnd_var, 9, 2) # Create a clean datetime object for plotting based on time of Geopotential heights vtime = ds.time.data[0].astype('datetime64[ms]').astype('O') Use MetPy to compute the baroclinic potential vorticity on all isobaric levels and other variables .. code-block:: default # Compute dx and dy spacing for use in vorticity calculation dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats) # Comput the PV on all isobaric surfaces pv = mpcalc.potential_vorticity_baroclinic(thta, pres[:, None, None], uwnd, vwnd, dx[None, :, :], dy[None, :, :], lats[None, :, None] * units('degrees')) # Use MetPy to compute the divergence on the pressure surfaces div = mpcalc.divergence(uwnd, vwnd, dx[None, :, :], dy[None, :, :], dim_order='yx') # Find the index value for the 250-hPa surface i250 = list(pres.m).index(((250 * units('hPa')).to(pres.units)).m) Map Creation ------------ This next set of code creates the plot and draws contours on a Lambert Conformal map centered on -100 E longitude. The main view is over the CONUS with isobaric PV map with PV contoured every 1 PVU and divergence colorshaded. .. code-block:: default # Set up the projection that will be used for plotting mapcrs = ccrs.LambertConformal(central_longitude=-100, central_latitude=35, standard_parallels=(30, 60)) # Set up the projection of the data; if lat/lon then PlateCarree is what you want datacrs = ccrs.PlateCarree() # Start the figure and create plot axes with proper projection fig = plt.figure(1, figsize=(14, 12)) ax = plt.subplot(111, projection=mapcrs) ax.set_extent([-130, -72, 20, 55], ccrs.PlateCarree()) # Add geopolitical boundaries for map reference ax.add_feature(cfeature.COASTLINE.with_scale('50m')) ax.add_feature(cfeature.STATES.with_scale('50m')) # Plot the contours of PV at 250 hPa, scaling 10^6 every 1 PVU clevs_pv = np.arange(0, 25, 1) cs1 = ax.contour(lons, lats, pv[i250]*1e6, clevs_pv, colors='black', transform=datacrs) plt.clabel(cs1, fmt='%d', fontsize='large') # Plot the colorfill of divergence, scaled 10^5 every 1 s^1 clevs_div = np.arange(-15, 16, 1) cs1 = ax.contourf(lons, lats, div[i250]*1e5, clevs_div, cmap=plt.cm.PuOr, extend='both', transform=datacrs) plt.colorbar(cs1, orientation='horizontal', pad=0, aspect=50, extendrect=True) # Plot the wind barbs at 250 hPa wind_slice = slice(None, None, 6) ax.barbs(lons[wind_slice], lats[wind_slice], uwnd[i250][wind_slice, wind_slice].to('kt').m, vwnd[i250][wind_slice, wind_slice].to('kt').m, pivot='middle', color='black', length=6.5, transform=datacrs) # Plot some titles to tell people what is on the map plt.title('250-hPa GFS PV (PVU), Divergence ($10^5$ s$^{-1}$),' ' and Wind Barbs (kt)', loc='left') plt.title('Valid Time: {}'.format(vtime), loc='right') plt.show() .. image:: /examples/images/sphx_glr_PV_baroclinic_isobaric_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 4.004 seconds) .. _sphx_glr_download_examples_PV_baroclinic_isobaric.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: PV_baroclinic_isobaric.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: PV_baroclinic_isobaric.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_