from owslib.wms import WebMapService
#We just need a WMS url from one TDS dataset...
serverurl ='http://thredds.ucar.edu/thredds/wms/grib/NCEP/NAM/CONUS_12km/best'
wms = WebMapService( serverurl, version='1.1.1')
The WebMapService object gets all the information available about the service through a GetCapabilities request:
#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
#Listing all available layers...
layers = list(wms.contents)
for l in layers:
print 'Layer title: '+wms[l].title +', name:'+wms[l].name
Layer title: Vertical velocity (pressure) @ Isobaric surface, name:Vertical_velocity_pressure_isobaric Layer title: Thunderstorm probability (3_Hour Accumulation) @ Ground or water surface, name:Thunderstorm_probability_surface_3_Hour_Accumulation 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: Geopotential height @ Level of 0C isotherm, name:Geopotential_height_zeroDegC_isotherm Layer title: Temperature @ Tropopause, name:Temperature_tropopause Layer title: Categorical Snow @ Ground or water surface, name:Categorical_Snow_surface Layer title: Convective precipitation (3_Hour Accumulation) @ Ground or water surface, name:Convective_precipitation_surface_3_Hour_Accumulation Layer title: Relative humidity @ Specified height level above ground, name:Relative_humidity_height_above_ground Layer title: Pressure @ Tropopause, name:Pressure_tropopause Layer title: Snow depth @ Ground or water surface, name:Snow_depth_surface Layer title: u-component of wind @ Specified height level above ground, name:u-component_of_wind_height_above_ground Layer title: Minimum temperature @ Specified height level above ground, name:Minimum_temperature_height_above_ground Layer title: Temperature @ Level at specified pressure difference from ground to level layer, name:Temperature_pressure_difference_layer Layer title: wind @ Tropopause, name:wind @ Tropopause 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: Convective available potential energy @ Ground or water surface, name:Convective_available_potential_energy_surface Layer title: Maximum temperature @ Specified height level above ground, name:Maximum_temperature_height_above_ground Layer title: Pressure reduced to MSL @ Mean sea level, name:Pressure_reduced_to_MSL_msl Layer title: u-component of wind @ Maximum wind level, name:u-component_of_wind_maximum_wind Layer title: Convective inhibition @ Level at specified pressure difference from ground to level layer, name:Convective_inhibition_pressure_difference_layer Layer title: MSLP (Eta model reduction) @ Mean sea level, name:MSLP_Eta_model_reduction_msl Layer title: Convective inhibition @ Ground or water surface, name:Convective_inhibition_surface Layer title: Relative humidity @ Isobaric surface, name:Relative_humidity_isobaric Layer title: Geopotential height @ Isobaric surface, name:Geopotential_height_isobaric Layer title: Composite reflectivity @ Entire atmosphere layer, name:Composite_reflectivity_entire_atmosphere Layer title: v-component of wind @ Tropopause, name:v-component_of_wind_tropopause Layer title: Parcel lifted index (to 500 hPa) @ Level at specified pressure difference from ground to level layer, name:Parcel_lifted_index_to_500_hPa_pressure_difference_layer Layer title: Wind speed @ Specified height level above ground, name:Wind_speed_height_above_ground Layer title: wind @ Specified height level above ground, name:wind @ Specified height level above ground Layer title: Water equivalent of accumulated snow depth @ Ground or water surface, name:Water_equivalent_of_accumulated_snow_depth_surface Layer title: v-component of wind @ Specified height level above ground, name:v-component_of_wind_height_above_ground Layer title: u-component of wind @ Level at specified pressure difference from ground to level layer, name:u-component_of_wind_pressure_difference_layer Layer title: Probability of Freezing Precipitation (3_Hour Accumulation) @ Ground or water surface, name:Probability_of_Freezing_Precipitation_surface_3_Hour_Accumulation Layer title: v-component of wind @ Maximum wind level, name:v-component_of_wind_maximum_wind Layer title: Relative humidity @ Sigma level layer, name:Relative_humidity_sigma_layer Layer title: Probability of 0.01 inch of precipitation (POP) (3_Hour Accumulation) @ Ground or water surface, name:Probability_of_0p01_inch_of_precipitation_POP_surface_3_Hour_Accumulation Layer title: Relative humidity @ Level of 0C isotherm, name:Relative_humidity_zeroDegC_isotherm Layer title: Categorical Ice Pellets @ Ground or water surface, name:Categorical_Ice_Pellets_surface Layer title: Total precipitation (3_Hour Accumulation) @ Ground or water surface, name:Total_precipitation_surface_3_Hour_Accumulation Layer title: Dewpoint temperature @ Specified height level above ground, name:Dewpoint_temperature_height_above_ground Layer title: v-component of wind @ Level at specified pressure difference from ground to level layer, name:v-component_of_wind_pressure_difference_layer Layer title: u-component of wind @ Tropopause, name:u-component_of_wind_tropopause Layer title: Categorical Freezing Rain @ Ground or water surface, name:Categorical_Freezing_Rain_surface Layer title: v-component of wind @ Isobaric surface, name:v-component_of_wind_isobaric Layer title: Pressure @ Ground or water surface, name:Pressure_surface Layer title: wind @ Level at specified pressure difference from ground to level layer, name:wind @ Level at specified pressure difference from ground to level layer Layer title: u-component of wind @ Isobaric surface, name:u-component_of_wind_isobaric Layer title: Total cloud cover @ Entire atmosphere layer, name:Total_cloud_cover_entire_atmosphere Layer title: Reflectivity @ Hybrid level, name:Reflectivity_hybrid Layer title: wind @ Isobaric surface, name:wind @ Isobaric surface Layer title: Wind direction (from which blowing) @ Specified height level above ground, name:Wind_direction_from_which_blowing_height_above_ground Layer title: Reflectivity @ Specified height level above ground, name:Reflectivity_height_above_ground Layer title: Pressure @ Maximum wind level, name:Pressure_maximum_wind Layer title: Absolute vorticity @ Isobaric surface, name:Absolute_vorticity_isobaric Layer title: Surface Lifted Index @ Isobaric surface layer, name:Surface_Lifted_Index_isobaric_layer Layer title: Precipitable water @ Entire atmosphere layer, name:Precipitable_water_entire_atmosphere Layer title: Storm relative helicity @ Specified height level above ground layer, name:Storm_relative_helicity_height_above_ground_layer Layer title: Probability of Frozen Precipitation (3_Hour Accumulation) @ Ground or water surface, name:Probability_of_Frozen_Precipitation_surface_3_Hour_Accumulation Layer title: Visibility @ Ground or water surface, name:Visibility_surface Layer title: Temperature @ Specified height level above ground, name:Temperature_height_above_ground Layer title: Categorical Rain @ Ground or water surface, name:Categorical_Rain_surface Layer title: Relative humidity @ Level at specified pressure difference from ground to level layer, name:Relative_humidity_pressure_difference_layer Layer title: wind @ Maximum wind level, name:wind @ Maximum wind level Layer title: Temperature @ Isobaric surface, name:Temperature_isobaric
#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
['image/png', 'image/png;mode=32bit', 'image/gif', 'image/jpeg', 'application/vnd.google-earth.kmz'] {'Get': {'url': 'http://thredds.ucar.edu/thredds/wms/grib/NCEP/NAM/CONUS_12km/best'}} (-152.9859023956905, 12.123672938563633, -49.308990314247126, 57.35625318852111, 'EPSG:4326') ['EPSG:3857', 'EPSG:32761', 'EPSG:4326', 'EPSG:32661', 'EPSG:41001', 'CRS:84', 'EPSG:3409', 'EPSG:3408', 'EPSG:27700']
#Function that saves the layer as an image
def saveLayerAsImage(layer, inname):
out = open(inname, 'wb')
out.write(layer.read())
out.close()
#let's get the image...
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'
)
#Save it..
saveLayerAsImage(img_wind, 'test_wind.png')
#Display the image we've just saved...
from IPython.core.display import Image
Image(filename='test_wind.png')
#Times are available in the timepositions property of the layer
times= [time.strip() for time in wind.timepositions]
print times
['2013-11-27T00:00:00.000Z', '2013-11-27T03:00:00.000Z', '2013-11-27T06:00:00.000Z', '2013-11-27T09:00:00.000Z', '2013-11-27T12:00:00.000Z', '2013-11-27T15:00:00.000Z', '2013-11-27T18:00:00.000Z', '2013-11-27T21:00:00.000Z', '2013-11-28T00:00:00.000Z', '2013-11-28T03:00:00.000Z', '2013-11-28T06:00:00.000Z', '2013-11-28T09:00:00.000Z', '2013-11-28T12:00:00.000Z', '2013-11-28T15:00:00.000Z', '2013-11-28T18:00:00.000Z', '2013-11-28T21:00:00.000Z', '2013-11-29T00:00:00.000Z', '2013-11-29T03:00:00.000Z', '2013-11-29T06:00:00.000Z', '2013-11-29T09:00:00.000Z', '2013-11-29T12:00:00.000Z', '2013-11-29T15:00:00.000Z', '2013-11-29T18:00:00.000Z', '2013-11-29T21:00:00.000Z', '2013-11-30T00:00:00.000Z', '2013-11-30T03:00:00.000Z', '2013-11-30T06:00:00.000Z', '2013-11-30T09:00:00.000Z', '2013-11-30T12:00:00.000Z', '2013-11-30T15:00:00.000Z', '2013-11-30T18:00:00.000Z', '2013-11-30T21:00:00.000Z', '2013-12-01T00:00:00.000Z', '2013-12-01T03:00:00.000Z', '2013-12-01T06:00:00.000Z', '2013-12-01T09:00:00.000Z', '2013-12-01T12:00:00.000Z', '2013-12-01T15:00:00.000Z', '2013-12-01T18:00:00.000Z', '2013-12-01T21:00:00.000Z', '2013-12-02T00:00:00.000Z', '2013-12-02T03:00:00.000Z', '2013-12-02T06:00:00.000Z', '2013-12-02T09:00:00.000Z', '2013-12-02T12:00:00.000Z', '2013-12-02T15:00:00.000Z', '2013-12-02T18:00:00.000Z', '2013-12-02T21:00:00.000Z', '2013-12-03T00:00:00.000Z', '2013-12-03T03:00:00.000Z', '2013-12-03T06:00:00.000Z', '2013-12-03T09:00:00.000Z', '2013-12-03T12:00:00.000Z', '2013-12-03T15:00:00.000Z', '2013-12-03T18:00:00.000Z', '2013-12-03T21:00:00.000Z', '2013-12-04T00:00:00.000Z', '2013-12-04T03:00:00.000Z', '2013-12-04T06:00:00.000Z', '2013-12-04T09:00:00.000Z', '2013-12-04T12:00:00.000Z', '2013-12-04T15:00:00.000Z', '2013-12-04T18:00:00.000Z', '2013-12-04T21:00:00.000Z', '2013-12-05T00:00:00.000Z', '2013-12-05T03:00:00.000Z', '2013-12-05T06:00:00.000Z', '2013-12-05T09:00:00.000Z', '2013-12-05T12:00:00.000Z', '2013-12-05T15:00:00.000Z', '2013-12-05T18:00:00.000Z', '2013-12-05T21:00:00.000Z', '2013-12-06T00:00:00.000Z', '2013-12-06T03:00:00.000Z', '2013-12-06T06:00:00.000Z', '2013-12-06T09:00:00.000Z', '2013-12-06T12:00:00.000Z', '2013-12-06T15:00:00.000Z', '2013-12-06T18:00:00.000Z', '2013-12-06T21:00:00.000Z', '2013-12-07T00:00:00.000Z', '2013-12-07T03:00:00.000Z', '2013-12-07T06:00:00.000Z', '2013-12-07T09:00:00.000Z', '2013-12-07T12:00:00.000Z', '2013-12-07T15:00:00.000Z', '2013-12-07T18:00:00.000Z', '2013-12-07T21:00:00.000Z', '2013-12-08T00:00:00.000Z', '2013-12-08T03:00:00.000Z', '2013-12-08T06:00:00.000Z', '2013-12-08T09:00:00.000Z', '2013-12-08T12:00:00.000Z', '2013-12-08T15:00:00.000Z', '2013-12-08T18:00:00.000Z', '2013-12-08T21:00:00.000Z', '2013-12-09T00:00:00.000Z', '2013-12-09T03:00:00.000Z', '2013-12-09T06:00:00.000Z', '2013-12-09T09:00:00.000Z', '2013-12-09T12:00:00.000Z', '2013-12-09T15:00:00.000Z', '2013-12-09T18:00:00.000Z', '2013-12-09T21:00:00.000Z', '2013-12-10T00:00:00.000Z', '2013-12-10T03:00:00.000Z', '2013-12-10T06:00:00.000Z', '2013-12-10T09:00:00.000Z', '2013-12-10T12:00:00.000Z', '2013-12-10T15:00:00.000Z', '2013-12-10T18:00:00.000Z', '2013-12-10T21:00:00.000Z', '2013-12-11T00:00:00.000Z', '2013-12-11T03:00:00.000Z', '2013-12-11T06:00:00.000Z', '2013-12-11T09:00:00.000Z', '2013-12-11T12:00:00.000Z', '2013-12-11T15:00:00.000Z', '2013-12-11T18:00:00.000Z', '2013-12-11T21:00:00.000Z', '2013-12-12T00:00:00.000Z', '2013-12-12T03:00:00.000Z', '2013-12-12T06:00:00.000Z', '2013-12-12T09:00:00.000Z', '2013-12-12T12:00:00.000Z', '2013-12-12T15:00:00.000Z', '2013-12-12T18:00:00.000Z', '2013-12-12T21:00:00.000Z', '2013-12-13T00:00:00.000Z', '2013-12-13T03:00:00.000Z', '2013-12-13T06:00:00.000Z', '2013-12-13T09:00:00.000Z', '2013-12-13T12:00:00.000Z', '2013-12-13T15:00:00.000Z', '2013-12-13T18:00:00.000Z', '2013-12-13T21:00:00.000Z', '2013-12-14T00:00:00.000Z', '2013-12-14T03:00:00.000Z', '2013-12-14T06:00:00.000Z', '2013-12-14T09:00:00.000Z', '2013-12-14T12:00:00.000Z', '2013-12-14T15:00:00.000Z', '2013-12-14T18:00:00.000Z']
#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
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',
time= times[len(times)-1]
)
saveLayerAsImage(img_wind, 'test_wind.png')
Image(filename='test_wind.png')
#We can also specify a time interval to get an animated gif
#Format must be image/gif
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',
time= times[len(times)-4]+'/'+times[len(times)-1]
)
#Image(url='http://python.org/images/python-logo.gif')
#saveLayerAsImage(img_wind, 'test_anim_wind.gif')
Image(url=img_wind.url)
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
#Next version of OWSLib will support this...
#elevations = [el.strip() for el in wind.elevations]
#print elevations
#In the meantime...
def find_elevations_for_layer(wms, layer_name):
"""
parses the wms capabilities document searching
the elevation dimension for the layer
"""
#Get all the layers
levels =None;
layers = wms._capabilities.findall(".//Layer")
layer_tag = None
for el in layers:
name = el.find("Name")
if name is not None and name.text.strip() == layer_name:
layer_tag = el
break
if layer_tag is not None:
elevation_tag = layer_tag.find("Extent[@name='elevation']")
if elevation_tag is not None:
levels = elevation_tag.text.strip().split(',')
return levels;
elevations = find_elevations_for_layer(wms, wind.name)
print elevations
['100000.0', '97500.0', '95000.0', '92500.0', '90000.0', '87500.0', '85000.0', '82500.0', '80000.0', '77500.0', '75000.0', '72500.0', '70000.0', '67500.0', '65000.0', '62500.0', '60000.0', '57500.0', '55000.0', '52500.0', '50000.0', '45000.0', '40000.0', '35000.0', '30000.0', '25000.0', '20000.0', '15000.0', '10000.0']
#now we can change our vertical level with the parameter elevation
#If no elevation parameter is provided the default is the first vertical level in the dimension.
img_wind = wms.getmap( layers=['wind @ Isobaric surface'], #only takes one layer
srs='EPSG:4326',
bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
size=(600, 600),
format='image/png',
time= times[0],
elevation=elevations[len(elevations)-1 ]
)
saveLayerAsImage(img_wind, 'test_wind.png')
Image(filename='test_wind.png')
#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',
time= times[0]
)
saveLayerAsImage(img_wind, 'test_wind_barb.png')
Image(filename='test_wind_barb.png')
#Reproject the bounding box to a global mercator (EPSG:3875, projection used by Google Maps, OSM...) using pyproj
from mpl_toolkits.basemap import pyproj
epsg = '3857'
psproj = pyproj.Proj(init="epsg:%s" % epsg)
xmin, ymin = psproj(wind.boundingBox[0], wind.boundingBox[1])
xmax, ymax = psproj(wind.boundingBox[2], wind.boundingBox[3])
img_wind = wms.getmap( layers=[wind.name],
srs='EPSG:'+ epsg,
bbox=(xmin, ymin, xmax, ymax),
size=(600, 600),
format='image/png',
time= times[0]
)
saveLayerAsImage(img_wind, 'test_wind_3857.png')
Image(filename='test_wind_3857.png')
Cool, we already know how to make get map requests. Let's change our layer...
temp =wms['Temperature_isobaric']
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]
)
saveLayerAsImage(img_temp, 'test_temp.png')
Image(filename='test_temp.png')
...well not that cool.
img_temp = wms.getmap( layers=[temp.name],
styles=['boxfill/rainbow'],
srs='EPSG:4326',
bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
size=(600, 600),
format='image/png',
time= times[0],
colorscalerange='250,320'
)
saveLayerAsImage(img_temp, 'test_temp.png')
Image(filename='test_temp.png')
colorscalerange='290,310'
img_temp = wms.getmap( layers=[temp.name],
styles=['boxfill/rainbow'],
srs='EPSG:4326',
bbox=(wind.boundingBox[0],wind.boundingBox[1], wind.boundingBox[2], wind.boundingBox[3]),
size=(600, 600),
format='image/png',
time= times[0],
colorscalerange=colorscalerange,
abovemaxcolor='transparent',
belowmincolor='transparent'
)
saveLayerAsImage(img_temp, 'test_temp.png')
Image(filename='test_temp.png')
The GetLegendGraphic request gives us a legend for the map, but the request is not supported by OWSLib.
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)
Image(url=legendUrl)
We can use basemap to overlay the layer with a coastline...
import os
import urllib2
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnnotationBbox, OffsetImage
from matplotlib._png import read_png
m = Basemap(llcrnrlon=temp.boundingBox[0], llcrnrlat=temp.boundingBox[1],
urcrnrlon=temp.boundingBox[2], urcrnrlat=temp.boundingBox[3]+5.0,
resolution='l',epsg=4326)
plt.figure(1, figsize=(16,12))
plt.title(temp.title +' '+times[0] )
m.wmsimage(serverurl,xpixels=600, ypixels=600, verbose=False,
layers=[temp.name],
styles=['boxfill/rainbow'],
time= times[0],
colorscalerange=colorscalerange,
abovemaxcolor='extend',
belowmincolor='transparent'
)
m.drawcoastlines(linewidth=0.25)
#Annotating the map with the legend
#Save the legend as image
cwd = os.getcwd()
legend = urllib2.urlopen(legendUrl)
saveLayerAsImage(legend, 'legend_temp.png')
#read the image as an array
arr = read_png('legend_temp.png')
imagebox = OffsetImage(arr, zoom=0.7)
xy =[ temp.boundingBox[2], temp.boundingBox[1] ]
#Gets the current axis
ax = plt.gca()
#Creates the annotation
ab = AnnotationBbox(imagebox, xy,
xybox=(-46.,100.),
xycoords='data',
boxcoords="offset points",
pad=0.)
#Adds the legend image as an AnnotationBbox to the map
ax.add_artist(ab)
plt.show()