How do you get a GIS feature from a Python collection/layer/provider/thingy? Let's look at 4 different popular GIS scripting environments.

ESRI (ArcGIS 9.3):

row = rows.next()
while row:
    # work with row
    ...
    row = rows.next()

FME:

self.my_collection_count = len(self.my_feature_collection)
for i in range(self.my_collection_count):
    self.my_feature_part = self.my_feature_collection[i]

    # work with self.my_feature_part
    ...

Note: don't use instance attributes like that. Use local variables.

QGIS:

feat = QgsFeature()
while provider.getNextFeature(feat):
    # work with feat
    ...

The oddball.

OGR (earlier than 1.5):

feature = layer.GetNextFeature()
while feature:
    # work with feature
    ...
    feature.Destroy()
    feature = layer.GetNextFeature()

Destroy!

4 different environments, 4 different ways, and none of them the natural Python way. There's one obviously right way to do it for Python, and that's the way that it's done in ogr.py versions >= 1.5, and how it's done in WorldMill. GeoDjango, too.

for f in layer:
    # work with f
    ...

where layer is among other things a generator that provides the iterator protocol just like Python strings, lists, and files do. It has a next method that yields a value, or raises a StopIteration exception when there are no more values. The advantages:

• Clarity: it's agonizingly clear. More clear to a non-programmer, in my opinion, than the other alternatives. For each feature in the set: do something. And then forget about the feature and move on to the next.

• Less error-prone: even a non-programmer can't screw up that one line of code any worse than to get a standard, understandable Python NameError, TypeError, or SyntaxError.

• Standardization: core Python modules such as itertools and many other useful add-on packages reward you for using the iterator protocol.

for f in feature_collection:
  coords = f.getCoordinates()
  for m in coords:
    x = m[0]
    y = m[1]

for feature in layer:
    # COLUMN_NAME must be normalized to all caps due to some
    # data sources not caring.
    print feature.COLUMN_NAME

I'd like to see GIS students taught to program in Python using Python idioms, not Avenue idioms. I hate to pick on Utah State's GIS Programming with Python just because of its popularity, but it contains some good introductory code that can be easily tuned up to teach even better Python GIS programming skills. For example, let's look at Lesson 5: Writing/Reading Text Files, Python Lists & Dictionaries, part 5a:

# Author:  John Lowry
# Date: Dec. 21, 2007
# Purpose: Lesson 5a: Use split and Write to a text file
#############################################################

#Import modules
import string

# Open a new text file to write to
outFile = open(r"C:\john\TeachingGIS\WILD6900_ArcGISPython\Lesson5_results\write_example.txt", "w")

# Make a string variable of featureclass names
fcString = "nests1990.shp,nests1995.shp,nests2000.shp"

# Make a list variable from the string variable using the split method, then print
fcList = fcString.split(",")

# Write each item in the list to a separate line the output file
outFile.write(fcList[0]+ "\n")
outFile.write(fcList[1]+ "\n")
outFile.write(fcList[2]+ "\n")

# Close the files
outFile.close()

There are 3 defects in that code:

• It imports the string module but never uses anything from it. You should almost always use string object methods anyway.

• It presumes knowledge of the number of items in the comma-separated string, specifically that it is at least 3.

• It needlessly concatenates newlines to items before writing.

The equivalent code, with none of those defects, looks like this:

outFile = open('/tmp/out.txt', 'w')
for item in fcString.split(','):
    outFile.write(item)
    outFile.write('\n')
outFile.close()

In Python 2.6 you can use a with statement to make it even more compact. The file closes itself at the end of the block. And you can use file.writelines to write the data and newline in one go. It's more efficient to pass it a tuple than to pass it a list.

with open('/tmp/out.txt', 'w') as f:
    for item in fcString.split(','):
        f.writelines((item, '\n'))

In the next part of the lesson, 5b, we see:

# Author:  John Lowry
# Date: Dec. 21, 2007
# Purpose: Lesson 5:  Reading a and writing a textfile
#############################################################

# Open the text file in read mode
inFile = open(r"C:\john\TeachingGIS\WILD6900_ArcGISPython\Lesson5\nests2005_coords.csv", "r")

# Open a new text file to write to
outFile = open(r"C:\john\TeachingGIS\WILD6900_ArcGISPython\Lesson5_results\nests2005_format.txt", "w")

# Read entire file and print one line at a time
for line in inFile.readlines():
    nestList = line.split(",")
    id = nestList[0]
    cnd = nestList[1]
    x = nestList[2]
    y = nestList[3]
    outFile.write("Siteid: " + id + "\n")
    outFile.write("Condition: " + cnd + "\n")
    outFile.write("X Coordinate: " + x + "\n")
    outFile.write("Y Coordinate: " + y + "\n")
    outFile.write("\n")

# Close the files
inFile.close()
outFile.close()

The more effective version of the looping block is this:

for line in inFile:
    outFile.write(
      'Siteid: %s\nCondition: %s\nX Coordinate: %s\nY Coordinate: %s\n\n' \
      % tuple(line.split(','))
      )

String formatting is more efficient than string concatenation (or not -- see the update below) and you can avoid needless variable assignments by using the split results directly.

I blogged before about how smelly the ArcGIS scripting cursor syntax was. I hear it's better now, but you can still see the old style in the USU course code.

Update (2009-10-21): Here's my benchmark script. I'm isolating just the inner part of the loop and focusing just on the extra assignments and file writes.

import timeit

# Sample input line
line = '1,good,433207.8362,4518107.044'

# A file-like object
class MockFile(object):
    def write(self, line):
        pass

outFile = MockFile()

# GIS Style programming. Assignment to intermediate variables
# and each written separately.

s1 = """\
nestList = line.split(',')
id = nestList[0]
cnd = nestList[1]
x = nestList[2]
y = nestList[3]
outFile.write(id)
outFile.write(cnd)
outFile.write(x)
outFile.write(y)
"""

t1 = timeit.Timer(
    stmt=s1,
    setup='from __main__ import line, outFile'
    )
print "GIS style"
print "%.2f usec/pass" % t1.timeit()
print

# Idiomatic Python. No intermediate variables and all written as
# a group

s2 = """\
outFile.write(''.join(line.split(',')))
"""

t2 = timeit.Timer(
    stmt=s2,
    setup='from __main__ import line, outFile'
    )
print "Idiomatic Python"
print "%.2f usec/pass" % t2.timeit()
print

The results:

$ python benchmarks.py
GIS style
2.07 usec/pass

Idiomatic Python
1.29 usec/pass

Someone else has looked at string performance more closely than I, and it looks like I'm wrong. On my Python 2.6, too, concatenation wins over formatting. Use of join is faster than concatenation for lists longer than 1000 items or so.

template = '''\
Siteid: %s
Condition: %s
X Coordinate: %s
Y Coordinate: %s
'''

for line in inFile:
    outFile.write(template % tuple(line.split(',')))

I've learned that in German, there's a name for a disastrous string of characters, such as "REST endpoint", that insults human dignity: unwort. Yes, "REST endpoint" is two words in English, but it would be one in German. There's even an "Unword of the Year" website devoted to these entities (in German only). The level of negative meaning in "REST API" is almost as high.

Via @xrotwang.

Check this out if you're in New York City this winter. It looks fascinating:

The Lost World of Old Europe brings to the United States for the first time more than 160 objects recovered by archaeologists from the graves, towns, and villages of Old Europe, a cycle of related cultures that achieved a precocious peak of sophistication and creativity in what is now southeastern Europe between 5000 and 4000 BC, and then mysteriously collapsed by 3500 BC. Long before Egypt or Mesopotamia rose to an equivalent level of achievement, Old Europe was among the most sophisticated places that humans inhabited. Some of its towns grew to city-like sizes. Potters developed striking designs, and the ubiquitous goddess figurines found in houses and shrines have triggered intense debates about womens roles in Old European society. Old European copper-smiths were, in their day, the most advanced metal artisans in the world. Their intense interest in acquiring copper, gold, Aegean shells, and other rare valuables created networks of negotiation that reached surprisingly far, permitting some of their chiefs to be buried with pounds of gold and copper in funerals without parallel in the Near East or Egypt at the time. The exhibition, arranged through loan agreements with 20 museums in three countries (Romania, The Republic of Bulgaria and the Republic of Moldova), brings the exuberant art, enigmatic goddess cults, and precocious metal ornaments and weapons of Old Europe to American audiences.

This fall I'm working in office space graciously provided by the Centre de Biologie et de Gestion des Populations at the CIRAD International Campus de Baillarguet outside Montferrier-sur-Lez [OpenStreetMap]. My wife is a visiting scientist here. I'm baggage, but extremely well-treated baggage. I share a corner office with a population geneticist who is waiting for his space in the new wing to be completed. There are bioinformatics programmers diving into Python next door. The internet mostly (c'est le Sud) works, and everybody is friendly and tolerant of my mediocre (terrible, really) French.

That's an image capture from 2009-10-16. Yes, I've heard about the OSM static map API, just waiting till it moves to what looks like a more permanent URL.

CBGP is out in the countryside north of the Montpellier metro. I've alluded to how pleasant it is out here and recently took a camera on my run to try to capture the feeling. Here's a view out of my office window of the parking lot and the Pic St-Loup hill to the north. One of my favorite local wineries is tucked under the cliffs on the other side of it.

A semi-wild, wooded ridge runs for some kilometers to the east of CBGP. Here's the steep trail (which I've edited into OSM) up the hill from the back door of the Centre. I'm sorry that I can't identify the trees (update: there's a least two p. pinea, or stone pine, in the photo). There's an interesting mix, including what appear to be firs near the top. A little elevation makes a big difference around here.

At the top of the ridge runs a national bike route. The markers look like this:

Single-track and fire road on the crest. Sometimes I forget I'm in France and feel like I'm back on the Front Range.

You may have noticed the group of identical conifers on the right side of the last photo. There was a 50 hectare fire here in 1981. André Jeanjean, A young volunteer sapeur-pompier, died fighting the blaze. I found fresh flowers from the Mairie of Clapiers at his memorial.

The left fork there runs down to a vineyard and on into Clapiers. I desperately need to acquire a mountain bike because CBGP is a great trailhead for some nice long rides.

Not the Google MapReduce, that is, but Python's good old map() and reduce(). Waiting on something else this afternoon and just want to write this modestly neat code down for myself.

>>> from shapely.geometry import Point
>>> def union(a, b):
...     return a.union(b)
...
>>> print reduce(union, map(Point, zip(range(3), range(3)))).wkt
MULTIPOINT (0.0000000000000000 0.0000000000000000, 1.0000000000000000 1.0000000000000000, 2.0000000000000000 2.0000000000000000)

Update (2009-10-16): Could be made even shorter using lambda.

>>> from shapely.geometry import Point
>>> print reduce(lambda a, b: a.union(b), map(Point, zip(range(3), range(3)))).__geo_interface__
{'type': 'MultiPoint', 'coordinates': ((0.0, 0.0), (1.0, 1.0), (2.0, 2.0))}

Reading that Chris Garrard was open to suggestions for improving her Geoprocessing with Python using Open Source GIS section at Utah State, the first thing that came to mind was that there could be a lesson on working with object-oriented and functional programming paradigms and how to take advantage of Python's support for each of them. Not just how to use instances, but how to write a useful class. How to use iterators and generators. A lesson about how to start to move from simple scripting, and the inevitable folders of replicating, mutating files named processing.py, processing2.py, processing2b.py, and so on (believe me, I've been there, only it was EASI scripts), to modules of reusable classes and functions. I've tried to come up with a simple example using the conventional GIS geometry model, something that will be familiar to a GIS analyst.

Actually, I wasn't considering functional programming at all when I started this, but it's obviously very helpful to a GIS programmer. All those nested loops one uses to loop over the elements of rasters? There's a Python function for that. This has been an exercise for me, too. I've found some recipes in Python's itertools module that I've overlooked for too long and are going to immediately improve a few of my own projects. For better or worse, Python lets you mix the object-oriented and functional paradigms, and I do that in the following walk-through of yet another Point class for Python.

The C or C++ library under your favorite Python bindings already provides a point class, so why reinvent the wheel? One day (I'm not predicting when) geometry and topology operation code written entirely in Python will execute quickly enough on enough different platforms to be ready for serious work. Imagine we're coding toward that future. In the meanwhile we'll see how to get more out of today's Python, and dabble in doctest, json, array, and itertools.

I've used Python 2.6 for this exercise, and nothing from outside the standard library. The source for the completed pygis_101 module is at http://bitbucket.org/sgillies/python-gis-101/src/.

Let's start with module containing a do-nothing Point class and a test runner.

# Module: pygis_101.py
# We'll stick to plane geometry

class Point(object):
    """The building block of higher dimensional geometries."""

# and introduce a doctest test runner to the end of the module

def _test():
    import doctest
    doctest.testmod()

if __name__ == "__main__":
    _test()

You can run all doctests in the module if you execute it as a script.

$ python pygis_101.py -v
3 items had no tests:
    __main__
    __main__.Point
    __main__._test
0 tests in 3 items.
0 passed and 0 failed.
Test passed.

Now, let's add a docstring with tests.

class Point(object):
    """The building block of higher dimensional geometries.

    This point class has 2 coordinates, which we'll call `x` and `y`:

        >>> point = Point(0.0, 0.0)
        >>> point.x
        0.0
        >>> point.y
        0.0
    """

Run the tests again

$ python pygis_101.py -v
Trying:
    point = Point(0.0, 0.0)
Expecting nothing
**********************************************************************
File "pygis_101.py", line 8, in __main__.Point
Failed example:
    point = Point(0.0, 0.0)
Exception raised:
    Traceback (most recent call last):
      File "/Users/seang/code/python-macosx/parts/opt/lib/python2.6/doctest.py", line 1231, in __run
        compileflags, 1) in test.globs
      File "<doctest __main__.Point[0]>", line 1, in <module>
        point = Point(0.0, 0.0)
    TypeError: object.__new__() takes no parameters
Trying:
    point.x
Expecting:
    0.0
**********************************************************************
File "pygis_101.py", line 9, in __main__.Point
Failed example:
    point.x
Exception raised:
    Traceback (most recent call last):
      File "/Users/seang/code/python-macosx/parts/opt/lib/python2.6/doctest.py", line 1231, in __run
        compileflags, 1) in test.globs
      File "<doctest __main__.Point[1]>", line 1, in <module>
        point.x
    NameError: name 'point' is not defined
Trying:
    point.y
Expecting:
    0.0
**********************************************************************
File "pygis_101.py", line 11, in __main__.Point
Failed example:
    point.y
Exception raised:
    Traceback (most recent call last):
      File "/Users/seang/code/python-macosx/parts/opt/lib/python2.6/doctest.py", line 1231, in __run
        compileflags, 1) in test.globs
      File "<doctest __main__.Point[2]>", line 1, in <module>
        point.y
    NameError: name 'point' is not defined
2 items had no tests:
    __main__
    __main__._test
**********************************************************************
1 items had failures:
   3 of   3 in __main__.Point
3 tests in 3 items.
0 passed and 3 failed.
***Test Failed*** 3 failures.

We'll need to override object.__init__() in order to pass in coordinate values when we create point instances. The values will be stored in instance attributes.

class Point(object):
    """The building block of higher dimensional geometries.

    This point class has 2 coordinates, which we'll call `x` and `y`:

        >>> point = Point(0.0, 0.0)
        >>> point.x
        0.0
        >>> point.y
        0.0
    """

    def __init__(self, x, y):
        self.x = x
        self.y = y

Now, the 3 tests pass

$ python pygis_101.py -v
Trying:
    point = Point(0.0, 0.0)
Expecting nothing
ok
Trying:
    point.x
Expecting:
    0.0
ok
Trying:
    point.y
Expecting:
    0.0
ok
3 items had no tests:
    __main__
    __main__.Point.__init__
    __main__._test
1 items passed all tests:
   3 tests in __main__.Point
3 tests in 4 items.
3 passed and 0 failed.
Test passed.

With that brief introduction to Python classes, test-driven development, and doctest finished, let's move on to making the Point class more useful.

Often, one wants to serialize geometries for use with other programs. One useful format is "well-known text", or "WKT". Let's add the following to the Point class:

def toWKT(self):
    """Returns a WKT representation of a point.

        >>> point = Point(0.0, 0.0)
        >>> print point.toWKT()
        POINT (0.000000 0.000000)
    """
    return 'POINT (%f %f)' % (self.x, self.y)

JSON is a rising format. Python and JSON go together like sausage and beer. You might be tempted to produce JSON by string interpolation, as we did for WKT, but it's best to use the json module included with Python 2.6, or simplejson installed separately for earlier Python versions. Don't forget to import json at the top of the module.

def toJSON(self):
    """Returns a GeoJSON representation of a point.

        >>> point = Point(0.0, 0.0)
        >>> print point.toJSON()
        {"type": "Point", "coordinates": [0.0, 0.0]}
    """
    return json.dumps(
        dict(type='Point', coordinates=(self.x, self.y))
        )

If we add a __geo_interface__ attribute (implemented as a property) that returns a dict like the one above, we could even use the point objects with the Shapely and GeoJSON packages.

@property
def __geo_interface__(self):
    """Returns a GeoJSON-ish dict representation of a point for interop
    with some other Python + GIS software.

        >>> point = Point(0.0, 0.0)
        >>> point.__geo_interface__
        {'type': 'Point', 'coordinates': (0.0, 0.0)}
    """
    return dict(type='Point', coordinates=(self.x, self.y))

def toJSON(self):
    """Returns a GeoJSON representation of a point.

        >>> point = Point(0.0, 0.0)
        >>> print point.toJSON()
        {"type": "Point", "coordinates": [0.0, 0.0]}
    """
    return json.dumps(self.__geo_interface__)

Now we've got 9 passing tests.

$ python pygis_101.py -v
...
4 items passed all tests:
   3 tests in __main__.Point
   2 tests in __main__.Point.__geo_interface__
   2 tests in __main__.Point.toJSON
   2 tests in __main__.Point.toWKT
9 tests in 7 items.
9 passed and 0 failed.
Test passed.

And, btw, 100% code coverage:

$ coverage run pygis_101.py
$ coverage report
Name        Stmts   Exec  Cover
-------------------------------
pygis_101      16     16   100%

A geometry object ought to be able to provide its bounding box. In order to show how to do this for many more kinds of geometries, let's rewrite Point's coordinate storage using arrays, replace the x and y attributes with dynamic properties for convenience, and add a move method for shifting the point to new coordinates. Then we get the bounding box by checking the minimum and maximum values of the coordinate arrays. This technique will work for just about any type of geometry. Use from array import array at the top, of course.

def __init__(self, x, y):
    self.xs = array('d')
    self.xs.append(x)
    self.ys = array('d')
    self.ys.append(y)

@property
def x(self):
    return self.xs[0]

@property
def y(self):
    return self.ys[0]

def move(self, x, y):
    """Move point to new coordinates x, y.

        >>> point = Point(0.0, 0.0)
        >>> point.move(1.0, 1.0)
        >>> point.x
        1.0
        >>> point.y
        1.0
    """
    self.xs[0] = x
    self.ys[0] = y

@property
def bounds(self):
    """Returns the bounding box as a (minx, miny, maxx, maxy) tuple.

        >>> point = Point(0.0, 0.0)
        >>> point.bounds
        (0.0, 0.0, 0.0, 0.0)
    """
    return min(self.xs), min(self.ys), max(self.xs), max(self.ys)

All production geometry-for-GIS software I know implements these kinds of coordinate sequences. To get coordinate pairs from the sequences, operate on the coordinate arrays with Python's built in zip method or itertools.izip, which I mention because the itertools documentation is a motherlode of fast and efficient iteration recipes.

>>> from array import array
>>> xs = array('d', [0.0, 1.0, 2.0])
>>> ys = array('d', [0.0, 1.0, 2.0])
>>> from itertools import izip, islice
>>> list(izip(xs, ys))
[(0.0, 0.0), (1.0, 1.0), (2.0, 2.0)]
>>> list(islice(izip(xs, ys), 1, None))
[(1.0, 1.0), (2.0, 2.0)]

Zip, btw, is its own inverse.

>>> zip(*zip([1,2,3], ['a', 'b', 'c']))
[(1, 2, 3), ('a', 'b', 'c')]

Going further, we can consolidate the toWKT and toJSON methods in a way that also works very efficiently for lines or paths (using list() instead of next() in that case):

@property
def coords(self):
    """Iterate over x,y pairs.

        >>> point = Point(0.0, 0.0)
        >>> next(point.coords)
        (0.0, 0.0)
    """
    return izip(self.xs, self.ys)

def toWKT(self):
    """Returns a WKT representation of a point.

        >>> point = Point(0.0, 0.0)
        >>> print point.toWKT()
        POINT (0.000000 0.000000)
    """
    return 'POINT (%f %f)' % next(self.coords)

@property
def __geo_interface__(self):
    """Returns a GeoJSON-ish dict representation of a point for interop
    with some other Python + GIS software.

        >>> point = Point(0.0, 0.0)
        >>> point.__geo_interface__
        {'type': 'Point', 'coordinates': (0.0, 0.0)}
    """
    return dict(type='Point', coordinates=next(self.coords))

Now we've got 17 passing tests and still 100% coverage.

$ coverage run pygis_101.py
$ coverage report
Name        Stmts   Exec  Cover
-------------------------------
pygis_101      31     31   100%

Next chance I get, I'll dig into some of my old blogged code to see what might be better done in a more functional style.

This release fixes a bug in writing out WKB on 64-bit systems. It was solved by a MapFish user. If you are one, or are writing WKB to PostGIS or another application, an upgrade is recommended. I've uploaded an sdist to http://pypi.python.org/pypi/Shapely/1.0.14.

HTML is the "HTML of the GeoWeb". From Mano Marks:

For a basic Google Maps mashup (or Bing, or Yahoo! Maps too), you need at least HTML, CSS, and JavaScript. More complicated mashups will use KML, GeoRSS, or GeoJSON, maybe Flash instead of JS, and wait. Notice the first one of those technologies, HTML. Still fundamentally, the GeoWeb relies on HTML. Sure, with a GeoWeb browser like Google Earth, you don't need HTML, though you can use it in the description balloons. But most mashups require HTML as the carrier.

Actually, you do still need HTML with Google Earth. Lacking forms or other intrinsic publishing mechanisms, KML is a format for read-only applications. Theoretically, one might put an HTML form in a placemark bubble and use that to modify the state of the resource behind the KML representation. Or javascript to modify the resource via XHR like many of us do with OpenLayers. Unfortunately, Google Earth doesn't surface forms in popups or fully support javascript. We'll have to wait for more convergence between Earth and Chrome, perhaps. Meanwhile, I'm making users hop between Earth and their browser. It's easy enough to cut and paste the placemark's KML text, but imagine what you could do with javascript and DOM-like access to the nodes in a KML doc.

body_html_template