MapServer has very basic geometry support in its mapscript module, and is adding some capabilities from GEOS, but the MapServer shapeObj suffers from a lack of separation between geometry and feature, and lacks support for geometry collections. MapServer features can have only one geometry. MultiPoint, MultiLineString, and MultiPolygon geometry types are not fully supported.
The case against using geometries from ogr.py is harder to make. OGR has none of the limitations of MapServer listed above (except that of single geometry, though this is more a limitation of the feature model), and also draws upon GEOS for the full range of predicates and operations. It is complete and powerful, no doubt about it. The only thing it lacks is usability. It has an un-Pythonic API, a single Geometry class (no Point, LineString, etc), and doesn't install itself properly. The GDAL and OGR modules would be much better organized as a real Python package, gdal for example, with gdal.ogr, and gdal.osr. As it is, the modules are turned loose in the top of the Python module namespace, a practice that is rightly discouraged. My last minor beef about ogr.py geometries is that they are lumped into the same module as the data drivers and features. I prefer some separation, so that I can deploy in smaller chunks, building up exactly just what I need. These are all minor warts, but they add up to the last 10-15% in missing usability for a Python programmer. Howard Butler has been working with Frank to chip away at this number, but I don't see how they can go any farther without compromising OGR's commitment to backwards compatibility.
A while back I committed a SWIG interface file to the GEOS project that allowed one to generate a Python geos.py module, but over time I have become a bit disenchanted with SWIG. The Python C API is no harder to understand than SWIG, and allows me to get right to the objective: a Python extension module with all the power of the underlying C/C++ library and the right Python feel. There's also a Windows platform problem unique to Python which makes the GEOS SWIG bindings unusable for many applications. The Python community is still straddling the 2.3/2.4 fence right now, particularly because Zope and Plone are not yet ready to go to Python 2.4. GEOS cannot be built using MSVC 6, which is the compiler for Python 2.3.x. However, we can build GEOS using MSVC 7+ and then link against the GEOS DLL using MSVC 6. This mandates that the new Python geometry packages use the GEOS C API rather than the C++/SWIG bindings.
Howard recently showed me a pure Python implementation of the OGC Simple Features for SQL spec where geometry classes were derived from Python's list, and imported ogr.py or geos.py to get spatial predicates and geometry operations. I liked it a lot, but the dependencies on ogr.py and geos.py, for all the above reasons, were non-starters for me. So, I did a rewrite using the Python and GEOS C APIs instead of geos.py [source].
The module contains GEOSGeometry and GEOSCoordinateSequence classes, and a GEOSError exception that encapsulates errors raised by the GEOS library. An accompanying Python cartography.spatial.geometry module implements
- Point
- LineString
- LinearRing
- Polygon
- MultiPoint
- MultiLineString
- MultiPolygon
classes that extend GEOSGeometry and provide the Python feel that I'm looking for.
The PCL-Spatial component depends only on the GEOS library (2.2.2) and PROJ.4. There's no requirement for any data or cartographic packages, which should make it more useful for simple and specialized applications. Here's an excerpt from the new usage example:
# ============================================================================
# Geometries
#
# from cartography.spatial.geometry import Geometry, Point, LinearRing, Polygon
#
# Create two triangular geometries, one from a well-known text string, the
# other from points. Next, get their intersection and union as two new
# geometries.
# ============================================================================
wkt = 'POLYGON ((-1.0 50.5, -0.5 51.2, 0.3 50.9, -1 50.5))'
tri_a = Geometry.fromWKT(wkt)
shell = LinearRing([Point(-0.7, 50.3), Point(0.1, 51.0), Point(0.6, 50.1),
Point(-0.7, 50.3)])
tri_b = Polygon(shell)
# the intersection and union of these geoms
inter = tri_a.intersection(tri_b)
union = tri_a.union(tri_b)
The new component is currently in PCL's featuretypes branch, and will be merged into the trunk first thing next week. We'll be basing our work at the Dublin Sprint on it.
