.. 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_MSLP_temp_winds.py: ============================= Model Surface Output ============================= Plot an surface map with mean sea level pressure (MSLP), 2m Temperature (F), and Wind Barbs (kt). Imports .. code-block:: default from datetime import datetime import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt from metpy.units import units from netCDF4 import num2date import numpy as np import scipy.ndimage as ndimage from siphon.ncss import NCSS Helper functions .. code-block:: default # Helper function for finding proper time variable def find_time_var(var, time_basename='time'): for coord_name in var.coordinates.split(): if coord_name.startswith(time_basename): return coord_name raise ValueError('No time variable found for ' + var.name) Create NCSS object to access the NetcdfSubset --------------------------------------------- Data from NCEI GFS 0.5 deg Analysis Archive .. code-block:: default base_url = 'https://www.ncei.noaa.gov/thredds/ncss/grid/gfs-g4-anl-files/' dt = datetime(2018, 1, 4, 12) ncss = NCSS('{}{dt:%Y%m}/{dt:%Y%m%d}/gfsanl_4_{dt:%Y%m%d}' '_{dt:%H}00_000.grb2'.format(base_url, dt=dt)) # Create lat/lon box for location you want to get data for query = ncss.query().time(dt) query.lonlat_box(north=65, south=15, east=310, west=220) query.accept('netcdf') # Request data for model "surface" data query.variables('Pressure_reduced_to_MSL_msl', 'Apparent_temperature_height_above_ground', 'u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground') data = ncss.get_data(query) Begin data maipulation ----------------------- Data for the surface from a model is a bit complicated. The variables come from different levels and may have different data array shapes. MSLP: Pressure_reduced_to_MSL_msl (time, lat, lon) 2m Temp: Apparent_temperature_height_above_ground (time, level, lat, lon) 10m Wind: u/v-component_of_wind_height_above_ground (time, level, lat, lon) Height above ground Temp from GFS has one level (2m) Height above ground Wind from GFS has three levels (10m, 80m, 100m) .. code-block:: default # Pull out variables you want to use mslp = data.variables['Pressure_reduced_to_MSL_msl'][:].squeeze() temp = units.K * data.variables['Apparent_temperature_height_above_ground'][:].squeeze() u_wind = units('m/s') * data.variables['u-component_of_wind_height_above_ground'][:].squeeze() v_wind = units('m/s') * data.variables['v-component_of_wind_height_above_ground'][:].squeeze() lat = data.variables['lat'][:].squeeze() lon = data.variables['lon'][:].squeeze() time_var = data.variables[find_time_var(data.variables['Pressure_reduced_to_MSL_msl'])] # Convert winds to knots u_wind.ito('kt') v_wind.ito('kt') # Convert number of hours since the reference time into an actual date time = num2date(time_var[:].squeeze(), time_var.units) lev_10m = np.where(data.variables['height_above_ground3'][:] == 10)[0][0] u_wind_10m = u_wind[lev_10m] v_wind_10m = v_wind[lev_10m] # Combine 1D latitude and longitudes into a 2D grid of locations lon_2d, lat_2d = np.meshgrid(lon, lat) # Smooth MSLP a little # Be sure to only put in a 2D lat/lon or Y/X array for smoothing smooth_mslp = ndimage.gaussian_filter(mslp, sigma=3, order=0) * units.Pa smooth_mslp.ito('hPa') Begin map creation ------------------ .. code-block:: default # Set Projection of Data datacrs = ccrs.PlateCarree() # Set Projection of Plot plotcrs = ccrs.LambertConformal(central_latitude=[30, 60], central_longitude=-100) # Create new figure fig = plt.figure(figsize=(11, 8.5)) # Add the map and set the extent ax = plt.subplot(111, projection=plotcrs) plt.title('GFS Analysis MSLP, 2m Temperature (F), Wind Barbs (kt)' ' {0:%d %B %Y %H:%MZ}'.format(time), fontsize=16) ax.set_extent([235., 290., 20., 55.]) # Add state boundaries to plot states_provinces = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lakes', scale='50m', facecolor='none') ax.add_feature(states_provinces, edgecolor='black', linewidth=1) # Add country borders to plot country_borders = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m', facecolor='none') ax.add_feature(country_borders, edgecolor='black', linewidth=1) # Plot MSLP Contours clev_mslp = np.arange(0, 1200, 4) cs = ax.contour(lon_2d, lat_2d, smooth_mslp, clev_mslp, colors='black', linewidths=1.5, linestyles='solid', transform=datacrs) plt.clabel(cs, fontsize=10, inline=1, inline_spacing=10, fmt='%i', rightside_up=True, use_clabeltext=True) # Plot 2m Temperature Contours clevtemp = np.arange(-60, 101, 10) cs2 = ax.contour(lon_2d, lat_2d, temp.to(units('degF')), clevtemp, colors='tab:red', linewidths=1.25, linestyles='dotted', transform=datacrs) plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=10, fmt='%i', rightside_up=True, use_clabeltext=True) # Plot 10m Wind Barbs ax.barbs(lon_2d, lat_2d, u_wind_10m.magnitude, v_wind_10m.magnitude, length=6, regrid_shape=20, pivot='middle', transform=datacrs) plt.show() .. image:: /examples/images/sphx_glr_MSLP_temp_winds_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 4.080 seconds) .. _sphx_glr_download_examples_MSLP_temp_winds.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: MSLP_temp_winds.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: MSLP_temp_winds.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_