# Import essential packages and setup inline plotting.
import matplotlib.pyplot as plt
from metpy.units import units
import numpy as np
%matplotlib inline
# Start setting up thredds (Thematic Real-time Environmental Distributed Data Services) access
# using siphon ncss (netCDF subset service) method.
from siphon.catalog import get_latest_access_url
gfs_catalog = "http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.html"
latest_gfs_ncss = get_latest_access_url(gfs_catalog, "NetcdfSubset")
# Set up access via NCSS
from siphon.ncss import NCSS
ncss = NCSS(latest_gfs_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 = 42
lon = -92.5
width = 25
height = 20
/usr/local/miniconda3/lib/python3.5/site-packages/siphon/catalog.py:185: UserWarning: URL http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.html returned HTML. Changing to: http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.xml new_url))
ncss.variables
{'5-Wave_Geopotential_Height_isobaric', 'Absolute_vorticity_isobaric', 'Albedo_surface_Mixed_intervals_Average', 'Apparent_temperature_height_above_ground', 'Best_4_layer_Lifted_Index_surface', 'Categorical_Freezing_Rain_surface_Mixed_intervals_Average', 'Categorical_Ice_Pellets_surface_Mixed_intervals_Average', 'Categorical_Rain_surface_Mixed_intervals_Average', 'Categorical_Snow_surface_Mixed_intervals_Average', 'Cloud_Work_Function_entire_atmosphere_single_layer_Mixed_intervals_Average', 'Cloud_mixing_ratio_isobaric', 'Cloud_water_entire_atmosphere_single_layer', 'Convective_Precipitation_Rate_surface_Mixed_intervals_Average', 'Convective_available_potential_energy_pressure_difference_layer', 'Convective_available_potential_energy_surface', 'Convective_inhibition_pressure_difference_layer', 'Convective_inhibition_surface', 'Convective_precipitation_surface_Mixed_intervals_Accumulation', 'Dewpoint_temperature_height_above_ground', 'Downward_Long-Wave_Radp_Flux_surface_Mixed_intervals_Average', 'Downward_Short-Wave_Radiation_Flux_surface_Mixed_intervals_Average', 'Field_Capacity_surface', 'Geopotential_height_highest_tropospheric_freezing', 'Geopotential_height_isobaric', 'Geopotential_height_maximum_wind', 'Geopotential_height_potential_vorticity_surface', 'Geopotential_height_surface', 'Geopotential_height_tropopause', 'Geopotential_height_zeroDegC_isotherm', 'Ground_Heat_Flux_surface_Mixed_intervals_Average', 'Haines_index_surface', 'ICAO_Standard_Atmosphere_Reference_Height_maximum_wind', 'ICAO_Standard_Atmosphere_Reference_Height_tropopause', 'Ice_cover_surface', 'Icing_severity_isobaric', 'Land-sea_coverage_nearest_neighbor_land1sea0_surface', 'Land_cover_0__sea_1__land_surface', 'Latent_heat_net_flux_surface_Mixed_intervals_Average', 'MSLP_Eta_model_reduction_msl', 'Maximum_temperature_height_above_ground_Mixed_intervals_Maximum', 'Meridional_Flux_of_Gravity_Wave_Stress_surface_Mixed_intervals_Average', 'Minimum_temperature_height_above_ground_Mixed_intervals_Minimum', 'Momentum_flux_u-component_surface_Mixed_intervals_Average', 'Momentum_flux_v-component_surface_Mixed_intervals_Average', 'Ozone_Mixing_Ratio_isobaric', 'Per_cent_frozen_precipitation_surface', 'Planetary_Boundary_Layer_Height_surface', 'Potential_Evaporation_Rate_surface', 'Potential_temperature_sigma', 'Precipitable_water_entire_atmosphere_single_layer', 'Precipitation_rate_surface_Mixed_intervals_Average', 'Pressure_convective_cloud_bottom', 'Pressure_convective_cloud_top', 'Pressure_height_above_ground', 'Pressure_high_cloud_bottom_Mixed_intervals_Average', 'Pressure_high_cloud_top_Mixed_intervals_Average', 'Pressure_low_cloud_bottom_Mixed_intervals_Average', 'Pressure_low_cloud_top_Mixed_intervals_Average', 'Pressure_maximum_wind', 'Pressure_middle_cloud_bottom_Mixed_intervals_Average', 'Pressure_middle_cloud_top_Mixed_intervals_Average', 'Pressure_of_level_from_which_parcel_was_lifted_pressure_difference_layer', 'Pressure_potential_vorticity_surface', 'Pressure_reduced_to_MSL_msl', 'Pressure_surface', 'Pressure_tropopause', 'Relative_humidity_entire_atmosphere_single_layer', 'Relative_humidity_height_above_ground', 'Relative_humidity_highest_tropospheric_freezing', 'Relative_humidity_isobaric', 'Relative_humidity_pressure_difference_layer', 'Relative_humidity_sigma', 'Relative_humidity_sigma_layer', 'Relative_humidity_zeroDegC_isotherm', 'Sensible_heat_net_flux_surface_Mixed_intervals_Average', 'Snow_depth_surface', 'Soil_temperature_depth_below_surface_layer', 'Specific_humidity_height_above_ground', 'Specific_humidity_pressure_difference_layer', 'Storm_relative_helicity_height_above_ground_layer', 'Sunshine_Duration_surface', 'Surface_Lifted_Index_surface', 'Temperature_altitude_above_msl', 'Temperature_height_above_ground', 'Temperature_high_cloud_top_Mixed_intervals_Average', 'Temperature_isobaric', 'Temperature_low_cloud_top_Mixed_intervals_Average', 'Temperature_maximum_wind', 'Temperature_middle_cloud_top_Mixed_intervals_Average', 'Temperature_potential_vorticity_surface', 'Temperature_pressure_difference_layer', 'Temperature_sigma', 'Temperature_surface', 'Temperature_tropopause', 'Total_cloud_cover_boundary_layer_cloud_Mixed_intervals_Average', 'Total_cloud_cover_convective_cloud', 'Total_cloud_cover_entire_atmosphere_Mixed_intervals_Average', 'Total_cloud_cover_high_cloud_Mixed_intervals_Average', 'Total_cloud_cover_low_cloud_Mixed_intervals_Average', 'Total_cloud_cover_middle_cloud_Mixed_intervals_Average', 'Total_ozone_entire_atmosphere_single_layer', 'Total_precipitation_surface_Mixed_intervals_Accumulation', 'U-Component_Storm_Motion_height_above_ground_layer', 'Upward_Long-Wave_Radp_Flux_atmosphere_top_Mixed_intervals_Average', 'Upward_Long-Wave_Radp_Flux_surface_Mixed_intervals_Average', 'Upward_Short-Wave_Radiation_Flux_atmosphere_top_Mixed_intervals_Average', 'Upward_Short-Wave_Radiation_Flux_surface_Mixed_intervals_Average', 'V-Component_Storm_Motion_height_above_ground_layer', 'Ventilation_Rate_planetary_boundary', 'Vertical_Speed_Shear_potential_vorticity_surface', 'Vertical_Speed_Shear_tropopause', 'Vertical_velocity_pressure_isobaric', 'Vertical_velocity_pressure_sigma', 'Visibility_surface', 'Volumetric_Soil_Moisture_Content_depth_below_surface_layer', 'Water_equivalent_of_accumulated_snow_depth_surface', 'Water_runoff_surface_Mixed_intervals_Accumulation', 'Wilting_Point_surface', 'Wind_speed_gust_surface', 'Zonal_Flux_of_Gravity_Wave_Stress_surface_Mixed_intervals_Average', 'u-component_of_wind_altitude_above_msl', 'u-component_of_wind_height_above_ground', 'u-component_of_wind_isobaric', 'u-component_of_wind_maximum_wind', 'u-component_of_wind_planetary_boundary', 'u-component_of_wind_potential_vorticity_surface', 'u-component_of_wind_pressure_difference_layer', 'u-component_of_wind_sigma', 'u-component_of_wind_tropopause', 'v-component_of_wind_altitude_above_msl', 'v-component_of_wind_height_above_ground', 'v-component_of_wind_isobaric', 'v-component_of_wind_maximum_wind', 'v-component_of_wind_planetary_boundary', 'v-component_of_wind_potential_vorticity_surface', 'v-component_of_wind_pressure_difference_layer', 'v-component_of_wind_sigma', 'v-component_of_wind_tropopause'}
query.all_times().accept('netcdf4').variables('u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground', 'Wind_speed_gust_surface', '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!
nc = ncss.get_data(query)
nc
<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: Analysis from GFS (Global Forecast System) 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/GFS/Global_0p25deg/GFS_Global_0p25deg_20171211_1200.grib2.ncx3#LatLon_721X1440-p125S-180p0E; Translation Date = 2017-12-11T20:34:33.839Z geospatial_lat_min: 32.0 geospatial_lat_max: 52.0 geospatial_lon_min: -105.0 geospatial_lon_max: -80.0 dimensions(sizes): time2(93), height_above_ground3(3), lat(81), lon(101), height_above_ground4(3) variables(dimensions): float32 u-component_of_wind_height_above_ground(time2,height_above_ground3,lat,lon), float64 time2(time2), float32 height_above_ground3(height_above_ground3), float32 lat(lat), float32 lon(lon), float32 v-component_of_wind_height_above_ground(time2,height_above_ground3,lat,lon), float32 Temperature_height_above_ground(time2,height_above_ground4,lat,lon), float32 height_above_ground4(height_above_ground4), float32 Wind_speed_gust_surface(time2,lat,lon) groups:
temp = nc['Temperature_height_above_ground']
u_wind = nc['u-component_of_wind_height_above_ground']
v_wind = nc['v-component_of_wind_height_above_ground']
wind_gust = nc['Wind_speed_gust_surface']
height_temp = nc['height_above_ground4']
height_wind = nc['height_above_ground3']
print(temp.units,u_wind.units,v_wind.units,wind_gust.units,height_temp.units,height_wind.units)
K m/s m/s m/s m m
temp, u_wind, v_wind, wind_gust, height_temp, height_wind
(<class 'netCDF4._netCDF4.Variable'> float32 Temperature_height_above_ground(time2, height_above_ground4, lat, lon) long_name: Temperature @ Specified height level above ground units: K abbreviation: TMP missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 height_above_ground4 lat lon 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: 103 Grib2_Level_Desc: Specified height level above ground Grib2_Generating_Process_Type: Forecast unlimited dimensions: current shape = (93, 3, 81, 101) filling on, default _FillValue of 9.969209968386869e+36 used, <class 'netCDF4._netCDF4.Variable'> float32 u-component_of_wind_height_above_ground(time2, height_above_ground3, lat, lon) long_name: u-component of wind @ Specified height level above ground units: m/s abbreviation: UGRD missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 height_above_ground3 lat lon Grib_Variable_Id: VAR_0-2-2_L103 Grib2_Parameter: [0 2 2] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Momentum Grib2_Parameter_Name: u-component of wind Grib2_Level_Type: 103 Grib2_Level_Desc: Specified height level above ground Grib2_Generating_Process_Type: Forecast unlimited dimensions: current shape = (93, 3, 81, 101) filling on, default _FillValue of 9.969209968386869e+36 used, <class 'netCDF4._netCDF4.Variable'> float32 v-component_of_wind_height_above_ground(time2, height_above_ground3, lat, lon) long_name: v-component of wind @ Specified height level above ground units: m/s abbreviation: VGRD missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 height_above_ground3 lat lon Grib_Variable_Id: VAR_0-2-3_L103 Grib2_Parameter: [0 2 3] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Momentum Grib2_Parameter_Name: v-component of wind Grib2_Level_Type: 103 Grib2_Level_Desc: Specified height level above ground Grib2_Generating_Process_Type: Forecast unlimited dimensions: current shape = (93, 3, 81, 101) filling on, default _FillValue of 9.969209968386869e+36 used, <class 'netCDF4._netCDF4.Variable'> float32 Wind_speed_gust_surface(time2, lat, lon) long_name: Wind speed (gust) @ Ground or water surface units: m/s abbreviation: GUST missing_value: nan grid_mapping: LatLon_Projection coordinates: time2 lat lon Grib_Variable_Id: VAR_0-2-22_L1 Grib2_Parameter: [ 0 2 22] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Momentum Grib2_Parameter_Name: Wind speed (gust) Grib2_Level_Type: 1 Grib2_Level_Desc: Ground or water surface Grib2_Generating_Process_Type: Forecast unlimited dimensions: current shape = (93, 81, 101) filling on, default _FillValue of 9.969209968386869e+36 used, <class 'netCDF4._netCDF4.Variable'> float32 height_above_ground4(height_above_ground4) units: m long_name: Specified height level above ground positive: up Grib_level_type: 103 datum: ground _CoordinateAxisType: Height _CoordinateZisPositive: up unlimited dimensions: current shape = (3,) filling on, default _FillValue of 9.969209968386869e+36 used, <class 'netCDF4._netCDF4.Variable'> float32 height_above_ground3(height_above_ground3) units: m long_name: Specified height level above ground positive: up Grib_level_type: 103 datum: ground _CoordinateAxisType: Height _CoordinateZisPositive: up unlimited dimensions: current shape = (3,) filling on, default _FillValue of 9.969209968386869e+36 used)
height_temp[:], height_wind[:]
(array([ 2., 80., 100.], dtype=float32), array([ 10., 80., 100.], dtype=float32))
u_wind_surf = u_wind[:,0,:,:] * units('m/s').to('mph')
v_wind_surf = v_wind[:,0,:,:] * units('m/s').to('mph')
t_surface = temp[:,0,:,:] * units('K').to('degF')
surf_wind = np.sqrt((u_wind_surf[:,:,:]*u_wind_surf[:,:,:])+(v_wind_surf[:,:,:]*v_wind_surf[:,:,:]))
u_wind_surf.shape, v_wind_surf.shape, surf_wind.shape, t_surface.shape
((93, 81, 101), (93, 81, 101), (93, 81, 101), (93, 81, 101))
import metpy.calc
surf_wind_metpy = metpy.calc.get_wind_speed(u_wind_surf,v_wind_surf)
surf_wind_metpy.shape
(93, 81, 101)
wind_chill_metpy = metpy.calc.windchill(t_surface,surf_wind,mask_undefined=False)
wind_chill_calc = 35.74* (0.6215*t_surface[:]) - (35.75*(surf_wind[:])**0.16) + (0.4275*t_surface[:]*(surf_wind[:])**0.16)
--------------------------------------------------------------------------- DimensionalityError Traceback (most recent call last) <ipython-input-61-924fc757dc42> in <module>() ----> 1 wind_chill_calc = 35.74* (0.6215*t_surface[:]) - (35.75*(surf_wind[:])**0.16) + (0.4275*t_surface[:]*(surf_wind[:])**0.16) /usr/local/miniconda3/lib/python3.5/site-packages/pint/quantity.py in __sub__(self, other) 678 679 def __sub__(self, other): --> 680 return self._add_sub(other, operator.sub) 681 682 def __rsub__(self, other): /usr/local/miniconda3/lib/python3.5/site-packages/pint/quantity.py in _add_sub(self, other, op) 584 raise DimensionalityError(self._units, other._units, 585 self.dimensionality, --> 586 other.dimensionality) 587 588 # Next we define some variables to make if-clauses more readable. DimensionalityError: Cannot convert from 'degF' ([temperature]) to 'mph ** 0.16' ([length] ** 0.16 / [time] ** 0.16)
import cartopy
import cartopy.crs as ccrs
from netCDF4 import num2date
# find the correct time dimension name
for d in temp.dimensions:
if "time" in d:
timevar = d
time = num2date(nc[timevar][:],nc[timevar].units)
time[6]
datetime.datetime(2017, 12, 12, 6, 0)
# Open and read netCDF variables
lat = nc['lat']
lon = nc['lon']
data = wind_chill_metpy[92, :, :]
# Set up a lat/lon projection
proj = ccrs.PlateCarree()
# Construct figure
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Setting the plot size and text
plt.figtext(.3, .15,
'Figure 1. This is not what I am concerned about',
fontsize=12, ha='center')
# define color map
cmap=plt.get_cmap("Blues")
#cmap = plt.cm.RdBu_r
# Nice high-level, human-readable abstractions for dealing with maps.
# add some common geographic features
import cartopy.feature as cfeature
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
ax.add_feature(states_provinces, edgecolor='black', zorder=2)
ax.add_feature(cartopy.feature.BORDERS,edgecolor='black',zorder=2)
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(zorder=3)
ax.gridlines()
ax.set_title(time[92].isoformat());
# Color-filled contour plot
cs = ax.contourf(lon[:], lat[:], data[:], 10, cmap=cmap,transform=ccrs.PlateCarree(), zorder=2)
# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('Wind Chill')