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)
In [2]:
from owslib.wms import WebMapService
#We just need a WMS url from one TDS dataset...
serverurl ='http://thredds.ucar.edu/thredds/wms/grib/NCEP/RAP/CONUS_13km/Best'
wms = WebMapService( serverurl, version='1.1.1')
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
In [4]:
#Listing all available layers...
layers = list(wms.contents)
for l in layers:
print 'Layer title: '+wms[l].title +', name:'+wms[l].name
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
In [6]:
#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..
saveLayerAsTmpImage(img_wind, 'test_wind.png')
#Display the image we've just saved...
loadTmpImage('test_wind.png')
Out[6]:
In [7]:
#Times are available in the timepositions property of the layer
times= [time.strip() for time in wind.timepositions]
print times
In [8]:
#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]
)
saveLayerAsTmpImage(img_wind, 'test_wind.png')
loadTmpImage('test_wind.png')
Out[8]:
In [9]:
#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]
)
from IPython.core.display import Image
loadRemoteImage(img_wind.url)
Out[9]:
In [10]:
#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
In [11]:
#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 ]
)
saveLayerAsTmpImage(img_wind, 'test_wind.png')
loadTmpImage('test_wind.png')
Out[11]:
In [12]:
#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]
)
saveLayerAsTmpImage(img_wind, 'test_wind_barb.png')
loadTmpImage('test_wind_barb.png')
Out[12]:
In [13]:
#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]
)
saveLayerAsTmpImage(img_wind, 'test_wind_3857.png')
loadTmpImage('test_wind_3857.png')
Out[13]:
Cool, we already know how to make get map requests. Let's change our layer...
In [14]:
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]
)
saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')
Out[14]:
...well not that cool.
ncWMS/THREDDS provides some non-standard WMS parameters that allow clients some control on the styling.
Change the scale range:
In [15]:
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'
)
saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')
Out[15]:
In [16]:
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'
)
saveLayerAsTmpImage(img_temp, 'test_temp.png')
loadTmpImage('test_temp.png')
Out[16]:
The GetLegendGraphic request gives us a legend for the map, but the request is not supported by OWSLib.
In [17]:
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[17]:
We can use basemap to overlay the layer with a coastline...
In [18]:
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
%matplotlib inline
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)
saveLayerAsTmpImage(legend, 'legend_temp.png')
#read the image as an array
arr = read_png(os.path.join(getTmpImgDir(),'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()