Siphon is a python package that makes downloading data from Unidata data technologies a breeze! In our examples, we'll focus on interacting with the netCDF Subset Service (NCSS) as well as the radar server to retrieve grid data and radar data.
But first! Bookmark these resources for when you want to use Siphon later!
First, we'll import the TDSCatalog class from Siphon and put the special 'matplotlib' line in so our map will show up later in the notebook. Let's construct an instance of TDSCatalog pointing to our dataset of interest. In this case, I've chosen the TDS' "Best" virtual dataset for the GFS global 0.25 degree collection of GRIB files. This will give us a good resolution for our map. This catalog contains a single dataset.
%matplotlib inline
from siphon.catalog import TDSCatalog
best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
'Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
best_gfs.datasets
OrderedDict([('Best GFS Quarter Degree Forecast Time Series', <siphon.catalog.Dataset at 0x10eb8f240>)])
We pull out this dataset and look at the access urls. Note there are many ways to access the data.
best_ds = list(best_gfs.datasets.values())[0]
best_ds.access_urls
{'CdmRemote': 'http://thredds.ucar.edu/thredds/cdmremote/grib/NCEP/GFS/Global_0p25deg/Best', 'ISO': 'http://thredds.ucar.edu/thredds/iso/grib/NCEP/GFS/Global_0p25deg/Best', 'NCML': 'http://thredds.ucar.edu/thredds/ncml/grib/NCEP/GFS/Global_0p25deg/Best', 'NetcdfSubset': 'http://thredds.ucar.edu/thredds/ncss/grib/NCEP/GFS/Global_0p25deg/Best', 'OPENDAP': 'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p25deg/Best', 'UDDC': 'http://thredds.ucar.edu/thredds/uddc/grib/NCEP/GFS/Global_0p25deg/Best', 'WCS': 'http://thredds.ucar.edu/thredds/wcs/grib/NCEP/GFS/Global_0p25deg/Best', 'WMS': 'http://thredds.ucar.edu/thredds/wms/grib/NCEP/GFS/Global_0p25deg/Best'}
The NetcdfSubset
entry is what we're after...we'll use this in our NCSS class. Let's import the NCSS class from Siphon and then pass in the NetcdfSubset access url.
from siphon.ncss import NCSS
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])
We can then use the ncss
object to create a new query object, which
facilitates asking for data from the server.
query = ncss.query()
We can look at the ncss.variables
object to see what variables are available from the dataset:
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_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', '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'}
We construct a query asking for data corresponding to a latitude and longitude box where 43 lat is the northern extent, 35 lat is the southern extent, 260 long is the western extent and 249 is the eastern extent. Note that longitude values are the longitude distance from the prime meridian. We request the data for the current time. This request will return all surface temperatures for points in our bounding box for a single time. Note the string representation of the query is a properly encoded query string.
from datetime import datetime
query.lonlat_box(north=43, south=35, east=260, west=249).time(datetime.utcnow())
query.accept('netcdf4')
query.variables('Temperature_surface')
var=Temperature_surface&time=2016-10-17T03%3A29%3A36.065934&north=43&east=260&south=35&west=249&accept=netcdf4
We now request data from the server using this query. The NCSS
class handles parsing this NetCDF data (using the netCDF4
module). If we print out the variable names, we see our requested variables, as well as a few others (more metadata information)
data = ncss.get_data(query)
list(data.variables.keys())
['Temperature_surface', 'reftime', 'time', 'lat', 'lon']
We'll pull out the temperature variable.
temp_3d = data.variables['Temperature_surface']
We need to look at the coordinates
attribute to determine the name of the proper time variable for this variable, since the GRIB collection can have multiple time variables.
temp_3d.coordinates
'reftime time lat lon '
# 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)
time_coord_name = find_time_var(temp_3d)
We'll pull out the useful variables for latitude, and longitude, and time (which is the time, in hours since the forecast run). Notice the variable names are labeled to show how many dimensions each variable is. This will come in to play soon when we prepare to plot. Try printing one of the variables to see some info on the data!
time_1d = data.variables[time_coord_name]
lat_1d = data.variables['lat']
lon_1d = data.variables['lon']
time_1d
<class 'netCDF4._netCDF4.Variable'> float64 time(time) units: Hour since 2016-10-14T00:00:00Z standard_name: time long_name: GRIB forecast or observation time calendar: proleptic_gregorian _CoordinateAxisType: Time unlimited dimensions: current shape = (1,) filling on, default _FillValue of 9.969209968386869e+36 used
Now we make our data suitable for plotting. We'll import numpy so we can combine lat/longs (meshgrid) and remove one-dimensional entities from our arrays (squeeze). Also we'll use netCDF4's num2date to change the time since the model run to an actual date.
import numpy as np
from scipy.constants import K2F
from netCDF4 import num2date
# Reduce the dimensions of the data
temp_2d = temp_3d[:].squeeze()
lat_1d = lat_1d[:].squeeze()
lon_1d = lon_1d[:].squeeze()
# Convert the number of hours since the reference time to an actual date
time_val = num2date(time_1d[:].squeeze(), time_1d.units)
# Convert temps to Fahrenheit from Kelvin
temp_2d = K2F(temp_2d)
# Combine latitude and longitudes
lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
Now we can plot these up using matplotlib. We import cartopy and matplotlib classes, create our figure, add a map, then add the temperature data and grid points.
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.plots import ctables
# Create a new figure
fig = plt.figure(figsize=(15, 12))
# Add the map and set the extent
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-100.03, -111.03, 35, 43])
# Retrieve the state boundaries using cFeature and add to plot
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
ax.add_feature(states_provinces, edgecolor='gray')
# Contour temperature at each lat/long
contours = ax.contourf(lon_2d, lat_2d, temp_2d, 200, transform=ccrs.PlateCarree(),
cmap='RdBu_r')
#Plot a colorbar to show temperature and reduce the size of it
fig.colorbar(contours)
# Make a title with the time value
ax.set_title(u'Temperature forecast (\u00b0F) for {0}Z'.format(time_val), fontsize=20)
# Plot markers for each lat/long to show grid points for 0.25 deg GFS
ax.plot(lon_2d.flatten(), lat_2d.flatten(), linestyle='none', marker='o',
color='black', markersize=2, alpha=0.3, transform=ccrs.PlateCarree());
Create your own map using the same projection as above but plot different data variables such as dewpoint or relative humidity.