Sorting features on spatial relations

I'm working on a map that shows all locations associated with an ancient place. The map is part of the HTML representation of the place (or more specifically: the resource that describes the place) and the data that goes into the map comes from a JSON representation of the place. I've written before about this hypertext-based design.

For this particular application, the ancient place is modeled as a GeoJSON feature collection and locations as GeoJSON features. A feature collection has a list of features and, GeoJSON being under-specified in comparison to other GIS standards, clients can't expect any particular ordering of the features. OpenLayers renders them in the order they are listed. We were running into trouble with feature collections that listed large feature patches after smaller feature patches or point features. The big features would be rendered last by OpenLayers and would pave over the smaller features, making it impossible to select them.

OpenLayers should continue to render features in the given order; predictable behavior is a very good thing. The answer, for us, is to get the features properly sorted in our JSON representation: a feature should be listed before the features it contains. Solved. We don't even need fancy geoprocessing to do it. This is easily accomplished using the sorting capability built into Python as sorted(), with a little help from Shapely.

Let's say we have 4 stereotypic features: a point that is contained by a polygon which is itself contained by another polygon, and a free spirited point contained by none.

>>> a = {'type': 'Point', 'coordinates': [2, 2]}
>>> b = {'type': 'Polygon', 'coordinates': [[[1, 1], [1, 3], [3, 3], [3, 1]]]}
>>> c = {'type': 'Polygon', 'coordinates': [[[0, 0], [0, 4], [4, 4], [4, 0]]]}
>>> d = {'type': 'Point', 'coordinates': [-1, -1]}

And that copies of these are collected into a list like

>>> features = [c, a, d, b, c]

Render this with OpenLayers and the big polygon (c) covers everything except point (d). Rendering them in reverse doesn't help because there's another big polygon at the other end of the list.

Sorting using the default key (comparison by characters) seems to almost produce the right order for the map:

>>> from pprint import pprint
>>> pprint(sorted(features, reverse=True))
[{'type': 'Polygon', 'coordinates': [[[1, 1], [1, 3], [3, 3], [3, 1]]]},
 {'type': 'Polygon', 'coordinates': [[[0, 0], [0, 4], [4, 4], [4, 0]]]},
 {'type': 'Polygon', 'coordinates': [[[0, 0], [0, 4], [4, 4], [4, 0]]]},
 {'type': 'Point', 'coordinates': [2, 2]},
 {'type': 'Point', 'coordinates': [-1, -1]}]

The points are ordered after the polygons, which is good, but the bigger polygons still cover the contained one, not what we want. We'll need a different sorting key, one that considers whether features are spatially within another.

As explained in the Python Sorting HowTo, we can define a key function that operates on each list element and returns a value for comparison. Our key function will be a wrapper class that implements __lt__() using Shapely's binary within() predicate and asShape() adapter.

from shapely.geometry import asShape

class Within(object):
    def __init__(self, o):
        self.o = o
    def __lt__(self, other):
        return asShape(self.o).within(asShape(other.o))

As the howto says, the less than comparison is guaranteed to be used in sorting. That's what we'll rely on to spatially sort, and is why we're using within() in reverse instead of the contains() predicate. Trying it out on features b and c, we see that it works:

>>> b < c
False
>>> Within(b) < Within(c)
True

Now, the wrapper is used to sort by spatial containment like so

>>> pprint(sorted(features, key=Within, reverse=True))
[{'type': 'Point', 'coordinates': [-1, -1]},
 {'type': 'Polygon', 'coordinates': [[[0, 0], [0, 4], [4, 4], [4, 0]]]},
 {'type': 'Polygon', 'coordinates': [[[0, 0], [0, 4], [4, 4], [4, 0]]]},
 {'type': 'Polygon', 'coordinates': [[[1, 1], [1, 3], [3, 3], [3, 1]]]},
 {'type': 'Point', 'coordinates': [2, 2]}]

Voila: feature d is within no other feature and is listed first; the other features are listed as we wish. This little trick is now bound for the Shapely manual.