Linfiniti Geo Blog

GIS for Open Source People

Browsing Posts in General FOSSGIS

Last week I had to cope with one of the (few!) limitations of Django ORM model about inheritance. For the GeoDjango website we are developing, we created a model that stores user-layer pairs, so that a logged user can restore the state of map and legend as [s]he left it at previous login, and anonymous users receive a default list of layers. Here is the code for the model:

class UserWmsLayer( models.Model ):
  """ Stores custom user preferences for a layer. """
  wmslayer = models.ForeignKey( WmsLayer )
  user = models.ForeignKey( User )
  is_visible = models.NullBooleanField( null=True, blank=True )
  is_deleted =  models.NullBooleanField( null=True, blank=True, default=False )

The Layers in Legend can be of WmsLayer class or any of its subclasses, at the moment the only subclass is DateQueryLayer, a WMS layer with a filter on the date attribute. We also have a Layer class, that is an ABC (Abstract Base Class), but it’s not possible to create a ForeignKey to an ABC as it has no table in the database. So we decided to create the ForeignKey to WmsLayer as a working, half-hack solution.

That worked nicely until I relied on inheritance and method overriding. I expected Django to store the ForeignKey to the actual class of the object, but the ORM stores it as a reference to the superclass, in our case of WmsLayer. This is very well explained in enlightening answer #4 of this stackoverflow QA.

In few words, I can actually have a User-DateQueryLayer pair, but in the model it will be stored as WmsLayer, and all methods I call will be WmsLayer‘s, not DateQueryLayer‘s as I would expect from OO programming. DateQueryLayer redefines asOpenLayer method, but due to this ORM limitation I could never call it.

The solutions were multiple, more or less Pythonic, Django-ish, brittle and scalable. I took two days to google and get through them, and the only two applicable solutions seemed to be:

  1. Run Time Type identification via a ifelse cascade: that would be the less scalable and most brittle hack.
  2. Django’s Generic Relations: better, but still overcomplicated.

I was about to implement Generic Relations when I realised that asOpenLayers() simply returns a JavaScript string. A simple solution came to me in a flash. Instead of creating the JavaScript code on the fly every time the layer is used somewhere, I decided to store that string into an attribute owed by Layer, the abstract superclass, so that all subclasses will simply inherit it and populate it at instantiation time using their own asOpenLayers method. The User-*Layer pair will therefore rely on an attribute in WmsLayer table without any inheritance issue.

Here is what the relevant code looks like:

class Layer( models.Model ):
  """This is an ABC (Abstract Base Class) for all models that are layers.
  It provides common api so that all layers can be treated in a similar way."""
  name = models.CharField( max_length=256 )
  owner = models.ForeignKey( User, related_name = 'owner' )
  # store the javascript as *attribute*, instead of generating it on-the-fly
  as_open_layer = models.TextField( )

class WmsLayer( Layer ):
  url = models.URLField( max_length=1024, verify_exists=True )
  layers = models.CharField( max_length=256 )
  # link to user-layer model to keep status of legend
  users = models.ManyToManyField(User, through='UserWmsLayer')

  def asOpenLayer( self ):
    """Return a string representation of this model as an open layers
    layer definition. The created layer def will be added
    to the openlayers map of name theMap (which defaults to "map". """
    return "the JavaScript string" #omitted for brevity

  def save( self, *args, **kwargs ):
    """ Overrides standard save, generating the asOpenLayers javascript
    and storing it into as_open_layer attribute. """
    self.as_open_layer = self.asOpenLayer()
    super(WmsLayer, self).save( *args, **kwargs)

class DateQueryLayer( WmsLayer ):
  """A layer model for storing user date range queries persistently"""
  date_query_type = models.ForeignKey( DateQueryType )
  sensor = models.ForeignKey( Sensor )
  start_date = models.DateTimeField( null=True, blank=True )
  end_date = models.DateTimeField( null=True, blank=True )

  def asOpenLayer( self ):
  """ Overrides WmsLayer's method. """
  return "the JavaScript string" #omitted for brevity

  def save( self, *args, **kwargs ):
    self.as_open_layer = self.asOpenLayer()
    super(DateQueryLayer, self).save( *args, **kwargs)

That’s it.

Once you filter the UserWmsLayer table by user, you won’t call anymore the asOpenLayer method…

myObjects = UserWmsLayer.objects.filter(user__username = "anonymous")
for myObject in myObjects:
  myObjects.asOpenLayer()

but fetch the as_open_layer attribute instead!

myObjects = UserWmsLayer.objects.filter(user__username = "anonymous")
for myObject in myObjects:
  myObjects.as_open_layer

This allows you to pick up the correct string, without worrying about which type myObject is. This is good. Clean. Scalable. Pythonic. Yes, the asOpenLayer method simply returns a string according to some Layer properties, this solution could not work with more complex methods. But if you face a similar inheritance issue, and your method returns something that can be stored in an attribute of the superclass, that’s one of the cleanest solutions I’m aware of.

Hope this helps, and better solutions are most welcome :)
Happy coding,
–anne

Are you looking for giving your OpenLayers map controls a cool appearance, smoothly integrated with the site’s theme, without writing a papyrus and scatter code among lot of files?

Then have a look at jQuery UI CSS framework, a system of classes developed for jQuery UI widgets.

This is the map toolbar of the webGIS site we are busy developing, rendered with UI-Darkness theme:
Fancy OpenLayers control toolbar

The controls (pan, measure and zoom) are OpenLayers’ controls. They are all created in the map’s init() (see first js snippet below). The binding with the buttons is made by name – therefore be sure that the names of the OpenLayers controls match exactly the name properties of the respective buttons. The activation of the selected control is done by the toggleControl() function, further below in the js snippet. That way you can add as many control-button pairs as you need.

Let’s see what the code looks like. It is not so much indeed.. My tribute to the proverbial programmer’s laziness and to the koan of Master Foo’s and the Ten Thousand Lines.

First be sure to include all necessary scripts in the html head:

<link type="text/css" href="/static/css/jquery.fancybox-1.2.6.css" rel="stylesheet" media="screen" />
<script type="text/javascript" src="/static/js/jquery-1.3.2.min.js"></script>
<script type="text/javascript" src="/static/js/jquery-ui-1.7.2.custom.min.js"></script>
<script type="text/javascript" src="/static/js/jquery.fancybox-1.2.6.pack.js"></script>
<script type="text/javascript" src="/static/js/jquery.easing.1.3.js"></script>

then add the buttons – I took inspiration from Filament Group excellent post:

<div id="mapcontrols" class="fg-buttonset fg-buttonset-single ui-helper-clearfix">
<button name='navigate'class="fg-button ui-state-default ui-state-active ui-priority-primary ui-corner-left" >Navigate</button>
<button name='line' class="fg-button ui-state-default ui-priority-primary">Measure line</button>
<button name='polygon' class="fg-button ui-state-default ui-priority-primary">Measure area</button>
<a href="#" name='zoomin' class="fg-button ui-state-default fg-button-icon-solo" title="Zoom in"><span class="ui-icon ui-icon-circle-zoomin"></span> Zoom in</a>
<a href="#" name='zoomout' class="fg-button ui-state-default fg-button-icon-solo ui-corner-right" title="Zoom out"><span class="ui-icon ui-icon-circle-zoomout"></span> Zoom out</a>
</div>

And in the end the JavaScript – add the following snippet in the map’s init() function:

    mapControls = {
	line: new OpenLayers.Control.Measure(
	    OpenLayers.Handler.Path, {
		persist: true
	    }
	),
	polygon: new OpenLayers.Control.Measure(
	    OpenLayers.Handler.Polygon, {
		persist: true
	    }
	),
	zoomin: new OpenLayers.Control.ZoomBox(
	    {title:"Zoom in box", out: false}
	),
	zoomout: new OpenLayers.Control.ZoomBox(
	    {title:"Zoom out box", out: true}
	)
    };

    var control;
    for(var key in mapControls) {
	control = mapControls[key];
	control.events.on({
	    "measure": handleMeasurements,
	    "measurepartial": handleMeasurements
	});
	map.addControl(control);
    }

and these functions at the bottom of your js file:

function handleMeasurements(event) {
    var geometry = event.geometry;
    var units = event.units;
    var order = event.order;
    var measure = event.measure;
    var element = document.getElementById('output'); //TODO redirect to other area?
    var out = "";
    if(order == 1) {
	out += "Measure: " + measure.toFixed(3) + " " + units;
    } else {
	out += "Measure: " + measure.toFixed(3) + " " + units + "2";
    }
    element.innerHTML = out;
}

function toggleControl(element) {
    for(key in mapControls) {
	var control = mapControls[key];
	//alert ($(element).is('.ui-state-active'));
	if(element.name == key && $(element).is('.ui-state-active')) {
	    control.activate();
	} else {
	    control.deactivate();
	}
    }
}

$(function(){
    //all hover and click logic for buttons
    $(".fg-button:not(.ui-state-disabled)")
    .hover(
	function(){
	    $(this).addClass("ui-state-hover");
	},
	function(){
	    $(this).removeClass("ui-state-hover");
	}
    )
    .mousedown(function(){
	$(this).parents('.fg-buttonset-single:first').find\
            (".fg-button.ui-state-active").removeClass("ui-state-active");
	if( $(this).is('.ui-state-active.fg-button-toggleable, \
            .fg-buttonset-multi .ui-state-active') )
	    { $(this).removeClass("ui-state-active"); }
	else { $(this).addClass("ui-state-active"); }
    })
    .mouseup(function(){
	if(! $(this).is('.fg-button-toggleable, .fg-buttonset-single .fg-button,  \
            .fg-buttonset-multi .fg-button') ){
	    $(this).removeClass("ui-state-active");
	}
	//TODO use this else only for measure/pan toggle.
	else {toggleControl(this);}
    });
});

Ok, should be all you have to know to set up the toolbar! Feel free to reuse the code and improve it :)
Oh, and don’t forget to tweak the CSS to get the perfect look and feel ;)

–anne

Another busy week…and not much time for blogging, so I will just put this presentation online that I gave to a Global Change thinktank meeting I attended in Pretoria.

Also a few weeks back someone I met from INPE, Brasil mentioned Spring as a potential point of collaboration with QGIS. Tonight I finally got a chance to take a look at their site only to be dissapointed when looking at their license to find that it is free as in beer but not open source or Free in the ‘libre’ sense of the word :-(

Intro

A client recently asked me to prepare some aeronautical data for use with QGIS. The data was supplied in about twenty 8 bit paletted GeoTiff images. The client wanted to produce a seamless mosaic so that operators only needed to open a single file to view the rasters. Also he wanted good read access times for the data (hey who doesn’t :-) ).

Trying things out…

Compressed Tiffs

The first thing I tried was simply resaving the imagery as compressed tiffs, and then building a virtual raster dataset using gdalbuildvrt. Compressing the images was easy with a few lines of bash:

#!/bin/bash

for FILE in *.tif
do
  BASE=`basename $FILE .tif`
  NEWFILE=${BASE}_c.tif
  gdal_translate -of GTiff -co COMPRESS=DEFLATE -co TILED=YES $FILE $NEWFILE
done

The script creates a new version of each tif file with a _c.tif suffix. One caveat is that last time I checked, ESRI GIS products didn’t support compressed tiff images (though it was a while ago that I last checked). The compressed tiff format works especially well on imagery with large areas of contiguous colour.

Creating a mosaic

More recent versions of gdal include the excellent gdalbuildvrt application that will build a seamless mosaic from a collection of images. The output is a ‘virtual raster’ – the original files are left untouched and the virtual raster contains a catalogue of the original files. This .vrt file can then be treated as a single file from within your GIS application. Here is an example:

gdalbuildvrt mosaic.vrt *_c.tif

Pop around to QGIS and use the add raster function to add mosaic.vrt and you should see a nice looking mosaic. However, there was one big problem with this approach – using gdalbuildvrt on paletted imagery will result in one master palette of 256 colours being assigned to all images. This results in the images having a speckled, unusable appearance.

So the next thing I did was to use gdal to expand the imagery out to 24bit rgba imagery. Doing this loses most of the gains we made by compressing the tifs, but resolves the palette incompatibility issue.

First I removed all the _c.tif files:

rm *_c.tif

Next I adapted my script to look like this:


#!/bin/bash

for FILE in *.tif
do
  BASE=`basename $FILE .tif`
  NEWFILE=${BASE}_c.tif
  gdal_translate -expand rgb -of GTiff -co COMPRESS=DEFLATE -co TILED=YES $FILE $NEWFILE
done

The noteworthy bit here is the -expand rgb option I have added which will convert the image from a paletted image to a true colour 24 bit image.

Now we can rerun our mosaic build to create a true-colour virtual raster:

gdalbuildvrt mosaic.vrt *_c.tif

Once again you can open the dataset in QGIS to test it. There are some things to bear in mind when building image mosaics. Firstly the imagery should all be in the same Coordinate Reference System. Secondly, they should all have the same spatial resolution. Lastly the band configuration should be the same for all images.

Building pyramids

If you have a particularly large collection of images, you may still find access times a little slow. You can speed things up (at the expense of some disk space) by building pyramids (also know as overviews). Pyramids presample the imagery at various scales so that it does not to resample as much data as you zoom and out of the image. Building pyramids is easy using the gdaladdo command (there is also an option within the QGIS raster properties dialog to do this).

gdaladdo -ro --config INTERLEAVE_OVERVIEW PIXEL --config \
           COMPRESS_OVERVIEW JPEG mosaic.vrt 2 4 8 16

Performance optimisation – Take 1

Because of the nature of my dataset, even with the above tweaking, data access was still too slow for me to be comfortable with passing it along to my client as a solution. To for my next party trick, I tried chopping up the mosaic into little 1000×1000 pixel tiles. I wrote a simple script to do it:

#!/bin/bash
mkdir tiles
XDIM=79011
YDIM=88977
BLOCKSIZE=1000
XPOS=0
YPOS=0
BLOCKNO=0
while [ $YPOS -le $YDIM ]
do
  while [ $XPOS -le $XDIM ]
  do
    echo "$XPOS $YPOS : ${BLOCKNO}.tif"
    gdal_translate -of GTiff -srcwin $XPOS $YPOS $BLOCKSIZE $BLOCKSIZE mosaic.vrt \
       tiles/${BLOCKNO}.tif
    BLOCKNO=`echo "$BLOCKNO + 1" | bc`
    XPOS=`echo "$XPOS + $BLOCKSIZE" | bc`
  done
  YPOS=`echo "$YPOS + $BLOCKSIZE" | bc`
  XPOS=0
done

This is just a quick and dirty approach to creating tiles and I didn’t spend more than about ten minutes writing the script, so it no doubt can be improved in all sorts of ways. One thing you should note is that I hard-coded the image dimensions.

Next I opened the newly created tiles directory and manually removed all unneeded tiles (e.g. all black ones on the periphery of the map.

Once again I built a virtual raster from this dataset:

cd tiles
gdalbuildvrt mosaic_for_tiles.vrt *_c.tif

I also made a pyramid. However performance was still too slow for my liking so I thought I would try the same approach with bigger tiles and compressing the tiles as ECW.

Performance optimisation – Take 2 (ecw)

GDAL includes support for the ECW wavelet compression format. There is a limitation in that if you don’t possess a commercial license, you may not compress a dataset larger than 500mb. Ubuntu does not ship with ecw support so you need to manually build gdal with it enabled. Here is how I built gdal from the source trunk. First I installed all of the build dependencies with some extra packages so that I could get hdf support:

sudo apt-get install build-essential libjpeg62-dev libtiff4-dev libhdf5-serial-dev \
libhdf5-serial-1.6.6-0 libhdf4g-dev libhdf4g libhdf4g-doc

Next I go hold of a copy of the ecw developer kit. Go to the erdas ecw sdk web page (link below) or just search their site for ‘ecw sdk’.

http://www.erdas.com/Products/ERDASProductInformation/tabid/84/CurrentID/1142/Default.aspx

Get the ECW JPEG2000 Codec SDK Source Code (second item listed – first if win version).

After the download you should have:

ecw_jpeg_2000_sdk_3_3_source.zip

Now unzip it:

cd /tmp
unzip ecw_jpeg_2000_sdk_3_3_source.zip
unzip ImageCompressionSDKSourceCode3.3Setup_20070509.zip
sudo mv libecwj2-3.3 /usr/local/src/
cd /usr/local/src/
sudo chown -R timlinux libecwj2-3.3/
cd libecwj2-3.3

Now build the ecw lib:

./configure
make
sudo make install

Next up we will build gdal. First install the subversion client if you don’t already have it:

sudo apt-get install subversion

Nowcheck it out gdal from svn (perhaps you prefer to use one of the stable release versions – I like to live on the edge…).

cd ~
mkdir -p dev/cpp
cd dev/cpp
svn co https://svn.osgeo.org/gdal/trunk/gdal gdal_svn
cd gdal_svn

Now configure and make gdal. Our build will go into /usr/local by default so it can exist side by side with any package installed version you may have.

make clean
./configure --with-ecw=/usr/local
make -j8
sudo make install
sudo echo "/usr/local/lib" | /etc/ld.so.conf
sudo ldconfig

It might take a little while to build but when it is done you can verify that you now have gdal built with ecw support like this:

/usr/local/bin/gdalinfo --formats | grep ECW

The response should look something like this:

  ECW (rw): ERMapper Compressed Wavelets
  JP2ECW (rw+): ERMapper JPEG2000

Ok now we can make ecw compressed images! Lets head back to our data directory and rechop up our mosaic into larger, ecw compressed pieces. My client has an ermapper license so I believe it would be ok to create a single large image for them using gdal, but for our purposes its better to stick under the 500mb limit. Here is a revision of my bash script from earlier:

#!/bin/bash
mkdir ecwtiles
XDIM=79011
YDIM=88977
#Notice I used a bigger block size now
BLOCKSIZE=10000
XPOS=0
YPOS=0
BLOCKNO=0
while [ $YPOS -le $YDIM ]
do
  while [ $XPOS -le $XDIM ]
  do
    echo "$XPOS $YPOS : ${BLOCKNO}.tif"
   # Notice I am using my hand built gdal now!
    /usr/local/bin/gdal_translate -of ECW -co LARGE_OK=NO -srcwin $XPOS $YPOS \
      $BLOCKSIZE $BLOCKSIZE mosaic.vrt ecwtiles/${BLOCKNO}.ecw
    BLOCKNO=`echo "$BLOCKNO + 1" | bc`
    XPOS=`echo "$XPOS + $BLOCKSIZE" | bc`
  done
  YPOS=`echo "$YPOS + $BLOCKSIZE" | bc`
  XPOS=0
done

When I ran that script, it created a bunch of image tiles for me in ecwtiles/ directory (it took a while!). Adding individual tiles to a QGIS project I found they loaded very quickly and pan and zoom operations were satisfyingly responsive. So how would this work out for a mosaic comprised of these tiles? Its easy to find out:

cd ecwtiles
gdalbuildvrt mosaic_for_tiles.vrt *.ecw

Afterwards I was able to open my mosaic in QGIS. It loaded in a few seconds and pan and soom operations take a second or two each. Quite a good result considering the dataset is quite large (~4GB).

Some reflections

GDAL and Linux provide an incredible toolbox to do geospatial dataprocessing. Its easy to experiment with the aid of writing bash scripts. Many people are command line averse, but I love the repeatability that using the command line gives you. Once I work out a procedure, it is easy to document and / or rerun on a different dataset. If I need to tweak it, I just create a copy of the script containing my tweaks. Compare this to a gui environment for dataprocessing where you typically have to click through sequences of dialogs in a very un-repeatable manner.

Sometimes you don’t know what will work best when you start out processing some data, but having a plethora of free tools to use makes it simple to mix and match, chop and change, until you find just the right combination that works for you.

An aside: dealing with libtiff incompatibilities

I did encounter one issue with building pyramids which (along with its solution) is described here:

http://trac.osgeo.org/gdal/ticket/3139

Here is the symptom:

gdaladdo --config COMPRESS_OVERVIEW JPEG --config INTERLEAVE_OVERVIEW PIXEL mosaic.vrt 2 4 8 16
0gdaladdo: tif_jpeg.c:1299: JPEGPreEncode: Assertion `!sp->cinfo.comm.is_decompressor' failed.
Aborted (core dumped)

Eric Raymond and others worked hard to produce and excellent guide on how to ask smart questions. I often give talks about QGIS or open source development in general and I usually end it with a slide or two about how to get help, in which I cover some of the same topics listed in the aforementioned article.

Unsurprisingly, I have been ‘victem’ to many of the scenarios they paint. I use ‘victem’ in a gentle sense as it’s often a source of amusement or triggering mild disinterest rather then getting my hackles up. Asking questions and asking for help is an art that many people just havent grasped. As an open source developer, our most important currency is time. Any question that wastes time is likely to have unsatisfactory outcome. Here then is my bullet list to problem resolution that will hopefully help you get a good result in trying to solve your problems.

  • Try every solution you can think of first
  • Google (and learn to use google advanced search) for a solution
  • Decide on the nature of your question – is it a general user question (how do I do xyz?) or a technical question (feature xyz should do this but does something else, how come?).
  • Post your question to an appropriate list based on the step above. Use your head and pick a list where the people most likely to be able to answer you are most active and responsive
  • Provide just enough detail to make the issue clear
  • Ask only one question at a time / in a thread
  • Dont hijack an existing thread for a new question
  • If someone offers you a suggestion, respond with the outcome of their suggestion when you tried it.
  • If you get the problem resolved, always write back to the list summarising what it was and how it was resolved
  • Always thank the person who helped you – it will make them more inclined to help you in the future
  • Sometimes people don’t know the answer to your question. If you get no answer do some more independent research, reformulate your question and try to repost it.
  • Dont take the fact that someone has tried to help you as an invitation to start an offlist conversation with them. Keeping your discourse on the list makes it searchable and the developer won’t have to answer the same thing another time for someone else.

Sometimes you just won’t get an answer to your question. Thats life, don’t take it personally – you could be asking when others are busy, occupied elsewhere having an actual real life etc. If it happens take another crack at it yourself , try and identify someone who can help you commercially, or start thinking out of the box to get it solved.

A common activity when building a web application is the ‘tweak code -> reload page -> view page source’ cycle. Obviously this gets a bit time consuming and boring after a while. Under Firefox you can use tools like firebug to dynamically inspect and alter the source for your web page, but it adds quite a bit of overhead to the load times of web pages.

Googles chromium (the open source, ubuntu incarnation of their Chrome web browser) is a nice addition to the browser line up on Linux. You can easily add it by dropping a couple of lines into your /etc/apt/sources.list. Visit the Ubuntu Chromium PPA for details.

One of the features I love about chromium is its view-source url prefix e.g.:


view-source:http://yourhost/somepage/?page=4

The above is displayed in the url bar of chromium when viewing the source of a page. The great part is for example above where I am working with a django paginated dataset, I can simply change the url and press enter e.g.:


view-source:http://yourhost/somepage/?page=5

When I do that the source is immediately fetched for the next page and loaded into the page source view. No need to wait for the page to render out and then switch to source view.

Chromium is great to try out and use for cross browser testing. Its very stable in my experience, but note that it has no flash support at this stage.

I’ve been researching options for thin clients over the past few months and trying to get my Fit PCs to boot over LTSP / etherboot (still no joy so far). I do have them working as thin clients via XDMCP thought. Along the way I have found some handy websites which I thought I would make a note of here:

VM from GDM – This one explains how to set up a custom session type so that it logs you directly into a virtual machine instance from the Gnome Display Manager. So you get a full screened windows or whatever session with no Gnome running in the background.

Thin Client Howto – This one is a slightly out of date but still useful run through on how to get clients booting via etherboot.

disklessworkstations.com – This is a commercial vendor of tin client workstations, some of which support dual screen which is useful for GIS workers and software developers. I havent purchased any of their products & would love to hear from those who have.

Ubuntu LTSP quick install guide – contains a useful hint when setting up a 64bit server with 32 bit clients (sudo ltsp-build-client –arch i386 )

This guide has a few handy nuggets like reminding you that every time you update the server you should rerun

ltsp-update-kernels
ltsp-update-image
ltsp-update-sshkeys

One other thing to note is that I believe LTSP is now bundled with the Ubuntu alternate CD so I’m thinking I might do a complete system re-install when 9.10 comes out.

I really want to get this LTSP stuff working nicely – I think its an enviromentally responsible way to go (my FIT-PC’s only consume around 7 Watts of electricity) and makes a lot of sense in terms of configuring and securing one server and all your clients automatically benifiting from a single admin action.

Hey hey hey! Our openModeller article has finally been published in the GeoInformatica Journal. Unfortunately it’s not an open access article so contact me if you want a copy (assuming I am allowed to redistribute it). Here is the abstract:


Species’ potential distribution modelling is the process of building a representation of the fundamental ecological requirements for a species and extrapolating these requirements into a geographical region. The importance of being able to predict the distribution of species is currently highlighted by issues like global climate change, public health problems caused by disease vectors, anthropogenic impacts that can lead to massive species extinction, among other challenges. There are several computational approaches that can be used to generate potential distribution models, each achieving optimal results under different conditions. However, the existing software packages available for this purpose typically implement a single algorithm, and each software package presents a new learning curve to the user. Whenever new software is developed for species’ potential distribution modelling, significant duplication of effort results because many feature requirements are shared between the different packages. Additionally, data preparation and comparison between algorithms becomes difficult when using separate software applications, since each application has different data input and output capabilities. This paper describes a generic approach for building a single computing framework capable of handling different data formats and multiple algorithms that can be used in potential distribution modelling. The ideas described in this paper have been implemented in a free and open source software package called openModeller. The main concepts of species’ potential distribution modelling are also explained and an example use case illustrates potential distribution maps generated by the framework.

University of Witwatersrand

Today Graeme McFerren and I gave a presentation about the FOSSGIS software stack to a bunch of honours students at the University of Witwatersrand. The group was receptive and since they had pretty much never heard of FOSSGIS, it was a great opportunity to show them that there are alternatives out there! The presentation we gave is attached to this post…

TheFossGisStack (odp)
TheFossGisStack (pdf)

Mapnik is an open source, high quality map rendering engine. From the front page of their site:


"Mapnik is a Free Toolkit for developing mapping applications. Above all Mapnik is about making beautiful maps. It is easily extensible and suitable for both desktop and web development."

For some time I have been planning to take a closer look at it but I’ve been to busy to sit down and pick my way through. Over the last two evenings I finally made a go of it and built it on my system and took it for a test drive. First step to installing (on Ubuntu Jaunty 64 bit) was to get the build dependencies from apt.

Mapnik Library Setup

sudo apt-get install g++ cpp libboost-dev libxml2 libxml2-dev \
libfreetype6 libfreetype6-dev libjpeg62 libjpeg62-dev libltdl7 \
libltdl7-dev libpng12-0 libpng12-dev libgeotiff-dev libtiff4 libtiff4-dev \
libcairo2 libcairo2-dev python-cairo python-cairo-dev libcairomm-1.0-1 \
libcairomm-1.0-dev ttf-dejavu ttf-dejavu-core ttf-dejavu-extra libgdal1-dev \
python-gdal postgresql-8.3-postgis postgresql-8.3 postgresql-server-dev-8.3 \
postgresql-contrib-8.3 libsqlite3-dev  subversion build-essential

Note:Above updated for Karmic, 29 Nov 2009

I am installing the most current state of the software from the mapnik subversion trunk. To use mapnik with the Quantumnik QGIS plugin (which we cover further down), you need at least version 0.6.1 of Mapnik. Using the trunk covers that base, but you could also install from the official release tarballs. The mapnik in apt is however too old so we need to build from source.

Next I went to my development dir:

cd ~/dev/cpp

And then checked out the mapnik sources:

svn co http://svn.mapnik.org/trunk

Then build mapnik:

cd mapnik
python scons/scons.py configure INPUT_PLUGINS=all \
  OPTIMIZATION=3 SYSTEM_FONTS=/usr/share/fonts/truetype/ttf-dejavu/
python scons/scons.py
sudo python scons/scons.py install

It didnt take a long time to build and install. Since I am running on the 64bit version of Jaunty, I had to make a small tweak to my system library search path:

sudo vim /etc/ld.so.conf

To which I added the following line:

/usr/local/lib64

And then I updated the library search path:

sudo ldconfig

Next I did the obligatory ‘hello world’ test:

python
import mapnik

Installing the Quantumnik plugin

So the mapnik library is a high quality map rendering engine. To define the maps, you need to create xml definition files (mapfiles). The Quantumnik plugin for QGIS will take an existing QGIS project and generate a mapnik mapfile for it – with some limitations. The limitations stem from limitations in the QGIS symbology infrastructure which does not allow for named styles, multiple symbols per feature and scale based symbol switching in a single layer. Some of these issues should be resolved in QGIS 1.3 with the introduction of some symbology work that Martin Dobias is busy with. In the mean time you can still create a fairly pleasing mapfile using QGIS.

To start I added the mapnik repository to QGIS using the python plugin manager. The repo url is:


http://qgis.dbsgeo.com/

After adding the repo I went ahead and installed the Quantumnik plugin, enabled it in the Plugin Manager in QGIS and then restarted QGIS.

Next I compiled a simple project in QGIS using vector layers stored in a PostGIS database.

Note that the layers need to be listed in the geometry_columns table in PostGIS in order for mapnik to recognise them.

Once my project was compiled, I could use the Plugins -> Quantumnik -> Quantumnik menu to switch the QGIS map canvas to instead use Mapnik to render the view. Once you are happy with the product, use the Plugins -> Quantumnik -> Export XML menu to write my mapnik xml mapfile out to disk.

Here is what my project looked like using the standard QGIS render engine.

Same scene rendered using QGIS native engine

And here is what the mapnik renderer looks like:

Quantumnik in action

Testing the tilelite standalone server

You can test your mapnik mapfile out using a lightweight web serverakin to the webserver used for testing django projects. First we need to install the mercurial code revision system:

sudo apt-get install mercurial

Then install the tilelite server:

cd ~/dev/python
hg clone http://bitbucket.org/springmeyer/tilelite/
cd tilelite
sudo python setup.py install

Now go to where you exported the mapnik mapfile using Quantumnik – I put mine into /tmp/ while testing. Then run the tilelite ‘liteserv’ server:

liteserv.py mapnik.xml

Finally point your web browser at http://localhost:8000 and you should see your map.

Publishing your map

From here you have basically three options for publishing your maps to the greater world:

- use the apache mod_tile module to serve the tiles in google mercator
- use tilecache (see my previous blog post on setting up tilecache). Tilecache has special support to caching mapnik data. You just need to remember to specify the projection = <proj4literal> and extension=png256 options. The latter will reduce the size of the rendered tiles.
- use OGCServer (also a mapnik project) to publish the map as wms. Publishing the map as WMS means you should be able to use it from a WMS client like QGIS since it does not render to tile boundaries.

You can get an idea of how simple it is to overlay your mapnik data onto google maps from this code example: http://bitbucket.org/springmeyer/tilelite/src/tip/utils/gmap.html#cl-28
In a follow up post I plan to walk through the configuration I do to publish a pleasing looking map via mapnik and TileCache.

Thanks

Many thanks to Dane Springmeyer (mapnik developer) who walked me through the whole process outlined above – much appreciated!