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
>>> Within(b) < Within(c)

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.


Re: Sorting features on spatial relations

Author: Tim Schaub

While this wouldn't provide a better solution than you've come up with, note that the graphicZIndex property of a symbolizer can be used to determine the z-ordering of features rendered by OpenLayers (for all feature types).

Relevant examples:

Re: Sorting features on spatial relations

Author: Sean

Thanks for the pointers, Tim, and for dealing with my crappy commenting system. Entering links is no fun here. Z level via styling is neat, but I think that in this particular case I'd still need to sort to produce an attribute that could be used in a styling rule, yes?

Re: Sorting features on spatial relations

Author: Vivien Deparday

Hi Sean,

I think you are right about producing an attribute to use in the styling rule.

I had a similar issue as you but I was using geoserver to produce the GeoJSON. So I used the graphicZIndez as mentionned by Tim but I had to add a rank field to my PostGIS table that I then used in the styling rule. I calculated the rank field based on the geometry type and the polygon areas, not as elegant as the within relationship.

It worked well but it is not ideal if the data is dynamic and I guess in this case a trigger would be needed to update the rank field so it is nice to see another solution.


Re: Sorting features on spatial relations

Author: Tim Schaub

Right, to use rules you need an attribute, and you'd only be determining symbolizer properties based on a single feature at a time.

If it would work to determine z-index based on geometry type (point on top of polygon) and area (small polygon on top of big polygon), you could provide a custom "context" to an OpenLayers.Style object [1]. This is a bit forced, but your symbolizer zIndex property value could reference a function that returned an integer z-index based on geometry type and area (assuming you know the range of poly area).

Anyway, sorting on the server makes good sense (as long as you can). Another downside of the custom context function is that we don't serialize/deserialize this (if you ever need to persist your style).

[1] styles-context.html example

Re: Sorting features on spatial relations

Author: Sean

Indeed, ordering features by decreasing area would work equally well in my application since a feature can't completely pave over one with a larger area. I got caught up in the excitement of finding something that I could do easily with Python and can't do easily in a RDBMS with SQL. I wouldn't be surprised if it's possible to sort rows of a table based on mutual relationships, but I don't have that degree of SQL fu.

This post is really about a Python programming technique; OpenLayers is just the point of departure.