2005 (old posts, page 2)

Python Cartographic Library 0.8 Source Release

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

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

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

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://sgillies.net:8080/PCL/api

Celestial Cartographic Library

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


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://sgillies.net/blog/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]))

# Text symbolizer to label brightest stars
font = mapping.Font(family='Vera', size=8.0)
text_symb = mapping.TextSymbolizer(label='FB_DES', font=font,
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')



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

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


Re: Celestial Cartographic Library

Author: Rob

Hi Sean, You say "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", could you possibly elaborate on this? Why do we multiply by 15? Is there any way we can obtain the altitude of the star relative to the earths surface? Rob

Re: Celestial Cartographic Library

Author: Sean Gillies

I'm projecting the stars onto a celestial sphere, where right ascension is analogous to a terrestrial longitude. Since my GIS software doesn't understand the celestial coordinate system, I have to convert RA to a longitude. Each hour is 15 degrees of arc, that's where the number comes from.

Coming Soon to a Bookstore Near You

Three books related to MapServer and/or Open Source GIS are due to be published this summer. It's fascinating how a community with no paper books can suddenly produce a trio. Almost like how Hollywood can burp up a pair of asteroid/volcano/virus disaster movies simultaneously -- except that these books won't suck.

Tyler Mitchell's Web Mapping Illustrated (oreilly) was the first that I noticed. My understanding is that Tyler is covering everything from data processing and desktop GIS with OpenEV to web applications using MapServer. I proof-read his chapter on mapscript and found it very good. He's not re-publishing the tired old community scripts; these are fresh and clearly coded Python, Perl, PHP, and Ruby examples. I'm definitely getting this one. Should be available in time for OSG05.

Next, i was clued-in to Mapping Hacks (oreilly) by Schuyler Erle, Rich Gibson, and Jo Walsh. This book addresses not the mainstream analyst/developer, but hackers, inventors, artists, and hobbyists. Their interest in scripting MapServer seemed to grow with time and I'm very impressed with the sample code that I found on their companion site: http://mappinghacks.com. Has also been fun watching Schuyler latch onto Python. Or is it Python latching onto Schuyler? Being a fairly conventional GIS type, I think I have more to learn from this book than from Web Mapping Illustrated and am eager to pick it up. Should also be available in time for OSG05.

Lastly, there is Bill Kropla's MapServer: Open Source GIS Development (apress). I know very little about this book or its author, but do know the technical reviewer who assures me that it will be a comprehensive guide to MapServer installation, configuration, and templating. No publication date announced, but probably this summer post-OSG05.

Obligatory Google Entry

The initial release of Google Maps was news worthy. But this? Satellite or aerial imagery layers in web mapping are commonplace. Privacy concerns voiced in the media are overblown, unless you have been trying to hide an illegal Olympic-sized pool on your property. The imagery doesn't zoom in enough to reveal that you are stealing cable or violating your neighborhood's covenant against pets or air-drying laundry. There are real concerns for privacy today, but this is not one.

Graticule Hacking with OGR

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

# Create two fields for labeling
fd = ogr.FieldDefn('VALUE', ogr.OFTInteger)

fd = ogr.FieldDefn('ABS_VALUE', ogr.OFTInteger)

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

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

# destroying data source closes the output file

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

layer = mapscript.layerObj()
layer.status = mapscript.MS_DEFAULT
layer.data = 'graticule.shp'
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()

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.setExpression('([TYPE] EQ 0)')
co.label.color.setRGB(0, 255, 0)

co.setExpression('([TYPE] EQ 1)')
co.label.color.setRGB(255, 0, 0)

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

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')
im = mo.draw()

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


Quick MapServer Extents Calculation

MapServer's scripting module allows users access to the program's projection capabilities. One of these is particularly useful for quick calculation of map extents. Let's say you have a world-wide dataset with a lat/long WGS84 coordinate system, and want to render it onto a Mollweide projection map. Start your interpreter, import the mapscript module, and create output and input projection objects. For the input I use the EPSG shorthand for a lat/long WGS84 coordinate system. Initialize a rectangle to the full extent of the world in units of the input projection, and transform it to output units using its project() method:

>>> import mapscript
>>> proj_out = mapscript.projectionObj('+proj=moll +lon_0=0.0')
>>> proj_in = mapscript.projectionObj('+init=epsg:4326')
>>> e = mapscript.rectObj(-180,-90,180,90)
>>> e.project(proj_in, proj_out)
>>> print e.minx, e.miny, e.maxx, e.maxy
-18040095.6961 -9020047.84807 18040095.6961 9020047.84807

The output values are ready to be used within your map's extent definition.

Python GEOS Module

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

Update (2007-06-12): http://sgillies.net/blog/478/shapely-geometries-for-python/.

Update (2007-04-10): http://sgillies.net/blog/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},
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')

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

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

resulting image, union.png:



Re: Python GEOS Module

Author: Pablo Grigoletti

Is there a version of PyGEOS for linux? I found only a Windows Installer!

Goals for MapServer 4.6

There is no formal release plan yet, but MapServer users can depend on another release around the same time as this summer's OSG05 conference. In anticipation, I recently reviewed the MapServer bugs that I own to adjust priorities and set targets. There is not a lot that I intend to do for the 4.6 release: bug 1061 extends symbolObj::getImage() to vector and truetype font symbols as well as pixmap symbols, and bugs 1128, 1203, 1209, 1266 address makefile and threading issues for the Java mapscript module.

There's a roughly 50% chance that I will also tackle the issue of non-square pixels for mapscript in April or May. This is a feature that is currently only available through the mapserv program's WMS dispatcher, but needed by a few in-development scripts and modules such as the Python Cartographic Library.


Re: Goals for MapServer 4.6

Author: Bob Almodovar

I am in the process of converting an ArcIMS site to MapServer. I have accomplised about 80% of this. The one tool that I have not been able to create is a buffer tool. This tool would need to create a buffer around a feature and query any features that fall within the buffer. It would be nice to have this kind of capability built in, as it is in ArcIMS. If anyone knows how to do this, any help would be appreciated.

MapServer Anti-Patterns

Reviewing recent threads on the mapserver-users list during a chat with a new MapServer user on IRC, I had a modest eureka moment: these particular struggles aren't just isolated cases. These are anti-patterns, patterns of thinking about MapServer that lead to inefficient or broken applications.

Death by a Thousand Classifications is my name for the common overuse of MapServer's feature classification mechanism. The classic case arises when users begin with TIGER lines (perhaps an anti-pattern all its own) and attempt to create separate freeway, highway, and road map layers by classifying on the "CFCC" attribute and throwing away all the features that do not match. The more of these single-class layers the user adds, the more expensive the map. The proper solution to the problem of creating multiple layers from a large dataset is to pre-process the data into separate files or feature tables.

Data Mutation is the anti-pattern that I perceived in this morning's chat. The user wants to have a mutable classification scheme; depending on parameters of a map request some set of counties (in this case) are rendered by the styles of class A and the rest by the styles of class B. The user's solution was to keep the classification rules fixed and instead change the attribute values of data features. Data Mutation is never the right solution to a tricky MapServer problem. First of all, it's unwise to give the MapServer user account write privileges on data. One could work around this by rendering from a temporary copy of the master data, but then there's extra overhead. Lastly, if the Data Mutation anti-pattern is followed, the application will not be able to easily switch to new data sources such as replacing shapefiles on disk with PostGIS or WFS. These features come over the wire and can't be mutated.