New ZCO and PCL releases

ZCO 0.5.1 and PCL 0.8.2 releases are available from http://sourceforge.net/project/showfiles.php?group_id=71074. The major new feature (see issue12 in the new issue tracker) will help new users properly configure their data stores.

The source for these projects has been moved to a Subversion repository hosted by Fort Collins' very own tummy.com, who do a lot of other great things for the Linux and Python communities. I've also moved from Sourceforge's issue tracker to a new tracker at http://sgillies.net:8900/primagis/ that is shared with PrimaGIS. It's filling up much faster than the old one.

GDAL-Warping MOLA Topography

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):

http://sgillies.net/images/mars_topo.jpg

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.

Maps of Il Giro and Le Tour

It's too wet to go out for a run or ride after work today, so I'm inside comparing the cartographic qualities of the official maps of the 88th Giro d'Italia and the 92nd Tour de France.

In 2002, the year in which I went to France to see the Tour, the official map changed to an apparently permanent yellow color. Maillot jaune, fine, but I can't stand big blobs of yellow in my browser window. This year's Tour de France has botched the terrain as well, and is stingy about the size of the map.

Wisely choosing not to go with pink, the Giro d'Italia map has natural color hypsometric terrain within Italy and a greyscale shaded relief (re-used from 2004?) without. That and the slightly more grande version (still tiny by my standards) give the 2005 Giro's map the edge at the tape.

Controversy at the Giro today! Check out this series of photos in which Baden Cooke gets chopped down at the finish.

MapServer 4.6 Beta 2

Daniel Morrisette has uploaded mapserver-4.6.0-beta2.tar.gz. Changes since the first beta:

  • add a small buffer around the cliping rectangle to avoid lines around the edges (bug 179).

  • Finished code to convert back and forth between GEOS geometries. Buffer and convex hull operations are exposed in mapscript.

  • fontset.fonts hash now exposed in mapscript (bug 1345).

  • Retreive distance value for DWithin filter request done with line and polygon shapes (bug 1336).

  • Don't render raster layers as classified if none of the classes has an expression set (bugs 985, 1015).

  • Fixed several issues in writing of inline SYMBOLS when saving mapfile (missing quotes around CHARACTER and other string members of SYMBOL object, check for NULLs, and write correct identifiers for POSITION, LINECAP and LINEJOIN) (bug 1344).

Sensor Webs Snag The Times

Nice article yesterday in the NY Times about sensor webs. We need more environmental applications of location technology, and fewer silly authoritarian applications.

By coincidence, I've just finished reading the 21st annual Year's Best Science Fiction anthology. The last story in this collection is Terry Bisson's "Dear Abbey", featuring a planetary environmental sensor web grown omniscient and omnipotent. A fantastic story. Charles Stross's "Rogue Farm" was the only one I enjoyed more.

Putting the Discount WFS to Work

Here's a more interesting use of the cheapo WFS from last week: mining the next generation MapServer website's portal membership for the latitude and longitude properties (just added by Hobu) and publishing them as a WFS.

Here's the layer definition for MapServer users who'd like to try it out:

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

The service root is

http://sgillies.net:9001/mapserver/members/

and the online resource URL for capabilities is

http://sgillies.net:9001/mapserver/members/description.rpy

So far it's just Hobu, Yves, and me:

http://sgillies.net/images/members.png

As other members provide a location, the service will update (every 5 minutes). The picture above will not. I've tried no other clients save uDig, which likes my capabilities, but gags on my feature response. Your mileage may vary.

Draft Program for OSG05

Steve Lime has cut a draft of the OSG05 program, and done a great job as far as I'm concerned. No conflicts for me. I'm scheduled to present in session 11 during the first slot of day 3 along with Schuyler, and some guys who offer "Talking to MapServer with PERL". I hope that's just a typo of Steve's. Are there still people who go around shouting "PERL"?

Discount WFS Source for MapServer

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.

Comments

Re: Discount WFS Source for MapServer

Author: Myunghwa Hwang

Sean, would you mind letting me know about writing a general python cgi script which can serve as features.rpy? This example is what I looked for, except features.rpy. from Myunghwa Hwang

Re: Discount WFS Source for MapServer

Author: Sean

The quick and dirty CGI would be something like this: ... print 'Content-type: text/html\n\n' print '' print csvstore.CSVFeatureCollection().toxml()