mapscript recipe - python mapscript and feature annotation with PostGIS

The situation

Given a set of lines which represent road segments, derived from GPS traces or another source, i want to be able to annotate them with properties and have those properties affect the display in mapserver.

As this is my first foray into python mapscript, i found a few 'tricks and traps', mostly due to the literal nature of the API docs. I found Sean Gillies mapscript API reference invaluable and still look at it daily.

Annotating existing features

I started building this with CGI mapserv, as it provides a convenient wrapper interface to querying points corresponding to a mouse click, this saved a lot of setup time.

In the class stanza of each layer i want to be able to annoate, i add a template reference to my script. The layer also needs a tolerance which is the area (in pixels?) around each queryable object - vital for point query.

    layer
	...
	connectiontype  postgis
        connection      dbname=freemap
        data            "geom from roadseg"
        filter          "type = 'http://openstreetmap.org/schema/Road'"
        class
	    ...
            template "http://uo.space.frot.org/freemap/metamap.cgi?feature=[seg_id]"
	end
	tolerance 10
    end 

The template line passes the value in seg_id in the database table roadseg through to the CGI script.

Now the script can use the seg_id to get the geometry out of the database table, display it using mapscript and use it to make annotations in the db.

        if q.has_key('feature'):

            seg_id = q['feature']
            vars['seg_id'] = seg_id.value
            r = self.get_feature_metadata(seg_id.value)
	    ...

    def get_feature_metadata(self,seg_id):
	db = pgdb.connect( database = "freemap" )
	c = db.cursor()
	c.execute("SELECT name, type, X(StartPoint(geom)) as x, Y(StartPoint(geom)), GeometryFromText(geom) as y, X(EndPoint(geom)), Y(EndPoint(geom)), geom from roadseg where seg_id ="+seg_id)

	""" Returns an arrayref of name, type, x1, y1, x1, y2, geom"""
	return c.fetchone()

What winds up in r is a reference to a seven-member list, the name and type of the geometry (in this case always a line and usually only made up out of two points), its start and end x and y coordinates, and the whole geometry for good measure.

From here it's easy to check the values of HTML form parameters, and send database UPDATEs if the metadata has changed. The main property i'm changing is the road type.

Drawing layers on-the-fly

I want to be able to highlight, on the map, the feature being annotated. I have basically two choices; i can create a layer with a filter, or i can create a shapeObj and insert it into a new layer 'by hand'

Having tried both approaches, i'm still having trouble with labelling objects.

	    map = mapscript.mapObj('bristol.map')

	    # use our mapOb to create a new layer
            layer = mapscript.layerObj(map)
	    layer.type = mapscript.MS_LAYER_LINE
            layer.status = mapscript.MS_ON
            layer.connectiontype = mapscript.MS_POSTGIS 
            layer.connection = 'dbname=freemap'
            layer.data = "geom from roadseg"
            filter = "seg_id = "+seg_id.value
            status = layer.setFilter(filter) 

            # create a new display class using this layer
	    c = mapscript.classObj(layer)
            c.name = 'highlighted'
    
	    # the symbol set is defined in our mapfile
            sset = map.symbolset
            sym = sset.index('solid')

            style = mapscript.styleObj(c)
            style.color.setRGB(255,255,0)
	    style.size = 10
            style.symbol = sym

This is enough to add a new layer with the feature matching the database filter, and colour it a nice yellow so the user can see what they're doing.

The from-scratch approach is a little messier, but is certainly necessary for future applications in which what we're drawing on the map doesn't live in the database (yet).

First, we construct a lineObj from two pointObj objects, and add the line to a wrapper shapeObj. Then we add the shape to the layer that we prepared earlier. The steps for allocating a class and style to the layer are the same as illustrated above.

	    line = mapscript.lineObj()
            p1 = mapscript.pointObj( int(x1), int(y1) )
            p2 = mapscript.pointObj( int(x2), int(y2) )
            line.add(p1)
            line.add(p2)
            shape = mapscript.shapeObj(mapscript.MS_SHAPE_LINE)
            shape.add(line)
                    
	    added = layer.addFeature(shape)

Finding contiguous lines

Initial feedback from the demo suggested that it took a long time to traverse a road of many segments, annotating each one with the same properties, and that on a tiny portion of a tiny map of one city, this quickly became disheartening. So i started looking for ways to be able to suggest, given one line segment, which of those adjacent to it were probably part of the same underlying road.

Recall we got an array r back from our database query, featuring the x and y values, of the start and end points, of each piece of line geometry. Using these to work out the tangent angle of the line, we can compare this with the tangents of the lines that touch it, and if the angles are 'close' enough, suggest that the touching lines should also be annotated with the same metadata properties:

	    contiguous = []
	    [x1,y1,x2,y2] = [r[2],r[3],r[4],r[5]]
	    run = x2 - x1
	    rise = y2 - y1
	    angle = math.atan2(float(run),float(rise))

	    geom = r[6]
	    adjacent = self.get_adjacent_segments(geom)
	
	    for a in adjacent:
                [x1,y1,x2,y2] = [a[2],a[3],a[4],a[5]]
                run = a[4] - a[2]
                rise = a[5] - a[3]

                a_angle = math.atan2(float(rise), float(run))
		
		# the difference in radians between the angles -
		# marginally bigger or smaller 

		diff = angle - a_angle
		if diff < 0: diff = -diff

		# 0.75 radians was an experiential effective value
		if diff < 0.75:

		    # add this to another layer using the approach below 
		    ...
		    # add this to a list of points for the cgi
		    contiguous.append(a[0])

    def get_adjacent_segments(self,geom):
        db = pgdb.connect( database = "freemap" )
	c = db.cursor()

	c.execute("SELECT seg_id, name, X(StartPoint(geom)), Y(StartPoint(geom)), X(EndPoint(geom)), Y(EndPoint(geom)) from roadseg where Touches(GeometryFromText('"+geom+"'),geom)")
	return c.fetchall()

Labelling shapes

Here's an outline of the approach to adding labels on shapes in a layer, which involves starting a class object with a layer, then assigning properties directly to the classObj.label property (which is pre-populated with a labelObj inside).

            class = mapscript.classObj(layer)
            class.label.type = mapscript.MS_TRUETYPE
            class.label.antialias = mapscript.MS_TRUE
            class.label.position = mapscript.MS_AUTO
            class.label.font = 'Vera'
            class.label.color.setRGB(255,255,0)
            class.label.size =  12 

Navigation UI

At this point i was getting pretty frustrated with my inability to program the mapserv CGI interface and making it respond to things that happened in previous transactions. So i started trying to make a pure-mapscript interface to replace it. Unfortunately this involves replacing the neat navigation that mapserv provides.

In my adventures with perl mapscript i figured out how to do a simple 'pan' function with left, right, up, down buttons. Mapserv provides click-to-pan, and since the advent of google maps everyone seems to expect that now.

Looking at cgi-mapserv, i discovered that it's possible to use an image as a form INPUT:

<input type="image" name="img" src="/tmp/image.jpg" width="640" height="480"/>
This returns img.x and img.y parameters to the CGI interface, indicating the X and Y of where has been clicked. The following code is what i am using to derive a 'pan' action from the image click:
    def point_selected(self):
        [imgx,imgy] = [self.q['img.x'].value,self.q['img.y'].value]
        
	x_span = self.vars['maxx'] - self.vars['minx']
        y_span = self.vars['maxy'] - self.vars['miny']
        
	[height, width] = [480,640];

        imgy = height-int(imgy)
        x_ratio = float(imgx)/width
        y_ratio = float(imgy)/height

        x = x_span * x_ratio
        y = y_span * y_ratio

        x = self.vars['minx']+x
        y = self.vars['miny']+y
        z = x_span / 2
        return [x,y,z]

This converts an x-y point click into a UTM grid reference. The 'z' parameter, a misnomer for 'zoom', indicates half the width of the image in metres. To generate an extent for this, to centre the map on the clicked point, later on, simply:

[self.vars['minx'],self.vars['miny'],self.vars['maxx'],self.vars['maxy']] = [x-z,y-z,x+z,y+z]

(self.vars here, i am using as a searchList for a cheetah template.

Now i can convert point clicks on an image, to UTM coordinates which i can perform actions on, my head fills with visions of drawing on maps with python mapscript, and i spend the next three days in the code hole.

Next recipe: Drawing on a PostGIS database with mapscript

gotchas