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/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
NIWA Open Web Services 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
stations:stations_360 Stations Catalogue - public stations (0 - 360) nemo:fchem_chemistry_metadata fchem_chemistry_metadata rec:rivers REC2 River network lakespi:lakespi LakeSPI Survey Results niwa:urban_area urban_area nemo:fbis_sample_details Freshwater Biodiversity - Sample Details nemo:mbio_sample_details Marine Biodiversity - Sample Details lakespi:lakespi_lake_types LakeSPI Geomorphological lake classifications gn:ne_50m_boundary_da Admin 0 – Breakaway, disputed areas stations:organisations Stations Catalogue - Organisations nemo:mbio_sampling_events Marine Biodiversity - Sampling Events niwa:regional_council_bbox regional_council_bbox rec:catchments NZ WONI Catchments with names stations:stations_180 Stations Catalogue - stations lakespi:lakespi_latest LakeSPI Summary of most recent surveys stations:stations_publish_to_sos Stations Catalogue - Public SOS Stations stations:timeseries_channels timeseries_channels niwa:urban_area_gen urban_area_gen nemo:mbio_datasources Marine Biodiversity - Data Sources niwa:regional_council_gen regional_council_gen stations:timeseries_types Stations Catalogue - Time Series Types lakespi:lakespi_lakes LakeSPI Lake Names lakespi:lakespi_regions LakeSPI Region Names niwa:territorial_authority territorial_authority stations:station_types Stations Catalogue - Station Types niwa:regional_council regional_council nemo:fbis_sampling_events Freshwater Biodiversity - Sampling Events nemo:fbis_datasources Freshwater Biodiversity - Data Sources nemo:fchem_chemistry fchem_chemistry gn:world Blue Marble world image gn:ne_50m_boundary_lines_land Admin 0 – Boundary Lines niwa:nz-topo-50-map-sheet nz-topo-50-map-sheet niwa:territorial_authority_gen territorial_authority_gen nemo:fbis_fish_assistant_report Freshwater Biodiversity - Fish Assistant Report gn:ne_50m_coastline Coastline nemo:fbis_subjects Freshwater Biodiversity - Subjects nemo:fchem_locations fchem_locations
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: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 Title Freshwater Biodiversity - Sampling Events Queryable 1 Opaque 0 BBOX (-180.0, -90.0, 180.0, 90.0, 'EPSG:4326') Styles: small_yellow_dot medium_blue_dot small_red_circle
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/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'
)
parallels = np.arange(0.,81,10.)
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians,labels=[True,False,False,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()
plt.title('Orthographic projection of NIWA Stations' )
m.wmsimage(niwa_ows,
xpixels=600,
layers=['stations:stations_180'],
styles=['StationType'],
transparent='true'
)
plt.show()