[CentralOH] Getting the elevation of a point on the Earth in Python

[Eric Floehr]

All,

Thought I'd share this with you.  I needed elevation data as part of a
project I'm working on with climate and weather.  Specifically, I
needed the elevation of a zip code.  I have the latitude, longitude of
the zip code centroid, but not the elevation at that latitude,
longitude pair.

I thought I might use the Google Elevation API:
http://code.google.com/apis/maps/documentation/elevation/

However, the TOS states it must be used with a map display.  I could
not ensure that, so I could not use that service.  Thankfully, through
some searching, I discovered the USGS has an elevation web service
interface:

http://gisdata.usgs.net/XMLWebServices/TNM_Elevation_Service.php

Interacting with a web service with Python is super-easy.  I used the
"suds" package (https://fedorahosted.org/suds/wiki/Documentation).
Just "pip install suds".  To create a client, just pass in the
web-service description file:

>>> from suds.client import Client
>>> USGS_Client = Client('http://gisdata.usgs.net/XMLWebServices/TNM_Elevation_service.asmx?WSDL')

Python's dynamic nature makes using a web service easy.  The client
object (USGS_Client) now has methods attached to it that match what is
described in the WSDL file.  And if you print the USGS_Client object,
it's representation is the description of the methods available.

In this case, I want to use the "getElevation" method
(http://gisdata.usgs.net/XMLWebServices/TNM_Elevation_Service_Methods.php).
 According to the doc, it takes a longitude, latitude, measurement
unit, and which elevation data set to use.  A -1 means use the best
available elevation.  So I can just call:

>>> result = USGS_Client.service.getElevation(longitude, latitude, 'FEET', '-1')

The result that is returned is also dynamically generated, based on
the WSDL file.  In this case, it's a number of things:

>>> result.USGS_Elevation_Web_Service_Query.Elevation_Query
(Elevation_Query){
   _y = "40"
   _x = "-80"
   Data_Source = "NED 1/9th arc-second: Southwest Pennsylvania"
   Data_ID = "Elev_PA_Southwest"
   Elevation = "910.211698276791"
   Units = "FEET"
 }

So to get elevation, I just had to:

>>> elevation = result.USGS_Elevation_Web_Service_Query.Elevation_Query.Elevation

Piece of cake!

Best Regards,
Eric