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') #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 #Listing all available layers... layers = list(wms.contents) for l in layers: print 'Layer title: '+wms[l].title +', name:'+wms[l].name #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 #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 #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) #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 #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') 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') 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') 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) 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()