# 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


#Listing all available layers...
layers = list(wms.contents)
for l in layers:
    print  wms[l].name, wms[l].title


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')    

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

#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

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()

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()


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 ? 

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


wms.items()