import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# Resolve the latest HRRR dataset
from siphon.catalog import get_latest_access_url
hrrr_catalog = "http://thredds.ucar.edu/thredds/catalog/grib/NCEP/HRRR/CONUS_2p5km/catalog.xml"
latest_hrrr_ncss = get_latest_access_url(hrrr_catalog, "NetcdfSubset")
# Set up access via NCSS
from siphon.ncss import NCSS
ncss = NCSS(latest_hrrr_ncss)
# Create a query to ask for all times in netcdf4 format for
# the Temperature_surface variable, with a bounding box centered
# on lat,lon with a height and width speficied in degrees.
query = ncss.query()
lat = 40
lon = -105
width = 25
height = 20
We can see what variables are available from ncss as well:
ncss.variables
{'Best_4_layer_Lifted_Index_pressure_difference_layer', 'Categorical_freezing_rain_surface', 'Categorical_ice_pellets_surface', 'Categorical_rain_surface', 'Categorical_snow_surface', 'Composite_reflectivity_entire_atmosphere', 'Convective_available_potential_energy_pressure_difference_layer', 'Convective_available_potential_energy_surface', 'Convective_inhibition_pressure_difference_layer', 'Convective_inhibition_surface', 'Dewpoint_temperature_height_above_ground', 'Dewpoint_temperature_isobaric', 'Echo_top_cloud_tops', 'Geopotential_height_adiabatic_condensation_lifted', 'Geopotential_height_cloud_ceiling', 'Geopotential_height_cloud_tops', 'Geopotential_height_isobaric', 'Geopotential_height_surface', 'High_cloud_cover_high_cloud', 'Hourly_Maximum_of_Downward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum', 'Hourly_Maximum_of_Simulated_Reflectivity_at_1_km_AGL_height_above_ground_1_Hour_Maximum', 'Hourly_Maximum_of_Updraft_Helicity_over_Layer_2km_to_5_km_AGL_height_above_ground_layer_1_Hour_Maximum', 'Hourly_Maximum_of_Upward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum', 'Lightning_entire_atmosphere', 'Low_cloud_cover_low_cloud', 'Medium_cloud_cover_middle_cloud', 'Per_cent_frozen_precipitation_surface', 'Planetary_boundary_layer_height_surface', 'Precipitable_water_entire_atmosphere_single_layer', 'Pressure_of_level_from_which_parcel_was_lifted_pressure_difference_layer', 'Pressure_reduced_to_MSL_msl', 'Pressure_surface', 'Reflectivity_height_above_ground', 'Snow_depth_surface', 'Storm_relative_helicity_height_above_ground_layer', 'Surface_lifted_index_isobaric_layer', 'Temperature_height_above_ground', 'Temperature_isobaric', 'Total_cloud_cover_entire_atmosphere', 'Total_column_integrated_graupel_entire_atmosphere_single_layer_1_Hour_Maximum', 'Total_precipitation_surface_1_Hour_Accumulation', 'Vertical_u-component_shear_height_above_ground_layer', 'Vertical_v-component_shear_height_above_ground_layer', 'Vertical_velocity_geometric_sigma_layer_1_Hour_Average', 'Vertically_integrated_liquid_water_VIL_entire_atmosphere', 'Visibility_surface', 'Water_equivalent_of_accumulated_snow_depth_surface_1_Hour_Accumulation', 'Wind_speed_gust_surface', 'Wind_speed_height_above_ground_1_Hour_Maximum', 'u-component_of_wind_height_above_ground', 'u-component_of_wind_isobaric', 'u-component_storm_motion_height_above_ground_layer', 'v-component_of_wind_height_above_ground', 'v-component_of_wind_isobaric', 'v-component_storm_motion_height_above_ground_layer'}
query.all_times().accept('netcdf4').variables('Temperature_height_above_ground')
query.lonlat_box(north=lat + height/2., south=lat - height/2.,
east=lon + width/2., west=lon - width/2.)
# get the data!
data = ncss.get_data(query)
data
<class 'netCDF4._netCDF4.Dataset'> root group (NETCDF4 data model, file format HDF5): Originating_or_generating_Center: US National Weather Service, National Centres for Environmental Prediction (NCEP) Originating_or_generating_Subcenter: 0 GRIB_table_version: 2,1 Type_of_generating_process: Forecast Analysis_or_forecast_generating_process_identifier_defined_by_originating_centre: High Resolution Rapid Refresh (HRRR) Conventions: CF-1.6 history: Read using CDM IOSP GribCollection v3 featureType: GRID History: Translated to CF-1.0 Conventions by Netcdf-Java CDM (CFGridWriter2) Original Dataset = /data/ldm/pub/native/grid/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_20170407_0200.grib2.ncx3#LambertConformal_1377X2145-38p22N-95p43W; Translation Date = 2017-04-07T04:04:46.736Z geospatial_lat_min: 28.4129707527 geospatial_lat_max: 51.1859678224 geospatial_lon_min: -122.261533359 geospatial_lon_max: -91.9086814311 dimensions(sizes): time1(19), height_above_ground1(1), y(973), x(949) variables(dimensions): float32 Temperature_height_above_ground(time1,height_above_ground1,y,x), float64 time1(time1), float32 height_above_ground1(height_above_ground1), float32 y(y), float32 x(x), int32 LambertConformal_Projection() groups:
var = 'Temperature_height_above_ground'
ncvar = data[var]
ncvar
<class 'netCDF4._netCDF4.Variable'> float32 Temperature_height_above_ground(time1, height_above_ground1, y, x) long_name: Temperature @ Specified height level above ground units: K abbreviation: TMP missing_value: nan grid_mapping: LambertConformal_Projection coordinates: time1 height_above_ground1 y x Grib_Variable_Id: VAR_0-0-0_L103 Grib2_Parameter: [0 0 0] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Temperature Grib2_Parameter_Name: Temperature Grib2_Level_Type: Specified height level above ground Grib2_Generating_Process_Type: Forecast unlimited dimensions: current shape = (19, 1, 973, 949) filling on, default _FillValue of 9.969209968386869e+36 used
grid = data[ncvar.grid_mapping]
grid
<class 'netCDF4._netCDF4.Variable'> int32 LambertConformal_Projection() grid_mapping_name: lambert_conformal_conic latitude_of_projection_origin: 25.0 longitude_of_central_meridian: 265.0 standard_parallel: 25.0 earth_radius: 6371229.0 _CoordinateTransformType: Projection _CoordinateAxisTypes: GeoX GeoY unlimited dimensions: current shape = () filling on, default _FillValue of -2147483647 used
lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin
lat1 = grid.standard_parallel
earth_radius = grid.earth_radius
import cartopy
import cartopy.crs as ccrs
x, y
grid from the dataset. We will use cartopy to compute the lat/lon values which corespond to the x,y
grid. Let's grab the x and y data, and use metpy to assign the arrays a unit.¶from metpy.units import units
x = data.variables['x'][:] * units(data.variables['x'].units)
y = data.variables['y'][:] * units(data.variables['x'].units)
print(x.units, y.units)
kilometer kilometer
# cartopy wants x and y in meters, so let's make sure we are good!
x = x.to('m')
y = y.to('m')
print(x.units, y.units)
meter meter
globe = ccrs.Globe(ellipse='WGS84') # default
hrrr_proj = ccrs.LambertConformal(central_longitude=lon0, central_latitude=lat0,
standard_parallels=(lat0, lat1), globe=globe)
print(x.shape)
print(y.shape)
print(ncvar.shape)
(949,) (973,) (19, 1, 973, 949)
datetime
object¶from netCDF4 import num2date
# find the correct time dimension name
for coord in ncvar.coordinates.split():
if 'time' in coord:
timevar = data.variables[coord]
break
time = num2date(timevar[:], timevar.units)
time[6]
datetime.datetime(2017, 4, 7, 8, 0)
import cartopy.feature as cfeat
states_provinces = cfeat.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
t_step = 0
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())
cmap = plt.get_cmap("coolwarm")
color_range = plt.Normalize(243, 303)
mesh = ax.pcolormesh(x, y, ncvar[t_step].squeeze(),
transform=hrrr_proj, zorder=0,
cmap=cmap, norm=color_range)
# add some common geographic features
ax.coastlines(resolution='10m', color='black', zorder=1)
ax.add_feature(states_provinces, edgecolor='black', zorder=1)
ax.add_feature(cfeat.BORDERS)
# add some lat/lon gridlines
ax.gridlines()
ax.set_title(time[t_step].isoformat())
# add a colorbar
cax = fig.colorbar(mesh)
cax.set_label(ncvar.units)
point_query = ncss.query()
point_query.all_times()
point_query.accept('netcdf4')
point_query.variables('Temperature_isobaric', 'Dewpoint_temperature_isobaric')
point_query.variables('u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
point_query.lonlat_point(lon, lat)
# get the data!
point_data = ncss.get_data(point_query)
# get netCDF variables
pressure = point_data["isobaric"]
temp = point_data["Temperature_isobaric"]
dewpt = point_data["Dewpoint_temperature_isobaric"]
u_cmp = point_data["u-component_of_wind_isobaric"]
v_cmp = point_data["v-component_of_wind_isobaric"]
# download data and assign the units based on the variables metadata
p = pressure[:].squeeze()[t_step, :] * units(pressure.units)
T = temp[:].squeeze()[t_step, :] * units(temp.units)
Td = dewpt[:].squeeze()[t_step, :] * units(dewpt.units)
u = u_cmp[:].squeeze()[t_step, :] * units(u_cmp.units)
v = v_cmp[:].squeeze()[t_step, :] * units(v_cmp.units)
p = p.to(units.millibar)
T = T.to(units.degC)
Td = Td.to(units.degC)
u = u.to(units.knot)
v = v.to(units.knot)
from metpy.calc import lcl, dry_lapse, parcel_profile
from metpy.plots import SkewT
from metpy.units import concatenate
# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=45)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
<matplotlib.collections.LineCollection at 0x120292eb8>