Python Mapscript Tricks

I mentioned previously a site which uses google maps with a WMS. The problem with using points (or labels) in a tiled application such as google maps, or OpenLayers is that each tile can only draw its own contents. So if you draw a point with a radius of 10 pixels whose center is 2 px away from the edge of the tile, then 8px of the entire 20px will be chopped. By default, those lost 8 pixels will not be drawn in the adjacent tile because the center of the point does not fall in that tile. so that gives something that looks like this for 2 adjacent tiles:



It's hard to even tell that those tiles belong together! Using some python mapscript and PIL, it's pretty simple to make those look like this( except I dont know how to tell blogger not to add the spacin...) :



That's actually the same WMS request(s), just changing the call from a simple WMS CGI script that calls mapserver to a WSGI script that uses python mapscript and PIL.
The script:
1. takes the current bounding box, width and height, and expands them all proportionally.
2. has mapscript draw the expanded image which is
3. saved into a stringIO object which is
4. sent that to PIL to be cropped down to size and extent.
5. returned to the browser.


#!/usr/bin/python

import mapscript
from cgi import parse_qsl
import os
import Image
from cStringIO import StringIO

BUFFER_PCT = 0.1
MAPFILE = 'soil.map'

########################################################
MAPFILE = os.path.join(os.path.dirname(__file__), MAPFILE)

def application(environ, start_response):
winput = dict(parse_qsl(environ['QUERY_STRING']))
format = winput.get('FORMAT', 'image/png')
start_response('200 OK', [('Content-Type', format)])
extension = format[format.find('/') + 1:]


wms = mapscript.mapObj(MAPFILE)
req = mapscript.OWSRequest()

for k,v in winput.items():
req.setParameter(k, v)

buffer = BUFFER_PCT/2

bbox = map(float, winput['BBOX'].split(","))
rangex = bbox[2] - bbox[0]
rangey = bbox[3] - bbox[1]
w = int(winput['WIDTH'])
h = int(winput['HEIGHT'])
xdelta = int(round(w * buffer)) # e.g add 13px on each side
ydelta = int(round(h * buffer)) # for 256x256 image

bbox[0] -= rangex * buffer # extend the
bbox[1] -= rangey * buffer # bbox in
bbox[2] += rangex * buffer # all
bbox[3] += rangey * buffer # directions

# http://trac.osgeo.org/mapserver/ticket/2299
req.setParameter('WIDTH', str(w + 2 * xdelta)) # and adjust the
req.setParameter('HEIGHT', str(h + 2 * ydelta)) # h, w by the same
req.setParameter('BBOX', ",".join(map(str,bbox))) # amount
req.setParameter("STYLES", "")
req.setParameter("REQUEST", "GetMap")

# PIL doesnt like interlace.
wms.outputformat.setOption('FORMATOPTIONS', 'INTERLACE=OFF')
wms.loadOWSParameters(req)

im = Image.open(StringIO(wms.draw().getBytes()))
if im is None: return ['']

# crop the image back to the requested w, h
im = im.crop((xdelta, ydelta, w + xdelta, h + ydelta))

buffer = StringIO()
im.save(buffer, extension)
buffer.seek(0)
data = buffer.read()
return [ data ]


That's pretty simple, but pretty cool. And it actually runs slightly faster than the CGI version. It runs even more quickly if "wms = mapscript.mapObj(MAPFILE)" is moved into the global scope, so the mapfile doesnt have to be re-parsed each time, but I didnt test enough to make sure that's ok.

Comments

Megan D said…
Wow that is so perfect. And exactly the same "mashcache.wsgi" file for our michigan site! good work friend, and now i understand it more :)

Popular posts from this blog

filtering paired end reads (high throughput sequencing)

python interval tree

needleman-wunsch global sequence alignment -- updates and optimizations to nwalign