Tobin Bradley writes:

We found out a few days before the conference that we had won the G. Herbert Stout Award for Visionary Use of GIS for our REST Web Services Framework. Mrs. Stout was there to present the award, and it was a very humbling and happy experience for us. Congratulations to the Mecklenburg County GIS staff, and congratulations to Asheville for winning in the city category - it was definitely well deserved!

On the one hand, it's neat that anything "REST" wins an award in geospatial. On the other hand, REST can't take any credit because there actually isn't any of it in these services. And that's fine (except for the misleading label) because Mecklenburg County GIS doesn't need REST for this. REST is a style for huge, distributed information architectures that need to last for decades, like the Web (or even a national spatial data infrastructure).

I'd like to see more GIS developers follow the lead of CloudMade and tout HTTP APIs. Not only would it be in almost all cases more truthful, there's the advantage to being able to point users to the HTTP specification, HTTP libraries and tools, and not having to explain why there's no "REST specification".

Update (2009-02-23): Where are my manners? Congratulations on the award. Making an HTTP API that users like is no small thing, and neither is open sourcing a well-documented implementation to a GIS community that still tends to be locked up by proprietary vendors.

I haven't started working on any pure Python (for App Engine) geometry predicates or operators, or any other alternatives to GEOS, but I've begun to make it possible. All GEOS dependency is being moved to its own shapely.geos module, which will be the default plugin provider for Shapely's new plugin framework.

The framework defines a few interfaces (geometry compiling or checking provider, area and length provider, bounds provider, etc) and entry points. When a geometry is asked to compute its bounds, for example, it delegates to a provider that has been wired up to the geometry factory on import. Here's a peek at how the default providers are registered in Shapely's setup.py:

...
setup(name          = 'Shapely',
      version       = '1.1.0',
      ...
      entry_points = """
          [shapely.geometry]
          geometryChecker=shapely.geos.geometry:geometryChecker
          metricProperties=shapely.geos.metrics:metricProperties
          geometryProperties=shapely.geos.topology:geometryProperties
          """,
)

And in shapely.geometry, the providers are activated with a function:

...
use('Shapely>=1.1.0')

from geo import shape, asShape
from point import Point asPoint
from linestring import LineString, asLineString
from polygon import Polygon, asPolygon
from multipoint import MultiPoint, asMultiPoint
from multilinestring import MultiLineString, asMultiLineString
from multipolygon import MultiPolygon, asMultiPolygon
from collection import GeometryCollection

Users can override the default GEOS providers by writing and installing (with setuptools) new packages that provide the same named end points, and "using" those at run time:

from shapely import geometry
geometry.use('MyGeometry')
...

point = geometry.Point(0.0, 0.0)
x = point.buffer(10.0) # calls on provider defined by MyGeometry package

In theory, this makes it possible to a write an application using Shapely that can run on either C Python, Jython, or IronPython using the appropriate backend for the platform, as long as setuptools and pkg_resources work. I've read that they will in the next Jython release. IronPython seems to be a little farther behind. In practice, the plugin framework is helping to improve the testability and quality of Shapely's code.

I think this is the first time I've been translated. Thank you, René-luc D'Hont. I read French just well enough to verify that it's a faithful translation. I disagree with the comment at the end about GML's fitness: GML is XML, which is right for the Web; GML has a unique content type (not just text/xml), which is also right for the Web. The problems with GML are XML Schema and that it's not RDF.

Shorter Kurt Schwehr: your data needs a cool URI.

Bryan Lawrence is another scientist that blogs about this issue, including the messy details of citing data in traditional hypertext-less paper journals. Of course, getting credit for publishing digital data is a completely different, social issue, which probably has to wait for a generational change in the sciences.

So, I'm going to be in Montpellier, France, for one year starting 1 June. I'm completely new to the neighborhood. I've got the First Open Source GIS UK Conference, EuroPython 2009, and OGRS 2009 on my calendar. What else is going on?

The answer to:

@sgillies What's the beef with OGC WMS and WFS?

requires a few more than 140 characters.

First, let me review the good about the OGC service architecture and its W*S specs. The OGC has made interoperability a top priority in GIS. Everybody recognizes this is a huge accomplishment. I do too. My favorite byproduct is the increasing priority of open access. It's no accident. The OGC intended that interoperability would lead to more open access to data, and it has. It's a wonderful thing. My other favorite, and perhaps more accidental, byproduct is that thinking of GIS services as interchangeable commodity components leads rather quickly to considering open source implementations. I also think the OGC has done a fine job identifying and standardizing the parameters of our common processes, and a generally good job on message formats. So much good, I must be in heaven, right?

My beef with W*S is that its architects didn't do their Web homework. Despite the "Web" in the name, service design isn't informed by Web architecture and the understanding of HTTP (the Hypertext Transfer Protocol) begins and ends with CGI (the Common Gateway Interface). W*S understands and uses the Web as an alchemist understood and used the elements. We bear the cost of needless reinvention: "Update sequence" instead of HTTP Expiration and Validation, "Web Geolinking Service" instead of standard HTTP interaction, "GeoDDS" instead of Atom. Despite the idea that W*S are designed to be transport-neutral, HTTP is the only significant "distributed computing type" (what the architects call "transport"). The USGS Framework WFS uses no other transport than HTTP. GeoBase uses no other transport than HTTP. Still, our "Web" services remain things that are not really of the Web.

Another minor beef is that in our interoperability fervor we have made standardization holy. The GIS community largely believes that standards should come before implementation, should be built in clean rooms by an elite group of standards scientists, and this stifles innovation. I depend on standards as much as anyone, but I feel we should be standardizing on best practices more than we currently do.

I'm going to use the announcement of Nanaimo's "authentic Web" GIS as an occasion to debunk some myths about REST and the Web, and their fitness for designing alternatives to the OGC's service architecture, that surfaced on Twitter last week.

Myth: RESTful Web services aren't based on standards.

Indeed, there are APIs on the programmable web touting "REST" which are very unlike each other. Not all of them are even RESTful when you get right down to it. They come from different and varying domains. It's understandable that a quick glance leaves some with the mistaken impression that interoperability is not a property of RESTful Web services.

First, and this can't be said enough, because it still isn't really understood in the GIS community: REST is a particularly constrained style of architecture which just so happens to be the style of the World Wide Web. It is not a standard, but needs and shapes standards. Interoperability depends on standards, whether your architecture is RESTful or not. Paul Prescod, who I'm quoting often these days, enumerates the necessary kinds of standards:

In application-level networking, there are three basic things to be standardized

• Addressing -- how do we locate objects
• Methods or Verbs -- what can we ask objects to do
• Message payloads or Nouns -- what data can we pass to the objects to ask them to accomplish their goals

For RESTful Web services, the first two are standardized by HTTP/1.1:

The Hypertext Transfer Protocol (HTTP) is an application-level protocol for distributed, collaborative, hypermedia information systems. It is a generic, stateless, protocol which can be used for many tasks beyond its use for hypertext, such as name servers and distributed object management systems, through extension of its request methods, error codes and headers. A feature of HTTP is the typing and negotiation of data representation, allowing systems to be built independently of the data being transferred.

HTTP has been in use by the World-Wide Web global information initiative since 1990. This specification defines the protocol referred to as "HTTP/1.1", and is an update to RFC 2068.

You wanted standards? Even OGC service specifications acknowledge HTTP/1.1 (though not without misusing it, more on that in a future post). And HTTP/1.1 has been shaped by REST.

As Prescod points out, RESTful Web services push all interoperability problems into the third standardization category: message payload. This is a conscious decision. Interoperability may not be perfectly solvable, but you can isolate its problems, and RESTful services should do so. Do read all of Prescod's article on standardization. We spend way too much time talking about the "which" of standards in GIS without really thinking about the "what".

Myth: RESTful Web services aren't "lights out" accessible.

In fact, a properly RESTful service has better accessibility than an OGC service. To use the same analogy, what if you drop your special OGC service client in the dark and can't find it? How do you access your OGC service? Pardon me, but you're screwed: standing knee-deep in other web programming tools, and none of them can make sense of the OGC's unique addressing schemes and unique methods of interaction. With a RESTful service you can poke at it in a standardized way (HTTP/1.1 again) with curl, or XHR, or whatever, and get an actionable description. Of course, you're a GIS geek, and you keep your OGC client firmly attached to a reeling key chain, but consider the other agents on the Web that don't have an OGC service client at all. The OGC's special addressing scheme and methods make it very hard for those agents to get even a partial understanding of the service.

Myth: REST is too immature for GIS.

In 2000, Roy Fielding wrote:

Since 1994, the REST architectural style has been used to guide the design and development of the architecture for the modern Web. This chapter describes the experience and lessons learned from applying REST while authoring the Internet standards for the Hypertext Transfer Protocol (HTTP) and Uniform Resource Identifiers (URI), the two specifications that define the generic interface used by all component interactions on the Web, as well as from the deployment of these technologies in the form of the libwww-perl client library, the Apache HTTP Server Project, and other implementations of the protocol standards.

Web architecture has been done in the REST style since 1994. HTTP/1.1 began to be adopted in 1997. The Web itself goes back to 1990. The OGC's service architecture is not more mature than this.

Myth: RESTful Web services are fine for small solutions, not for large, inter-agency solutions.

The World Wide Web is our largest network application. It's global. Interoperability is promoted through the use of common hypertext formats. One would think that geographic information systems using the same architecture, in the same style, could scale equally well, right? I think the onus is instead on the naysayers to make the case that their favored architectures can scale like the Web and bridge domains like the Web.

On to cascaded unions for Shapely ...

>>> from osgeo import ogr
>>> from shapely.wkb import loads
>>> ds = ogr.Open('/Users/seang/data/census/co99_d00.shp')
>>> co99 = ds.GetLayer(0)
>>> counties = []
>>> while 1:
...    f = co99.GetNextFeature()
...    if f is None: break
...    g = f.geometry()
...    counties.append(loads(g.ExportToWkb()))
...
>>> len(counties)
3489

Matplotlib makes a pretty picture of those 3489 polygons:

>>> import pylab
>>> from numpy import asarray
>>> fig = pylab.figure(1, figsize=(5.5, 3.0), dpi=150, facecolor='#f8f8f8')
>>> for co in counties:
...    a = asarray(co.exterior)
...    pylab.plot(a[:,0], a[:,1], aa=True, color='#666666', lw=0.5)
...
>>> pylab.show()

Shapely never had the power to dissolve adjacent polygons in a collection before, or at least not over large collections of real-world data. GEOS 3.1's cascaded unions are a big help:

>>> from shapely.ops import cascaded_union
>>> u = cascaded_union(counties)
>>> len(u.geoms)
219
>>> for part in u.geoms:
...     a = asarray(part.exterior)
...     pylab.plot(a[:,0], a[:,1], aa=True, color='#666666', lw=0.5)
...
>>> pylab.show()

There's user interest in leveraging the new reentrant API in GEOS 3.1, and releasing the GIL when calling GEOS functions to improve performance in multithreaded apps. I'm all for it.

I began exposing some of the new features of GEOS 3.1 in Shapely today. Geometries can be tuned up, or "prepared" ala SQL, for efficient batch operations. For example consider this set of random polygons scattered around (0.5, 0.0) and their intersection with a triangle:

>>> from shapely.geometry import Point, Polygon
>>> from random import random
>>> spots = [Point(random()*2.0-0.5, random()*2.0-1.0).buffer(0.1) for i in xrange(200)]
>>> triangle = Polygon(((0.0, 0.0), (1.0, 1.0), (1.0, -1.0)))
>>> x = [s for s in spots if triangle.intersects(s)]
>>> len(x)
67

Without preparing the triangle for batch intersection (or using an index), it takes about 12 ms to get the 67 intersecting spots:

>>> import timeit
>>> t = timeit.Timer('x = [s for s in spots if triangle.intersects(s)]',
...                  setup='from __main__ import spots, triangle')
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100)/100)
11917.51 usec/pass

A prepared triangle finds the same 67 in just a little more than 1/4 of the time:

>>> from shapely.prepared import prep
>>> pt = prep(triangle)
>>> t = timeit.Timer('x = [s for s in spots if pt.intersects(s)]',
...                  setup='from __main__ import spots, pt')
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100)/100)
3145.49 usec/pass

Update (2009-01-26): let's try an example from the real world. The boundary of Larimer County, Colorado, as represented in 2000 Census data, is a polygon with a 370 vertex exterior ring:

>>> from osgeo import ogr
>>> ds = ogr.Open('/Users/seang/data/census/co99_d00.shp')
>>> counties = ds.GetLayer(0)
>>> co069 = counties.GetFeature(1151)
>>> g = co069.geometry()
>>> larimer = loads(g.ExportToWkt())
>>> # nb: OGR properties and methods usually can't be safely chained
>>> len(larimer.exterior.coords)
370

Scatter 200 spots around the neighborhood of Larimer County:

>>> b = larimer.bounds
>>> w = b[2] - b[0]
>>> h = b[3] - b[1]
>>> cx = b[0] + w/2.0
>>> cy = b[1] + h/2.0
>>> from random import random
>>> from shapely.geometry import Point
>>> def spotter():
...     x = (random()*2.0-1.0)*w + cx
...     y = (random()*2.0-1.0)*h + cy
...     return Point(x, y).buffer(0.05)
...
>>> spots = [spotter() for i in xrange(200)]
>>> x = [s for s in spots if larimer.intersects(s)]
>>> len(x)
48

48 spots intersect the county, and it takes about 22 ms to find them using the unoptimized method:

>>> import timeit
>>> t = timeit.Timer('x = [s for s in spots if larimer.intersects(s)]',
...                  setup='from __main__ import spots, larimer')
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100)/100)
21728.65 usec/pass

It only takes about 1.5 ms using a prepared geometry, a speed-up of about 14X:

>>> from shapely.prepared import prep
>>> pt = prep(larimer)
>>> t = timeit.Timer('x = [s for s in spots if pt.intersects(s)]',
...                  setup='from __main__ import spots, pt')
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100)/100)
1570.41 usec/pass