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 :)
9:57 AM
Post a Comment

Popular posts from this blog

python interval tree

EDIT: added a couple points inline. I'm obsessed with trees lately -- of the CS variety, not the plant variety. Although we are studying poplar, so I'll be using trees to study trees . I'd tried a couple times to implement an interval tree from scratch following the Wikipedia entry . Today I more or less did that in python. It's the simplest possible form. There's no insertion (though that's trivial to add), it just takes a list of 'things' with start and stop attributes and creates a tree with a .find() method. The wikipedia entry that baffled me was about storing 2 copies of each node's intervals --one sorted by start and the other by stop. I didn't do that as I think in most cases it won't improve lookup time. I figure if you have 1 million elements and a tree of depth 16, then you have on average 15 intervals per node (actually fewer since there are the non-leaf nodes). So I just brute force each of those nodes and move to the next. I th...
Read more

levenshtein in cython

EDIT(2): fix markup for <char *> casts ... fix malloc (see comments. thanks Bao). NOTE: using a kwarg for limit slows things down. setting that to a required arg and using calloc for m2 speed things up to nearly as fast as the pylevenshtein. Well, it seems to be popular to code up the levenshtein . I actually have a use for this and wanted to practice some Cython , so I've written a version. I used bearophile's recipe , wikibooks , this (from k4st on reddit) and this for reference. It follows bearophile's code closely, using only O(m) space instead of O(mn). cdef extern from "stdlib.h": ctypedef unsigned int size_t size_t strlen(char *s) void *malloc(size_t size) void free(void *ptr) int strcmp(char *a, char *b) cdef inline size_t imin(int a, int b, int c): if a if c return c return a if c return c return b cpdef int levenshtein(char *a, char *b, int limit=100): cdef int m = strlen(a),...
Read more

synteny mapping

Image
Living in Synteny I've been working on automating synteny mapping between any pairs of genomes. Synteny is where there's a stretch of DNA or genes in some order on chromosomeA of organismX and due to a shared evolutionary history, you can find a similar stretch of genes in order on chromosomeB of organismY. Often there are small losses and inversions, but between closely related organisms like man and mouse , there's still a lot of synteny. Plants can undergo polyploidy , following which, a species can have 2 entire copies of its genome. Over time, much of the duplicated cruft is lost, and the homologous chromosomes diverge, but if the divergence is not too great, it's still possible (actually common) to find synteny within the genome of a single organism--as well as between organisms. I've written my own algorithm to find synteny which uses python sets, and numpy array slicing to do the heavy lifting. It is quite clever [wink]. And it _almost_ works but it'...
Read more