.. 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_Isentropic_Interpolation.py: =================== Isentropic Analysis =================== The MetPy function `metpy.calc.isentropic_interpolation` allows for isentropic analysis from model analysis data in isobaric coordinates. .. code-block:: default from datetime import datetime, timedelta import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt import metpy.calc from metpy.units import units from netCDF4 import num2date import numpy as np from siphon.catalog import TDSCatalog Getting the data In this example, the latest GFS forecasts data from the National Centers for Environmental Information (https://nomads.ncdc.noaa.gov) will be used, courtesy of the Univeristy Corporation for Atmospheric Research Thredds Data Server. .. code-block:: default # Latest GFS Dataset cat = TDSCatalog('http://thredds-jetstream.unidata.ucar.edu/thredds/catalog/grib/' 'NCEP/GFS/Global_0p5deg/catalog.xml') ncss = cat.latest.subset() # Find the start of the model run and define time range start_time = ncss.metadata.time_span['begin'] start = datetime.strptime(start_time, '%Y-%m-%dT%H:%M:%Sz') end = start + timedelta(hours=9) # Query for Latest GFS Run gfsdata = ncss.query().time_range(start, end).accept('netcdf4') gfsdata.variables('Temperature_isobaric', 'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric', 'Relative_humidity_isobaric').add_lonlat() # Set the lat/lon box for the data you want to pull in. # lonlat_box(north_lat,south_lat,east_lon,west_lon) gfsdata.lonlat_box(-150, -50, 15, 65) # Actually getting the data data = ncss.get_data(gfsdata) dtime = data.variables['Temperature_isobaric'].dimensions[0] dlev_hght = data.variables['Temperature_isobaric'].dimensions[1] dlev_uwnd = data.variables['u-component_of_wind_isobaric'].dimensions[1] lat = data.variables['lat'][:] lon = data.variables['lon'][:] lev_hght = data.variables[dlev_hght][:] * units.Pa lev_uwnd = data.variables[dlev_uwnd][:] * units.Pa # Due to a different number of vertical levels find where they are common _, _, common_ind = np.intersect1d(lev_uwnd, lev_hght, return_indices=True) times = data.variables[dtime] vtimes = num2date(times[:], times.units) temps = data.variables['Temperature_isobaric'] tmp = temps[:, common_ind, :, :] * units.kelvin uwnd = data.variables['u-component_of_wind_isobaric'][:] * units.meter / units.second vwnd = data.variables['v-component_of_wind_isobaric'][:] * units.meter / units.second relh = data.variables['Relative_humidity_isobaric'][:] To properly interpolate to isentropic coordinates, the function must know the desired output isentropic levels. An array with these levels will be created below. .. code-block:: default isentlevs = np.arange(310, 316, 5) * units.kelvin Conversion to Isentropic Coordinates Once model data in isobaric coordinates has been pulled and the desired isentropic levels created, the conversion to isentropic coordinates can begin. Data will be passed to the function as below. The function requires that isentropic levels, isobaric levels, and temperature be input. Any additional inputs (in this case relative humidity, u, and v wind components) will be linearly interpolated to isentropic space. .. code-block:: default isent_anal = metpy.calc.isentropic_interpolation(isentlevs, lev_uwnd, tmp, relh, uwnd, vwnd, axis=1) The output is a list, so now we will separate the variables to different names before plotting. .. code-block:: default isentprs, isentrh, isentu, isentv = isent_anal A quick look at the shape of these variables will show that the data is now in isentropic coordinates, with the number of vertical levels as specified above. .. code-block:: default print(isentprs.shape) print(isentrh.shape) print(isentu.shape) print(isentv.shape) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none (4, 2, 101, 201) (4, 2, 101, 201) (4, 2, 101, 201) (4, 2, 101, 201) **Plotting the Isentropic Analysis** Set up our projection .. code-block:: default crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0) # Set up our array of latitude and longitude values and transform to # the desired projection. clons, clats = np.meshgrid(lon, lat) # Get data to plot state and province boundaries states_provinces = cfeature.NaturalEarthFeature( category='cultural', name='admin_1_states_provinces_lakes', scale='50m', facecolor='none') level = 0 FH = 0 fig = plt.figure(1, figsize=(14., 12.)) ax = plt.subplot(111, projection=crs) # Set plot extent ax.set_extent((-121., -74., 25., 50.), crs=ccrs.PlateCarree()) ax.coastlines('50m', edgecolor='black', linewidth=0.75) ax.add_feature(states_provinces, edgecolor='black', linewidth=0.5) # Plot the 300K surface clevisent = np.arange(0, 1000, 25) cs = ax.contour(clons, clats, isentprs[FH, level, :, :], clevisent, transform=ccrs.PlateCarree(), colors='k', linewidths=1.0, linestyles='solid') plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7, fmt='%i', rightside_up=True, use_clabeltext=True) cf = ax.contourf(clons, clats, isentrh[FH, level, :, :], range(10, 106, 5), transform=ccrs.PlateCarree(), cmap=plt.cm.gist_earth_r) plt.colorbar(cf, orientation='horizontal', extend=max, aspect=65, pad=0, extendrect='True') wind_slice = [FH, level, slice(None, None, 5), slice(None, None, 5)] ax.barbs(clons[wind_slice[2:]], clats[wind_slice[2:]], isentu[wind_slice].m, isentv[wind_slice].m, length=6, transform=ccrs.PlateCarree()) # Make some titles plt.title('{:.0f} K Isentropic Level'.format(isentlevs[level].m), loc='left') plt.title('VALID: {:s} UTC'.format(str(vtimes[FH])), loc='right') plt.show() .. image:: /examples/images/sphx_glr_Isentropic_Interpolation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none /home/travis/miniconda/envs/gallery/lib/python3.7/site-packages/numpy/ma/core.py:3172: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. dout = self.data[indx] /home/travis/miniconda/envs/gallery/lib/python3.7/site-packages/numpy/ma/core.py:3204: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. mout = _mask[indx] .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 41.007 seconds) .. _sphx_glr_download_examples_Isentropic_Interpolation.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: Isentropic_Interpolation.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: Isentropic_Interpolation.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_