Notebook

Exploring Web Map Service (WMS)¶

Image Helper Funcitons
WMS and OWSLib
Getting some information about the service
Getting the basic information we need to perform a GetMap request
More on GetMap request
TDS-ncWMS styles and extensions
WMS and Cartopy

1. Image Helper Functions ¶

At a very basic level, a WMS service returns an image
Let's define a few helper functions to save and load images from a temporairy directory or url

In [1]:

# Function to get tmp image dir. If it does not exist, 
#  create it
def getTmpImgDir():
    from os import makedirs
    from os.path import exists, join
    tmp_img_dir = "tmp_img"
    if (not exists(tmp_img_dir)):
        makedirs(tmp_img_dir)
    return tmp_img_dir

#Function that saves the layer as an image
def saveLayerAsTmpImage(layer, inname):
    from os.path import join
    tmp_img_dir = getTmpImgDir()
    full_img_path = join(tmp_img_dir, inname)
    out = open(full_img_path, 'wb')
    out.write(layer.read())
    out.close()

# Function to load image
def loadTmpImage(image_name):
    from IPython.core.display import Image
    from os.path import join
    tmp_img_dir = getTmpImgDir()
    filename = join(tmp_img_dir, image_name)
    return Image(filename)

# Function to display image from a url
def loadRemoteImage(imgUrl):
    from IPython.core.display import Image
    return Image(url = imgUrl)

2. About WMS and OWSLib ¶

WMS is the Open Geospatial Consortium (OGC) standard interface for requesting georeferenced images through HTTP.
OWSLib is part of geopython, a GitHub organization comprised of Python projects related to geospatial.
OWSLib is a Python package for client programming with OGC Web Services (OWS) developed by Tom Kralidis.
OWSLib supports several OGC standards: WFS, WCS, SOS...and of course WMS 1.1.1. More.

3. Getting some information about the service ¶

We will use OWSLib package and in particular the owslib.wms module.
Within the TDS context, if WMS is enabled and set up in the catalogs, each dataset in the catalog will have a WMS access url
Let's use Siphon to get WMS access url for the latest run of the NCEP HRRR model.

In [2]:

from owslib.wms import WebMapService
from siphon.catalog import get_latest_access_url

# just need a WMS url from one TDS dataset...NCEP HRRR 2.5km forecast model
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')
wms.url

Out[2]:

'http://thredds-jumbo.unidata.ucar.edu/thredds/wms/grib/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_20150722_1200.grib2'

The WebMapService object gets all the information available about the service through a GetCapabilities request:¶

In [3]:

#This is general information, common to all datasets in a TDS server
operations =[ op.name for op in  wms.operations ]
print('Available operations: ')
print(operations)

print('General information (common to all datasets):')
print(wms.identification.type)
print(wms.identification.abstract)
print(wms.identification.keywords)
print(wms.identification.version)
print(wms.identification.title)

Available operations: 
['GetCapabilities', 'GetMap', 'GetFeatureInfo']
General information (common to all datasets):
OGC:WMS
Scientific Data
['meteorology', 'atmosphere', 'climate', 'ocean', 'earth science']
1.1.1
THREDDS Data Server

What's in the dataset?¶

Each variable in a dataset translates into a layer in the WMS service.
Bounding boxes, styles and dimensions are specific to each layer.
The server creates virtual layers if it finds vector components in CF-1 or Grib conventions.

In [4]:

#Listing all available layers...
layers = list(wms.contents)
for l in layers:
    print('Layer title: '+wms[l].title +', name:'+wms[l].name)

Layer title: Composite reflectivity @ Entire atmosphere, name:Composite_reflectivity_entire_atmosphere
Layer title: High cloud cover @ High cloud layer, name:High_cloud_cover_high_cloud
Layer title: Pressure @ Ground or water surface, name:Pressure_surface
Layer title: Geopotential height @ Level of cloud tops, name:Geopotential_height_cloud_tops
Layer title: Pressure reduced to MSL @ Mean sea level, name:Pressure_reduced_to_MSL_msl
Layer title: u-component storm motion @ Specified height level above ground layer, name:u-component_storm_motion_height_above_ground_layer
Layer title: Best (4-layer) lifted index @ Level at specified pressure difference from ground to level layer, name:Best_4-layer_lifted_index_pressure_difference_layer
Layer title: u-component of wind @ Isobaric surface, name:u-component_of_wind_isobaric
Layer title: Total cloud cover @ Entire atmosphere, name:Total_cloud_cover_entire_atmosphere
Layer title: Geopotential height @ Level of adiabatic condensation lifted from the surface, name:Geopotential_height_adiabatic_condensation_lifted
Layer title: v-component storm motion @ Specified height level above ground layer, name:v-component_storm_motion_height_above_ground_layer
Layer title: Wind speed (gust) @ Ground or water surface, name:Wind_speed_gust_surface
Layer title: Categorical snow @ Ground or water surface, name:Categorical_snow_surface
Layer title: Vertical u-component shear @ Specified height level above ground layer, name:Vertical_u-component_shear_height_above_ground_layer
Layer title: Vertically integrated liquid water (VIL) @ Entire atmosphere, name:Vertically_integrated_liquid_water_VIL_entire_atmosphere
Layer title: Water equivalent of accumulated snow depth (1_Hour Accumulation) @ Ground or water surface, name:Water_equivalent_of_accumulated_snow_depth_surface_1_Hour_Accumulation
Layer title: Categorical freezing rain @ Ground or water surface, name:Categorical_freezing_rain_surface
Layer title: Wind speed (1_Hour Maximum) @ Specified height level above ground, name:Wind_speed_height_above_ground_1_Hour_Maximum
Layer title: Storm relative helicity @ Specified height level above ground layer, name:Storm_relative_helicity_height_above_ground_layer
Layer title: Lightning @ Ground or water surface, name:Lightning_surface
Layer title: Low cloud cover @ Low cloud layer, name:Low_cloud_cover_low_cloud
Layer title: Visibility @ Ground or water surface, name:Visibility_surface
Layer title: Surface lifted index @ Isobaric surface layer, name:Surface_lifted_index_isobaric_layer
Layer title: Hourly Maximum of Downward Vertical Velocity in the lowest 400hPa (1_Hour Maximum) @ Level at specified pressure difference from ground to level layer, name:Hourly_Maximum_of_Downward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum
Layer title: Convective available potential energy @ Ground or water surface, name:Convective_available_potential_energy_surface
Layer title: Convective inhibition @ Ground or water surface, name:Convective_inhibition_surface
Layer title: Precipitable water @ Entire atmosphere layer, name:Precipitable_water_entire_atmosphere_single_layer
Layer title: Medium cloud cover @ Middle cloud layer, name:Medium_cloud_cover_middle_cloud
Layer title: Temperature @ Isobaric surface, name:Temperature_isobaric
Layer title: Snow depth @ Ground or water surface, name:Snow_depth_surface
Layer title: Categorical rain @ Ground or water surface, name:Categorical_rain_surface
Layer title: Total column integrated graupel (1_Hour Maximum) @ Entire atmosphere layer, name:Total_column_integrated_graupel_entire_atmosphere_single_layer_1_Hour_Maximum
Layer title: Reflectivity @ Specified height level above ground, name:Reflectivity_height_above_ground
Layer title: Geopotential height @ Cloud ceiling, name:Geopotential_height_cloud_ceiling
Layer title: Echo top @ Level of cloud tops, name:Echo_top_cloud_tops
Layer title: Pressure of level from which parcel was lifted @ Level at specified pressure difference from ground to level layer, name:Pressure_of_level_from_which_parcel_was_lifted_pressure_difference_layer
Layer title: wind @ Isobaric surface, name:wind @ Isobaric surface
Layer title: Geopotential height @ Ground or water surface, name:Geopotential_height_surface
Layer title: Dewpoint temperature @ Isobaric surface, name:Dewpoint_temperature_isobaric
Layer title: Vertical velocity (geometric) (1_Hour Average) @ Sigma level layer, name:Vertical_velocity_geometric_sigma_layer_1_Hour_Average
Layer title: v-component of wind @ Specified height level above ground, name:v-component_of_wind_height_above_ground
Layer title: v-component of wind @ Isobaric surface, name:v-component_of_wind_isobaric
Layer title: Hourly Maximum of Updraft Helicity over Layer 2km to 5 km AGL (1_Hour Maximum) @ Specified height level above ground layer, name:Hourly_Maximum_of_Updraft_Helicity_over_Layer_2km_to_5_km_AGL_height_above_ground_layer_1_Hour_Maximum
Layer title: Dewpoint temperature @ Specified height level above ground, name:Dewpoint_temperature_height_above_ground
Layer title: Geopotential height @ Isobaric surface, name:Geopotential_height_isobaric
Layer title: Planetary boundary layer height @ Ground or water surface, name:Planetary_boundary_layer_height_surface
Layer title: Temperature @ Specified height level above ground, name:Temperature_height_above_ground
Layer title: Categorical ice pellets @ Ground or water surface, name:Categorical_ice_pellets_surface
Layer title: Hourly Maximum of Upward Vertical Velocity in the lowest 400hPa (1_Hour Maximum) @ Level at specified pressure difference from ground to level layer, name:Hourly_Maximum_of_Upward_Vertical_Velocity_in_the_lowest_400hPa_pressure_difference_layer_1_Hour_Maximum
Layer title: wind @ Specified height level above ground, name:wind @ Specified height level above ground
Layer title: Convective inhibition @ Level at specified pressure difference from ground to level layer, name:Convective_inhibition_pressure_difference_layer
Layer title: Convective available potential energy @ Level at specified pressure difference from ground to level layer, name:Convective_available_potential_energy_pressure_difference_layer
Layer title: Total precipitation (1_Hour Accumulation) @ Ground or water surface, name:Total_precipitation_surface_1_Hour_Accumulation
Layer title: Vertical v-component shear @ Specified height level above ground layer, name:Vertical_v-component_shear_height_above_ground_layer
Layer title: Per cent frozen precipitation @ Ground or water surface, name:Per_cent_frozen_precipitation_surface
Layer title: u-component of wind @ Specified height level above ground, name:u-component_of_wind_height_above_ground
Layer title: Hourly Maximum of Simulated Reflectivity at 1 km AGL (1_Hour Maximum) @ Specified height level above ground, name:Hourly_Maximum_of_Simulated_Reflectivity_at_1_km_AGL_height_above_ground_1_Hour_Maximum

4. Getting the basic information we need to perform a GetMap request ¶

All the information that clients need is available in the capabilities document, which is stored in the WebMapService object.
The WMS included in the TDS (ncWMS) only supports GetMap requests on one layer (variable).
We need to choose our layer, bounding box, spatial reference system (SRS), the size, and format of the image.

What is available?¶

In [5]:

#Values common to all GetMap requests: formats and http methods:
print(wms.getOperationByName('GetMap').formatOptions)
print(wms.getOperationByName('GetMap').methods)

#Let's choose: 'wind @ Isobaric surface' (the value in the parameter must be name of the layer)
wind = wms['wind @ Isobaric surface']

#What is its bounding box?
print(wind.boundingBox)

#available CRS
print(wind.crsOptions)
# --> NOT ALL THE AVAILABLE CRS OPTIONS ARE LISTED
# print elevations at which the wind is avaliable
elevations = [elevation.strip() for elevation in wind.elevations]
print(elevations)
# units on the elevation are not exposed. 50,000 for a height -> probably Pa

['image/png', 'image/png;mode=32bit', 'image/gif', 'image/jpeg', 'application/vnd.google-earth.kmz']
[{'type': 'Get', 'url': 'http://thredds-jumbo.unidata.ucar.edu/thredds/wms/grib/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_20150722_1200.grib2'}]
(-130.12289354741372, 20.17862215527375, -60.86645064452432, 50.11284676676232, 'EPSG:4326')
['EPSG:3857', 'CRS:84', 'EPSG:32661', 'EPSG:3409', 'EPSG:32761', 'EPSG:4326', 'EPSG:3408', 'EPSG:41001', 'EPSG:27700']
['50000.0', '70000.0', '85000.0', '92500.0', '100000.0']

Ok, let's request an image from the layer!¶

In [6]:

img_wind = wms.getmap( layers=[wind.name], #only takes one layer
                  srs='EPSG:4326',
                  bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
                  size=(512, 512),
                  format='image/png',
                  elevation=elevations[0])

# Save it..
saveLayerAsTmpImage(img_wind, 'test_wind.png')

# Display the image we've just saved...
loadTmpImage('test_wind.png')

Out[6]:

5. More on GetMap requests ¶

Handling time and vertical dimensions
Changing styles
Changing the spatial reference system (SRS)

Handling time and vertical dimensions¶

Getting available times and elevations for a layer
We can choose any of the available times and make a request for it with the parameter time
If no time is provided the default in TDS is the closest available time to the current time

In [7]:

times = [time.strip() for time in wind.timepositions]
print(times)
elevations = [elevation.strip() for elevation in wind.elevations]
print("")
print(elevations)

['2015-07-22T12:00:00.000Z', '2015-07-22T13:00:00.000Z', '2015-07-22T14:00:00.000Z', '2015-07-22T15:00:00.000Z', '2015-07-22T16:00:00.000Z', '2015-07-22T17:00:00.000Z', '2015-07-22T18:00:00.000Z', '2015-07-22T19:00:00.000Z', '2015-07-22T20:00:00.000Z', '2015-07-22T21:00:00.000Z', '2015-07-22T22:00:00.000Z', '2015-07-22T23:00:00.000Z', '2015-07-23T00:00:00.000Z', '2015-07-23T01:00:00.000Z', '2015-07-23T02:00:00.000Z', '2015-07-23T03:00:00.000Z']

['50000.0', '70000.0', '85000.0', '92500.0', '100000.0']

In [8]:

img_wind = wms.getmap( layers=[wind.name], 
                  srs='EPSG:4326',
                  bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
                  size=(600, 600),
                  format='image/png',
                  elevation=elevations[0],
                  time= times[-1]
)

saveLayerAsTmpImage(img_wind, 'test_wind.png')
loadTmpImage('test_wind.png')

Out[8]:

Animations¶

We can also specify a time interval to get an animated gif
Format must be image/gif

In [9]:

img_wind = wms.getmap( layers=[wind.name], 
                  srs='EPSG:4326',
                  bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
                  size=(600, 600),
                  format='image/gif',
                  elevation=elevations[0],
                  time= "{}/{}".format(times[0],times[4])
)

from IPython.core.display import Image
loadRemoteImage(img_wind.geturl())

Out[9]:

Getting the available vertical levels:

OWSLib does not support vertical levels, meaning the layer objects do not have a property "elevations" with the vertical levels. So, we need a little extra work to get the available vertical levels for a layer

Changing styles¶

We can specify the style (any from the available styles for a layer) in the param styles

In [10]:

#available styles: 
#print wind.styles
#Change the style of our layer
img_wind = wms.getmap( layers=[wind.name], #only takes one layer
                  styles=['barb/rainbow'], #one style per layer    
                  srs='EPSG:4326',
                  bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
                  size=(600, 600),
                  format='image/png',
                  elevation=elevations[0],
                  time= times[0]
)

saveLayerAsTmpImage(img_wind, 'test_wind_barb.png')
loadTmpImage('test_wind_barb.png')

Out[10]:

Changing the spatial reference system (SRS)¶

We can reproject to any of the available SRS.

In [11]:

#Reproject the bounding box to a global mercator (EPSG:3875, projection used by Google Maps, OSM...) using cartopy
import cartopy.crs as ccrs
epsg = 3857
psproj = ccrs.epsg(epsg)
xmin, ymin = psproj.transform_point(wind.boundingBox[0], wind.boundingBox[1], ccrs.Geodetic())
xmax, ymax = psproj.transform_point(wind.boundingBox[2], wind.boundingBox[3], ccrs.Geodetic())
img_wind = wms.getmap( layers=[wind.name],
                  srs='EPSG:%d' % epsg,
                  bbox=(xmin, ymin, xmax, ymax),
                  size=(600, 600),
                  format='image/png',
                  elevation=elevations[0],
                  time= times[0]
)

saveLayerAsTmpImage(img_wind, 'test_wind_3857.png')
loadTmpImage('test_wind_3857.png')

Out[11]:

Finally, let's see how much work it is to change our layer¶

In [12]:

temp = wms['Temperature_height_above_ground']
img_temp = wms.getmap( layers=[temp.name], 
                  styles=['boxfill/rainbow'], 
                  srs='EPSG:4326',
                  bbox=(temp.boundingBox[0],temp.boundingBox[1], temp.boundingBox[2], temp.boundingBox[3]),
                  size=(600, 600),
                  format='image/png',
                  elevation=2,
                  time= times[0]
)

saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')

Out[12]:

6. TDS-ncWMS styles and extensions

ncWMS/THREDDS provides some [non-standard WMS parameters](http://www.resc.rdg.ac.uk/trac/ncWMS/wiki/WmsExtensions) that allow clients some control on the styling.
Change the scale range:
Default is -50,50. Parameter colorscalerange allows us to use a different scale

In [13]:

img_temp = wms.getmap( layers=[temp.name], 
                  styles=['boxfill/rainbow'], 
                  srs='EPSG:4326',
                  bbox=(temp.boundingBox[0],temp.boundingBox[1], temp.boundingBox[2], temp.boundingBox[3]),
                  size=(600, 600),
                  format='image/png',
                  time= times[0],
                  colorscalerange='250,320'
)

saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')

Out[13]:

More control over colors ¶

abovemaxcolor, belowmincolor params give us control on how we want the values out of range to be displayed.
valid values for those params are: extend (will use the highest/lowest value of the palette for values larger/smaller than the maximun/minimun), transparent and a color in 0xRRGGBB format
sytles: defines the colorbar used

In [14]:

colorscalerange='290,310'

img_temp = wms.getmap( layers=[temp.name], 
                  styles=['boxfill/rainbow'], 
                  srs='EPSG:4326',
                  bbox=(temp.boundingBox[0],temp.boundingBox[1], temp.boundingBox[2], temp.boundingBox[3]),
                  size=(600, 600),
                  format='image/png',
                  time= times[0],
                  colorscalerange=colorscalerange,
                  abovemaxcolor='transparent',
                  belowmincolor='transparent'
)
saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')

Out[14]:

The GetLegendGraphic request gives us a legend for the map, but the request is not supported by OWSLib.¶

In [15]:

params ={'request': 'GetLegendGraphic',
         'colorbaronly':'False', #want the text in the legend
         'layer':temp.name,
         'colorscalerange':colorscalerange}

legendUrl=serverurl+'?REQUEST={request:s}&COLORBARONLY={colorbaronly:s}&LAYER={layer:s}&COLORSCALERANGE={colorscalerange:s}'.format(**params)

loadRemoteImage(legendUrl)

Out[15]:

7. ncWMS and Cartopy ¶

We can use Cartopy geolocate a layer and overlay it with political borders and coastlines
Note that Cartopy directly use the return from the WMS service, so no need to use the helper functions to save an image to a temporary location
Cartopy is *not* limited to using the colorbars from ncWMS.
Cartopy is *not* limited to using the projections enabled by ncWMS.

In [16]:

%matplotlib inline
import cartopy
import matplotlib as mpl
import matplotlib.pyplot as plt

# get the first time - Forecast Hour 0
times = [time.strip() for time in temp.timepositions]
time = times[0]
# only one elevation, so use it
elevations = [elevation.strip() for elevation in temp.elevations]
elevation = elevations[0]
# have to guess the range
color_max = 315  # K
color_min = 273  # K
colorscalerange = '{},{}'.format(color_min,color_max)

# 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
cmap.set_under('None') # Colors below minimum should be transparent
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
im = ax.add_wms(wms=serverurl, layers=[temp.name], wms_kwargs=wms_kwargs, cmap=cmap, zorder=2)

# Range for image data should above pixels with 0 value. This makes the 0 values
# transparent
im.norm.vmin = 1

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

Out[16]:

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

Exercise (Ideas):¶

What happens if you remove all of the "zorder" keyword args in cell 16?
Plot a new layer (look at cell 4 for the available layers).
How would you request a specific layer for a variable with multiple vertical levels (for example, "Temperature_isobaric")
Use metpy to change the units from Kelvin to degC or degF in cell 16