Simple Python Interface to USGS TNM Elevation Service

Premise
Wanted a simpler way to access the USGS seamless elevation look-up service. Python seemed like a logical start. Note that the response from the USGS webservice is not correctly identified as valid XML by the python XML-parser. Therefore there is a small amount of scrubbing used to coerce the response into valid XML. Comments on why this is, or is not, a good idea are welcome.

Update It looks like the USGS service does not accept POST-style requests. I have made some small changes to the script below.

Example Listing

# nice CLI interface

import sys

from optparse import OptionParser

# url access

import urllib, urllib2

# simple XML parsing

from xml.dom import minidom

# CSV file parsing

import csv

# define service URL

base_url = 'http://gisdata.usgs.gov/XMLWebServices/TNM_Elevation_Service.asmx/getElevation?X_Value=%f&Y_Value=%f&Elevation_Units=meters&Source_Layer=-1&Elevation_Only=1'

#process command line (optparse)

parser = OptionParser()

parser.add_option("-f", "--file", dest="infile", help="input csv file containing WGS84 (lon,lat,site_id) record", metavar="FILE")

# process args

(options, args) = parser.parse_args()

#require an input file

if not options.infile:

        print "ERROR: must supply an input file!"

        sys.exit(1)

#open input file

try:

        infile = options.infile

        reader = csv.reader(open(infile, "rb"))

except:

        print "ERROR: Cannot open: " + infile

        sys.exit(1)

# read the csv file

for row in reader:

        lon = float(row[0])

        lat = float(row[1])

        site_id = row[2]

        # the base URL of the web service

        url = 'http://gisdata.usgs.gov/XMLWebServices/TNM_Elevation_Service.asmx/getElevation'

        # url GET args

        values = {'X_Value' : lon,

        'Y_Value' : lat,

        'Elevation_Units' : 'meters',

        'Source_Layer' : '-1',

        'Elevation_Only' : '1', }

        # make some fake headers, with a user-agent that will 

        # not be rejected by bone-headed servers

        user_agent = 'Mozilla/4.0 (compatible; MSIE 5.5; Windows NT)'

        headers = {'User-Agent' : user_agent}

        # encode the GET arguments

        data = urllib.urlencode(values)

        # make the URL into a qualified GET statement:

        get_url = url + '?' + data

        # make the request: note that by ommitting the url arguments

        # we force a GET request, instead of a POST

        req = urllib2.Request(url=get_url, headers=headers)

        response = urllib2.urlopen(req)

        the_page = response.read()

        # convert the HTML back into plain XML

        for entity, char in (('lt', '<'), ('gt', '>'), ('amp', '&')):

                the_page = the_page.replace('&%s;' % entity, char)

        # clean some cruft... XML won't parse with this stuff in there...

        the_page = the_page.replace('<string xmlns="http://gisdata.usgs.gov/XMLWebServices/">', '')

        the_page = the_page.replace('<?xml version="1.0" encoding="utf-8"?>\r\n', '')

        the_page = the_page.replace('</string>', '')

        the_page = the_page.replace('<!-- Elevation Values of -1.79769313486231E+308 (Negative Exponential Value) may mean the data source does not have values at that point.  --> <USGS_Elevation_Web_Service_Query>', '')

        # parse the cleaned XML

        dom = minidom.parseString(the_page)

        children = dom.getElementsByTagName('Elevation_Query')[0]

        # extract the interesting parts

        elev = float(children.getElementsByTagName('Elevation')[0].firstChild.data)

        data_source = children.getElementsByTagName('Data_Source')[0].firstChild.data

        # print to stdout

        print "%f,%f,%s,%f,%s" % (lon, lat, site_id, elev, data_source)

Incantation

python elev_lookup.py -f ca_cities.txt

Input File CA cities in WGS84 (lon,lat,id) records

-122.32,37.78,Alameda NAS
-120.53,41.48,Alturas
-124.1,40.98,Arcata
-119.05,35.43,Bakersfield
-121.45,39.13,Beale AFB
-116.95,33.93,Beaumont
-116.62,35.28,Bicycle Lk
-116.68,34.27,Big Bear Apt
-118.6,37.60,Bishop
-120.7,39.28,Blue Canyon
-114.72,33.62,Blythe
-118.37,34.20,Burbank

Output

-122.320000,37.780000,Alameda NAS,2.123548,NED Contiguous U. S. 1/3W arc second elevation data
-120.530000,41.480000,Alturas,1331.056396,NED Contiguous U. S. 1/3W arc second elevation data
-124.100000,40.980000,Arcata,63.749836,NED Contiguous U. S. 1/3W arc second elevation data
-119.050000,35.430000,Bakersfield,148.473618,NED Contiguous U. S. 1/3W arc second elevation data
-121.450000,39.130000,Beale AFB,29.950783,NED Contiguous U. S. 1/3W arc second elevation data
-116.950000,33.930000,Beaumont,792.840881,NED Contiguous U. S. 1/3W arc second elevation data
-116.620000,35.280000,Bicycle Lk,711.765869,NED Contiguous U. S. 1/3W arc second elevation data
-116.680000,34.270000,Big Bear Apt,1697.037720,NED Contiguous U. S. 1/3W arc second elevation data
-118.600000,37.600000,Bishop,2209.536865,NED Contiguous U. S. 1/3W arc second elevation data
-120.700000,39.280000,Blue Canyon,1603.887573,NED Contiguous U. S. 1/3W arc second elevation data
-114.720000,33.620000,Blythe,120.626007,NED Contiguous U. S. 1/3W arc second elevation data
-118.370000,34.200000,Burbank,224.287033,NED Contiguous U. S. 1/3W arc second elevation data
-117.350000,33.300000,Camp Pendlet,21.895588,NED Contiguous U. S. 1/3W arc second elevation data