General Purpose Programming with Scripting Languages

Some initial hints of use of various scripting languages.

Python

Ruby

Interactive Ruby tutorial

Python Image Module Example: How much ink on that page?

Premise:
Thought it would be fun to compute how much ink a given poster requires, per unit area of paper, when sending to the department large-format printer. The Python Imaging Library provides several modules suitable for low-level operation on image data. A simple (and probably very inefficient) script was developed to compute the white/black percentage of an image. A script like this could be used to adjust a per-poster "ink cost", which would hopefully prevent people from wasting ink. Obviously, this computation is scale-dependent, so standardized rasterization parameters would have to be set in order for the "ink cost" calculation to be fair. More generalized or efficient approaches are always welcomed.

Implementation: (when copying/pasting note whitespace in blocks)

#!/usr/bin/env python

# invoke like this: ./white.py image.png

# idea from http://www.halfbakery.com/user/lurch

import sys

import Image

# cheap and non-robust way to access first argument

file = sys.argv[1]

# open the image, load into memory

pic = Image.open( file )

imgdata = pic.load()

# extract the dimensions of the image

xsize, ysize = pic.size

# init some counters

black = 0

white = 0

# define some colors (R,G,B)

b_color = (0,0,0)

w_color = (255,255,255)

# loop over columns

for x in xrange(xsize):

        # loop over rows

        for y in xrange(ysize):

                if imgdata[x,y] == w_color:

                        white+=1

                if imgdata[x,y] == b_color:

                        black+=1

# compute percentages

total_pixels = xsize * ysize

pct_white = (float(white) / float(total_pixels)) * 100

pct_black = (float(black) / float(total_pixels)) * 100

# formatted printing

print "%.3f%% of image is white!" % (pct_white)

print "%.3f%% of image is black!" % (pct_black)

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. Here is a link to an interactive version.

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

Text processing with the UNIX toolbox

Working with DBF files

Extract the text contents of a DBF file with dbfxtrct, select out the 3rd and 2nd columns with awk, keep only unique rows with sort and uniq, and finally saving to a new file called 'ca_fips-name_list.txt' with shell redirection.

# dbfxtrct utility
dbfxtrct -ica_counties_aea83.dbf | awk -F"," '{print $3","$2}' | sort -g | uniq > ca_fips-name_list.txt
Other DBF tools include:
- dbf2mysql
- dbf2pg
- dbfadd
- dbfcreate
- dbfdump
- dbfutil1
- dbfxtrct