Skip to content

Instantly share code, notes, and snippets.

@rbejar

rbejar/IDEE_Web_Services_Example.py

Created September 29, 2011 15:27
Show Gist options
  • Save rbejar/1250982 to your computer and use it in GitHub Desktop.
Save rbejar/1250982 to your computer and use it in GitHub Desktop.
Download ZIP
A 'quick and dirty' experiment to show how to create a service chain based on Spain SDI (http://www.idee.es) OGC web services.
Raw
IDEE_Web_Services_Example.py
"""
A 'quick and dirty' experiment to show how to create a service chain
based on Spain SDI (http://www.idee.es) OGC web services.
Most of the code is boilerplate code: an OGC WPS client library
could be useful.
Many things should be parameterized etc. The WPS used have more
operations, and the used operations have more options: this code
is just an example.
The services should be normally up and running. If that is not
the case, you can go to http://www.idee.es, find the
contact e-mail and ask there if they can solve the problem.
This code has been tested with Python 2.6 in a Linux environment.
It requires Shapely (https://github.com/sgillies/shapely) and
ElementTree (http://effbot.org/zone/element-index.htm).
Tested with ElementTree 1.2.7 and Shapely 1.2.1.
Permission is hereby granted, free of charge, to any
person obtaining a copy of this software (the "Software"),
to deal in the Software without restriction, including without
limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the
Software, and to permit persons to whom the Software is
furnished to do so, subject to the following condition:
This permission notice shall be included in all copies or
substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
import urllib2
import shapely.geometry
import elementtree.ElementTree as ET
_domain = 'http://www.idee.es'
_url = _domain+'/WPS/services'
def _request_Post_XML(url, xml_request):
""" Posts xml_request to url and returns its results.
Url must be willing to accept requests in XML.
"""
req = urllib2.Request(url, xml_request)
req.add_header('Content-Type', 'text/xml')
response = urllib2.urlopen(req)
results = response.read()
return results
def raster_statistics_WPS(boundary):
""" Boundary must be a shapely.geometry.Polygon.
This function wraps a WPS service in the IDEE. It returns raster information
about a given boundary in XML. The boundary must be in the Iberian Peninsula
and described in UTM 30N, ED50 (EPSG:23030). It should not be a very large area.
"""
xmin, ymin, xmax, ymax = boundary.bounds
# To format boundary in the GRASS ASCII format the WPS expects (see http://grass.fbk.eu/grass64/manuals/html64_user/v.in.ascii.html)
# this steps are taken:
# list of tuples (x1, y1),(x2, y2)... where x and y are numbers
boundary_list = list(boundary.exterior.coords)
# list of strings 'x y','x2 y2'... where x and y were numbers
boundary_str_list = [str(x)+' '+str(y) for x,y in boundary_list]
# string of 'x1 y1\nx2 y2 ...'
boundary_str = reduce(lambda x, y: x+'\n '+y, boundary_str_list)
wcs_request = _domain + '/wcs/IDEE-WCS-UTM30N/wcsServlet?SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.0.0&CRS=EPSG:23030&BBOX={xmin},{ymin},{xmax},{ymax}&COVERAGE=MDT25_peninsula_ZIP&RESX=10&RESY=10&FORMAT=FloatGridExtended_Zip&EXCEPTIONS=XML'.format(xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax)
# The WPS is very picky with the GRASS ASCII parameter, so any change in spaces, tabs etc. can result
# in difficult to understand service exceptions.
# The xml_request is formatted to include the wcs_request and the adapted boundary
xml_request="""<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<wps:Execute service="WPS" version="0.4.0" store="false" status="false" xmlns:wps="http://www.opengeospatial.net/wps" xmlns:ows="http://www.opengis.net/ows" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengeospatial.net/wps..\wpsExecute.xsd">
<ows:Identifier>RasterStatistics</ows:Identifier>
<wps:DataInputs>
<wps:Input><ows:Identifier>ASCIIPolygon</ows:Identifier><ows:Title>ASCIIPolygon</ows:Title><ows:Abstract>ASCIIPolygon</ows:Abstract>
<wps:ComplexValue format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">ORGANIZATION:
DIGIT DATE:
DIGIT NAME: user_name
MAP NAME:
MAP DATE: Mon Feb 18 10:30:48 2008
MAP SCALE: 1
OTHER INFO:
ZONE: 0
MAP THRESH: 0.000000
VERTI:
B {boundary_len}
{boundary_str}</wps:ComplexValue>
</wps:Input>
<wps:Input><ows:Identifier>URLCoverageServer</ows:Identifier><ows:Title>URLCoverageServer</ows:Title><ows:Abstract>URLCoverageServer</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">{wcs_request}</wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>URLFeatureServer</ows:Identifier><ows:Title>URLFeatureServer</ows:Title><ows:Abstract>URLFeatureServer</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String"></wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>FeatureRequest</ows:Identifier><ows:Title>FeatureRequest</ows:Title><ows:Abstract>FeatureRequest</ows:Abstract>
<wps:ComplexValue format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"></wps:ComplexValue>
</wps:Input>
<wps:Input><ows:Identifier>CRS</ows:Identifier><ows:Title>CRS</ows:Title><ows:Abstract>CRS</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">EPSG:23030</wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>CRSResponse</ows:Identifier><ows:Title>CRSResponse</ows:Title><ows:Abstract>CRSResponse</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">EPSG:23030</wps:LiteralValue>
</wps:Input>
</wps:DataInputs>
<wps:OutputDefinitions>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/basicTypes.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>CoverageResponse</ows:Identifier><ows:Title>CoverageResponse</ows:Title><ows:Abstract>CoverageResponse</ows:Abstract>
</wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>Maximo</ows:Identifier><ows:Title>Maximo</ows:Title><ows:Abstract>Maximo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>xMaximo</ows:Identifier><ows:Title>xMaximo</ows:Title><ows:Abstract>xMaximo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>yMaximo</ows:Identifier><ows:Title>yMaximo</ows:Title><ows:Abstract>yMaximo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>Minimo</ows:Identifier><ows:Title>Minimo</ows:Title><ows:Abstract>Minimo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>xMinimo</ows:Identifier><ows:Title>xMinimo</ows:Title><ows:Abstract>xMinimo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>yMinimo</ows:Identifier><ows:Title>yMinimo</ows:Title><ows:Abstract>yMinimo</ows:Abstract></wps:Output>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Double">
<ows:Identifier>Media</ows:Identifier><ows:Title>Media</ows:Title><ows:Abstract>Media</ows:Abstract>
</wps:Output></wps:OutputDefinitions></wps:Execute>""".format(wcs_request = wcs_request, boundary_len=len(boundary_list), boundary_str=boundary_str)
return _request_Post_XML(_url, xml_request)
def relief_profile_WPS(point1, point2, format='Png'):
"""Point1 and point2 are shapely.geometry.Points.
Format must be 'Txt' or 'Png' (default), for the answer to be a text file with coordinates and heights,
or an image in Png format.
This function wraps a WPS service in the IDEE. It returns a relief profile
between two points. These points must be in the Iberian Peninsula
and in UTM 30N, ED50 (EPSG:23030). They should not be very far.
"""
# Calculate bounds for the WCS request
tolerance = 500
xmin, ymin, xmax, ymax = shapely.geometry.MultiPoint([point1, point2]).bounds
xmin = xmin - tolerance
ymin = ymin - tolerance
xmax = xmax + tolerance
ymax = ymax + tolerance
wcs_request = _domain + '/wcs/IDEE-WCS-UTM30N/wcsServlet?SERVICE=WCS&amp;REQUEST=GetCoverage&amp;VERSION=1.0.0&amp;CRS=EPSG:23030&amp;BBOX={xmin},{ymin},{xmax},{ymax}&amp;COVERAGE=MDT25_peninsula_ZIP&amp;RESX=10&amp;RESY=10&amp;FORMAT=AsciiGrid&amp;EXCEPTIONS=XML'.format(xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax)
xml_request="""<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<wps:Execute service="WPS" version="0.4.0" store="false" status="false" xmlns:wps="http://www.opengeospatial.net/wps"
xmlns:ows="http://www.opengis.net/ows" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://www.opengeospatial.net/wps..\wpsExecute.xsd"> <ows:Identifier>GetProfile{format}</ows:Identifier>
<wps:DataInputs>
<wps:Input><ows:Identifier>Resolution</ows:Identifier><ows:Title>Resolution</ows:Title><ows:Abstract>Resolution</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:Double" uom="urn:ogc:def:dataType:OGC:0.0:Double">10</wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>Profile</ows:Identifier><ows:Title>Profile</ows:Title><ows:Abstract>Profile</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">{x1},{y1},{x2},{y2}</wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>URLCoverageServer</ows:Identifier><ows:Title>URLCoverageServer</ows:Title>
<ows:Abstract>URLCoverageServer</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">{wcs_request}</wps:LiteralValue>
</wps:Input>
<wps:Input><ows:Identifier>CRS</ows:Identifier><ows:Title>CRS</ows:Title><ows:Abstract>CRS</ows:Abstract>
<wps:LiteralValue dataType="urn:ogc:def:dataType:OGC:0.0:String" uom="urn:ogc:def:dataType:OGC:0.0:String">EPSG:23030</wps:LiteralValue> </wps:Input>
</wps:DataInputs><wps:OutputDefinitions>
<wps:Output format="text/xml" encoding="UTF-8" schema="http://schemas.opengis.net/gml/3.0.0/base/gml.xsd" uom="urn:ogc:def:dataType:OGC:0.0:Integer">
<ows:Identifier>Profile</ows:Identifier> <ows:Title>Profile</ows:Title> <ows:Abstract>Profile</ows:Abstract>
</wps:Output> </wps:OutputDefinitions></wps:Execute>""".format(format = format, x1 = point1.x, y1 = point1.y, x2 = point2.x, y2 = point2.y, wcs_request = wcs_request)
return _request_Post_XML(_url, xml_request)
def POI_elevation_info(point, radius):
"""Point is a shapely.geometry.Point.
Radius is in meters.
A simple example of OGC service chain that shows how WPS and WCS services from
the Spain SDI (http://www.idee.es) can be chained together to create higher
level, reusable geospatial functions.
Given the coordinates of a Point of Interest in Spain (for now it must be located
in the Iberian Peninsula and be in UTM 30N, ED50 (EPSG:23030)) and a radius in
meters, this function returns the highest and lowest points in that radius, the
average height and an image with a relief profile between the points of maximum and minimum
height.
"""
point_buffer = point.buffer(radius)
# Obtain elevation info from a WPS
elevation = raster_statistics_WPS(point_buffer)
# parse the XML answer of the WPS
element = ET.fromstring(elevation)
process_outputs = element.find('{http://www.opengeospatial.net/wps}ProcessOutputs')
outputs = process_outputs.findall('{http://www.opengeospatial.net/wps}Output')
for output in outputs:
identifier = output.find('{http://www.opengis.net/ows}Identifier')
# Not every output has a LiteralValue element, so we find it only for the
# outputs which identifier we need
if identifier.text == 'Maximo':
max_height = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'xMaximo':
xMax = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'yMaximo':
yMax = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'Minimo':
min_height = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'xMinimo':
xMin = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'yMinimo':
yMin = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
elif identifier.text == 'Media':
avg_height = float(output.find('{http://www.opengeospatial.net/wps}LiteralValue').text)
# Obtain relief profile from a WPS
profile = relief_profile_WPS(shapely.geometry.Point(xMax, yMax), shapely.geometry.Point(xMin, yMin), format='Png')
# parse the XML answer of the WPS
element = ET.fromstring(profile)
profileURL = element.find('.//URL').text
return (max_height, xMax, yMax, min_height, xMin, yMin, avg_height, profileURL)
# Mulhacen mountain coordinates in UTM 30N, ED50: 472400E 4101046N
mulhacen = shapely.geometry.Point(472400, 4101046)
distance = 500
#max_height, xMax, yMax, min_height, xMin, yMin, avg_height, profileURL = POI_elevation_info(mulhacen, distance)
elevation_info = POI_elevation_info(mulhacen, distance)
print("In a distance of {0} meters around {1},{2}:".format(distance,
mulhacen.x,mulhacen.y))
print("Average height is {0}.".format(elevation_info[6]))
print("Maximum height is {0} in {1},{2}.".format(elevation_info[0],
elevation_info[1],
elevation_info[2]))
print("Minimum height is {0} in {1},{2}.".format(elevation_info[3],
elevation_info[4],
elevation_info[5]))
print("A relief profile between points of minimum and maximum heights is in {0}.".format(elevation_info[7]))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment