Skip to content

Update for metpy1 #113

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dcamron merged 22 commits into from
Mar 29, 2021
Merged

Update for metpy1 #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
c7142c4
update hemispheric plot for metpy 1.0
kgoebber Nov 28, 2020
46aaec4
update 500 abs vort for metpy 1.0
kgoebber Nov 28, 2020
2e58a22
update basic 500 plot for metpy 1.0
kgoebber Nov 28, 2020
d595f03
update 500 vort adv for metpy 1.0
kgoebber Nov 28, 2020
a49fbd7
update 700 basic plot for metpy 1.0
kgoebber Nov 28, 2020
bae5d6f
update 850 fronto for metpy 1.0
kgoebber Nov 28, 2020
9bee40f
update 850 qvectors for metpy 1.0
kgoebber Nov 28, 2020
777829a
update 850 basic plot for metpy 1.0
kgoebber Nov 28, 2020
2996d1e
update 850 temp adv for metpy 1.0
kgoebber Nov 28, 2020
bc69622
update ageo wind for metpy 1.0
kgoebber Nov 28, 2020
0267475
update declarative 300 for metpy 1.0
kgoebber Nov 28, 2020
7f19a6f
update diff temp adv for metpy 1.0
kgoebber Nov 28, 2020
35bac4a
update HILO symbol for metpy 1.0
kgoebber Nov 28, 2020
f5d6adc
update hovmoller for metpy1.0
kgoebber Nov 28, 2020
5525983
update isentropic for metpy 1.0
kgoebber Nov 28, 2020
9258a9d
update MSLP example for metpy 1.0
kgoebber Nov 28, 2020
8add306
update PV example for metpy 1.0
kgoebber Nov 28, 2020
6fff437
update smooth contour example metpy 1.0
kgoebber Nov 28, 2020
1b2f712
update wind shear example for metpy 1.0
kgoebber Nov 28, 2020
e8b0f3e
update declarative 500 for metpy 1.0
kgoebber Nov 28, 2020
3ce1fdc
update DT example for metpy 1.0
kgoebber Nov 28, 2020
ad81875
update notebook for metpy 1.0
kgoebber Nov 28, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 68 additions & 36 deletions pages/gallery/250hPa_Hemispheric_Plot.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,13 @@
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"import cartopy.util as cutil\n",
"import matplotlib.gridspec as gridspec\n",
"import matplotlib.pyplot as plt\n",
"import metpy.calc as mpcalc\n",
"from metpy.units import units\n",
"from netCDF4 import num2date\n",
"import numpy as np\n",
"import scipy.ndimage as ndimage\n",
"from siphon.catalog import TDSCatalog\n",
"from siphon.ncss import NCSS\n",
"from xarray.backends import NetCDF4DataStore\n",
"import xarray as xr\n",
"\n",
"# Latest GFS Dataset\n",
"cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/'\n",
Expand All @@ -58,6 +56,9 @@
"# Actually getting the data\n",
"data_hght = ncss.get_data(gfsdata_hght)\n",
"\n",
"# Make into an xarray Dataset object\n",
"ds_hght = xr.open_dataset(NetCDF4DataStore(data_hght)).metpy.parse_cf()\n",
"\n",
"# Query for Latest GFS Run\n",
"gfsdata_wind = ncss.query().time(now).accept('netcdf4')\n",
"gfsdata_wind.variables('u-component_of_wind_isobaric',\n",
Expand All @@ -71,7 +72,10 @@
"gfsdata_wind.vertical_level(25000)\n",
"\n",
"# Actually getting the data\n",
"data_wind = ncss.get_data(gfsdata_wind)"
"data_wind = ncss.get_data(gfsdata_wind)\n",
"\n",
"# Make into an xarray Dataset object\n",
"ds_wind = xr.open_dataset(NetCDF4DataStore(data_wind)).metpy.parse_cf()"
]
},
{
Expand All @@ -92,28 +96,28 @@
"metadata": {},
"outputs": [],
"source": [
"dtime = data_hght.variables['Geopotential_height_isobaric'].dimensions[0]\n",
"dlat = data_hght.variables['Geopotential_height_isobaric'].dimensions[2]\n",
"dlon = data_hght.variables['Geopotential_height_isobaric'].dimensions[3]\n",
"lat = data_hght.variables[dlat][:]\n",
"lon = data_hght.variables[dlon][:]\n",
"dtime = ds_hght.Geopotential_height_isobaric.dims[0]\n",
"lat = ds_hght.lat.values\n",
"lon = ds_hght.lon.values\n",
"\n",
"# Converting times using the num2date function available through netCDF4\n",
"times = data_hght.variables[dtime]\n",
"vtimes = num2date(times[:], times.units)\n",
"vtimes = ds_hght[dtime].values.astype('datetime64[ms]').astype('O')\n",
"\n",
"# Smooth the 250-hPa heights using a gaussian filter from scipy.ndimage\n",
"hgt_250, lon = cutil.add_cyclic_point(data_hght.variables['Geopotential_height_isobaric'][:],\n",
" coord=lon)\n",
"Z_250 = ndimage.gaussian_filter(hgt_250[0, 0, :, :], sigma=3, order=0)\n",
"\n",
"u250 = (units(data_wind.variables['u-component_of_wind_isobaric'].units) *\n",
" data_wind.variables['u-component_of_wind_isobaric'][0, 0, :, :])\n",
"v250 = (units(data_wind.variables['v-component_of_wind_isobaric'].units) *\n",
" data_wind.variables['v-component_of_wind_isobaric'][0, 0, :, :])\n",
"u250 = u250.units * cutil.add_cyclic_point(u250)\n",
"v250 = v250.units * cutil.add_cyclic_point(v250)\n",
"wspd250 = mpcalc.wind_speed(u250, v250).to('knots')"
"hgt = ds_hght.Geopotential_height_isobaric.squeeze()\n",
"hgt_250, lon = cutil.add_cyclic_point(hgt, coord=lon)\n",
"\n",
"Z_250 = mpcalc.smooth_n_point(hgt_250, 9, 2)\n",
"\n",
"# Calculate wind speed from u and v components, add cyclic point,\n",
"# and smooth slightly\n",
"u250 = ds_wind['u-component_of_wind_isobaric'].squeeze()\n",
"v250 = ds_wind['v-component_of_wind_isobaric'].squeeze()\n",
"\n",
"wspd250 = mpcalc.wind_speed(u250, v250).metpy.convert_units('knots')\n",
"wspd250 = cutil.add_cyclic_point(wspd250)\n",
"\n",
"smooth_wspd250 = mpcalc.smooth_n_point(wspd250, 9, 2)"
]
},
{
Expand All @@ -137,43 +141,71 @@
"datacrs = ccrs.PlateCarree()\n",
"plotcrs = ccrs.NorthPolarStereo(central_longitude=-100.0)\n",
"\n",
"# Make a grid of lat/lon values to use for plotting with Basemap.\n",
"# Make a grid of lat/lon values to use for plotting.\n",
"lons, lats = np.meshgrid(lon, lat)\n",
"\n",
"fig = plt.figure(1, figsize=(12., 13.))\n",
"gs = gridspec.GridSpec(2, 1, height_ratios=[1, .02],\n",
" bottom=.07, top=.99, hspace=0.01, wspace=0.01)\n",
"fig = plt.figure(1, figsize=(12., 15.))\n",
"ax = plt.subplot(111, projection=plotcrs)\n",
"\n",
"ax = plt.subplot(gs[0], projection=plotcrs)\n",
"# Set some titles for the plots\n",
"ax.set_title('250-hPa Geopotential Heights (m)', loc='left')\n",
"ax.set_title('VALID: {}'.format(vtimes[0]), loc='right')\n",
"ax.set_title(f'VALID: {vtimes[0]}', loc='right')\n",
"\n",
"# Set the extent of the image for the NH and add\n",
"# ax.set_extent([west long, east long, south lat, north lat])\n",
"ax.set_extent([-180, 180, 10, 90], ccrs.PlateCarree())\n",
"ax.coastlines('50m', edgecolor='black', linewidth=0.5)\n",
"ax.add_feature(cfeature.STATES, linewidth=0.5)\n",
"ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='black',\n",
" linewidth=0.5)\n",
"ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=0.5)\n",
"\n",
"# Add geopotential height contours every 120 m\n",
"clev250 = np.arange(9000, 12000, 120)\n",
"cs = ax.contour(lons, lats, Z_250, clev250, colors='k',\n",
" linewidths=1.0, linestyles='solid', transform=datacrs)\n",
"plt.clabel(cs, fontsize=8, inline=1, inline_spacing=10, fmt='%i',\n",
" rightside_up=True, use_clabeltext=True)\n",
"\n",
"clevsped250 = np.arange(50, 230, 20)\n",
"# Add colorfilled wind speed in knots every 20 kts\n",
"clevsped250 = np.arange(50, 200, 20)\n",
"cmap = plt.cm.get_cmap('BuPu')\n",
"cf = ax.contourf(lons, lats, wspd250, clevsped250, cmap=cmap, transform=datacrs)\n",
"cax = plt.subplot(gs[1])\n",
"cbar = plt.colorbar(cf, cax=cax, orientation='horizontal', extend='max', extendrect=True)"
"cf = ax.contourf(lons, lats, smooth_wspd250, clevsped250, cmap=cmap,\n",
" extend='max', transform=datacrs)\n",
"cbar = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50,\n",
" extendrect=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"jupytext": {
"cell_metadata_filter": "-all",
"main_language": "python",
"notebook_metadata_filter": "-all"
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
68 changes: 43 additions & 25 deletions pages/gallery/500hPa_Absolute_Vorticity_winds.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,9 @@
"outputs": [],
"source": [
"def earth_relative_wind_components(ugrd, vgrd):\n",
" \"\"\"Calculate the north-relative components of the wind from the grid-relative\n",
" components using Cartopy transform_vectors.\n",
" \"\"\"Calculate the north-relative components of the\n",
" wind from the grid-relative components using Cartopy\n",
" transform_vectors.\n",
"\n",
" Parameters\n",
" ----------\n",
Expand All @@ -84,11 +85,12 @@
" Returns\n",
" -------\n",
" unr, vnr : tuple of array-like Quantity\n",
" The north-relative wind components in the X (East-West) and Y (North-South)\n",
" directions, respectively.\n",
" The north-relative wind components in the X (East-West)\n",
" and Y (North-South) directions, respectively.\n",
" \"\"\"\n",
" if 'crs' not in list(ugrd.coords):\n",
" raise ValueError('No CRS in coordinate, be sure to use the MetPy accessor parse_cf()')\n",
" if 'metpy_crs' not in ugrd.coords:\n",
" raise ValueError('No CRS in coordinate, be sure to'\n",
" 'use the MetPy accessor parse_cf()')\n",
"\n",
" data_crs = ugrd.metpy.cartopy_crs\n",
"\n",
Expand All @@ -97,10 +99,16 @@
"\n",
" xx, yy = np.meshgrid(x, y)\n",
"\n",
" ut, vt = ccrs.PlateCarree().transform_vectors(data_crs, xx, yy, ugrd.values, vgrd.values)\n",
" ut, vt = ccrs.PlateCarree().transform_vectors(data_crs, xx, yy,\n",
" ugrd.values, vgrd.values)\n",
"\n",
" uer = ut * units(ugrd.units)\n",
" ver = vt * units(vgrd.units)\n",
" # Make a copy of u and v component DataArrays\n",
" uer = ugrd.copy()\n",
" ver = vgrd.copy()\n",
"\n",
" # Update values with transformed winds\n",
" uer.values = ut\n",
" ver.values = vt\n",
"\n",
" return uer, ver\n"
]
Expand Down Expand Up @@ -167,8 +175,10 @@
" vertical=level).squeeze(), 9, 50)\n",
"\n",
"# Select and grab 500-hPa wind components\n",
"uwnd_500 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=level).squeeze()\n",
"vwnd_500 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=level).squeeze()\n",
"uwnd_500 = ds['u-component_of_wind_isobaric'].metpy.sel(\n",
" vertical=level).squeeze().metpy.assign_latitude_longitude()\n",
"vwnd_500 = ds['v-component_of_wind_isobaric'].metpy.sel(\n",
" vertical=level).squeeze().metpy.assign_latitude_longitude()\n",
"\n",
"# Compute north-relative wind components for plotting purposes\n",
"uwnd_er, vwnd_er = earth_relative_wind_components(uwnd_500, vwnd_500)\n",
Expand All @@ -177,7 +187,8 @@
"uwnd_er = mpcalc.smooth_n_point(uwnd_er, 9, 50)\n",
"vwnd_er = mpcalc.smooth_n_point(vwnd_er, 9, 50)\n",
"\n",
"# Create a clean datetime object for plotting based on time of Geopotential heights\n",
"# Create a clean datetime object for plotting based on time\n",
"# of Geopotential heights\n",
"vtime = ds.time.data[0].astype('datetime64[ms]').astype('O')"
]
},
Expand All @@ -202,12 +213,8 @@
"metadata": {},
"outputs": [],
"source": [
"# Calculate grid spacing that is sign aware to use in absolute vorticity calculation\n",
"dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)\n",
"\n",
"# Calculate absolute vorticity from MetPy function\n",
"avor_500 = mpcalc.absolute_vorticity(uwnd_er, vwnd_er, dx, dy, lats * units.degrees,\n",
" dim_order='yx')"
"avor_500 = mpcalc.absolute_vorticity(uwnd_er, vwnd_er)"
]
},
{
Expand Down Expand Up @@ -235,7 +242,8 @@
"mapcrs = ccrs.LambertConformal(central_longitude=-100, central_latitude=35,\n",
" standard_parallels=(30, 60))\n",
"\n",
"# Set up the projection of the data; if lat/lon then PlateCarree is what you want\n",
"# Set up the projection of the data;\n",
"# if lat/lon then PlateCarree is what you want\n",
"datacrs = ccrs.PlateCarree()\n",
"\n",
"# Start the figure and create plot axes with proper projection\n",
Expand All @@ -255,27 +263,37 @@
"colors = np.vstack((colors2, (1, 1, 1, 1), colors1))\n",
"\n",
"# Plot absolute vorticity values (multiplying by 10^5 to scale appropriately)\n",
"cf = ax.contourf(lons, lats, avor_500*1e5, clevs_500_avor, colors=colors, extend='max',\n",
" transform=datacrs)\n",
"cb = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50, extendrect=True)\n",
"cf = ax.contourf(lons, lats, avor_500*1e5, clevs_500_avor, colors=colors,\n",
" extend='max', transform=datacrs)\n",
"cb = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50,\n",
" extendrect=True)\n",
"cb.set_label('Abs. Vorticity ($s^{-1}$)')\n",
"\n",
"# Plot 500-hPa Geopotential Heights in meters\n",
"clevs_500_hght = np.arange(0, 8000, 60)\n",
"cs = ax.contour(lons, lats, hght_500, clevs_500_hght, colors='black', transform=datacrs)\n",
"cs = ax.contour(lons, lats, hght_500, clevs_500_hght, colors='black',\n",
" transform=datacrs)\n",
"plt.clabel(cs, fmt='%d')\n",
"\n",
"# Set up a 2D slice to reduce the number of wind barbs plotted (every 20th)\n",
"wind_slice = (slice(None, None, 20), slice(None, None, 20))\n",
"ax.barbs(lons[wind_slice], lats[wind_slice],\n",
" uwnd_er.to('kt')[wind_slice].m, vwnd_er[wind_slice].to('kt').m,\n",
" uwnd_er.metpy.convert_units('kt')[wind_slice].values,\n",
" vwnd_er[wind_slice].metpy.convert_units('kt').values,\n",
" pivot='middle', color='black', transform=datacrs)\n",
"\n",
"# Plot two titles, one on right and left side\n",
"plt.title('500-hPa NAM Geopotential Heights (m)'\n",
" ' and Wind Barbs (kt)', loc='left')\n",
"plt.title('Valid Time: {}'.format(vtime), loc='right')"
"plt.title(f'Valid Time: {vtime}', loc='right')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -299,7 +317,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
"version": "3.8.6"
}
},
"nbformat": 4,
Expand Down
Loading