This example demonstrates connecting to a NIWA Web Map Service from a Python Script. The WMS Service allow the selection of layers filtered by location and parameters and returns styled image layers that can be overplotted on a map.
This example uses the stations catalogue feed to show the location of climate and hydrology stations.
OWSLib is a Python package for client programming with Open Geospatial Consortium (OGC) web service (hence OWS) interface standards, and their related content models.
http://geopython.github.io/OWSLib/
easy_install OWSLib
see http://nbviewer.ipython.org/github/rsignell-usgs/notebook/blob/master/wms_sample.ipynb for more examples
# NIWA WMS and WFS Server endpoints are at
niwa_ows = "http://gs.niwa.co.nz/geoserver/ows"
niwa_ows = "http://gs.niwa.co.nz/geoserver/stations/ows"
%matplotlib inline
from owslib.wms import WebMapService
wms = WebMapService(niwa_ows, version='1.1.1')
print wms.identification.title
print wms.identification.type
print wms.identification.abstract
print wms.identification.keywords
print wms.identification.version
GeoServer Web Map Service OGC:WMS A compliant implementation of WMS plus most of the SLD extension (dynamic styling). Can also generate PDF, SVG, KML, GeoRSS ['WFS', 'WMS', 'GEOSERVER'] 1.1.1
#Listing all available layers...
layers = list(wms.contents)
for l in layers:
print wms[l].name, wms[l].title
timeseries_types Stations Catalogue - Time Series Types organisations Stations Catalogue - Organisations stations_360 Stations Catalogue - public stations (0 - 360) timeseries_channels timeseries_channels stations_publish_to_sos Stations Catalogue - Public SOS Stations stations_180 Stations Catalogue - stations station_types Stations Catalogue - Station Types
def printlayer(layer) :
print "Title\t", wms[layer].title
print "Queryable\t", wms[layer].queryable
print "Opaque\t", wms[layer].opaque
print "BBOX\t", wms[layer].boundingBox
# print "BBOX\t", wms[layer].boundingBoxWGS84 # the same
# print "CRS\t", wms[layer].crsOptions # long list
print "Styles:"
for s in wms[layer].styles: print "\t", s
print
printlayer('stations_360')
# printlayer('nemo:fbis_sampling_events')
Title Stations Catalogue - public stations (0 - 360) Queryable 1 Opaque 0 BBOX (0.0, -90.0, 360.0, 30.0, 'EPSG:4326') Styles: small_yellow_dot small_red_circle StationType
print "Operations\t", [op.name for op in wms.operations]
print "GetMap Methods\t", wms.getOperationByName('GetMap').methods
print "GetMap Formats\t", wms.getOperationByName('GetMap').formatOptions
Operations ['GetCapabilities', 'GetMap', 'GetFeatureInfo', 'DescribeLayer', 'GetLegendGraphic', 'GetStyles'] GetMap Methods [{'url': 'http://gs.niwa.co.nz:80/geoserver/stations/wms?SERVICE=WMS&', 'type': 'Get'}] GetMap Formats ['image/png', 'application/atom xml', 'application/atom+xml', 'application/openlayers', 'application/pdf', 'application/rss xml', 'application/rss+xml', 'application/vnd.google-earth.kml', 'application/vnd.google-earth.kml xml', 'application/vnd.google-earth.kml+xml', 'application/vnd.google-earth.kml+xml;mode=networklink', 'application/vnd.google-earth.kmz', 'application/vnd.google-earth.kmz xml', 'application/vnd.google-earth.kmz+xml', 'application/vnd.google-earth.kmz;mode=networklink', 'atom', 'image/geotiff', 'image/geotiff8', 'image/gif', 'image/gif;subtype=animated', 'image/jpeg', 'image/png8', 'image/png; mode=8bit', 'image/svg', 'image/svg xml', 'image/svg+xml', 'image/tiff', 'image/tiff8', 'kml', 'kmz', 'openlayers', 'rss', 'text/html; subtype=openlayers']
#Function that saves the layer as an image
def saveLayerAsImage(layer, inname):
out = open(inname, 'wb')
out.write(layer.read())
out.close()
# mark out a bounding box
llong = 160
rlong = 180
tlat = -30
blat = -50
layer = wms.getmap( layers=['stations:stations_360'],
styles=['StationType'],
srs='EPSG:4326',
bbox=(llong, blat, rlong, tlat),
size=(512,512),
format='image/png',
transparent=True
)
from IPython.core.display import Image
stationsmap = Image(layer.read())
stationsmap
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
import numpy as np
llong = 172
rlong = 179
tlat = -35
blat = -40
m = Basemap(llcrnrlon=llong, llcrnrlat=blat,
urcrnrlon=rlong, urcrnrlat=tlat,
resolution='l',epsg=4326)
# plot the image in Basemap
plt.figure(figsize=(12,10),frameon=True)
m.shadedrelief()
m.wmsimage(niwa_ows,
xpixels=600,
# ypixels=600,
verbose=False,
layers=['stations:stations_360'],
styles=['StationType'],
transparent='true'
)
plt.show()
We can zoom in - but then we need to find a better base map.
TODO use google or open street map as the basemap.
llong = 174.5
rlong = 175.5
tlat = -36.5
blat = -37.5
m = Basemap(llcrnrlon=llong, llcrnrlat=blat,
urcrnrlon=rlong, urcrnrlat=tlat,
resolution='l',epsg=4326)
# plot the image in Basemap
plt.figure(figsize=(12,10),frameon=True)
m.shadedrelief()
m.drawlsmask(resolution='f', grid=1.25, land_color='darkgreen', ocean_color="aqua")
m.drawcoastlines()
m.wmsimage(niwa_ows,
xpixels=600,
layers=['stations:stations_360'],
styles=['StationType'],
transparent='true'
)
plt.show()
plt.figure(figsize=(12,10),frameon=True)
m = Basemap(projection='cyl',resolution='l')
m.bluemarble()
m.drawmapboundary()
m.drawcoastlines()
m.wmsimage(niwa_ows,
xpixels=600,
layers=['stations:stations_180'],
styles=['StationType'],
transparent='true'
)
plt.show()
To build a filter we need to use the OGC Filter Encoding System (FES). This is a consistant method for describing filters across WMS/WFS and CSW services. Filters can be spatial, temporal and property based.
Using the NIWA Environmental Information Browser (http://ei.niwa.co.nz) we can quickly select a subset of stations. We can then obtain the filter matching that request: e.g.
<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs"
service="WFS" version="1.1.0"
xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
resultType="results">
<wfs:Query
typeName="feature:stations_180"
srsName="EPSG:4326"
xmlns:feature="http://niwa.co.nz/ns/niwa">
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
<ogc:And>
<ogc:And>
<ogc:PropertyIsLike matchCase="false" wildCard="*" singleChar="." escapeChar="!">
<ogc:PropertyName>types</ogc:PropertyName><ogc:Literal>*Climate*</ogc:Literal>
</ogc:PropertyIsLike><ogc:PropertyIsLike matchCase="false" wildCard="*" singleChar="." escapeChar="!">
<ogc:PropertyName>organisation</ogc:PropertyName>
<ogc:Literal>*National Rural Fire Authority *</ogc:Literal>
</ogc:PropertyIsLike>
</ogc:And>
<ogc:BBOX>
<ogc:PropertyName>geom</ogc:PropertyName>
<gml:Envelope xmlns:gml="http://www.opengis.net/gml">
<gml:lowerCorner>144.2123945312524 -53.366243075652555</gml:lowerCorner>
<gml:upperCorner>205.5161054687471 -26.09026816479101</gml:upperCorner>
</gml:Envelope>
</ogc:BBOX>
</ogc:And>
</ogc:Filter>
</wfs:Query>
</wfs:GetFeature>
fes = """
<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs"
service="WFS" version="1.1.0"
xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
resultType="results">
<wfs:Query
typeName="feature:stations_180"
srsName="EPSG:4326"
xmlns:feature="http://niwa.co.nz/ns/niwa">
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
<ogc:And>
<ogc:And>
<ogc:PropertyIsLike matchCase="false" wildCard="*" singleChar="." escapeChar="!">
<ogc:PropertyName>types</ogc:PropertyName><ogc:Literal>*Climate*</ogc:Literal>
</ogc:PropertyIsLike><ogc:PropertyIsLike matchCase="false" wildCard="*" singleChar="." escapeChar="!">
<ogc:PropertyName>organisation</ogc:PropertyName>
<ogc:Literal>*National Rural Fire Authority *</ogc:Literal>
</ogc:PropertyIsLike>
</ogc:And>
<ogc:BBOX>
<ogc:PropertyName>geom</ogc:PropertyName>
<gml:Envelope xmlns:gml="http://www.opengis.net/gml">
<gml:lowerCorner>144.2123945312524 -53.366243075652555</gml:lowerCorner>
<gml:upperCorner>205.5161054687471 -26.09026816479101</gml:upperCorner>
</gml:Envelope>
</ogc:BBOX>
</ogc:And>
</ogc:Filter>
</wfs:Query>
</wfs:GetFeature>
"""
# FES is the Filter Encoding Standard.
from owslib.fes import PropertyIsEqualTo, PropertyIsLike, BBox
from owslib.etree import etree
station_type = "climate"
stationTypeFes = PropertyIsLike(propertyname='types', literal=('*%s*' % station_type))
stationTypeFes.toXML()
# How to put the FES in the WMS ?
<Element {http://www.opengis.net/ogc}PropertyIsLike at 0x108d90d40>
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
llong = 172
rlong = 179
tlat = -35
blat = -40
# plot the image in Basemap
plt.figure(figsize=(12,10),frameon=True)
layer = wms.getmap( layers=['stations:stations_360'],
styles=['StationType'],
srs='EPSG:4326',
bbox=(llong, blat, rlong, tlat),
size=(512,512),
format='image/png',
transparent=True,
constraints=[stationTypeFes]
)
from IPython.core.display import Image
stationsmap = Image(layer.read())
stationsmap
<matplotlib.figure.Figure at 0x108fa8c90>
wms.items()
[('timeseries_types', <owslib.wms.ContentMetadata instance at 0x108d66b90>), ('organisations', <owslib.wms.ContentMetadata instance at 0x108d61a70>), ('stations_360', <owslib.wms.ContentMetadata instance at 0x108d66758>), ('timeseries_channels', <owslib.wms.ContentMetadata instance at 0x108d66ab8>), ('stations_publish_to_sos', <owslib.wms.ContentMetadata instance at 0x108d664d0>), ('stations_180', <owslib.wms.ContentMetadata instance at 0x108d61fc8>), ('station_types', <owslib.wms.ContentMetadata instance at 0x108d61908>)]