Sean Gillies Blog

A while back my colleague Hugh Cayless wrote a Djatoka layer for use with OpenLayers. It is deployed at UNC's library, allowing you to view over 3000 North Carolina Maps. View the source on any map and you'll quickly get a sense of how it's rigged up. You pass a URL encoding requests for image regions and a URL encoding requests for image metadata to the OpenURL layer constructor, it figures out the grid and zoom levels for you, and you add the layer to a map that's configured using the layer's resolution, maximum extent, and tile size parameters. It works much like a TMS layer, but in pixel units.

I began with this code for a different project and modified it to satisfy some of my project's special needs:

• The imagery of a physical object and its annotations will be shown in their natural units – centimeters.
• Square tiles will be used instead of rectangular tiles with the same aspect ratio as the physical object (approximately 13:1).
• Browser and image server are in different domains, complicating discovery of image metadata.

I decided to go with Djatoka instead of an open source GIS map server like GeoServer or MapServer because the project, while spatial in nature, doesn't require cartographic projections or spatial indexing; I was interested in trying something new; and NYU's already got one running (Thanks, Hugh!). I ended up modifying Hugh and Cliff's OpenURL.js considerably. Now, it's more like a hybrid between OpenLayers.Layer.Image and OpenLayers.Layer.WMS. The code is at http://github.com/sgillies/tpeutb/blob/master/OpenURL.js.

Unlike a GIS map server, Djatoka doesn't pad requests: it fails if you ask for regions outside the image. To have square pixels, I've had to pad my imagery a bit. It's now 43008 pixels wide by 3072 pixels high, exactly 14:1. This makes the right tile size 192 x 192 pixels. As things stand, my new OpenURL layer can't be used in the same map as Hugh's, so I haven't bothered to rename it yet. It's used like this:

var RES = 1.0/60.3; // Centimeters per pixel
var size = new OpenLayers.Size(43008, 3072);
var origin = new OpenLayers.LonLat(-16.6222, -1.2174);
var extent = new OpenLayers.Bounds(
    origin.lon, origin.lat, origin.lon + size.w*RES, origin.lat + size.h*RES);

// Zoom factors of 1/8 to 4
var resolutions = [8.0*RES, 4.0*RES, 2.0*RES, 1.0*RES, 0.5*RES, 0.25*RES];

// At zoom level 0, the image is covered by 2 x 28 tiles
var tileSize = new OpenLayers.Size(192, 192);

var mosaic = new OpenLayers.Layer.OpenURL(
  'All section, mosaicked',
  'http://dl-img.home.nyu.edu/',
  'http://pipsqueak.atlantides.org/mosaic.tif', // for example, 404 in reality
  extent,
  size,
  { format: 'image/jpeg',
    units: 'cm',
    tileSize: tileSize,
    resolutions: resolutions }
  );

var map = new OpenLayers.Map(
  'map',
  { units: 'cm',
    tileSize: tileSize,
    resolutions: resolutions }
  );

A very simple application of this layer can be seen at http://pipsqueak.atlantides.org/tpeutv/djatoka.html.

Shapely's had a shape() function that makes geometric objects from GeoJSON-like mappings for some time. Now, by popular request, is the inverse: mapping():

>>> from shapely.geometry import Point, mapping
>>> m = mapping(Point(0, 0))
>>> m['type']
'Point'
>>> m['coordinates']
(0.0, 0.0)

The new function operates on any object that provides __geo_interface__ [link]. This release also fixes a major bug involving GEOS versions < 3.1 introduced in 1.2.2. You can download an sdist from http://pypi.python.org/pypi/Shapely/. Windows installers should be coming soon.

As I like to remind people every few releases, major portions of this work were supported by a grant (for Pleiades) from the U.S. National Endowment for the Humanities (http://www.neh.gov).

Update (2010-08-22): Win32 and Win-AMD64 installers have been uploaded to the locations above.

We're back in Fort Fun. After an initial snafu the travel gods smiled on us and the flights and connections went more smoothly than we had expected. There was a bit of airsickness and our poussette went missing, but we got some sleep on the long CDG to ATL leg, Customs didn't balk at the number of wine bottles we brought in our luggage, and we arrived not too burned out to celebrate with take-out pizza and a St. Lupulin. The kids even slept until five o'clock this morning. I'm declaring it a win.

The help and hospitality of many people made our séjour an amazing experience. Thank you, everyone, and particularly our nounou, Célia Mourrut; our hosts and friends at CBGP, Denis Bourguet, Arnaud Estoup, Benoît Facon, Renaud Vitalis, and Miguel Navascues; our friend across the street at the USDA, René Sforza; in Heidelberg the Epigraphische Datenbank Heidelberg team: Francisca Feraudi-Gruénais, Brigitte Graef, and James Cowey; in Montpellier, René-Luc D'Hont (3LIZ) and Gerald Fenoy (Zoo Project), and Vincent Heurteaux, Adrian Custer, and Martin Deruisseaux from Geomatys; in Torino, Giorgio Borelli and Silvio Tomatis (collective.geo and more), and Stefano Costa (Open Archaeology); in Isère (the Department of), Eric Lemoine. I'm looking forward to returning the favor some day.

Next steps are making some fresh salsa (I could never find ripe tomatoes, onion, cilantro, and chiles simultaneously at Montpellier's Arceaux market), unpacking, and getting some fall season vegetables planted.

Waiting at the Louis Blanc tram stop after probably our last family dinner in the centre ville – at Le Roule ma Poule in the relatively sleepy north side of the Ecusson – we saw a guy using the tramway as bike lane.

If you keep an eye on the board at every stop and time it right, you can roll along unimpeded.

Into the sunset.

Au revoir, Montpellier.

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.

Waouh (scroll down, no internal anchors), this is a big deal to me:

• OGC TC directs the OGC-NA that all new OGC identifiers issued for persistent public OGC resources shall be http URIs, instead of URNs
• New standards and new major versions of existing standards shall use http URIs for persistent public OGC resources to replace OGC URN identifiers defined in previous standards and versions, unless OGC-NA approves an exception

Use of HTTP URIs is axiomatic for "linked data". This begins to make OGC stuff relevant to the web.

Update (2010-07-08): More on Allan Doyle's blog.

Here is ISAW's Associate Director for Digital Programs (my boss), Tom Elliot, on collaboration between Pleiades and DARMC [Ramping up Pleiades 2]:

Meanwhile, we've been in dialog with Michael McCormick, Guoping Huang and Kelly Gibson at Harvard. They're the driving force behind the Digital Atlas of Roman and Medieval Civilization, with whom we're collaborating under the new grant. Our aim is to collate and share the datasets assembled by both projects and to cross-link our web applications. This will bring more accurate coordinates for many features into Pleiades, as well as a number of new features that will expand our time horizon into the middle ages. You'll get a choice of display and map interaction modes and, eventually, the ability to move back and forth between both resources. We'll keep you posted as the timeline for this portion of the work is refined.

GeoScript now has something like descartes [plotting geometries]:

The plotting functionality makes use of the JFreeChart library, a popular open source Java framework for creating diagrams and charts. Another example of one of the benefits of the marriage of Java and Python that is Jython.

It reminds me of the run of such posts I went on this past spring. JFreeChart might be great, but Java platform developers shouldn't underestimate what matplotlib brings to the C Python platform: loads of tools for earth science visualization and excellent integration with tried and true numerical libraries. That said, I'll freely admit to some JTS envy.

BTW, Shapely 1.2.1 is now accepted to Debian unstable.