Sean Gillies Blog

I'm disatisfied with Sourceforge's project news capabilities and starting this next week will be publishing news of the Python Cartographic Library and Cartographic Objects for Zope projects on this site.

Artem is the developer of mapnik, a C++ mapping toolkit, and we have been talking about collaborating on a mapnik-based mapping engine for PCL that might, at a user's preference, be used in place of the current MapServer-based engine. Artem is keen to use Anti-Grain Geometry for rendering. I do not know how its performance compares to GD -- but after a look at the AGG demo page I can not deny that AGG is far more slick and featureful. Mapnik will also include a C++ implementation of OGC-style filters comparable to that of GeoTools.

Depending on how quickly Artem works, and my impression is that he's more productive than I, we may have some results in time for OSG05.

I'm taking a short vacation, and have just wrapped up the last things I wanted to finish before leaving: putting the ZCO demo application online, and setting up a limited role for guests to see the application from inside Zope's management interface. The URL is http://zcologia.com/zco.

The most interesting combination of demo layers are global_mosaic, zpoints/x, and cpoints/dot. The first demonstrates ZCO's capability to consume WMS, while the others show layers derived from data maintained in the ZODB rather than on disk in the form of WKT geometries in a Zope File (zpoints) and points from a Zope catalog query (cpoints). See demo/zdata/points_wkt and demo/zdata/points_catalog within the Zope management interface.

Good news: heard from Steve Lime this morning that my proposal to present a talk about PCL and ZCO at OSG05 has been accepted.

Update (2008-01-24): See http://pypi.python.org/pypi/Shapely/1.0.

Update (2007-06-12): http://zcologia.com/news/478/shapely-geometries-for-python/.

Update (2007-04-10): http://zcologia.com/news/429/geometries-for-python-update/.

Update (2006-06-01): I'm no longer using SWIG to make an interface for GEOS. Instead, I'm working on a better GEOS-based geometry module for Python.

GEOS is the Geometry Engine Open Source, a library for programming with 2D geometries and spatial predicates. Just before the 2.0 release, I committed a SWIG interface file and Python distutils setup script for generating a Python 'geos' module on Linux. Hobu has put together a geos (from GEOS 2.1.1) module for Python 2.4 and is distributing it from http://hobu.biz/software/PyGEOS. This is the module that he and I will be using in our OSG05 workshop.

It's easy to get started with the Python geos module. There's no specific documentation for it, but you can readily adapt the GEOS docs to use with Python; the API is almost exactly the same. Below are examples of using the geos GeometryFactory, WKTReader, and WKTWriter classes, and also of the Union and intersection methods of the geos Geometry class.

For kicks, I'll render the results using PCL's cartography module. First, an example of two triangle polygons rendered separately using the default polygon symbolizer and a simple text symbolizer:

from cartography import data, mapping, rendering, spatial
import geos

# Setup the renderer and styling
# ----------------------------------------------------------------------
FONTSET = '/home/sean/projects/ms_45/mapserver/tests/fonts.txt'
mapper = rendering.MapRenderer('MAPSERVER', fontset=FONTSET)

# Text and Polygon symbolizers
text_symb = mapping.TextSymbolizer('FNAME', {'font-family': 'Vera',
                                             'font-size': 10.0},
                                   placement='AUTO')
poly_symb = mapping.PolygonSymbolizer()

# Style and two rules
pe_style = mapping.ExtrinsicStyle('pe_style')
pe_style.rules.append(mapping.Rule('catch-all', text_symb))
pe_style.rules.append(mapping.Rule('catch-all', poly_symb))

# Create two intersecting triangles using PCL
# ----------------------------------------------------------------------
# We're treating everything as WGS 1984
srs = spatial.SpatialReference(epsg=4326)

t1wkt = 'POLYGON ((-1.0 50.5, -0.5 51.2, 0.3 50.9, -1 50.5))'
triangle1 = spatial.Polygon(wkt=t1wkt)
atts1 = ('A Polygon',)

t2wkt = 'POLYGON((-0.7 50.3, 0.1 51.0, 0.6 50.1, -0.7 50.3))'
triangle2 = spatial.Polygon(wkt=t2wkt)
atts2 = ('B Polygon',)

# Render the two shapes separately
# -----------------------------------------------------------------------
store = data.MemoryFeatureStore([('FNAME', 'string')],
            zip([triangle1, triangle2], [atts1, atts2]), srs=srs)
layer = mapping.Layer(store, 'POLYGON')
view = mapping.View(srs, spatial.Point(-1.2, 51.2), spatial.Point(1, 50))
size = [300, 300]
image = mapper.render([(layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('separate.png', 'wb')
f.write(image.getdata())
f.close()

the resulting image, separate.png:

Use the read() method of WKTWriter to convert the PCL polygons to geos Geometries and then call the intersection method of one geometry with the other as the single argument. A new geometry is returned, and can be rendered by stuffing it into a PCL data store and wrapping with a layer:

# Next, render intersection of the polygons
# -----------------------------------------------------------------------
geom_factory = geos.GeometryFactory()
wkt_reader = geos.WKTReader(geom_factory)
wkt_writer = geos.WKTWriter()

geom1 = wkt_reader.read(triangle1.toWKT())
geom2 = wkt_reader.read(triangle2.toWKT())
g = geom1.intersection(geom2)

# make new PCL geometry
inter = spatial.Polygon(wkt=wkt_writer.write(g))

inter_store = data.MemoryFeatureStore([('FNAME', 'string')],
                                      [(inter, ('Intersection',))], srs=srs)
inter_layer = mapping.Layer(inter_store, 'POLYGON')
image = mapper.render([(inter_layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('intersection.png', 'wb')
f.write(image.getdata())
f.close()

the resulting image, intersection.png:

The union of the two triangle polygons is just as easily obtained:

# Now, render the union of the polygons
# -----------------------------------------------------------------------
geom_factory = geos.GeometryFactory()
wkt_reader = geos.WKTReader(geom_factory)
wkt_writer = geos.WKTWriter()

geom1 = wkt_reader.read(triangle1.toWKT())
geom2 = wkt_reader.read(triangle2.toWKT())
g = geom1.Union(geom2)  # <-- watch out for capitalized union

# make new PCL geometry
union = spatial.Polygon(wkt=wkt_writer.write(g))

union_store = data.MemoryFeatureStore([('FNAME', 'string')],
                                      [(union, ('Union',))], srs=srs)
union_layer = mapping.Layer(union_store, 'POLYGON')
image = mapper.render([(union_layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('union.png', 'wb')
f.write(image.getdata())
f.close()

resulting image, union.png:

MapServer has a built-in graticule object that can be configured in the same way as its scalebar and reference map. PCL's MapServer-based engine, however, is not going to use the built-in graticule. Users will instead be required to render a graticule layer from a set of parallel and meridian lines that I do not just happen to have lying around. In the past I have scripted such processes using the shapelib module, but for some practice with OGR I decided to write a simple script using ogr.py to generate my annotated graticule lines. OGR's pymod/samples/tigerpoly.py was my reference for creating and writing a shapefile of features. Right after I finished, it was pointed out to me that I'd overlooked the pymod/mkgraticule.py script recently committed to CVS. Still, mine has attributes for labeling. After a few more tweaks, I'll use the following script to generate parallel and meridian data to be distributed with PCL.

The graticule must consist of individual short line segments for use with MapServer, which throws out features that cannot be completely projected. This occurs in the case of Orthographic, and possibly other hemispheric type projections. The multitude of lines complicates labeling a bit, but that may be something I can overcome as the script is refined.

We begin by creating shapefile to contain our graticule line segments:

import math
import ogr

# Create an ESRI shapefile of parallels and meridians for a MapServer
# world map.

filename = 'graticule.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
driver.DeleteDataSource(filename)
ds = driver.CreateDataSource(filename)
layer = ds.CreateLayer('graticule', geom_type=ogr.wkbLineString)

Next, define several attribute fields. A 'TYPE' field will let me classify and render differently the parallels and meridians if I so choose. 'VALUE' and 'ABS_VALUE' fields will be useful for labeling the parallels and meridians:

# Create an integer field for classification
# 0: parallel
# 1: meridian
fd = ogr.FieldDefn('TYPE', ogr.OFTInteger)
fd.SetWidth(1)
fd.SetPrecision(0)
layer.CreateField(fd)

# Create two fields for labeling
fd = ogr.FieldDefn('VALUE', ogr.OFTInteger)
fd.SetWidth(4)
fd.SetPrecision(0)
layer.CreateField(fd)

fd = ogr.FieldDefn('ABS_VALUE', ogr.OFTInteger)
fd.SetWidth(4)
fd.SetPrecision(0)
layer.CreateField(fd)

Next, create parallels and insert them into the layer:

# First, the parallels at 10 degree intervals, with one degree resolution
for j in range(1,18):
    y = 10*float(j-9)
    for i in range(0, 360):
        line = ogr.Geometry(type=ogr.wkbLineString)

        # hack: MapServer has trouble within .1 decimal degrees of the
        # dateline
        if i == 0:
            x1 = -179.9
            x2 = -179.0
        elif i == 359:
            x1 = 179.0
            x2 = 179.9
        else:
            x1 = float(i-180)
            x2 = x1 + 1.0

        line.AddPoint(x1, y)
        line.AddPoint(x2, y)

        f = ogr.Feature(feature_def=layer.GetLayerDefn())
        f.SetField(0, 0)
        f.SetField(1, y)
        f.SetField(2, math.fabs(y))
        f.SetGeometryDirectly(line)
        layer.CreateFeature(f)
        f.Destroy()

And now the meridians:

# Next, the meridians at 10 degree intervals and one degree resolution

for i in range(0, 37):
    x = 10*float(i-18)

    # hack: MapServer has trouble within .1 decimal degrees of the
    # dateline
    if i == 0:
        x = -179.9
    if i == 36:
        x = 179.9

    for j in range(10, 170):
        line = ogr.Geometry(type=ogr.wkbLineString)
        y1 = float(j - 90)
        y2 = y1 + 1.0

        line.AddPoint(x, y1)
        line.AddPoint(x, y2)

        f = ogr.Feature(feature_def=layer.GetLayerDefn())
        f.SetField(0, 1)
        f.SetField(1, x)
        f.SetField(2, math.fabs(x))
        f.SetGeometryDirectly(line)
        layer.CreateFeature(f)
        f.Destroy()

# destroying data source closes the output file
ds.Destroy()

Because the Python Cartographic Library doesn't yet support non-EPSG coordinate systems (coming soon), I'm going to use MapServer's scripting module to draw the graticules:

import mapscript

mo = mapscript.mapObj()
mo.setFontSet('/home/sean/projects/ms_44/mapserver/tests/fonts.txt')
mo.selectOutputFormat('image/png')

layer = mapscript.layerObj()
layer.status = mapscript.MS_DEFAULT
layer.data = 'graticule.shp'
layer.setProjection('+init=epsg:4326')
layer.type = mapscript.MS_LAYER_LINE
layer.labelitem = 'VALUE'
layer.classitem = 'TYPE'
li = mo.insertLayer(layer)
layer = mo.getLayer(li)

# render parallels and meridians in same style
style = mapscript.styleObj()
style.color.setHex('#808080')

co = mapscript.classObj()
co.label.type = mapscript.MS_TRUETYPE
co.label.font = 'Vera'
co.label.size = 6
co.label.position = mapscript.MS_AUTO
co.label.mindistance = 200

co.insertStyle(style)
co.setExpression('([TYPE] EQ 0)')
co.label.color.setRGB(0, 255, 0)
layer.insertClass(co)

style.color.setHex('#0000FF')
co.insertStyle(style)
co.setExpression('([TYPE] EQ 1)')
co.label.color.setRGB(255, 0, 0)
layer.insertClass(co)

Next, draw the graticule using a Mollweide projection:

# Mollweide Projection
mo.setSize(600, 300)
mo.setProjection('+proj=moll +lon_0=0.0')
mo.setExtent(-18100000, -9050000, 18100000, 9050000)
im = mo.draw()
im.save('graticule-moll.png')

Looking good, graticule-moll.png:

And finally, an Orthographic projection:

mo.setSize(400, 400)
mo.setProjection('+proj=ortho +lon_0=0.0 +lat_0=40.0')
mo.setExtent(-6400000,-6400000,6400000,6400000)
im = mo.draw()
im.save('graticule-ortho.png')

Looks like we have an ortho-related bug in MapServer or PROJ.4, graticule-ortho.png:

Norman Vine (always a source of fun ideas) and I were speculating a few weeks ago about using MapServer to produce maps of the sky. I finally made the time to demonstrate that it is as straight-forward as we had thought. I'm using PCL with a MapServer-based renderer, but everything should be reproducible with more traditional methods. The star data was obtained from Steve Dick's Corrected Bright Star List, 1997 Astronomical Almanac H2-H31: http://aa.usno.navy.mil/publications/AstroAlmanac/bstar97.html.

To begin, I wrangled the list into an OpenOffice spreadsheet and saved it as a 15-field CSV-formatted file. The following Python script processed the spreedsheet into a PCL feature store. Believe it or not, we can use a lat/long coordinate system for objects in the sky. The right ascension is multiplied by 15 because it is in units of hh:mm:ss, and we take the opposite sign of this value because right ascension increases in the opposite sense of our longitude:

import csv
from cartography import data, spatial

properties = (['FB_DES', 'string'],        # Flamsteed/Bayer Designation
              ['MAGNITUDE', 'float'],      # Visual Magnitude
              ['SPECTRAL_TYPE', 'string']) # such as Sol's G2 V

star_srs = spatial.SpatialReference(['+proj=latlong'])
star_store = data.MemoryFeatureStore(properties)

fd = open('bstar97.csv', 'r')
reader = csv.reader(fd)
for row in reader:
    try:
        # values of our properties
        des = ' '.join([row[1].strip() or row[0].strip(), row[2].strip()])
        mag = float(row[11])
        spe = row[14].strip()

        # the location of the star as a Point
        # right ascension (longitude)
        x = -15.0*(float(row[4]) + float(row[5])/60 + float(row[6])/3600)
        x -= 180.0

        # declination (latitude)
        y = float(row[7]) + float(row[8])/60 + float(row[9])/3600
        p = spatial.Point(x, y)

        # add to the data store
        star_store.add(p, [des, mag, spe])
    except:
        pass

fd.close()

A graticule is as good a visual aid in a star map as it is in a map of the Earth. Set up a graticule (from http://zcologia.com/news/16) and styling rules for its presentation:

from cartography import mapping, query, rendering

grat_store = data.DiskFeatureStore('../grat/graticule.shp')
grat_style = mapping.Style()
line_symb = mapping.LineSymbolizer({'stroke': '#CCCCCC', 'opacity': 0.5})
grat_style.rules.append(mapping.Rule('catch-all', line_symb))

grat_layer = mapping.Layer(grat_store, 'LINE', star_srs)
graticule = (grat_layer, grat_style)

And now the stars themselves. First thing I do is define a sequence that associates filter expressions with marker sizes. The styling rules for different stellar brightness levels will be generated from this 'mags' sequence:

mags = \
[
['f.MAGNITUDE <= 0',       15],
['0 < f.MAGNITUDE <= 1',   11],
['1 < f.MAGNITUDE <= 2',    7],
['2 < f.MAGNITUDE <= 3',    5],
['3 < f.MAGNITUDE <= 4',    3],
['4 < f.MAGNITUDE <= 5',    1]
]

circle = mapping.Mark('circle')
star_style = mapping.ExtrinsicStyle()

for m in mags:
    symb = mapping.PointSymbolizer(mark=circle, size=m[1])
    rule = mapping.Rule(m[0], symb, query.PythonFilter(m[0]))
    star_style.rules.append(rule)

# Text symbolizer to label brightest stars
font = mapping.Font(family='Vera', size=8.0)
text_symb = mapping.TextSymbolizer(label='FB_DES', font=font,
                                       placement='AUTO')
star_style.rules.append(mapping.Rule('m0t', text_symb,
                  query.PythonFilter('f.MAGNITUDE <= 1')))

# Wrap the star data store up in a Layer and associate with the style
star_layer = mapping.Layer(star_store, 'POINT', star_srs)
stars = (star_layer, star_style)

Nothing fancy, just rendering the brighter stars as bigger filled circles. Now we can render a map of the stars of the Northern Hemisphere using a polar Lambert Azimuthal Equal Area projection:

FONTSET = '/home/sean/projects/ms_45/mapserver/tests/fonts.txt'
mapper = rendering.MapRenderer('MAPSERVER', fontset=FONTSET)

srs = spatial.SpatialReference(['+proj=laea', '+lat_0=90.0'])
ul = spatial.Point(-10000000, 10000000)
lr = spatial.Point(10000000, -10000000)
view = mapping.View(srs, ul, lr)
size = [600, 600]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_polar.png', 'wb')
f.write(image.getdata())
f.close()

stars_polar.png:

Polaris is right near the center, Cassiopeia above. Left of center is Vega (alpha Lyra), featured in Carl Sagan's novel Contact. Alpha Centauri appears in the lower left corner. Orion is on the right edge, greatly distorted. LAEA is not an ideal projection for a sky map. I also tried Polar Stereographic and Azimuthal Equidistant projections but ran into some issues with projecting the graticule.

Next we zoom in on Orion, my favorite constellation, using an LCC projection:

srs = spatial.SpatialReference(['+proj=lcc', '+lat_1=10.0', '+lon_0=100.0'])
ul = spatial.Point(-2251981, 2350184)
lr = spatial.Point(1611262, -1237944)
view = mapping.View(srs, ul, lr)
size = [400, 400]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_orion.png', 'wb')
f.write(image.getdata())
f.close()

The result, stars_orion.png:

From left to right: Betelgeuse (alpha Ori) the 10th brightess star in the sky, the Rigel (beta Ori) of many Star Trek references, and Aldebaran (alpha Tau).

Update: Frank Warmerdam educated me on PROJ4 and its use by MapServer during rendering. His suggestion was to "inline" the graticule as I had the stars. Simply change the grat_style to an instance of mapping.ExtrinsicStyle. It's costly, as any extrinsic style, but works. Here is the Polar Stereographic map that I wanted:

chart_srs = spatial.SpatialReference(['+proj=stere','+lat_0=90.0'])
ul = spatial.Point(-14000000, 14000000)
lr = spatial.Point( 14000000, -14000000)
view = mapping.View(chart_srs, ul, lr)
size = [600, 600]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_stere.png', 'wb')
f.write(image.getdata())
f.close()

It has been longer than I intended, but finally I have cut a 0.8 source release of the Python Cartographic Library: PCL-0.8.tar.gz (release notes -- sourceforge.net). There are a few new features, but the major work was significant refactoring that will better support future mapping engines. There are some API changes, hence the minor version number increment. This version of PCL is what I've used in the past week's entries. I'll try to have a Windows binary for Python 2.3.x up for download on Monday, and a new ZCO release later next week.

Here are the latest results of the coverage tests for all the fans of coverage.py:

running non-engine tests ...
 ----------------------------------------------------------------------
Ran 112 tests in 0.347s

OK
Name         Stmts   Exec  Cover
--------------------------------
data           127    103    81%
__init__         0      0   100%
interfaces      21      0     0%
mapping        105    103    98%
query           52     52   100%
rendering        7      0     0%
spatial        185    176    95%
--------------------------------
TOTAL          497    434    87%

running standard engine tests ...
----------------------------------------------------------------------
Ran 26 tests in 0.242s

OK
Name       Stmts   Exec  Cover
------------------------------
data          74     69    93%
__init__       0      0   100%
------------------------------
TOTAL         74     69    93%

running mapserver engine tests ...
----------------------------------------------------------------------
Ran 64 tests in 2.671s

OK
Name        Stmts   Exec  Cover
-------------------------------
data          149    110    73%
__init__        0      0   100%
msbridge       91     60    65%
rendering     460    385    83%
-------------------------------
TOTAL         700    555    79%

Not perfect, but not too bad. More on this release over the weekend.

Update: API documentation is now online at http://zcologia.com:8080/PCL/api

PCL-0.8-win32.zip (sourceforge) contains windows installer and 'dumb' binaries for Python 2.3. Do remember that this requires the MapServer mapscript module and Python modules from GDAL/OGR. I recommend the FWTools bundle as a source for those modules.

Matplotlib (http://matplotlib.sourceforge.net), a 2D plotting library, is attracting well deserved attention from Python users. It is a readily installed and rich environment for 2D plotting and visualization with a fat users guide and many nice examples. There even exists a cartographic module to generate background basemaps for plots of geophysical data. What intrigues me the most is its Agg backend and API that provides an easy way to learn about anti-grain geometry. It seemed like it wouldn't take long to find out if matplotlib and Agg were a good fit to PCL's mapping classes.

I began by revisiting the two slightly overlapping triangle polygons from a previous entry, then creating a polygon symbolizer, and defining map image output parameters. PCL polygon symbolizers have a default opacity of 50% -- we'll see the overlap effect in the output image:

from cartography import mapping, spatial
from cartography.spatial import Point, LinearRing, Polygon

t1wkt = 'POLYGON ((-1.0 50.5, -0.5 51.2, 0.3 50.9, -1 50.5))'
triangle1 = spatial.Polygon(wkt=t1wkt)

t2wkt = 'POLYGON((-0.7 50.3, 0.1 51.0, 0.6 50.1, -0.7 50.3))'
triangle2 = spatial.Polygon(wkt=t2wkt)

# Symbolizer for triangles
symb = mapping.PolygonSymbolizer(fill={'fill': '#FF8080'})

# Define output map parameters
epsg4326 = spatial.SpatialReference(epsg=4326)
view = mapping.View(epsg4326, spatial.Point(-1.2,51.25),
                    spatial.Point(1,50.05))

width = 550
height = 300

Define two utility functions: hex_to_rgb() converts PCL colors to matplotlib's normalized tuple, and transform_point() maps PCL points in world units to matplotlib image units. Image coordinates in matplotlib have their origin at the lower left corner, contrary to convention:

def hex_to_rgb(color):
    """convert a hex color such as #FFFFFF to a red, green, blue tuple"""
    c = eval('0x' + color[1:])
    r = (c >> 16) & 0xFF
    g = (c >> 8) & 0xFF
    b = c & 0xFF
    return (r, g, b)

def transform_point(point, view, width, height):
    """transform from world to image coordinates"""
    dx = (view.lr.x - view.ul.x)/width
    dy = (view.lr.y - view.ul.y)/height
    p = (point.x - view.ul.x)/dx
    l = (point.y - view.lr.y)/dy
    return (p, l)

Import the matplotlib module and render the two triangles. The matplotlib graphic context 'gc' is configured from properties of the PCL polygon symbolizer:

import matplotlib
matplotlib.use('Agg')
from matplotlib.backends.backend_agg import RendererAgg
from matplotlib.transforms import Value

dpi = Value(72.0)
o = RendererAgg(w, h, dpi)
gc = o.new_gc()

for polygon in [triangle1, triangle2]:
    points = [transform_point(p, view, width, height) for p in polygon[0]]

    # Fill
    gc.set_alpha(symb.fill['fill-opacity'])
    rgb = hex_to_rgb(symb.fill['fill'])
    face = (rgb[0]/255.0, rgb[1]/255.0, rgb[2]/255.0)
    o.draw_polygon(gc, face, points)

    # Stroke
    gc.set_alpha(symb.stroke['stroke-opacity'])
    gc.set_foreground(symb.stroke['stroke'])
    gc.set_linewidth(symb.stroke['stroke-width'])
    o.draw_polygon(gc, None, points)

o._renderer.write_png('agg_rendered.png')

the result, agg_rendered.png:

My impressions? So far, I like what I see. The matplotlib API is better suited to PCL's SLD-ish mapping objects than MapServer's mapscript module. That's a big plus for me. I don't yet appreciate why matplotlib excludes shape fill colors from the graphic context, or the quirky choice of lower left for image coordinate origin.

Clearly, the code within the loop over triangles needs to be written in C++ (like Agg) rather than Python to get decent performance when rendering hundreds or thousands of polygons at a time. I haven't written any C++ extensions for Python yet (C only), but this might be my first opportunity.

I couldn't resist writing a real matplotlib/Agg mapping engine for PCL to see how it could handle real-world features. In this example I am using the VMAP0-based world borders data available from http://mappinghacks.com. Usually available, that is -- Schuyler's site is offline today. The world borders shapefile has 3784 polygons and roughly 400,000 vertices in all.

My prototype engine is about 9x slower than the existing MapServer-based engine. I actually had expected it to be even slower than this because I am doing some expensive operations in Python rather than C. Turning PCL geometries into a Python list of vertex tuples for matplotlib and transforming from world to image coordinates are costly. This mere factor of 9 gives me hope that a mapping engine using Agg's C++ API could be at least as fast as MapServer's GD renderer.

The result is below. Its color scheme is from Color Brewer, courtesy of Cindy Brewer. The grey blobs around coastlines are polygon stroke noise; the data is too detailed for a map of this scale.

Below, for comparison, is a map using the same data rendered by MapServer (and GD). It renders small features such as the Maldives rather poorly compared to matplotlib and Agg. Indonesia and the Philippines are messy as well.

Have I mentioned how impressed I am by matplotlib? Swapping out the RendererAgg in my prototype mapping engine for a RendererSVG from matplotlib's SVG backend in no time yields a map as image/svg+xml: nwmed.svg (warning: 500K file). I downloaded an SVG-enabled Firefox CVS build and it displays a nice looking map of the Northwest Mediterranean.

The Cartographic Objects for Zope 0.5 release is now available. There are new features associated with styles and rules, improved spatial reference systems, and support for richer scripting with the Python Cartographic Library. ZCO-0.5.0.tar.gz (sourceforge) requires either the PCL-0.8.1.tar.gz (sourceforge) source release or the PCL-0.8.1-win32.zip (sourceforge) windows binary release.

This morning Kai Hänninen gave me a thorough demo of PrimaGIS, a collaborative mapping application based on ZCO. It's what geospatial Plone users have been looking for: a means to georeference their existing content and provide a mapping interface to their sites. I found it very easy to add a georeferenced document to Kai's demo. See the Python event at http://primagis.technocore.fi/pgdemo/map_view?centerx=-93.21&centery=44.93. If you're a Plone user, take a look and either write Kai for more info or drop in on freenode's #zco channel.

PrimaGIS-0.1.0.tar.gz. Details and demo at primagis.technocore.fi. If you like what you see, send Kai some fan mail and find out how you can get involved.

The preview requires ZCO 0.5 and PCL 0.8, which can be downloaded from links at zmapserver.sourceforge.net.