This example demonstrates connecting to a NIWA Web Feature Service from a Python Script.
We use the EI Browser to obtain a Filter and then show information about the resulting data.
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/stations/ows"
niwa_gis = "http://gis.niwa.co.nz/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer"
%matplotlib inline
from owslib.wfs import WebFeatureService
wfs = WebFeatureService(niwa_gis)
print "Title:\t", wfs.identification.title
print "Version:\t", wfs.identification.version
print "Abstract:\t", wfs.identification.abstract
print "Keywords:"
for key in wfs.identification.keywords: print "\t", key
# list operations
print "Operations:"
for op in wfs.operations: print "\t", op.name
Title: Virtual Climate Station Network Version: 1.0.0 Abstract: Virtual Climate Station Network Keywords: Virtual Climate Station Network Operations: GetCapabilities DescribeFeatureType GetFeature
#Listing all available layers...
layers = list(wfs.contents)
for l in layers: print l
CLIMATE_VCS_Network:VCS_Forcast_Sites CLIMATE_VCS_Network:VCS_Grid CLIMATE_VCS_Network:VCS_Network
def printlayer(layer) :
print "Title\t", layer.title
# print "BBOX\t", layer.boundingBox
print "BBOX\t", layer.boundingBoxWGS84
print "Abstract\t", layer.abstract
print "Keywords\t", layer.keywords
print "CRS Options\t", layer.crsOptions
print "Verb Options\t", layer.verbOptions
print
printlayer(wfs.contents['CLIMATE_VCS_Network:VCS_Forcast_Sites'])
printlayer(wfs.contents['CLIMATE_VCS_Network:VCS_Grid'])
printlayer(wfs.contents['CLIMATE_VCS_Network:VCS_Network'])
# printlayer(wfs.contents['nemo:fbis_sampling_events'])
Title VCS_Forcast_Sites BBOX (166.47504312700005, -47.273449178999954, 178.4752155010001, -34.47317700399998) Abstract None Keywords [] CRS Options [urn:ogc:def:crs:EPSG::4326] Verb Options ['{http://www.opengis.net/wfs}Query'] Title VCS_Grid BBOX (166.45003862400006, -47.29845010399998, 178.50021570400008, -34.39817650599997) Abstract None Keywords [] CRS Options [urn:ogc:def:crs:EPSG::4326] Verb Options ['{http://www.opengis.net/wfs}Query'] Title VCS_Network BBOX (166.47503985200012, -47.273449178999954, 178.4752155010001, -34.42317680399998) Abstract None Keywords [] CRS Options [urn:ogc:def:crs:EPSG::4326] Verb Options ['{http://www.opengis.net/wfs}Query']
print "Operations\t", [op.name for op in wfs.operations]
print "GetFeature Methods\t", wfs.getOperationByName('GetFeature').methods
print "GetFeature Formats\t", wfs.getOperationByName('GetFeature').formatOptions
Operations ['GetCapabilities', 'DescribeFeatureType', 'GetFeature'] GetFeature Methods [{'url': 'http://gis.niwa.co.nz/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer?', 'type': 'Get'}, {'url': 'http://gis.niwa.co.nz/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer', 'type': 'Post'}] GetFeature Formats ['{http://www.opengis.net/wfs}GML2']
Note that the ArcGIS server WFS only provides results in GML2. Not JSON
vcsn_forecast_sites = 'CLIMATE_VCS_Network:VCS_Forcast_Sites'
vcsn_grid = 'CLIMATE_VCS_Network:VCS_Grid'
vcsn_network = 'CLIMATE_VCS_Network:VCS_Network'
# mark out a bounding box - Whole of NZ.
llong = 160
rlong = 180
tlat = -30
blat = -50
get_stations = wfs.getfeature(
typename=[vcsn_forecast_sites],
bbox=(llong, blat, rlong, tlat),
maxfeatures = 5,
# outputFormat="application/json"
)
# the result returned is a stream on the request so we have to read or getValue it to get the actual data.
s = get_stations.read()
print s
<wfs:FeatureCollection xsi:schemaLocation='http://Beaker2.niwa.local:6080/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer http://gis.niwa.co.nz/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer?request=DescribeFeatureType%26version=1.0.0%26typename=VCS_Forcast_Sites http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd' xmlns:CLIMATE_VCS_Network='http://Beaker2.niwa.local:6080/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer' xmlns:gml='http://www.opengis.net/gml' xmlns:wfs='http://www.opengis.net/wfs' xmlns:xlink='http://www.w3.org/1999/xlink' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance'><gml:boundedBy><gml:Box srsName='EPSG:4326'><gml:coordinates>166.47504312700005,-47.273449178999954 178.47521550100009,-34.473177003999979</gml:coordinates></gml:Box></gml:boundedBy><gml:featureMember><CLIMATE_VCS_Network:VCS_Forcast_Sites fid='F422__710'><CLIMATE_VCS_Network:OBJECTID>710</CLIMATE_VCS_Network:OBJECTID><CLIMATE_VCS_Network:Shape><gml:Point srsName="EPSG:4326"><gml:coordinates>168.67507512900011,-45.473382317999949</gml:coordinates></gml:Point></CLIMATE_VCS_Network:Shape></CLIMATE_VCS_Network:VCS_Forcast_Sites></gml:featureMember><gml:featureMember><CLIMATE_VCS_Network:VCS_Forcast_Sites fid='F422__547'><CLIMATE_VCS_Network:OBJECTID>547</CLIMATE_VCS_Network:OBJECTID><CLIMATE_VCS_Network:Shape><gml:Point srsName="EPSG:4326"><gml:coordinates>170.07509738200008,-45.873384668999961</gml:coordinates></gml:Point></CLIMATE_VCS_Network:Shape></CLIMATE_VCS_Network:VCS_Forcast_Sites></gml:featureMember><gml:featureMember><CLIMATE_VCS_Network:VCS_Forcast_Sites fid='F422__1234'><CLIMATE_VCS_Network:OBJECTID>1234</CLIMATE_VCS_Network:OBJECTID><CLIMATE_VCS_Network:Shape><gml:Point srsName="EPSG:4326"><gml:coordinates>170.97509801900003,-43.573320896999974</gml:coordinates></gml:Point></CLIMATE_VCS_Network:Shape></CLIMATE_VCS_Network:VCS_Forcast_Sites></gml:featureMember><gml:featureMember><CLIMATE_VCS_Network:VCS_Forcast_Sites fid='F422__570'><CLIMATE_VCS_Network:OBJECTID>570</CLIMATE_VCS_Network:OBJECTID><CLIMATE_VCS_Network:Shape><gml:Point srsName="EPSG:4326"><gml:coordinates>169.97508378900011,-43.973337152999932</gml:coordinates></gml:Point></CLIMATE_VCS_Network:Shape></CLIMATE_VCS_Network:VCS_Forcast_Sites></gml:featureMember><gml:featureMember><CLIMATE_VCS_Network:VCS_Forcast_Sites fid='F422__902'><CLIMATE_VCS_Network:OBJECTID>902</CLIMATE_VCS_Network:OBJECTID><CLIMATE_VCS_Network:Shape><gml:Point srsName="EPSG:4326"><gml:coordinates>171.07509356000003,-44.673354605999975</gml:coordinates></gml:Point></CLIMATE_VCS_Network:Shape></CLIMATE_VCS_Network:VCS_Forcast_Sites></gml:featureMember></wfs:FeatureCollection>
This GML has fairly messy namespace definitions so we want to build the namespace map from the header element. ElementTree doesn't do this that well so we will use lxml
from lxml import etree
root = etree.fromstring(s)
nsmap = root.nsmap
print nsmap
#etree.tostring(root, method='html', pretty_print=True)
{'wfs': 'http://www.opengis.net/wfs', 'xsi': 'http://www.w3.org/2001/XMLSchema-instance', 'CLIMATE_VCS_Network': 'http://Beaker2.niwa.local:6080/arcgis/services/CLIMATE/VCS_Network/MapServer/WFSServer', 'xlink': 'http://www.w3.org/1999/xlink', 'gml': 'http://www.opengis.net/gml'}
Now we can find the features using tidy prefixes
features = root.findall('gml:featureMember/CLIMATE_VCS_Network:VCS_Forcast_Sites', root.nsmap)
print "Result %s features." % (len(features))
def print_feature(f):
print "ObjectID: ", f.find("CLIMATE_VCS_Network:OBJECTID", root.nsmap).text,
coords = f.find("CLIMATE_VCS_Network:Shape/gml:Point/gml:coordinates", root.nsmap).text
lon, lat = coords.split(',')
print "\tat:\t", lon, lat
for f in features:
print_feature(f)
Result 5 features. ObjectID: 710 at: 168.67507512900011 -45.473382317999949 ObjectID: 547 at: 170.07509738200008 -45.873384668999961 ObjectID: 1234 at: 170.97509801900003 -43.573320896999974 ObjectID: 570 at: 169.97508378900011 -43.973337152999932 ObjectID: 902 at: 171.07509356000003 -44.673354605999975
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
#Get all the point data
# mark out a bounding box - Whole of NZ.
llong = 165
rlong = 179
tlat = -33
blat = -48
get_stations = wfs.getfeature(
typename=[vcsn_forecast_sites],
bbox=(llong, blat, rlong, tlat),
maxfeatures = 12000
# outputFormat="application/json"
)
# the result returned is a stream on the request so we have to read or getValue it to get the actual data.
vcsn_stations = get_stations.read()
root = etree.fromstring(s)
features = root.findall('gml:featureMember/CLIMATE_VCS_Network:VCS_Forcast_Sites', root.nsmap)
print "Result %s features." % (len(features))
Result 2882 features.
# using list comprehensions separate out the lats and longs
def getcoords(f):
coords = f.find("CLIMATE_VCS_Network:Shape/gml:Point/gml:coordinates", root.nsmap).text
return coords.split(',')
points = [getcoords(f) for f in features]
lons = [float(p[0]) for p in points]
lats = [float(p[1]) for p in points]
# Note that we're drawing on a regular matplotlib figure, so we set the
# figure size just like we would any other.
plt.figure(figsize=(22,22))
# Create the Basemap
m = Basemap(llcrnrlon=llong, llcrnrlat=blat,
urcrnrlon=rlong, urcrnrlat=tlat,
resolution='h',
epsg=3994) # Mercator 41.
# epsg=4326) # WGS84.
# Draw important features
m.drawcoastlines(color='0.4')
m.drawcountries()
# m.fillcontinents(color='0.8') # Light gray
x, y = m(lons, lats) # Convert lat, long to y,x
m.plot(x,y, 'bo', markersize=4, alpha=0.7)
[<matplotlib.lines.Line2D at 0x10f51b410>]