Shapely recipes

I'm committed to keeping the Shapely API small. If it were a Swiss Army knife, it would be a slim one that fits comfortably in your hand and doesn't require you to wear a belt to keep your pants up when you pocket it. This means saying no to features. I recently read a person's opinion that Shapely should support SRIDs – like PostGIS does, I assume the reasoning to be. The GEOS library does in fact implement spatial reference system identification numbers for geometric objects, but I'm leaving this out of Shapely. Such identification numbers are relevant only to OGC standard (SFSQL) databases. On the web we identify coordinate systems by URIs. In a Python application, a spatial reference system should be a Python object not an integer pointing to some object outside the environment. A premise of Shapely is that if your application really needs to deal with transformations between coordinate systems, you'll know best how to add back this missing complexity. What's actually missing from Shapely is not the SRID concept, but a recipe or two for how to add in that concept.

Another requested feature left out is the splitting of lines. Inclusion of GEOS's linear referencing in the Shapely core barely made the cut (I think it could have been a good example of an extension), but makes it possible to implement line cutting algorithms in not very many lines of code. A naive one is shown below.

from shapely.geometry import LineString, Point

def cut(line, distance):
    # Cuts a line in two at a distance from its starting point
    if distance <= 0.0 or distance >= line.length:
        return [LineString(line)]
    coords = list(line.coords)
    for i, p in enumerate(coords):
        pd = line.project(Point(p))
        if pd == distance:
            return [
                LineString(coords[:i+1]),
                LineString(coords[i:])]
        if pd > distance:
            cp = line.interpolate(distance)
            return [
                LineString(coords[:i] + [(cp.x, cp.y)]),
                LineString([(cp.x, cp.y)] + coords[i:])]

The result (see cut.py):

This figure was made with matplotlib, of course. Control points of the line are shown as circles, and the boundary points are displayed in black.

I can imagine that different users and different applications will want to chop up lines in grossly or subtly different ways. The above returns a list of LineString, but another application might want a MultiLineString or an iterator that nibbles at the input while reading from another iterator of bite sizes. Who knows? Cutting up a line simply isn't a primitive operation. A recipe for doing this type of operation using Shapely seems better to me than including one implementation in the core.

Geometric object constructors that validate their input have been proposed from time to time. Validation is particular and expensive and better (in my opinion) left to a developer. Maybe your input has already been validated; maybe objects deemed "invalid" by JTS or GEOS are useful to you in some way. Shapely's constructors don't validate and aren't cluttered with validate=True arguments or the consequent compound if/else statement. Disappointed? The good news is that it's easy to add back the complexity of validation. The Shapely manual includes a recipe for a validating decorator.

from functools import wraps

from shapely.geos import TopologicalError

def validate(func):
    @wraps(func)
    def wrapper(*args, **kwargs):
        ob = func(*args, **kwargs)
        if not ob.is_valid:
            raise TopologicalError(
                "Given arguments do not determine a valid geometric object")
        return ob
    return wrapper

The decorator could be used like so:

>>> @validate
... def ring(coordinates):
...     return LinearRing(coordinates)
...
>>> coords = [(0, 0), (1, 1), (1, -1), (0, 1)]
>>> ring(coords)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 7, in wrapper
shapely.geos.TopologicalError: Given arguments do not determine a valid geometric object

Not everyone is fond of decorators and might want to do this differently. Subclass the geometry types, write some Zope-style adapters, or monkey-patch asShape – it's your call. If you will need validation, you will know (or discover) best how to exploit Shapely to implement it.

By keeping Shapely light and simple I expect we can keep it easy to use and extend. Indulge me in a food analogy, if you will. Shapely isn't ketchup, it's the ripe tomato: raw, but more useful in a wider range of recipes.