Archive for the ‘GIS’ Category

I’ve been playing with the much-enhanced raster support in PostGIS 2.1 in the last few days. Amongst other things I’ve been producing maps with ST_MapAlgebra, which has changed a little since 2.0.

The example callback function in the documentation is somewhat unhelpfully bare-bones: it always returns 0, touching none of its arguments. So here’s an equally useless — but rather more informative — example of a callback function for the simple one-raster, one-band, one-pixel-at-a-time case:

And here’s something that caught me out: unless you’re passing some userargs, don’t declare your callback function to be strict. If you do, your callback function will never be called because its userargs argument is always null.

Since light can affect happiness, two important pieces of environmental data I add to the Mappiness data set during analysis are: (1) whether a response was made in daylight; and (2) day length when and where the response was made.

To derive these, I need the date, time and location of the response, and sunrise and sunset times for that date and location. R’s StreamMetabolism library provides sunrise/sunset calculations based on NOAA routines. And since my data is in PostGIS, it’s handy to use PL/R to access these.

(Note: what I actually do with the Mappiness data is to use a single call to _r_daylight_period, and calculate the other quantities as a second step. This speeds things up a lot, because I have millions of rows to process and Postgres doesn’t appear to do as much caching of immutable function results as one would like here).

When converting coordinates between WGS84 (GPS) and OSGB36 (UK National Grid), using OSTN02 can gain us a few metres in accuracy over the basic parametric transformation provided by PostGIS’s st_transform via PROJ.4.

There are plenty of ways to get spatial data from a PostGIS database into a Processing sketch.

You can export to CSV or SVG and load it from there; you can query the database directly; or, depending on context, you might choose to generate Processing commands directly, which is the route I went to display a background map of the UK in a recent visualization project.

Imagine you want compare various locations in terms of the availability of a certain type of environment, such as fresh water.

You might want to use a measure of the proximity of that environment — such as the nearest neighbour distance.

You might want to use a measure of the quantity of that environment in the vicinity — such as the proportion of land within a specific radius that is of that type.

Or you might ideally like a measure that combines both of these: one that incorporates the quantity of that environment, but gives greater weight to areas that are nearer, and lesser weight to those that are further away.

In that third case, what you probably want is a kernel-weighted proportion.

Ever noticed how, in Google Earth, marker pins that overlap each other spring apart gracefully when you click them, so you can pick the one you meant?

And ever noticed how, when using the Google Maps API, the exact same thing doesn’t happen?

This code makes Google Maps API version 3 map markers behave in that Google Earth way. Small numbers of markers (up to 8, configurable) spiderfy into a circle. Larger numbers fan out into a (more space-efficient) spiral.

It follows the same basic idea of using series of enlarging search radii to restrict distance calculations to a manageable subset of things-that-might-be-near. The difference is that it uses a geometric progression of sizes (x, x * y, x * y^2, x * y^3, ...) instead of an arithmetic one (x, x + y, x + 2y, x + 3y, ...).

For some distributions of things-that-might-be-near, and tuned with the right parameters (x, y), this turns out substantially faster (I’ve used it to locate the nearest UK postcode to each mappiness response).

PostGIS has an st_askml function. This turns geometries into fragments of KML, and thus takes you most of the way to easy visualisation of spatial queries using Google Earth. But not all the way: these fragments have then to be assembled into a complete document.

I’ve written some wrapper and aggregate functions to automate this. They’re probably deeply inefficient — I wouldn’t advocate building your next web service on them — but for one-off eyeballing of query results I find them really useful.

The key functions are called as_kmldoc; you could see them as the missing aggregate versions of st_askml.