Barry Rowlingson's GeoSpatial Blog

Document Management System review - Mayan EDMS

2012-01-01T13:57:00.000Z

We're redoing all our web pages in Maths and Stats. This means getting rid of our old Plone 2.1 CMS and replacing it with a groovy new Django-based system using Arkestra - see previous blog entry for more on that.

There's a casualty though. We have an archive of old exam papers going back to 1997. These are sitting on our Plone server (technically as Zope File objects in the ZODB) and are accessed via some custom page templates I wrote so students can get them by calendar year, year of study, or search for a string.

Its a mess. I've been wanting to put them into a proper document management system for a while. Here's my requirements:

Users should be able to search by course code, year of study, and calendar year.
Users should be able to bulk download a zip file of selected PDFs
Admins should be able to bulk upload and enter metadata

I imagine Sharepoint can do this, but I'm no Sharepoint Wrangler and I don't really want to be one. There are a few open source solutions around, such as LogicalDoc and these all seem to be Java-based Enterprisey solutions. Often there's an open-source, or 'community', version, and then an Enterprise supported version with integration features such as drag n drop, OCR, versioning etc.

Then the other day I noticed Mayan EDMS. Its a Django-based DMS which looks like it might well do everything I want. And being open-source if it doesn't, I can have a go at making it do so.

Installation
Smooth as you like. The documentation was spot on, just download, get some dependencies with pip, and it all goes nicely in a virtual environment. The whole virtual env ends up about 90M big on my system, which is comparable to Java-based systems when they bundle the JVM with themselves (which they often seem to do).

Setup
The included settings.py file has options for development and production servers. For development you can run from the usual Django dev server, and for production you can collect all your static files into /static for serving from a web server. For test purposes I set the development mode and did the usual Django syncdb/migrate dance. I couldn't see what the default superuser name and password was, so I just created one with ./manage.py createsuperuser and I was in. I was uploading documents in no time.

Features
Can it do what I want? The version I downloaded currently couldn't. The big missing thing - anonymous read-only access. We don't want our students to need to login to download exam papers, but the permission system on Mayan EDMS doesn't allow anonymous access to anything.

I had a look at creating a new user with read-only permissions. The system has a very fine level of permission control, with lots of capabilities that can bet set on or off. The access to permission is via the usual Role/Group/User style of thing. You create a Role with permissions, then give that role to a User (or a Group of Users). So I created a 'Reader' role and a 'guest' user with that role, thinking we could use that for anonymous access. But the guest user could still create Folders. At that point I had a look at the code.

Code
The code is all on Github and it's incredibly well written, seemingly by one person. It uses lots of existing django apps for some of its functionality which is a good thing. It appears to comply with Python and Django best practices, is commented, and easy to find your way around. I found the Folder creation code and noticed there was no permission check like there was with, say, Document creation. I popped a quick issue on the issue tracker.

Author

The next morning I had a response from Roberto. And that was over New Year's Eve to New Year's Day. He's already added a permission check for Folder creation on a development branch. Coming out in the next version. Not sure if that will also have support for fully anonymous users without logging in (Django returns AnonymousUser if this is the case, rather than a real User object) but it should be sufficient for us. I don't think I've ever had such a quick response from an Open Source dev for a while - and certainly never over New Year celebrations! Credit also goes to Github for its fantastic tracking and notification system, but most of the credit goes to Roberto himself!

Metadata
You can define metadata sets in Mayan EDMS, so I defined a set for our exam papers and added the properties we want to classify our papers with - calendar year, year of study, and course code. Documents can also be tagged, so we could tag very old exams with 'obsolete' - since sometimes course codes are re-used for different courses altogether.

Problems
I'm not sure some of the features are working properly yet. In the document details view, there's a 'Contents' box that I thought would extract the text contents of the document (PDF, OCR'd image etc) but it's blank even for a text document. Searches for text that occurs in that text document don't return the document. Then I noticed one of my PDFs did have text in the box. But only one. Here the text was all run on together with no spaces, but the search engine could find 'scenario' in the text containing 'modelscenariowas', but it would also match on run-on strings like 'wehavealsofoundthat'. This might be a problem but I know it is hard to extract words from PDF files sometimes.

I also tried the OCR functionality on an image but that returned no text too - looking in the Django admin showed me the OCR Document queue - with an entry in it with an error message: "get_image_cache_name() takes exactly 3 non-keyword argument (2 given)". I can run tessaract from the command line successfully so something isn't quite right. I'll have a look at these issues before reporting them on Github (and checking to see if they are fixed in dev versions!)

The other missing function I need as far as I can see is bulk downloads. From a search I'd like to click a few boxes and have a 'Download All' option, that gets everything in a ZIP (or tar.gz) file.

More Features?

The documentation for some of the features are a bit thin. But clicking around and looking at the source code revealed that it has:

Themes: change the look with a config option, several themes supplied
Versioning: keep old docs around. Not sure we need this.
History: show all the changes to files and metadata
Multilingual: with English, Spanish, Portuguese and Russian.
Smart Links, Indexes: I don't know what these are and I couldn't make them do anything. I'm sure they are wonderful.
Signatures: documents can be signed by uploading a signature file. I don't know how this works.
Public Keys: In the setup, you can upload GPG public keys which you can easily get by querying a keyserver. I'm not sure what these keys then let you do.

So overall I'm very impressed with it - it does a lot of what the Enterprisey, Java-based solutions try to give you but in a light python package. Just about the only Enterprise feature I've seen that's not here is true Filesystem integration via SMB or WebDAV, even though there's "Filesystem Integration" on the features list. I suppose we need to see what that means. I'm sure it's good!

Reading AfriPop data

2011-12-07T10:56:00.001Z

The AfriPop project aims to produce high-resolution population density estimates for Africa. This could be very useful for some of our spatial epi projects. So I thought I'd have a look.

Their methods for computing population are explained. I can't find any licensing or usage restrictions on the data, so I'll have to clear that up before doing anything serious with it. Let's download something and read it in.

The download form asks for a few details, but there's no email check (or 'I have read your terms and conditions' checkbox) so fill in as honestly as your conscience requires. You'll get a zip file to download.

When unzipped you'll have three files. For Benin these are called apben10v3 with extensions hdr, flt and prj.

The flt file is the data, the hdr is header info, and prj is projection info. Use gdalinfo to get a summary - run it on the flt file:

$ gdalinfo apben10v3.flt

Driver: EHdr/ESRI .hdr Labelled

Files: apben10v3.flt

apben10v3.hdr

apben10v3.prj

[WGS84 coordinate system info here]

Origin = (0.773989957311480,12.409206711176100)

Pixel Size = (0.000833300000000,-0.000833300000000)

Corner Coordinates:

Upper Left (0.7739900,12.4092067) (0d46'26.36"E,12d24'33.14"N)

Lower Left  (0.7739900,6.2252874) (0d46'26.36"E,6d13'31.03"N)

Upper Right (3.8522002,12.4092067)(3d51' 7.92"E,12d24'33.14"N)

Lower Right (3.8522002,6.2252874) (3d51' 7.92"E,6d13'31.03"N)

Center      (2.3130951,9.3172471) (2d18'47.14"E,  9d19' 2.09"N)

Band 1 Block=3694x1 Type=Byte, ColorInterp=Undefined

  NoData Value=-9999

All pretty good. A standard ESRI file (except there's often nothing really standard about ESRI file formats). Let's read it and plot it using the raster package:

> bpop=raster("./apben10v3.flt")

> image(bpop)

Oh dear. Reading the docs about these files tells me that sometimes gdal just doesn't get it right. It said "Type=Byte", but the file is 109652696 bytes long, which is exactly 7421x3694x4 bytes. These are 4 byte values. In fact the AfriPop page says these files are "ESRI Float". Luckily there's an easy way to persuade gdal these are floats. Open up the hdr file (NOT the flt file. Don't open the flt file!) in a text editor (emacs, vi, notepad++, notepad) and add:

pixeltype float

to the items. Then try again:

> bpop=raster("./apben10v3.flt")

> plot(bpop,asp=1)

Looks good. Note we can get away with asp=1 since we're close to the equator. Or we can give it a projection:

projection(bpop)="+proj=tmerc"

plot(bpop)

and that's just as good. We can also zoom in using the very nifty zoom function from the raster package:

and inspect the data further.

Note that the raster package does a really good job of dealing with large rasters. When plotting these things it is sampling the raster file so not all of it is read into memory. I don't have a monitor capable of showing 7000x3000 pixels (yet) so raster just grabs enough to look good. But if you try and do anything with this size of raster, you will strangle your machine and it may die. Don't do this:

plot(log(bpop))

for example. Let's sample our raster down a bit for some more exploring. First I need to restart R after the last line crashed out. I'll sample down to something of the order of a few hundred pixels by a few hundred pixels by specifying the size as 500*500 - sampleRegular will work out the exact resolution, keeping the aspect ratio. It's smart like that:

> library(raster)

> bpop=raster("apben10v3.flt")

> projection(bpop)="+proj=tmerc"

> plot(bpop)

> bpop2 = sampleRegular(bpop, 500*500, asRaster=TRUE)

> bpop2

class : RasterLayer

dimensions : 708, 352, 249216 (nrow, ncol, ncell)

resolution : 0.008744915, 0.008734349 (x, y)

extent : 0.77399, 3.8522, 6.225287, 12.40921 (xmin, xmax, ymin, ymax)

coord. ref. : +proj=tmerc +ellps=WGS84

values : in memory

min value : 0

max value : 336.4539

> plot(bpop2)

> plot(log(bpop2))

[Perhaps we should use 'resample' here?] Now the log plot is doable:

BUT these values are now wrong. The units of the original data are persons per grid square. Our grid squares are now bigger then the original grid squares, so more people live in them. We can correct this by scaling by the ratio of the cell sizes:

> bpop2= bpop2 * prod(res(bpop2))/prod(res(bpop))

> plot(log(bpop2))

and now we get a more correct map:

I'm not sure how 'correct' or optimal this is of downsampling population density grids is - for example the total population (sum of grid cells) will be slightly different. Perhaps it is better to rescale the sampled grid to keep the total population the same? Thoughts and musings in comments for discussion please!

optifix - optim with fixed values

2011-10-18T14:24:00.000+01:00

So you've written your objective function, fr, and want to optimise it over the two parameters. Fire up optim and away you go:

> fr <- function(x) { ## Rosenbrock Banana function

x1 <- x[1]

x2 <- x[2]

100 * (x2 - x1 * x1)^2 + (1 - x1)^2

}

> optim(c(-1.2,1), fr)[c("par","value")]

$par

[1] 1.000260 1.000506

$value

[1] 8.825241e-08

Great. But sometimes you only want to optimise over some of your parameters, keeping the others fixed. Let's consider optimising the banana function for x2=0.3

> fr <- function(x) { ## Rosenbrock Banana function

x1 <- x[1]

x2 <- 0.3

100 * (x2 - x1 * x1)^2 + (1 - x1)^2

}

> optim(-1.2, fr,method="BFGS")[c("par","value")]

$par

[1] -0.5344572

$value

[1] 2.375167

(R gripes about doing 1-d optimisation with Nelder-Mead, so I switch to BFGS)

Okay, but we had to write a new objective function with 0.3 coded into it. Now repeat for a number of values. That's a pain.

The conventional solution would be to make x2 an extra parameter of fr:

> fr <- function(x,x2) {

x1 <- x[1]

100 * (x2 - x1 * x1)^2 + (1 - x1)^2

}

and call optim with this passed in the dot-dot-dots:

> optim(-1.2,fr,method="BFGS",x2=0.3)[c("par","value")]

$par

[1] -0.5344572

$value

[1] 2.375167

Okay great. But suppose now we want to fix x1 and optimise over x2. You're going to need another objective function. Suppose your function has twelve parameters... This is becoming a pain.

You need optifix. Go back to the original 2 parameter banana function, and do:

> optifix(c(-1.2,0.3),c(FALSE,TRUE), fr,method="BFGS")

$par

[1] -0.5344572

$value

[1] 2.375167

$counts

function gradient

27 9

$convergence

[1] 0

$message

NULL

$fullpars

[1] -0.5344572 0.3000000

$fixed

[1] FALSE TRUE

Here the function call is almost the same as the 2-parameter optimisation, except with a second argument that tells which parameters to fix - and they are fixed at the values given in the first argument. The values in that argument for non-fixed parameters are starting values for the optimiser as usual.

The output is the output from the reduced optimiser. In this case it is one parameter value, but there's also a couple of other components, $fullpars and $fixed, to tell you what the overall parameter values were and which ones were fixed.

Optifix takes your objective function and your specification of what's fixed to create a wrapper round your function that only takes the variable parameters. The wrapper inserts the fixed values in the right place and calls your original function. It then calls optim with a reduced vector of starting parameters and the wrapped objective function. Hence the wrapped function is only every called with the variable parameters, and the original function is called with the variables and the fixed values.

Optifix should work on optim calls that use a supplied gradient function - it wraps that as well to produce a reduced-dimensional gradient function for optim to work on.

Code? Oh you want the code? Okay. I've put a file in a folder on a server for you to download. Its very beta, needs testing and probably should find a package for its home. All things to do.

Mind The Gap

2011-07-26T18:41:00.000+01:00

Every so often Prof D watches Hans Rosling's TED talk or TV show and thinks to himself, "Barry could do that". And so the other week he said to me, "Could we do graphics like that in R?". Because he only really understands R, and Ben had recently cooked up some interactive R plots using tk widgets. Here's the kind of thing that www.gapminder.org can serve up with their Adobe Flash widget:

Well obviously you could do this in R. Or something like it. But why? On CRAN there's the googleVis package which pumps your data into google's visualisation library and you can do just this kind of chart. Doing it in vanilla R would be painful, and the whole chart would have to be redrawn every time. It's not really designed for it.

So I looked for another solution. Javascript in the browser is pretty swift now, and there's visualisation packages. I recently heard of D3, and thought I'd give it a go. Pretty soon I'd knocked out a space-time visualisation plot, and then spent another day or so mimicking gapminder's charts.

Here then is preview number 1.

This is created by a pretty tiny (but at the moment, untidy) chunk of javascript plus some R to write out the data in JSON format. The heavy lifting is done by D3, which has nice data-driven ways of controlling the visualisation. The user can click on the back/fwd links and the data advances a year, and it transitions smoothly between the data points. The colouring is by geographic region, and the circle size is by one of the variables. It's easy to switch the axis or size variables from the Firefox javascript console, and will be even easier when a proper GUI is wired up.

I'll never have the time to make it as polished as Gapminder or Googlevis, but here is the start of an open-source version of these charts.

Ooh did I say 'open-source'? Anyone want the code?

A Middle-aged Geek's Guide To The Arkestra : First Movement

2011-07-25T16:45:00.000+01:00

Intro

"Arkestra is an intelligent semantic Django-based CMS for organisations and institutions". Which is seemingly exactly what we want for our department web pages.

The showcase site, Cardiff University's School Of Medicine is rather impressive.

Arkestra is built on Django-CMS, and is open source. So I decided to play around with it. This blog entry will relate my progress. All this work is being done on an Ubuntu 10.04 box.

Installation Packages

First, create a new virtualenv and get into it. Install the required packages using pip, and the cutting-edge stuff from source. Some of these packages are necessary for django-cms and Arkestra, and some are just necessary for the Arkestra example. I'm not sure which is which, so I grab everything.

virtualenv --no-site-packages Test
. Test/bin/activate
cd Test
pip install PIL
pip install django
pip install BeautifulSoup
pip install django-typogrify
pip install pyquery
pip install -e git+http://github.com/stefanfoulis/django-widgetry#egg=django-widgetry
pip install django-filer # gets easy-thumbnails, django-mptt

pip uninstall -y django-mptt

pip install django-cms
pip install django-polymorphic
pip install django-appmedia
pip install -e hg+https://bitbucket.org/spookylukey/semanticeditor#egg=semanticeditor
pip install -e git+git://github.com/ojii/django-classy-tags.git#egg=classytags
git clone git://github.com/evildmp/Arkestra.git src/arkestra

I tried to do all this with a pip requirements file but failed - the problem seems to be that django-filer will pull in django-mptt, but django-cms can't work with it.

For reference, here's the output of pip freeze telling you what versions of what I have this working with:

BeautifulSoup==3.2.0
Django==1.3
PIL==1.1.7
South==0.7.3
distribute==0.6.10
django-appmedia==1.0.1
-e git://github.com/ojii/django-classy-tags.git@8ad613b3bc0310ba9afb206136cfadbdfa8e6b01#egg=django_classy_tags-dev
django-cms==2.1.3
django-filer==0.8.2
django-polymorphic==0.2
django-typogrify==1.2.2
-e git+http://github.com/stefanfoulis/django-widgetry@57bbf529509e805a74b3e8f36e6938bef591e505#egg=django_widgetry-dev
easy-thumbnails==1.0-alpha-17
lxml==2.3
parti-all==0.0.7.22
pyquery==1.0
-e hg+https://bitbucket.org/spookylukey/semanticeditor@3a40b4766f32262be2cfe131cf20a19703bfc339#egg=semanticeditor-dev
smartypants==1.6.0.3
textile==2.1.5
wsgiref==0.1.2

Config and Run

Head into the arkestra directory that git got for you, and try running the server:

cd src/arkestra/example
export PYTHONPATH=.. # to find the arkestra modules
python manage.py runserver

It should run, but you will fail to log in because there's no user accounts. Also, the example.db file isn't compatible with the version of django-cms and possibly other packages. So you can move it and run syncdb again:

mv example.db example.db.orig
python manage.py syncdb   # create db tables, create admin user
python manage.py runserver

Now you can log in and create Page objects, and view them on the site. Don't forget to tick the published and navigation boxes if needed. To add content, you add a Text/Layout plugin and put the content in there.

Some 404s

Watch your log output. At some point you may see some 404s caused by some magic autocomplete jQuery javascripts not being found. Try this from the example/ folder:

cd media
ln -s ../../../arkestra/arkestra_utilities/static/ javascript

That should stop the 404s.

Loading the Demo Site

Note that there's a JSON-encoded version of the database supplied in example.json. In theory all you need to do is:

python manage.py loaddata example.json

over an empty database and get it all going. But there's those pesky database differences to deal with. I did go through it editing the changes and got it mostly working, but you'll probably want to build your Arkestra site from scratch. Which is what the second movement in this piece will probably be about.

Thoughts on OSGIS 2011

2011-06-26T16:44:00.000+01:00

OSGIS is the UK Open-Source GIS conference. This is the third year it has run, and the first time I've made it there. Nottingham University has been its home every time, mostly due to the presence of the Centre for Geospatial Science on campus, and the enthusiasm of Suchith Anand and Jeremy Morley, who seem to be running on free energy all the time.

The meeting follows the standard template for scientific and technical conferences - there's some workshops, some invited keynote papers, and some general presentations. There's no poster session, and no novel ideas like lightning talks or group discussions - except perhaps in the Open Gov Workshop, but I wasn't there. Not that group discussions don't happen elsewhere, but only over coffee, great lunches, and beers.

Links to the agenda and presentations are all on the Web Site so I won't go into details here. The talks were also webcasted live and recorded, and should be available soon.

The workshop sessions on QGIS were run by Faunalia who let us take away their 8Gb memory sticks packed with a bootable Linux. Very handy. There were also lots of OSGEO Live DVD discs kicking around too, used for the workshop on Ordnance Survey data. Clearly open-sourcers are finding more ways of working in labs full of Windows boxes. I once ran a workshop where we used Virtual Box and an OSGEO Live DVD image to hide the horrors of Windows XP.

So what was hot or not at the conference? Here's my impressions:

On the desktop:

Hot hot hot: Lots of mentions of QGIS.
Also Hot: gvSIG (a day of workshops), GRASS here and there, OpenJump poked its nose in.
Not Hot: No-show GIS packages: uDig, SAGA.
Hot: OpenLayers and GeoExt.
Not hot: The Cloud. Its not the answer to everything.

On the server:

Hot: PostGIS continues its dominance.
Not Hot: Mentions of Oracle Spatial were usually accompanied by apologies. Usual excuse was 'we had the license anyway'.
Not Hot: No sign of NoSQL DB technology - where's the MongoDB Spatial people?
Hot: WFS, WMS, WCS and even WPS becoming almost commonplace.

In the hands of the programmer:

Warm: Only one statistics package mentioned - R of course - although the Orange Toolbox gets a look-in.
Hot Programming languages: Python, Java, JavaScript, C, R, SQL. I think I saw some Perl.
Cold Programming languages: C#, Ruby, Visual Basic/C.
Warm and fuzzy feelings from: Drupal, Django, ExtJS.
Cold and dead to us: Sharepoint. IIS. .net. Silverlight. Flash.

Hardware:

Hot: Android phones - nearly everywhere!
Cold: iPhones - hard to spot. Either nobody had them or were ashamed to show them.
Cold: Tablets. Think I saw one Samsung Galaxy Tab but nobody was showing off iPads. Or iPhone Maxis as they are sometimes known.

No great surprises there.

A few people were tweeting useful things during the conference, and so I should apologise to Jeremy Morley for not saying hello IRL and only communicating via twitter. I will try and submit a paper on doing GIS in R for next year. It will be more than 140 characters.

Also great to finally meet Suchith - the force behind lots of interesting things going on at Nottingham.

Tyler Mitchell rolled in from across the pond and did well against the jet lag to present a keynote speech and network like crazy for two days. OSGEO are very busy both globally and locally.

In all the conference did feel like a mini FOSS4G, but with a local feel - but not a restriction that cut down the quality of the talks or the participants or the coffee chat or the lunch discussions or the beer conversations. See you next year!

Reading Layers from a Qgis Project File.

2011-06-01T18:27:00.001+01:00

Some notes on reading individual layers out of a Qgis Project file. This should end up in a 'import layers from a project file' plugin.

1. get the xml into some text:

xml = file("project.qgs").read()

2. make a blank QDomDocument and set the content:

d = PyQt4.QtXml.QDomDocument()
d.setContent(xml)

3. Now find the maplayers:

maps = d.elementsByTagName("maplayer")

4. How many?

maps.length()

5. Now read the fourth layer (0-index) into running Qgis:

QgsProject.instance().read(maps.item(3))

Wrap this in a nice gui that tells the user what the map layers are, lets them check boxes, import the layers, marks the current project as dirty. Jobzagoodun.

OpenLayers 2.10 Beginners Guide Book Review

2011-04-12T16:24:00.000+01:00

A man is out driving in unfamiliar countryside. Inevitably he gets lost, so he pulls over to ask a local farmer for directions. "Oh no", says the farmer, "I wouldn't start from here".

It's an old joke (and not a very funny one), but one that springs to mind when people ask certain questions on mailing lists. Someone with very little computing experience will appear on the OpenLayers mailing list (or IRC) and say "How do I put a map on my web page?". My answer would be "Well, I wouldn't start from here".

First, I'd say, learn programming using a nice friendly language - python perhaps. Then get a basic understanding of how client-server systems work, and how the web is implemented via the HTTP protocol. Then follow the 20-something year evolution of HTML - from its humble beginnings through to its happy marriage with CSS and their successful menage-a-trois together with Javascript. Ignore its affair with Java - that was a moment of crazed passion, and Java is now out of the HTML/browser party and busily running the "enterprisey" affairs of corporations.

With a knowledge of HTTP, HTML, CSS and Javascript, you are ready to put maps on your web pages. And I'd choose OpenLayers to do it.

Erik Hazzard's book, "OpenLayers 2.10 Beginner's Guide", published by Packt Publishing, is designed to get people with very little computing experience all the way through that process to the point of quite complex map production. He has to cover everything, and does so pretty well (although I'm not sure about his exposition on Object-oriented programming). Although he doesn't teach programming, he encourages users to cut-and-paste code and then experiment.

The sample chapter will show you how the book works - simple explanations, example code, and open-ended questions and mini-assignments. The Firebug debugger is introduced early on in the book, and I think this is a great idea. Firebug didn't take long to repay me the time I took to grok what it could do. Erik sprinkles it liberally on the examples, since interactive inspection of a program is a great way to really see what it's doing.

Inevitably with a book like this, large chunks of its 372 pages will be similar to dumps of the API documentation from the OpenLayers web site. Erik has fleshed out a lot of the terse API information to include fuller examples to make the book more than just that - I can imagine hardcopies of this book with a multitude of flapping yellow sticky notes for their owners to get rapid access to useful pages. Packt Publishing will sell you the PDF version for those who like to search, and the paper version for those who like to stick.

By the end of the book you should have the knowledge to understand how the sample Flickr map application developed in the last chapter works, and you'll also have a handy reference to a large percentage of the features of OpenLayers. The one glaring omission is that there is no mention of the OpenLayers "Popup" class, typically used for "info-balloons" when clicking on map features. But if you can work through the book, then it should be no problem to work out how to use them from the API documentation.

There seem to be very few books on OpenLayers on the market at the moment, so if you want to make maps and would like a high-quality step-by-step guide to it, then this is the book for you.

Note that OpenLayers, like all web technologies is evolving fast, so watch out for updates and get the latest version of the book. A recent OpenLayers hackfest concentrated on developing functionality for mobile devices, so expect to see more OpenLayers applications for your smartphone soon. With Google insisting on presenting advertising along with their map images, it might be a good time to switch to OpenStreetMap for your map background and OpenLayers for your mapping client.

Barry

Note: I was provided with a PDF copy of the book for review, and I have no connection with the author or publisher.

Tolerance and bridges

2011-03-25T12:28:00.000Z

So I'm modelling the H1N1 Swine Flu outbreak in Glasgow. One of the questions is how different regions (specifically 'datazones') within the study area interact with each other. You'd expect zones next to each other to interact with each other more than those on opposite sides of the city, for example. So our first idea is to use the distance between centroids as an explanatory variable and a distance-decay parameterisation for the effect. That works quite nicely.

But the real world isn't quite like that. Glasgow has a river running through it, splitting our study region into north and south parts. I wanted a new distance metric to reflect this partition.

So I came up with this idea - that the distance between two datazones is the minimum distance between the centroids of the graph formed by adjacent datazones.

First problem of that is computing the adjacencies between datazones. In R, I can just use poly2nb from the spdep package and get that. But I found that the graph created wasn't connected over the study region. Further investigation showed that polygons that clearly looked adjacent weren't being considered as adjacent. Here's a screenshot that shows three regions that look adjacent:

Now why wasn't poly2nb saying they were connected? I zoomed in for the answer:

Check the scale bar. These polygons are about a hundredth of a millimetre apart. That's a very tiny no-man's land. More like a "no-fly zone".

I figured I needed some kind of buffer on the polygons, and started looking into the rgeos package, recently released onto CRAN. It'll do that, and all sorts of things, but then when I tried working out adjacencies using its gIntersects function I discovered it worked to those tolerances already. By some magic the apparently adjacent regions were adjacent again.

However the graph was still not fully connected - it was composed of two subgraphs - the river is wide enough, and runs fully from east to west, that no fuzzy tolerances could really overlap the polygons. And that would probably be the wrong thing to do anyway.

So I loaded the regions into Quantum GIS and then loaded some Google Maps as underlays using the marvellous Google Maps plugin. I created a new vector layer of lines and started adding extra connections where bridges cross the Clyde:

These lines, shown in red above, were saved as a shapefile, then loaded into R and their start and endpoints overlaid on the datazones. That gave me some extra edges to add to the adjacency graph, all the graph stuff being handled by R's igraph package (including the distance-weighted minimum distance algorithm along adjacent datazones).

So now I have a completely-connected graph of datazones, which I can convert into a full distance matrix that reflects the geometry of the study region. The university's High-End Cluster is currently churning out the number for the model fit, and the results will hopefully be published in a journal sometime...

Thanks to: R, Quantum GIS, Google Maps plugin authors Ivan Mincik and Nungyao Lin, R package developers for sp, spdep, rgeos, igraph.

Mapping The Kendal Mint Quake

2010-12-22T19:43:00.001Z

About 11pm last night my wall of photographs rattled. It's not unusual, it seems to happen whenever anyone is having a washing machine delivered and a truck crawls past my house. But a minute after that I was on the computer and someone local was tweeting about her house shaking. Earthquake? A few twitter searches revealed twitterers in Cumbria and North Lancashire had indeed experienced a small earthquake.

So I did what all statisticians would do - go looking for data. Tweets can be geo-coded, but getting reliable info might be hard. I checked the USGS earthquake site, but it's not real-time. But then I found the British Geological Society earthquakes page, with true real-time reports from a bunch of sensors around the country. Under the Monitoring menu there's the Real-Time Data option.

Here you can choose the station and a day and get a seismogram. Here's what the one in Keswick (KESW) was showing:

BIG spike. It went off the scale of the plot. I looked at a few other plots and discovered most of them had recorded it too. And I noticed the quake had hit at slightly different times. Maybe I could use that information to work out the location of the quake.

At this point the BGS site was getting slow as worried Cumbrians checked to see if it was an earthquake or if the local nuclear power station had exploded. Anyway, it was bedtime, so I left it to the morning, as something to do in my lunch hour.

I couldn't find a way to get the raw data, so I decided I'd get what I could from the images. A quick script pulled all the GIF images from the sites, except for a couple that seemed to be off-line. Another script pulled all the station information book HTML files and from there I got the locations.

I then stepped through each of the 24 images, identifying the point near 11pm where the tremors started, and entered the data.

Now I had the location and time of the points where the quake had been sensed. This was all converted to OSGB coordinates using rgdal in R.

I didn't have a zero-time for the quake or a location. A quick plot did show the quake had to be somewhere in Cumbria.

I don't know how professional seismologists work out earthquake epicentres, so I came up with something quick and dirty.

Assume the velocity of the quake waves are constant across the UK. Then if the quake was at location (x,y) then a plot of distance against time should be a straight line, with all the points on it. Some quick experiments showed I could fit a straight line then use the sum of the residuals squared as a measure to be minimised over (x,y). A few lines of code later and I had a coordinate!

> fit$par
[1] 307499.9 506035.7

By now the quake had hit the news and the USGS web site was showing it:

How did that compare to my point? I converted it to OSGB coordinates:

> coordinates(usgsos)
     coords.x1 coords.x2
[1,]  328994.6  500054.7

20km further east, and 6km south. Not bad considering I had to eyeball the time from the seismograms and various other assumptions. I don't know what kind of uncertainty the USGS put on their predictions, and I can't put an uncertainty on my location either since there's not a proper statistical model behind this. I can give you the data if you have a better idea - or you can scrape it again from the web site! Here's a comparison of the epicentre predictions on a map:

Here's what the distance/time plot looks like for my predicted epicentre:

And from the slope of this we can get the speed of travel of the quake waves - about 7,500 km/s - which is about what a quick google search on earthquakes tells me.

Footnote: I think on the bike home tonight I worked out a way to do this statistically. Suppose the earthquake is at (x,y), at time t0, and the waves move with speed v. Then suppose the time of arrival of the quake at a point at distance d is Normally distributed with mean time t0+d/v and variance k(d;theta), where k() is some increasing function of d. This allows us to consider the increasing noise in the system for points further away. Fit the whole thing by maximum likelihood and we get estimates (with errors) for location, time of origin, and wave speed. Something for next lunchtime...

Speeding up some R Code

2010-12-21T16:03:00.000Z

Part of my job involves looking at code written by research students and trying to make it do other things. Often a student is focused on their application area or particular dataset, and so the code and data are inextricably linked. Add another observation to the dataset and you'll find yourself having to change 3467 to 3468 in a dozen places. And maybe find some 6934s and change them to 6936s. The other problem is that research students are mostly new to serious programming, and so still need to develop both intuition and formal skills.

Today I've been working on such code for predicting epidemic outcomes. It was taking 14 seconds to run, which isn't too bad (especially considering some other work I'm doing takes 6 hours using OpenMP on an 8-core machine), but in my inspection I'd spotted some odd R code.

exp(var%*%t(coef[1,]))

The oddity here is that var is one-dimensional and so is coef[1,]. This matrix multiply is really just something like exp(sum(var*coef[,1])). So I tried that, and the code took 6 times longer to run. Matrix multiplications are efficient in R, even in one dimension.

So I switched back to the original code. Time to get profiling. This is easy in R.


Rprof("run.out")
dostuff()
Rprof("")
summaryProfile("run.out")

This told me most of the time was spent in the t function, within the function I was playing with earlier. Then I noticed that pretty much wherever coef was used, it was transposed with the t function.

I refactored the code to do t(coef) once, and then pass the transposed matrix around. That sped the code up from the original 14 seconds to 3 seconds. Sweet.

Note that at each stage I check that the new output is the same as the old output - it can be useful to write some little test functions at this point, but doing formal test-driven development can be a bit tricky for statistical applications in R.

There's still a bit of work to do to make this code a bit more general purpose. Currently it works on predictions for one particular country - we want to do predictions over half of Africa!

A queue in R

2010-09-07T23:27:00.000+01:00

Recently someone asked on the R-help mailing list if there was an implementation of a queue in R. I had a quick look (which the poster should have done anyway) and all I could find was an implementation in the filehash package that used database file storage.

I figured it couldn't be that hard to implement just by sticking things on one end of a list and taking them off the other, so I wrote one between having breakfast and setting off late for work.

This evening I polished it off a bit, and even added logging so you could see how big your queue had grown.

The code is available for download under a CC license from my Maths and Stats web pages.

If anyone wants to add this to a package or make a package of it, let me know and I'll help you write a stack and a deque too.

gvSIG, I want to love you...

2010-05-28T12:49:00.000+01:00

I've just been proofreading an article on gvSIG for the next OSgeo Journal (Volume 6, coming to an HTTP server near you soon). And I so want to love it. The author claims it does everything for him. I want it to do everything for me too.

Currently I use QGis for my day-to-day GIS operations. But the limitations are there - mostly in the geoprocessing arena. Vector and raster operations are a bit fragile and buggy at the moment, mostly because they are fairly new to QGis. A grand effort by Carson Farmer et al has greatly increased the analytical capabilities of QGis but there is still a way to go. Meanwhile gvSIG, with the Sextante toolbox, has it all. At least in theory. Time for me to try it out again.

So, off to the web site, and download the latest version. Its a Java app, and I keep getting bitten by incompatible Java environments so I downloaded the version for Linux with the included JRE. It's a big .bin file, so once downloaded you chmod it executable and run it. Choose an install directory and components, and there it goes. All ready to run.

First thing I notice is that the java environment doesn't resize properly. No matter how I stretch or refresh the window, the blue desktop remains fixed size. The grey area is unusable:

Oh well never mind, let's have a look at the interface:

yes, well peace to you too. This menu has one option concerning DBF file encodings. I went 'huh?' and tried something else. Let's load a shapefile.

First you have to create a new "View", so click on the View icon and there's a new view. Double-click or hit 'open'. Now, how to load that shapefile?

Ah, the menu bar has changed. Okay, let's hunt. Layer? No, not there. Shalom? Nope. File? Nope. Oh, it's under View then Add Layer. That's not a standard place for adding things. Usually the View menu is for stuff like zooming in and out. Oh well.

At least now I've got the option. Hit it!

Yeah. Blank. There's controls in there, since random clicking does things. But I see nothing. Oh well, let's quit and try something else:

Blank again. Most of the dialogs are blank, including Help. There's no obvious errors on the console either. I have to hit Ctrl-C to kill it.

Normally I blame these problems with java apps on JRE incompatibilities, but I have no idea what is going on here. I'm running on Ubuntu Karmic on a PC and everything else works nicely. Obviously something is messing up gvSIG since plenty of other people must be using it. Ideas welcome. I think my next step is to try it with a new clean username. But that shouldn't be necessary in this day and age...

Medical Research Council grant success!

2010-05-03T11:19:00.000+01:00

Back in October I talked about a grant proposal we'd put in to fund some research into disease surveillance. Well, a few weeks ago we heard that our application was successful! So that's nearly half a million pounds coming our way. It'll pay a large chunk of my salary for three years, a smaller chunk of the Prof, and a whole research associate. There's also money for new toys and travel to Africa, the USA and various conferences.

So if anyone with strong statistical skills in R is interested in a research project involving spatial and temporal statistics in epidemic prediction, get in touch and I'll send you some details. Probably starting in September.

Culture Colours and Infographics

2010-05-02T16:27:00.000+01:00

A graphic from the Information Is Beautiful blog has been doing the rounds lately. I first saw it on Apartment Therapy and decided to have a little grumble about it. Then it got some coverage on the Junk Charts blog. I decided to get the raw data and see what it looked like in a more informative presentation.

But how to get the raw data? There's no link to it I could find. So I decided to scrape it.

First, I found the highest resolution image of it I could find on the web. Wasn't massive, but it would do. Then I could load it into R using readGDAL from the rgdal package:

co=readGDAL("image.png")
image(co,red=1,green=2,blue=3)

and there it is in my R graphics window. Now I need to generate a whole bunch of circular points where I want to sample the colours. So out with the sines and cosines. Lots of trial and error getting the radii right for all ten sections, but eventually I had it. It was a seconds work with the overlay function from the sp package to get the 860 colour points.

860? There's only 84 traits? I found and dropped the two rows that corresponded to the A-J index and the J-A index. Now I had all the colours in a nice matrix which I could save to a file for safe keeping.

For a simple presentation, I just created an HTML table where the cell background was the colour. It's hard to do rotated column headers in HTML, so I added a fixed DIVas index, and multiple row headers so you can always see the column heading letters. I think the resulting HTML page makes it much easier to navigate and find the colour for the thing you are looking for, although nobody will want to buy large-format version posters of it.

If anyone wants the raw data, drop me an email.

Twitter Timeline Visualisations

2010-03-12T23:11:00.000Z

Ever wanted to make scrolling timeline visualisations of tweets?

Here's how.

You will need:

The R environment for statistics and graphics (www.r-project.org)
The twitteR and brew packages for R
Some web space

I'm using the MIT Simile Timeline widget to do all the hard work. All I need to do is get the tweets into the right format. To do that I wrote a little R program using the twitteR package to access tweets via the search API, and the brew package to reformat into XML. Here's my entire twitline function:

twitline <- function(q,outfile,num=15,...){
require(twitteR)
require(brew)
srch=searchTwitter(q,num=num,...)
brewSrc="

<% for(tweet in srch){ %>\
\"
title=\"<%=screenName(tweet) %>\">
<%=text(tweet) %>

<% } %>

"
brew(text=brewSrc,output=outfile)
}

Now all I do is: twitline("sea otters","twitter.xml",num=100) and I get an XML file in the right format. Read the timeline docs and you'll see that all you need is:

Timeline.loadXML("twitter.xml", function(xml, url) { eventSource.loadXML(xml, url); });

in the right place in your javascript.

There's a few things need doing here, like scaling the timeline to the events on startup, maybe styling things a bit better, having a way to get more than 100 tweets, linking to twitter.com and so on. Consider it a start.

Also, the grabbing and conversion could easily be done in Python or [IYFPLH]. Even cooler would be on-the-fly updates. Check out the SIMILE docs for more ideas.

Barry

Fun With OpenLayers and Flot

2010-03-11T16:45:00.000Z

OpenLayers is a javascript library for doing client-side maps. Flot is a javascript library for doing charts. Put them together and what have you got?

A system for dynamic interactive spatio-temporal graphics. It's something I've wanted to play with for a while, and now I had an excuse. Our meningitis data from Africa is exactly that - counts of meningitis cases in areas in time periods. Eventually we'll be getting it in real time, and predicting epidemics. But for now everything is simulated data.

What I wanted was a map of our country of interest (in this example, Niger). Clicking on the regions of Niger would show a time-series plot of the case count history for those regions. Flot can happily show multiple line charts, show dates on the X-axis, adjust its axes automatically, colour-code it's lines and so on.

So I started with an OpenLayers map and read in my Health District GML file. That's my spatial data. Then I read in a JSON file using JQuery's AJAX api. That's my time-series data (one time series per region). I create a new plot object with Flotr, initially empty.

Adding a SelectControl to the map lets you hook into clicks and hovers on the regions. Each click toggles a region on or off, and when toggled on the corresponding time-series is added to the Flotr plot. Toggle a region off and the time-series line is removed. Flot takes care of updating the legend and the axes.

The Flot charts get a bit cluttered with a dozen lines, so I decided to restrict the number of selected regions to ten. If the user tries to select an eleventh, the system doesn't let you. You have to deselect another region first. This was all done in the onSelect action of the controller. I should probably also make it display a message at this point.

I also decided to match the line colour on the chart with the area colour on the map. I generated a palette of colours from ColorBrewer and used those. When an area is selected a colour is popped from the palette and used to style the area and the line. When an area is deselected its colour is pushed back onto the array.

After fixing some dangling commas and a few other annoyances it even worked on Internet Explorer.

A plain demo is available, which possibly needs a bit of tweaking to make good. It doesn't have the styling of the above image, but does work in IE and lets you see how it all works.

You may use my javascript code in linkedmap.js freely. The map data is an approximate digitisation from a low-resolution image, so does not accurately represent the areas in Niger. All count data is simulated.

The great British Library DRM Flip-flop

2010-02-02T16:22:00.000Z

Our tech support helpdesk recently received requests from a PhD student and a member of staff to install something they needed to read documents requested from the British Library. A proprietary Adobe system called 'ADEPT' was being used to prevent copying and distribution of papers requested from the BL's loan system.

Problem number one was that our PhD students, like many scientists, are enlightened enough to use Linux. There is no Linux version of ADEPT. Apparently Adobe promised one a year or so ago, but that promise disappeared from their web page. Now there's no plans for it. Thanks.

I emailed BL customer service about this. Explaining both how DRM was a stupid idea and DRM with no Linux version was a dumb stupid idea. But not in those terms, of course. The customer services response was fun. First they said:

"It is unfortunate that your student has problems receiving documents in an elecronic format."

My response was that misfortune had nothing to do with it, but poor decisions on content-delivery did. They went on:

"We supply articles via Secure Electronic Delivery as customers want to receive them electronically."

Yup, we all love saving trees.

"The need for encrypted documents is due to the current agreement the British Library has with publishers"

Encrypted? Here's my GPG Public Key if you want to encrypt documents. I explained how this was not encryption, this was DRM, and was hence doomed to failure. At some point any DRM system has to let the user read the file, or hear the music, or see the movie, and at that point it is clear that any encryption has been decrypted. At that point the user has had, somehow, the decryption key. DRM is just obfuscation of encryption keys. Find the encryption keys and you can bypass the DRM.

And that has been done. There is published on the web a couple of quite short python scripts for bypassing ADEPT DRM. One gets the decryption key. From my reading of the code it seems to be stored in the Windows registry and mashed together with your PC's CPU identifier to tie it to your hardware. So you need an official ADEPT account and all that.

Once that script has got the key, you run the second script on any DRMd files, using the recovered key to produce an unencumbered document (PDF file, usually). This you can print, copy to your laptop, backup and know the backup will be readable in whatever PC you have in ten years time, run through ps2text to create plain text and so on. All the things you'd love to do, and are probably legal, if you had a real copy in the first place.

Customer service went on:

"Your student may find it helpful to contact your University Library as they may be able to help by receiving and printing the documents for her. We also still supply Xerox copies by post."

Tiiiiiimmmmmmberrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr!!!!!!

The great mystery is how the BL's policy has changed in the last four years. I found this quote from a news website:

[The] British Library when speaking to the All Parliamentary Internet Group in 2006, warned that the adoption of DRM technology would "fundamentally threaten the longstanding and accepted concepts of fair dealing and library privilege and undermine, or even prevent, legitimate public good access.

*sigh*

Generating Fake Geographies

2010-01-04T22:42:00.000Z

The other day I asked a question on the R-sig-geo mailing list. How can I generate a set of polygons that look a bit like a real set of districts or other subdivisions? Here's what I came up with:

Take a set of points
Generate the voronoi tiles of the points
Turn the straight-line edges into something wiggly
Clip or cut the resulting polygons

So I wrote some R code to do all this. 1 and 2 are pretty easy using the tripack package to generate the voronoi tiles. Making the edges wiggly is a bit trickier. I implemented a fractal algorithm. Split the edge into two pieces by making a new point somewhere near the centre of the edge. Then repeat for the two pieces to make four pieces. Do this three times. This should make a fairly straight wiggly line. For clipping and cutting the polygons I used the gpclib package (free for non-commercial use only - read the license).

Clipping and cutting is necessary because voronoi tiles can extend way beyond the initial points if the tile edges are nearly parallel. I implemented four clipping modes:

No clipping (so you get the full polygons)
Include polygons that overlap a given boundary polygon
Clip polygons to given boundary polygon
Include only polygons completely inside a boundary polygon

By default the boundary polygon is the convex hull of the input points, but you can specify another polygon for clipping.

Here's my sample data - some points and a blue polygon for clipping:

Here's a fake geography created using these points and clipmode zero - the blue polygon is ignored here:

Note that the polygons bottom left head off some distance south-west. The lack of polygons in the north and west is due to these areas being infinite. They don't form closed tiles. This means the number of polygons you get returned wont be the same as the number of input points.

The wiggling algorithm here may cause lines to cross and overlap. There's no test for that at the moment.

Clipmode 1 only includes polygons that are partly or wholly inside the clip polygon:

This has eliminated all the extending polygons in the south-west. Mode 2 clipping results in the polygons being clipped to the boundary polygon. This is useful if you are creating a set of polygon districts within an existing area such as a country boundary:

Note you can't see the red edges under the blue polygon boundary here - they are there though.

Mode 3 only returns polygons wholly inside the clip polygon:

This only leaves 8 polygons - it looks like one polygon on the inner corner has just been clipped out.

So there it is. I've not polished up the code yet, but it's pretty short and sweet. Interested? Let me know.

Addendum: I promise I'll bundle the code up shortly, but today I've been making tweaks to produce fake river networks:

It's pretty much the same basic procedure - generate a voronoi tiling and wiggle the edges. Then compute a minimum spanning tree (any tree will probably do) with distance as the weight. The resulting map looks a bit like a river network (although I wouldn't want to draw the contours or DEM it came from).

I tried various weighting schemes for the MST. Using distance means it uses up the small edges before the large edges, which tends to leave those big gaps that could be mountain ridges.

Other ideas welcome!

Making Interactive Plots with Flot

2009-12-27T19:50:00.000Z

I've been working on a demonstration web site of a meningitis case reporting system in Africa. The current demo uses OpenLayers to show a map of the spatial distribution of cases, and a couple of PNG plots generated from R to show the cases in time for a particular area. One plot showed the cases for the past year, and one was a plot of the previous four weeks with our prediction of the next four weeks cases (plus 95% confidence interval).

But I thought these were a bit static. So I thought I'd add a bit more functionality. Options? Well, I could use my imagemap package for R to create hotspot areas on the plot so users can hover over points and get information, or click on points to zoom to epidemic predictions and so on. But I decided to try a Javascript plotting library instead. There's a few out there, and I settled on using Flot - although jqPlot looks pretty capable too.

Flot lets you plot lines, points and bars, and you can make filled areas by constructing a data series that goes along the top and back along the bottom, completing a polygon. You can shade and style the points and lines, and you get a legend. If your X-axis is time, you can feed it milliseconds and get dates. Here's what I've got to start with:

The buttons at the bottom control the time period. By clicking on 'Prediction', the user can zoom in to see the most recent few weeks, plus our two- and four-week predictions with 50% and 95% confidence interval. Flot also lets you add hover and click events to points:

Here we're showing that our 2-week ahead prediction has a 92% chance of exceeding our 10-case threshold.

Flot also allows zooming in on plots. See that earlier peak around July? We can select that area:

and Flot can show that time period in detail:

All very nice. I've note tested this on Internet Explorer yet, but jQuery and Flot should hide browser-specifics from the author so it should work. This all worked out pretty smoothly, except for a few problems. The biggest one is that if the user resizes the browser window, the plot doesn't resize properly. Not a problem if your DIV element has a fixed width, but I had width: 80% in my CSS. The fix was to redo the plot on a resize event of the window. In jQuery, that's:

   $(window).bind("resize",function(eventObject){
                      plot = $.plot(plotDiv, d, options );
                });

If anyone wants some more code, or the full example, then just ask!

Barry

Polygon Shapefile Editing with QGIS

2009-12-17T15:40:00.000Z

This is a little excursion into editing polygons with some of the new features in Qgis. I'm using a fresh-out-of-SVN of Qgis 1.4. Thanks to whoever fixed the right-click crash bug one hour after I reported it!

I was given a bitmap file representing the new boundaries of regions we are supposed to be working on for a project. I couldn't get this in a spatial format no matter how hard I tried. I did have some similar boundaries but there were a few changes here and there. I loaded the bitmap into the Qgis georeferencer and using my old boundary shapefile matched some points up and created a world file. I could then load the bitmap into Qgis. Here it is:

There's no way I want to redigitize the boundaries when I have something pretty close, and all the overlapping labels on this mean it would be tricky to automatically redigitize the boundaries. Zoom in and see the resolution for yourself.

So I loaded the nearest shapefile and overlaid it, set the fill colour to something transparent (if you set a fill to 'no fill' then you can't see the selection highlight), and activated the magic super labelling from the labelling plugin:

They match up pretty well, but zooming around shows a few problems:

I want to fix the border between Loga and Dosso to better match the underlying raster. First up, activate the editing mode:

Now I hit a problem. Since polygon shapefiles are doubly-digitised, I couldn't see how to move the points on both polygonal boundaries together. And if I moved first one, and then the other, I would probably end up with slivers or polygons that didn't quite fit. I hit on another strategy.

First, select Loga region:

then use the 'Split Feature' tool to cut off a little bit of Loga in the corner away from the edge. This is just to keep the attributes of Loga somewhere, since in a second we're going to merge it with Dosso:

Now we have two features with the same attributes. It's possible that ID numbers have changed here, so don't rely on them. That might explain a bug to appear later... So now, select Dosso and the main part of Loga, and use the Merge Selected Features tool:

This pops up the dialog asking you where the merged feature should get its attributes from. In this case we want the feature to take the attributes of Dosso:

This leaves a few stray bits lying around on the border:

So use the Delete Ring tool to get rid of them:

Now we are going to split Dosso along the new border. First select it:

Then use the Split Feature tool to split by following the border on the underlying raster:

Here's the bug. It's labelled Dosso with 'Loga'. If you check the attribute table at this point it still says 'Dosso', so maybe everything is okay. We'll press on. What we do next is to select the two parts that will make up the new Loga in preparation for another 'Merge Selected Features':

We do the merge, and get the attributes from the little corner of Loga that we cut off earlier:

Now we have the geometry we want:

Except it's still labelled Dosso as Loga. The attribute table seems fine though - if I select Dosso the right line is highlighted in the attribute table:

Never mind, let's turn off the editing and see what happens:

Ooh! It's all looking good! Dosso is Dosso and Loga is Loga.

Now I hope there is an easier way of doing this - a way of editing points along common boundaries such that both polygons are modified. But I can't find it, and I've asked on #qgis and nobody seems to know. If not, then this might be the best way to do it!

Why Screencasts Can Suck

2009-12-11T22:48:00.000Z

A screencast is an instructional video where you see someone's computer screen as they do something, together with a voice-over of them telling you what they are doing. Google trends shows them starting in 2006 and increasing in search volume ever since. Everyone wants to do a screencast.

But are they being used for the right thing? They have their place, just as a screwdriver is the right tool for putting in screws, but there seems to me to be a lot of people using them in the wrong place.

Here's a few reasons why they suck, together with fixes for ameliorating the suckage:

Speech and graphics aren't searchable in a browser or indexable in a search engine

Fix: include a transcript or at least a number of keywords and keyphrases on the page.

Video is hard to jump to a precise point. If I didn't follow a point, or I want to see how you set something up near the start, I'd like to be able to jump back and forth to key points. This also makes it hard to go through a screencast at your own pace.

Fix: I think YouTube videos can have chapter points. It may also be possible that screencasting software can give you some links for this. Include them.
Fix: Alternatively, overlay your screencast with big chunky step numbers or text in the corner. Then as I scroll the fiddly little control back I can see it change from (3) to (2).

Speech and text can't be cut and pasted. This is a pain with command-line videos, less so with gui clickage.

Fix: supply a text transcript of commands and speech where appropriate.

I'm listening to music. The voiceover doesn't go with my tunes.

Fix: include subtitles. Also helps people with hearing difficulties, which will include me if I carry on listening to Muse at these volumes.

Those problems apply to all screencasts. Here are some things that make bad screencasts:

Loss of resolution. If you are squeezing your entire monitor onto a little YouTube video something is going to get lost. Zoom in if you need to show detail.
Two Many Misteaks. I don't want to see you mistyping things. Or going "oops". Or stuttering. Or forgetting something. Shout "cut" to yourself and do another take. Or do several takes and edit. Remember you only have to get this right once, thousands of people may have to watch it.
Nothing happening. Perhaps Andy Warhol would enjoy an unchanging screen with some rambling banter over the top. But I'm not sure everyone else does. This is supposed to be instructable, not minimalist.

Most of these problems can be worked round by creating a series of screenshots and annotating them. I've used Wink for this in the past. It creates a Flash file which you can either play or step through one stage at a time. You can add callout annotations and make the mouse cursor sparkle when it clicks something. The user can step back and forth.

It's extra work to do this, but then it's extra work to do a screencast properly with transcript, subtitles, or key point markers. Terse and precise text is much better than a rambling voice-over. Screencasts may be great for demonstrating anything that has lots of animation in it, but to demonstrate the operation of a piece of point-and-click software it has it's limitations.

Just sayin'.

Converting an mdb to spatialite

2009-12-01T19:03:00.000Z

There's loads of data on the European Environment Agency's web site. And I'm looking for data for my Open Source Geospatial course in January. That's next month now. A lot of the data is in shapefile or geoTiff form, which you can read right into Qgis, R, and with a bit of fiddling, OpenLayers. However there are some databases on there in .mdb format. That's Microsoft Access format. Oh great.

Luckily there's a bunch of command-line tools for extracting data from mdb files. Then I can import them into spatialite and add the geographic data to the tables.

Here's what I did. First get the mdb file from the EPER data set downloads. I'm not sure why the March 2008 version is considered later than the August version, but never mind. We unzip and now have EPER_dataset_27-03-2008.mdb to play with. I'm using the bash shell here for my unix commands.

Let's just set a variable to the filename because it's such a fiddle to type:

mdb=EPER_dataset_27-03-2008.mdb

Now fire up spatialite-gui and create a new database called eper.sqlite. It should have some tables in it, namely: geom_cols_ref_sys, geometry_columns, and spatial_ref_sys . Quit spatialite-gui.

The next job is to generate the tables in the spatialite file. You can dump the schema from the mdb file and read it into spatialite easy enough, and I found that the sybase version of the schema was accepted better by spatialite than the default access option, or the oracle version:

mdb-schema $mdb sybase | spatialite eper.sqlite

Ignore the errors, they occur because the schema script tries to drop all tables before creating them.

Next up, read all the table data from the mdb and load into sqlite. I use the -I option of mdb-export to dump SQL INSERT statements, and the -R option to add the semicolons that spatialite needs:

for t in `mdb-tables $mdb` ; do
echo $t
mdb-export -I -R ";\n" $mdb $t |spatialite eper.sqlite
done

All my table names and column names are single words, so I don't need any quotes or the -S option to sanitize names. That's not always the case.

Now we have to spatially enable the Facility table. Fire up the gui again on the eper.sqlite file. You should see the new tables.

The Facility table has Longitude and Latitude columns which look like WGS84 coords, and indeed there is a Geographic Coordinate System column. Let's just check:

SELECT * FROM "Facility" where GeographicCoordinateSystem <> 'WGS84'

Oh dear. 28 of them appear to be in ETRS89. We'll ignore that for now, but it's always good to check.

The table needs a geometry column, so we add it thus:

SELECT AddGeometryColumn('Facility', 'the_geom', 4326, 'POINT', 2);

And now we fill that column by making points from the Lat/Long columns:

UPDATE Facility SET the_geom = MakePoint(Longitude,Latitude,4326)

From that point we can load the database in Qgis 1.4.x using the spatialite layer and plot the points, query them etc etc.

The database is a complex beast if you wish to relate the other emissions data to the facility points. I've constructed what I think is an approximate picture of the relations between the tables, but I may be wrong!

The PDF/DBDesigner XML of this is available from me if you ask nicely! You could probably use this now to do stuff like creating a view in spatialite that mapped all emissions of a certain substance in various countries, or summarise emissions counts in countries and so on. There's a Flash Interface to the EPER data, but it could easily be done in OpenLayers... Anyone?

Anyway, that's how to get some mdb database into spatialite and plot it in Qgis!

Choropleth mapping challenge in Qgis

2009-11-16T19:01:00.000Z

Making choropleth maps - how hard can it be? It's just a bunch of data colouring another bunch of polygons, right?

This all started over at Flowing Data where Nathan used Python to hack an existing SVG file to colour the counties of the US according to data read in from a file of county unemployment levels. Dave Smith of Revolution Computing issued a challenge: do this in R. I'm always up for a challenge so I found a shapefile of the US counties with the FIPS county codes, and then using the rgdal and sp packages had a plot up in about twenty minutes. The next day Dave posted the challenge results.

Today on the R-Sig-Geo mailing list there has been a bit of traffic on using R or GIS for cartography. Various people arguing one way or the other about how what constitutes 'publication quality' and various opinions bouncing around. So I thought I'd revisit the challenge - this time using open-source GIS software.

Problem one is getting the unemployment data and map data - the unemployment data is linked as a CSV file from the Flowing Data blog entry, and shapefiles of US counties are downloadable after a quick google search.

There are various variations of the US counties depending on the time they were defined, whether Alaska and Hawaii are in their proper place or tucked neatly around the 'lower 48' to make mapping easier, and whether the data has the right FIPS codes. I found a shapefile with county data from 2000, full FIPS codes, and with AK and HI slotted in.

The next step is to create a geographic data set with the unemployment data as attributes. But first I had to fiddle with the CSV file. It had the FIPS codes in two parts - one for the state and one for the county. The full FIPS code is five digits, made up from the two parts with leading zeroes. I didn't want to use R or Python for this, so I loaded the CSV into Calc, the spreadsheet component of OpenOffice.

The first thing I did was to shift all the data down one row and add some headers. The only important ones are the two parts of the FIPS codes and the unemployment percentage.

To create the full FIPS code I have to concatenate the two parts (in columns B and C) with the leading zeroes - I put this formula into J2 and copy it down the column:

=CONCATENATE(TEXT(B2;"00");TEXT(C2;"000"))

then I add the word 'FIPS' as the header in J1 for this column.

Now I have a spreadsheet with a FIPS and an UNEMPLOYMENT column, and a shapefile with a FIPS attribute. How do I join them?

A couple of ways spring to mind:

Load everything into PostGIS and do a join or create a view
Create a shapefile with the Unemployment data in the DBF

I think doing (1) isn't too difficult, but because I had OpenOffice open I figured I'd save the data as a DBF and see what I could do with Quantum GIS

In the Tools menu of Qgis I found the Data Management options and the Join Attributes function. Exactly what I needed! I set my Join Data to be the DBF of the unemployment data and FIPS to be the field to join with. Hit OK, and in under a minute I had a shapefile with the DBF data in place.

The new layer was loaded. Now all that was left was a bit of clicking in the Layer.. Properties.. Symbology dialog to colour the map in six classes with the same intervals as in the challenge, and using the same colours from the Color Brewer system. I also overlaid a state map with transparent fill and a thick solid white outline to get the same effect as the challenge map.

Here's how it looks in Qgis:

Note the large chunk of Alaska missing - that's a non-matching FIPS code. Oh well, nothing's perfect!

Qgis also has a map composer where you can do basic page layout for cartography. Here's one rendering of the unemployment map:

Here you can add layers, legends, and text, and also play with fonts and spacing and sizes. The maps can be output to SVG and PDF.

I'd call these publication-quality - at least they're a better quality than many maps I've seen published (wacky handwriting font excepted).

This has taken me about 45 minutes while waiting for the rain to stop so I can go home. Any gvSIG orArcGIS people out there want to have a go?

Workshop Report

2009-11-05T22:37:00.000Z

I've just got back from London after a trip for a one-day workshop on mapping software in health research organised by The Infectious Disease Research Network. It was held at the Royal Geographical Society buildings - a wonderful place filled with history from RGS expeditions of the past.

The day was arranged by Mike Head of the IDRN and the meeting chaired by Prof. Graham Moon of Southampton University. The talks were varied enough to sustain interest for over a hundred people stuck in one room for the day. The applied talks before lunch combined the usual stories of fieldwork adventures with the gritty details of database servers and web mapping platforms.

During lunch I set up my e-presentation. It was a series of animated slides talking about the use of open-source geospatial software. I also took the opportunity to distribute some flyers for my course in January and the OSGeo group. I talked through my case-studies several times to interested people and was losing my voice by the time I was finishing off my little sticky chocolate cake.

After lunch two talks grabbed the audience - Tim Fendley's amusing tales of new map systems for pedestrians around London was illustrated by screenfuls of comic signage, many of which I encountered on my London wanderings the day after. Chris Phillips from MapAction talked about the work of getting maps and GIS technology out to disaster areas to help co-ordinate search, rescue, and aid deployment.

The final talk was from Mikaela Keller on the Healthmap.org project. This is a system that takes news feeds and health data and maps it in real time. She talked about how a combination of computer language processing and human scanning produce reliable maps. I'm still not sure how you could use this data to do rigorous statistical analysis but I don't think that's the point of it.

After a summing-up from Graham we adjourned for wine and juice outside the hall, and then as numbered dwindled we remaining few jumped into a couple of taxis and headed for the John Snow pub in Soho. Dr John Snow was perhaps the first person to combine health and mapping when he plotted cases of cholera in 1854 and concluded the source was a water pump near where the pub that bears his name now stands. We raised our glasses to the good doctor after a day of discussion of the sort of work made possible by standing on his shoulders.

Pushing OpenSource in a grant bid

2009-10-07T16:00:00.000+01:00

Today I've been writing a section on a research grant (as well as fixing the prof's typos and spellong misteaks). It's great to be able to put the following into a proposal:

"Software and accompanying documentation and training tools developed and released during the programme will be made available under an Open Source license. This will maximise its accessibility to users, especially those working in health institutions in developing countries. It also enables the project to accept improvements and fixes in the software from the wider community.

The project will also commit to using Free and Open Source software as widely as possible within the research team, so that software solutions are developed that do not have proprietary components as part of their architecture".

The computing infrastructure will include a database server (running PostGIS), a front-end (running Apache and serving maps via OpenLayers and maybe WFS/WMS services), and two compute engines churning out statistics with R. It'll be fun to put all these things together if we get the grant. I'm sure I'll blog the result when we get it.

Tools I Use

2009-10-03T17:52:00.000+01:00

This is just a list of software I use for work. Is there some way I can make these keywords on the blog?

Ubuntu Linux - the free operating system
R - statistics software
Python - flexible programming language
Quantum GIS - desktop GIS with Python extensibility
PostGIS - spatial database
OpenLayers - Javascript client-side web maps

All good free and open source stuff. With these things I can build desktop applications, web maps, spatial data infrastructure and all that good stuff.

GeoIntroduction

2009-10-03T09:31:00.000+01:00

Morning all. This is my new blog site for all things generally work-related.

I work in the School of Health and Medicine at Lancaster University within the 'CHICAS' group - Combining Health Information, Computation, and Statistics'. Basically a group of numerates doing data analysis of anything health or biomedicine related.

Most of my activity centres round geographic data, so here I'll be posting useful snippets, questions, comments, rants and opinions, which will not necessarily be those of my employer. I might also wander off into general computing and IT areas.