import cartopy
import matplotlib as mpl
import matplotlib.pyplot as plt
from siphon.catalog import get_latest_access_url
from siphon.ncss import NCSSGet the catalog
# get NCSS access url for the most recent GEFS run
gefs_catalog_url = "http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GEFS/Global_1p0deg_Ensemble/members/catalog.xml"
ncss_access_url = get_latest_access_url(gefs_catalog_url, "NetcdfSubset")
print(ncss_access_url)http://thredds.ucar.edu/thredds/ncss/grib/NCEP/GEFS/Global_1p0deg_Ensemble/members/GEFS_Global_1p0deg_Ensemble_20160120_1800.grib2
Generate the query
# create the dataset, make a subset query, and get the data
gefs_ds = NCSS(ncss_access_url)
# See what variables are available
gefs_ds.variables{'Categorical_Freezing_Rain_surface_6_Hour_Average_ens',
'Categorical_Ice_Pellets_surface_6_Hour_Average_ens',
'Categorical_Rain_surface_6_Hour_Average_ens',
'Categorical_Snow_surface_6_Hour_Average_ens',
'Convective_available_potential_energy_pressure_difference_layer_ens',
'Convective_inhibition_pressure_difference_layer_ens',
'Downward_Long-Wave_Radp_Flux_surface_6_Hour_Average_ens',
'Downward_Short-Wave_Radiation_Flux_surface_6_Hour_Average_ens',
'Geopotential_height_isobaric_ens',
'Geopotential_height_surface_ens',
'Latent_heat_net_flux_surface_6_Hour_Average_ens',
'Maximum_temperature_height_above_ground_6_Hour_Maximum_ens',
'Minimum_temperature_height_above_ground_6_Hour_Minimum_ens',
'Precipitable_water_entire_atmosphere_single_layer_ens',
'Pressure_reduced_to_MSL_msl_ens',
'Pressure_surface_ens',
'Relative_humidity_height_above_ground_ens',
'Relative_humidity_isobaric_ens',
'Sensible_heat_net_flux_surface_6_Hour_Average_ens',
'Snow_depth_surface_ens',
'Soil_temperature_depth_below_surface_layer_ens',
'Temperature_height_above_ground_ens',
'Temperature_isobaric_ens',
'Total_cloud_cover_entire_atmosphere_6_Hour_Average_ens',
'Total_precipitation_surface_6_Hour_Accumulation_ens',
'Upward_Long-Wave_Radp_Flux_atmosphere_top_6_Hour_Average_ens',
'Upward_Long-Wave_Radp_Flux_surface_6_Hour_Average_ens',
'Upward_Short-Wave_Radiation_Flux_surface_6_Hour_Average_ens',
'Vertical_velocity_pressure_isobaric_ens',
'Volumetric_Soil_Moisture_Content_depth_below_surface_layer_ens',
'Water_equivalent_of_accumulated_snow_depth_surface_ens',
'u-component_of_wind_height_above_ground_ens',
'u-component_of_wind_isobaric_ens',
'v-component_of_wind_height_above_ground_ens',
'v-component_of_wind_isobaric_ens'}
west = -130 # degrees east
east = -55 # degrees east
north = 55 # degrees north
south = 20 # degrees north
level = 50000.0 # Pa
query = gefs_ds.query()
query.lonlat_box(west, east, south, north)
query.all_times()
query.accept('netcdf4')
var_name = 'Geopotential_height_isobaric_ens'
query.variables(var_name)
query.vertical_level(level)
data = gefs_ds.get_data(query)Make a pretty plot
lon = data.variables['lon'][:]
lat = data.variables['lat'][:]
data_var = data.variables[var_name]
print(data_var)
geopotential_height = data_var[:].squeeze()float32 Geopotential_height_isobaric_ens(time2, ens, isobaric2, lat, lon) long_name: Geopotential height @ Isobaric surface units: gpm abbreviation: HGT missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 ens isobaric2 lat lon Grib_Variable_Id: VAR_0-3-5_L100 Grib2_Parameter: [0 3 5] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Mass Grib2_Parameter_Name: Geopotential height Grib2_Level_Type: Isobaric surface Grib2_Generating_Process_Type: Ensemble forecast unlimited dimensions: current shape = (65, 21, 1, 36, 76) filling on, default _FillValue of 9.969209968386869e+36 used
from netCDF4 import num2date
time_var_name = data_var.coordinates.split()[0]
time_var = data.variables[time_var_name]
datetimes = num2date(time_var[:], time_var.units)# plot projection - Albers Equal Area
plot_bounds = [-120, -70, 20, 50]
aee = cartopy.crs.AlbersEqualArea(central_longitude=-95, central_latitude=35)
# create image
fig = plt.figure(figsize=(18, 12))
for subfig, time_idx in enumerate([0, 4, 8, 12], 1):
ax = fig.add_subplot(2, 2, subfig, projection=aee)
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(resolution='50m')
ax.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines',
scale='50m', facecolor='none'))
ax.add_feature(cartopy.feature.BORDERS, linewidth='2', edgecolor='black')
ax.gridlines()
ax.set_extent(plot_bounds, cartopy.crs.PlateCarree())
for member_data in geopotential_height[time_idx]:
ax.set_title(datetimes[time_idx])
ax.contour(lon, lat, member_data,
levels=[5520],
transform=cartopy.crs.PlateCarree())
var_name = 'Snow_depth_surface_ens'
query.variables(var_name)
query.vertical_level(None)
data = gefs_ds.get_data(query)lon = data.variables['lon'][:]
lat = data.variables['lat'][:]
data_var = data.variables[var_name]
print(data_var)
snow_depth = data_var[:].squeeze()float32 Snow_depth_surface_ens(time2, ens, lat, lon) long_name: Snow depth @ Ground or water surface units: m abbreviation: SNOD missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 ens lat lon Grib_Variable_Id: VAR_0-1-11_L1 Grib2_Parameter: [ 0 1 11] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Moisture Grib2_Parameter_Name: Snow depth Grib2_Level_Type: Ground or water surface Grib2_Generating_Process_Type: Ensemble forecast unlimited dimensions: current shape = (65, 21, 36, 76) filling on, default _FillValue of 9.969209968386869e+36 used
time_var_name = data_var.coordinates.split()[0]
time_var = data.variables[time_var_name]
datetimes = num2date(time_var[:], time_var.units)
# create image
fig = plt.figure(figsize=(18, 12))
for subfig, time_idx in enumerate(range(8, 12), 1):
ax = fig.add_subplot(2, 2, subfig, projection=aee)
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(resolution='50m')
ax.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines',
scale='50m', facecolor='none'))
ax.add_feature(cartopy.feature.BORDERS, linewidth='2', edgecolor='black')
ax.gridlines()
ax.set_extent(plot_bounds, cartopy.crs.PlateCarree())
for member in range(len(data.dimensions['ens'])):
snow_inches = snow_depth[time_idx, member] * 39.3701
ax.set_title(datetimes[time_idx])
ax.contour(lon, lat, snow_inches,
levels=[6, 12], colors=['darkcyan', 'cyan'],
transform=cartopy.crs.PlateCarree())
ax.contour(lon, lat, geopotential_height[time_idx, member],
levels=[5520], colors='darkred',
transform=cartopy.crs.PlateCarree())
