Sean Gillies Blog

Guido van Rossum is thinking out loud about saying bye bye to lambda, reduce(), filter(), and map(). No comment on the first two -- but filter and map are pretty much obsolete thanks to generator expressions and conflict with the most natural names for two objects of the GIS problem domain. This proposed change won't guarantee that we never again see myMap in Python mapscript code, but it would be a start.

Fort Collins' own Linux and Python advocate Sean Reifschneider blogs PyCon 2005. Interesting photos and fun commentary. My suggestions for code sprints at the MapServer Users Meetings has become a running joke among MapServer developers, but these photos have inspired me to run the old flag up the pole one more time. Past two years, the MapServer project has scheduled releases to coincide with the meeting; this leaves most of the developers in no mood for programming. Any sprinting would probably have to be on some project other than MapServer itself.

update: Hobu just informed me that he proposed sprinting at the OSG05 organizing meeting and it was shot down immediately.

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.

How do I get my data into MapServer? Often, new MapServer users have a large legacy dataset that is not in a supported format. This might be a CSV file, a MySQL database, or a proprietary binary format. An overlooked solution to this problem involves WFS (Web Feature Service), the closest thing MapServer has to a generic data source. A Web Feature Service responds to requests for features with a XML document that contains features and attributes encoded using GML, the Geography Markup Language. There are a number of full-featured WFS, commercial and open source, but it is remarkably easy to roll your own that can serve as a lightweight MapServer data source.

Let's look at the simple case of publishing via WFS the data in a test CSV file. Granted, this is not the most compelling case, but is easy to reproduce and a good starting point. Consider a file with a header row and six point feature rows:

"FID:int","FNAME:string","SCORE:float","X:float","Y:float"
1,"Point A",0.2,-110.5,44.3
2,"Point B",0.4,-106.7,37.9
3,"Point C",0.15,-109.5,41.6
4,"Point D",0.6,-108.1,42.5
5,"Point E",0.8,-107.4,39.8
6,"Point F",0.55,-108.6,38.2

Each feature has an integer ID, a name, a float type score attribute, and, as is often the case with legacy or non-GIS data, coordinates as separate x and y values.

I will render this into XML using ElementTree, a flexible module for programming hierarchical data structures. The reading of data, and creation and arrangement of elements will be handled by a CSVFeatureCollection class. A toxml() method will return a WFS feature collection encoded as an XML string, and be used as the core of a response to a GetFeature request. To keep the example as simple as possible, I will ignore all WFS request parameters; this cheap-o WFS returns every feature from the data source. Because MapServer is a simple and forgiving WFS client, I can also neglect the GetCapabilities and DescribeFeatureType responses required of a full-featured WFS.

Since the WFS response makes much use of qualified names for element tags, I first define two functions for convenience:

# begin module csvsource.py

import csv
import elementtree
from elementtree.ElementTree import Element, QName, tostring

# Convenience wrappers for the Element factory
# --------------------------------------------

def WFSElement(tag):
    e = Element(tag)
    e.tag = QName('http://www.opengis.net/wfs', e.tag)
    return e

def GMLElement(tag):
    e = Element(tag)
    e.tag = QName('http://www.opengis.net/gml', e.tag)
    return e

Create the root element, defining the required WFS and GML namespaces, and creating the element which describes the bounding box for features in the response. One might use mapscript or the OGR Python module to compute the proper extents. In the interest of simplicity, I'm taking a shortcut and hard-coding the bounding box:

class CSVFeatureCollection:

    """Express the 'points.csv' file as a WFS Feature Collection"""

    def toxml(self):
        """Return an XML string, sans prolog"""
        # root element
        e_root = WFSElement('FeatureCollection')
        e_root.attrib['xmlns'] = 'http://localhost/myns'
        e_root.attrib[QName('http://localhost/myns', 'wfs')] \
                   = 'http://www.opengis.net/wfs'
        e_root.attrib[QName('http://localhost/myns', 'gml')] \
                   = 'http://www.opengis.net/gml'

        # bounds
        e_bounds = GMLElement('boundedBy')
        e_box = GMLElement('Box')
        e_box.attrib['srsName'] = 'EPSG:4326'
        e_coords = GMLElement('coordinates')
        e_coords.text = '-120,30 -100,50'
        e_box.append(e_coords)
        e_bounds.append(e_box)
        e_root.append(e_bounds)

Next, open the data file and iterate over its rows, producing and attaching new elements, and then returning the XML string:

# continued ...

        # open data file and iterate over rows
        f = file('points.csv', 'r')
        reader = csv.reader(f)
        headers = reader.next()
        self._fields = [(h.split(':')) for h in headers]

        for row in reader:
            if not row:
                break

            e_memb = GMLElement('featureMember')
            e_type = Element('points')

            # member bounds
            e_bounds = GMLElement('boundedBy')
            e_box = GMLElement('Box')
            e_box.attrib['srsName'] = 'EPSG:4326'
            e_coords = GMLElement('coordinates')
            e_coords.text = '%s,%s %s,%s' % (row[3], row[4], row[3], row[4])
            e_box.append(e_coords)
            e_bounds.append(e_box)
            e_type.append(e_bounds)

            # member properties
            e_ptprop = GMLElement('pointProperty')
            e_point = GMLElement('Point')
            e_point.attrib['srsName'] = 'EPSG:4326'
            e_coords = GMLElement('coordinates')
            e_coords.text = '%s,%s' % (row[3], row[4])
            e_point.append(e_coords)
            e_ptprop.append(e_point)
            e_type.append(e_ptprop)

            # other props
            for i in range(3):
                e = Element(self._fields[i][0])
                e.text = row[i]
                e_type.append(e)

            e_memb.append(e_type)
            e_root.append(e_memb)

        f.close()
        return tostring(e_root)

This module could be adapted to a number of different data sources. Instead of reading from a CSV file, one might make a database query, or fetch and parse RSS files.

The result of executing:

import csvsource
csvsource.CSVFeatureCollection().toxml()

is:

<?xml version="1.0" ?>
<ns0:FeatureCollection xmlns="http://localhost/myns"
    ns1:gml="http://www.opengis.net/gml"
    ns1:wfs="http://www.opengis.net/wfs"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:ns1="http://localhost/myns">
  <ns2:boundedBy xmlns:gml="http://www.opengis.net/gml">
    <ns2:Box srsName="EPSG:4326">
      <ns2:coordinates>-120,30 -100,50</ns2:coordinates>
    </ns2:Box>
  </ns2:boundedBy>
  <ns2:featureMember xmlns:gml="http://www.opengis.net/gml">
    <points>
      <ns2:boundedBy>
        <ns2:Box srsName="EPSG:4326">
          <ns2:coordinates>-110.5,44.3 -110.5,44.3</ns2:coordinates>
        </ns2:Box>
      </ns2:boundedBy>
      <ns2:pointProperty>
        <ns2:Point srsName="EPSG:4326">
          <ns2:coordinates>-110.5,44.3</ns2:coordinates>
        </ns2:Point>
        </ns2:pointProperty>
      <FID>1</FID>
      <FNAME>Point A</FNAME>
      <SCORE>0.2</SCORE>
    </points>
  </ns2:featureMember>
  ...
</ns0:FeatureCollection>

This is close, but it is not yet going to work with MapServer. The XML Namespace specification states that applications should use the element namespace URL not the tag prefix, but this detail is handled improperly. MapServer looks for tags with 'wfs' and 'gml' prefixes as it parses the response and does not recognize the prefixes generated by our root element's writer.

The solution, easy to do with Python, is to monkey-patch the ElementTree module. I simply add some namespace and prefix mappings to the module's well-known namespaces immediately after it is imported at the top of csvstore.py:

import csv
import elementtree
from elementtree.ElementTree import Element, QName, tostring

# Monkey Patch adds to the default well known namespaces
elementtree.ElementTree._namespace_map.update({
    'http://www.opengis.net/wfs': 'wfs',
    'http://www.opengis.net/gml': 'gml'})

This gives us what MapServer expects:

<?xml version="1.0" ?>
<wfs:FeatureCollection xmlns="http://localhost/myns"
    ns1:gml="http://www.opengis.net/gml"
    ns1:wfs="http://www.opengis.net/wfs"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:ns1="http://localhost/myns">
  <gml:boundedBy xmlns:gml="http://www.opengis.net/gml">
    <gml:Box srsName="EPSG:4326">
      <gml:coordinates>-120,30 -100,50</gml:coordinates>
    </gml:Box>
  </gml:boundedBy>
  <gml:featureMember xmlns:gml="http://www.opengis.net/gml">
    <points>
      <gml:boundedBy>
        <gml:Box srsName="EPSG:4326">
          <gml:coordinates>-110.5,44.3 -110.5,44.3</gml:coordinates>
        </gml:Box>
      </gml:boundedBy>
      <gml:pointProperty>
        <gml:Point srsName="EPSG:4326">
          <gml:coordinates>-110.5,44.3</gml:coordinates>
        </gml:Point>
      </gml:pointProperty>
      <FID>1</FID>
      <FNAME>Point A</FNAME>
      <SCORE>0.2</SCORE>
    </points>
  </gml:featureMember>
...
</wfs:FeatureCollection>

All that's left is a means to serve this data. One might use a Python CGI script and Apache, or use the code above in a Zope external method that is served to the web via Python Script. Instead, for extra flavor, I'm going to use twisted.web and write a Python resource script to answer GetFeatures requests:

# begin features.rpy

from twisted.web import server
from twisted.web.resource import Resource
from twisted.internet import reactor, defer

import csvstore

class Features(Resource):

    def get_data(self):
        d = defer.Deferred()
        reactor.callLater(0.25, d.callback,
                          csvstore.CSVFeatureCollection().toxml())
        return d

    def render_GET(self, request):
        """a default method called when the object is published"""
        request.setHeader('Content-type', 'text/xml')
        request.write('<?xml version="1.0" ?>')
        d = self.get_data()
        d.addCallback(lambda x:(request.write(x+'\n'), request.finish()))
        d.addErrback(lambda x:(request.write(x), request.finish()))
        return server.NOT_DONE_YET

    render_POST = render_GET

resource = Features()

The entire WFS and data consists of three files: points.csv, csvstore.py, and features.rpy, in a directory such as '/var/www/cheapowfs'. It can be put online using Twisted's mktap and twistd utilities:

mktap web -p 8080 --path /var/www/cheapowfs
twistd -f /var/www/cheapowfs/web.tap

All files: cheapowfs.tar.gz

This WFS can be accessed by MapServer using a layer definition like the one below:

LAYER
   NAME "cheapo_wfs_points"
   TYPE POINT
   STATUS DEFAULT
   CONNECTIONTYPE WFS
   CONNECTION "http://localhost:8080/features.rpy?"
   METADATA
     "wfs_service" "WFS"
     "wfs_typename" "points"
     "wfs_version" "1.0.0"
   END
   LABELITEM "FNAME"
   PROJECTION "+init=epsg:4326" END
END

Hobu and I will be expanding on this example in our Python workshop this summer, and I'll be back to it with a few more interesting examples in the next few weeks.

In response to a post on the Matplotlib-users list I almost automatically responded that, of course, we can use GDAL for that; GDAL can read anything. Being Friday the 13th, I decided that maybe it should actually be tried first, so I downloaded the coarse 4 pixels per degree MOLA (Mars Orbiter Laser Altimeter) topography raster from http://pds-geosciences.wustl.edu/missions/mgs/megdr.html. This is not a format directly supported by GDAL, but is a raw binary file with a header descriptive enough that I can wrap it with a GDAL virtual dataset:

<VRTDataset rasterXSize="1440" rasterYSize="720">
  <GeoTransform>-180.0,0.25,0.0,90.0,0.0,-0.25</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTRawRasterBand">
    <SourceFilename relativetoVRT="1">megt90n000cb.img</SourceFilename>
    <ImageOffset>0</ImageOffset>
    <PixelOffset>2</PixelOffset>
    <LineOffset>2880</LineOffset>
    <ByteOrder>MSB</ByteOrder>
  </VRTRasterBand>
 </VRTDataset>

and then warp it to an orthographic projection using the gdalwarp utility:

$ gdalwarp -tr 21875 21875 -te -7000000 -7000000 7000000 7000000 \
  -wo SOURCE_EXTRA=100 -wo SAMPLE_GRID=YES \
  -s_srs "+proj=latlong +ellps=sphere" \
  -t_srs "+proj=ortho +lon_0=0.0 +lat_0=0.0 +ellps=sphere" \
  -dstnodata -32700 topo.vrt topo_ortho.tif

The result (converted to a JPEG):

The functionality of gdalwarp is available through GDAL's gdal.py module as well, which means that we should be able to hook matplotlib's pylab up to GDAL, and thereby to almost any raster data source.

The April-June 2005 issue of ArcUser contains an article by Hobu on pages 34-36. It's an introduction to Python and GIS, and vector, grid, and raster processing with ogr.py and gdal.py for users of ArcGIS 9 -- which includes Python and makes a small subset of its functionality scriptable. Good writing, too bad he seems to have retired from blogging ;)

I ventured further into the magaine and found myself comparing the agility of ArcGIS and the GDAL/OGR Python modules while reading the tutorial on pages 46-49 ("Creating Accurate Footprints of Registered Aerial Images"). It's a rather involved and round-about manual process requiring an expensive ArcGIS add-on (Spatial Analyst). I'm sure that one could also manage to open a beer bottle using dentistry tools, but not as readily as with a bottle opener. Consider the following Python script:

# Create footprint shapefile
import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
footprints_shp = driver.CreateDataSource('footprints.shp')
footprints = footprints_shp.CreateLayer('footprints',
                 geom_type=ogr.wkbPolygon)
fd = ogr.FieldDefn('FILENAME', ogr.OFTString)
fd.SetWidth(30)
footprints.CreateField(fd)

# Loop over a number of georeferenced images
import glob
import gdal

files = glob.glob('*.jpg')
for file in files:
    # Get georeferencing and size of imagery
    dataset = gdal.Open(file)
    g = dataset.GetGeoTransform()
    pixels = dataset.RasterXSize
    lines = dataset.RasterYSize
    minx = g[0]
    maxx = minx + pixels * g[1]
    maxy = g[3]
    miny = maxy + lines * g[5]

    # append to the 'footprints' layer
    wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))' \
        % (minx, miny, minx, maxy, maxx, maxy, maxx, miny, minx, miny)
    g = ogr.CreateGeometryFromWkt(wkt)
    f = ogr.Feature(feature_def=footprints.GetLayerDefn())
    f.SetField(0, file)
    f.SetGeometryDirectly(g)
    footprints.CreateFeature(f)
    f.Destroy()

# destroy footprints_shp to flush and close
footprints_shp.Destroy()

No fooling about -- directly to the imagery bounds and write a polygon. Simple and straight-forward, the right tool for the job.

Howard has put online our workshop materials. I may be recycling some of these at GIS in the Rockies.

It's been a while, but I have finally made the time to upgrade my Cheapo WFS so that it can feed [screenshot] uDig, Refractions' new desktop GIS app.

I found the oXygen XML editor to be a handy program for finding bugs in my capabilities document. I'm only demo-ing oXygen now but might spring for it if it's equally useful with WMS and SLD. Getting my data into uDig was mostly a matter of tweaking the WFS response and schema into a form that GeoTools would accept. Evidentally, members is not a good typename no matter the namespace.

The service root is

http://zcologia.com:9001/mapserver/members/

and the online resource URL for capabilities, which you may enter in uDig's WFS data source wizard, is

http://zcologia.com:9001/mapserver/members/capabilities.rpy?request=GetCapabilities

MapServer users will need to adjust to the change in the typename, using the following layer definition:

LAYER
   NAME "members"
   TYPE POINT
   STATUS DEFAULT
   CONNECTIONTYPE WFS
   CONNECTION "http://zcologia.com:9001/mapserver/members/features.rpy?"
   METADATA
     "wfs_service" "WFS"
     "wfs_typename" "users"
     "wfs_version" "1.0.0"
   END
   LABELITEM "zco:fullname"
   PROJECTION "+init=epsg:4326" END
END

In the wake of Google's mapping API release, ESRI has opened up some of its SOAP-based web services for limited public use to see if they can lure in the map hackers. I played around with a trial subscription to ArcWeb services two years ago, and this seemed a good opportunity to dust off the old code and try again with the public services.

As before I'm using ZSI, the Zolera SOAP Infrastructure for Python. Happily, it's improved, and now features generation of service stubs and types from WSDL. Two years ago I was faced with implementing my own ArcWeb types. Now I can use ZSI's wsdl2py.py script to generate type and client stub modules from the ESRI authentication and map image WSDL URLs, which allows me to abandon my earlier code and start fresh:

$ wsdl2py.py -u https://arcweb.esri.com/services/v2/Authentication.wsdl
$ wsdl2py.py -u https://arcweb.esri.com/services/v2/MapImage.wsdl

resulting in the creation of three files:

$ ls *.py
Authentication_services.py  MapImage_services.py
MapImage_services_types.py

Following the example, I quickly coded up a request for an authenticated token using my existing ESRI global account:

import sys

# modules generated by wsdl2py
from Authentication_services import *

opts = {'tracefile': sys.stdout}

# Auth service
auth_loc = AuthenticationLocator()
auth_serv = auth_loc.getIAuthentication(**opts)

# Get a token
tok_req = getToken1InWrapper()
tok_req._username = '******'
tok_req._password = '******'
tok_resp = auth_serv.getToken(tok_req)
token_ustr = tok_resp._Result

and was informed that I could not use the public services with this account because I had an expired non-public subscription. Annoying, this was, and seems to support James Fee's contention that the public service release caught ESRI itself off-guard. So I marched through the ESRI account drill again, and continued on. The following code gets a token and uses it to make the simplest possible request for a map image and print the resulting map image URL:

## awmap.py

import sys

# modules generated by wsdl2py
from Authentication_services import *
from MapImage_services import *

opts = {'tracefile': sys.stdout}

def set_map_extent(o, e):
    o.Set_minx(e['minx'])
    o.Set_miny(e['miny'])
    o.Set_maxx(e['maxx'])
    o.Set_maxy(e['maxy'])

def set_map_size(o, s):
    o.Set_width(s[0])
    o.Set_height(s[1])

# extent and datasource
extent = {'minx': -105.14, 'miny': 40.53, 'maxx': -105.04, 'maxy': 40.63}
size = [400, 400]
datasource = 'GDT.Streets.US'

def main():
    # Auth service
    auth_loc = AuthenticationLocator()
    auth_serv = auth_loc.getIAuthentication()

    # get a token
    tok_req = getToken1InWrapper()
    tok_req._username = '******'
    tok_req._password = '******'
    tok_resp = auth_serv.getToken(tok_req)
    token_ustr = tok_resp._Result

    # map service implementation
    map_loc = MapImageLocator()
    map_serv = map_loc.getIMapImage()

    # get map image url
    map_req = getMap0InWrapper()
    map_req._token = token_ustr
    map_req._mapExtent = ns4.Envelope_Def()
    set_map_extent(map_req._mapExtent, extent)
    map_req._mapImageOptions = ns5.MapImageOptions_Def()
    map_req._mapImageOptions._dataSource = datasource
    map_resp = map_serv.getMap(map_req)

    print map_resp._Result._mapUrl

if __name__ == '__main__':
    main()

Running the script:

$ [sean@lenny arcweb]$ python awmap.py
http://redlandsarcweb.esri.com/out/maps/GDT_ArcWeb_US_2_blackmap2248241332.png

Produced the following image at the URL:

In my opinion, ESRI needs to make some cosmetic improvements to its default map if it is going to catch on with web designers and programmers. I know for a fact that the software behind the service allows for great looking maps, so let's see some fine cartography in this default map. This is the first look that many people will have of the public ESRI maps, so it must look sharp.

This exercise increases my fondness for SOAP very little, but ZSI makes using services such as these a fairly painless task for an intermediate Python programmer.

Next month I'll be at GIS in the Rockies, conducting a presentation/workshop entitled Data Processing and Mapping with Python: Open Source and Open Standards in the conference's Emerging Standards track. I'm pleased to have been invited to talk about the projects that we're working on in the Open Source GIS community, and show off some specific Python applications. Contents include:

• Overview of open source GIS projects.
• Vector and raster processing with GDAL, GEOS, and Python.
• Introduction to PostGIS and the Python DB API.
• Providing and consuming web map and feature services using MapServer and the Python mapscript module.

This won't be hands-on like the workshop Hobu and I ran at OSG05, but will be similarly aimed at programmer/analysts and their managers, with many scripting examples that can be brought back for immediate use in the shop.

While not strictly open, ArcWeb services are technology that some attendees will want to see covered, and so I'll work my previous ZSI and ArcWeb exercise into the mix. Using Python for ArcGIS geoprocessing is, however, outside the scope of this presentation.

One problem that new Python MapScript users immediately run into is the choice of web framework; Python is blessed (or afflicted) with a multitude. I think these users owe it to themselves to check out what Adrian Holovaty has been doing on chicagocrime.org using Django.

I'm still up to my ears in Zope 2 projects through the end of this summer, but will be trying out Django as soon as possible.

Frank Warmerdam's Geospatial Data Abstraction Library, GDAL, is a big part of the foundation of our open source GIS software stack, and the 1.3.0 release is now available for download.

This release corrects some GML parsing and spatial reference issues that I found in version 1.2.6. I'm only using a subset of GDAL's formats; there are improvements across the board. Hobu and Kevin Ruland have revised the GDAL and OGR scripting framework to make it possible to produce Java and C# bindings using SWIG. Perhaps even PHP, if SWIG support for that language improves.

Unfortunately for me, Python was the first scripting language for GDAL and its interface has a distinctly un-Pythonic feel. New language users have nothing but blue skies and can work with SWIG to make APIs with the right feel. Java users, for example, really want methods with a getFoo form rather than GetFoo for sake of consistency and clarity. This can be accomplished using SWIG's %rename directive:

%rename(getLayer) DataSource::GetLayer;

Perhaps it seems like a lot of work to some, but I guarantee that some Java user will make the time to do this. Hell, I would if I could, but there are too many legacy Python scripts out there to allow any such changes for the Python bindings. While I'm at it, I might as well wish for a pony ;)

I haven't looked at SWIG's online docs since the 1.3.23 release: when did they get so purty?

Through Planet Python I just saw that ESRI is hiring a Python programmer to, among other things

improve the Python language as required by ESRI and its thousands of users

As opposed to the mere dozens of existing Python users? ;)

I'm really curious about what improvements are required. Hopefully this would mean improvements to the standard library rather than big changes to the core language, though as Python is used as a scripting language for geo-processing chains they may be looking at coroutines, continuations, and closure.

At any rate, here's hoping that they hire somebody who really understands Python and open source. Like Hobu :)