How to use Siphon and Cartopy to visualize data served by a THREDDS Data Server (TDS) via ncWMS

In [1]:

    

import cartopy
import matplotlib as mpl
import matplotlib.pyplot as plt
from owslib.wms import WebMapService
from siphon.catalog import get_latest_access_url


    

In [2]:

    

catalog = 'http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/grib/NCEP/HRRR/CONUS_2p5km/catalog.xml'
serverurl = get_latest_access_url(catalog, 'WMS')

wms = WebMapService( serverurl, version='1.1.1')


    

In [3]:

    

#Listing all available layers...
layers = list(wms.contents)
for layer in layers:
    print('Layer name: {}'.format(wms[layer].name))


    

    




Layer name: Best_4-layer_lifted_index_pressure_difference_layer
Layer name: u-component_storm_motion_height_above_ground_layer
Layer name: wind @ Isobaric surface
Layer name: Categorical_ice_pellets_surface
Layer name: Temperature_isobaric
Layer name: u-component_of_wind_isobaric
Layer name: Convective_available_potential_energy_surface
Layer name: Echo_top_cloud_tops
Layer name: Snow_depth_surface
Layer name: Water_equivalent_of_accumulated_snow_depth_surface_1_Hour_Accumulation
Layer name: Wind_speed_gust_surface
Layer name: v-component_of_wind_isobaric
Layer name: Total_precipitation_surface_1_Hour_Accumulation
Layer name: Wind_speed_height_above_ground_1_Hour_Maximum
Layer name: Pressure_reduced_to_MSL_msl
Layer name: Dewpoint_temperature_height_above_ground
Layer name: Convective_inhibition_surface
Layer name: Geopotential_height_adiabatic_condensation_lifted
Layer name: Hourly_Maximum_of_Upward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum
Layer name: Visibility_surface
Layer name: Reflectivity_height_above_ground
Layer name: Convective_inhibition_pressure_difference_layer
Layer name: Geopotential_height_surface
Layer name: Convective_available_potential_energy_pressure_difference_layer
Layer name: Lightning_surface
Layer name: Per_cent_frozen_precipitation_surface
Layer name: v-component_storm_motion_height_above_ground_layer
Layer name: Geopotential_height_cloud_tops
Layer name: Hourly_Maximum_of_Downward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum
Layer name: Total_column_integrated_graupel_entire_atmosphere_single_layer_1_Hour_Maximum
Layer name: Low_cloud_cover_low_cloud
Layer name: Vertical_v-component_shear_height_above_ground_layer
Layer name: Precipitable_water_entire_atmosphere_single_layer
Layer name: Hourly_Maximum_of_Updraft_Helicity_over_Layer_2km_to_5_km_AGL_height_above_ground_layer_1_Hour_Maximum
Layer name: Categorical_freezing_rain_surface
Layer name: Dewpoint_temperature_isobaric
Layer name: Hourly_Maximum_of_Simulated_Reflectivity_at_1_km_AGL_height_above_ground_1_Hour_Maximum
Layer name: Geopotential_height_isobaric
Layer name: Pressure_surface
Layer name: Categorical_snow_surface
Layer name: Composite_reflectivity_entire_atmosphere
Layer name: u-component_of_wind_height_above_ground
Layer name: Vertical_velocity_geometric_sigma_layer_1_Hour_Average
Layer name: Surface_lifted_index_isobaric_layer
Layer name: Vertically_integrated_liquid_water_VIL_entire_atmosphere
Layer name: Pressure_of_level_from_which_parcel_was_lifted_pressure_difference_layer
Layer name: v-component_of_wind_height_above_ground
Layer name: wind @ Specified height level above ground
Layer name: Storm_relative_helicity_height_above_ground_layer
Layer name: High_cloud_cover_high_cloud
Layer name: Categorical_rain_surface
Layer name: Total_cloud_cover_entire_atmosphere
Layer name: Planetary_boundary_layer_height_surface
Layer name: Geopotential_height_cloud_ceiling
Layer name: Vertical_u-component_shear_height_above_ground_layer
Layer name: Temperature_height_above_ground
Layer name: Medium_cloud_cover_middle_cloud

In [4]:

    

temp = wms['Temperature_height_above_ground']

elevations = [elevation.strip() for elevation in temp.elevations]
print(elevations)
# only one elevation, so use it
elevation = elevations[0]
# have to guess the range
color_max = 315  # K
color_min = 273  # K
colorscalerange = '{},{}'.format(color_min,color_max)


    

    




['2.0']

In [5]:

    

times = [time.strip() for time in temp.timepositions]
print(times)
# get the first time - Forecast Hour 0
time = times[0]


    

    




['2015-07-03T15:00:00.000Z', '2015-07-03T16:00:00.000Z', '2015-07-03T17:00:00.000Z', '2015-07-03T18:00:00.000Z', '2015-07-03T19:00:00.000Z', '2015-07-03T20:00:00.000Z', '2015-07-03T21:00:00.000Z', '2015-07-03T22:00:00.000Z', '2015-07-03T23:00:00.000Z', '2015-07-04T00:00:00.000Z', '2015-07-04T01:00:00.000Z', '2015-07-04T02:00:00.000Z', '2015-07-04T03:00:00.000Z', '2015-07-04T04:00:00.000Z', '2015-07-04T05:00:00.000Z', '2015-07-04T06:00:00.000Z']

In [6]:

    

# pick a projection - going with Miller for this example
# note that with Cartopy you are NOT limited to the projections avaliable through ncWMS
plt_proj = cartopy.crs.Miller()

fig, ax = plt.subplots(figsize=(12,12), subplot_kw={'projection': plt_proj})

# Colorbar goodness
cax = fig.add_axes([0.95, 0.3, 0.02, 0.42])
norm = plt.Normalize(vmin=color_min, vmax=color_max)
cmap = plt.cm.gist_heat
cb = mpl.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm, spacing='proportional', orientation='vertical')
cb.set_label('Temperature [K]')

# use bounding box info obtained from the ncWMS service to frame the image
extent = (temp.boundingBox[0], temp.boundingBox[2], temp.boundingBox[1], temp.boundingBox[3])
ax.set_extent(extent)

# ncWMS keywords (which includes the WMS keywords as well)
wms_kwargs = {'colorscalerange': colorscalerange,
              'abovemaxcolor': 'transparent',
              'belowmincolor': 'transparent',
              'transparent': 'true',
              'elevation': elevation,
              'time': time}

# plot the layer using Cartopy's WMS interface
ax.add_wms(wms=serverurl, layers=[temp.name], wms_kwargs=wms_kwargs, cmap=cmap)

# add coastlines, country borders and state outlines
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none'), linestyle=':')


    

    
Out[6]:





<cartopy.mpl.feature_artist.FeatureArtist at 0x10aa5aba8>






    

In [ ]: