Building Anti-Grain Geometry without SDL

A short post about an annoying problem while building AGG 2.4:
If you don’t have SDL installed, the build fails because of

aclocal:configure.in:103: warning: macro `AM_PATH_SDL' not found in library
configure.in:103: error: possibly undefined macro: AM_PATH_SDL
      If this token and others are legitimate, please use m4_pattern_allow.
      See the Autoconf documentation.

Installing SDL solves this problem, but if you want to skip this, here’s a shorter workaround:

$ echo “AC_DEFUN([AM_PATH_SDL])” > acinclude.m4

Continue as usual with

$ source .autogen.sh
$ ./configure
$ make
$ sudo make install

Using Modcache with MapServer

Suppose you have a setup where you have a dynamic map image rendered by MapServer on a server, and multiple synchronized information displays on which the map should be displayed. Here’s how to use a shared web cache (or web accelerator) in order to reduce load on the MapServer application.

First, some preliminary remarks. If your map is fully or partly composed of static content, you should create a web accessible WMS-C tilecache with TileCache. The TileCache layer can be pre-seeded (pre-populated) through MapServer with tilecache_seed.py. This means that for static layers you don’t even need TileCache or MapServer running on the server at all (which is great, because in my case the server is a resource-constrained embedded platform where python has a footprint that is far too high). OpenLayers supports TileCache layers out of the box via OpenLayers.Layer.TileCache.

For dynamic content, you can use TileCache for caching as well if you are not afraid of running python, but here we are going to use a server-side reverse cache for the map images that are dynamically created by MapServer.

First, it is a good idea to include FastCGI support in MapServer. Head over to www.fastcgi.com and grab their FastCGI library. Building the FastCGI library is trivial:

$ ./configure --prefix=/usr/local
$ make
$ sudo make install

When building MapServer, be sure to enable FastCGI with the --with-fastcgi=/usr/local option.

Verify that your mapserv executable has FastCGI enabled:

$ ./mapserv -v 
MapServer version 6.0.0 OUTPUT=GIF OUTPUT=PNG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=ICONV SUPPORTS=WMS_SERVER SUPPORTS=FASTCGI INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

Check if the output of ./mapserv -v contains the item "SUPPORTS=FASTCGI”.

Now you have to tell your webserver how to access MapServer via FastCGI. If you are using Apache, this page should get you started. I’m using Lighttpd instead, so here are the instructions for editing your lighttpd.conf:

1. Append the FastCGI module mod_fastcgi to the server.modules directive:

server.modules += ( "mod_fastcgi" )

2. Add the following lines to your lighttpd.conf

fastcgi.server = ( 
    "/mapserver" =>
    (   "localhost" =>
        (   "socket" => "/tmp/mapserver-fastcgi.socket",
            "bin-path" => “/usr/lib/cgi-bin/mapserv",
            "max-procs" => 1,
            "check-local" => "disable"
        )
    )
)

Using FastCGI, the MapServer process is only created once. Database connections can be cached with the following PROCESSING directive:

PROCESSING “CLOSE_CONNECTION=DEFER”

Now, when sending a request like “http://localhost:81/mapserver?” with your browser, you should see the following message:

No query information to decode. QUERY_STRING is set, but empty.

(That’s OK, MapServer is working but we haven’t supplied it with a meaningful query string)

Now that MapServer is FastCGI-enabled, we can add mod_cache to the picture, assuming you are using lighttpd. If you are using Apache, the squid caching proxy provides similar functionality. mod_cache uses disk caching, but you can use mod_mem_cache for a memory based cache manager if you like.

Unfortunately, mod_cache for lighttpd is not included in the official lighttpd release and is only available as a patch. Download lighttpd with mod_cache support from lighttpd-improved.

Content is stored in the cache using URI based keys. In order to support dynamic generated FastCGI content, cache.dynamic-mode must be enabled. With cache.support-queries enabled as well, modcache will treat “/mapserver?querystr1” and “/mapserver?querystr2” as different resource and save them to different local cache files.

cache.support-queries   = "enable"      # try to cache query with '?'
cache.dynamic-mode      = "enable"      # modcache will treat "/uri?q1" and "/uri?q2" as different resource
cache.max-memory-size   = 2             # number of megabytes for modcache
cache.bases             = ( "/tmp" )    # make sure directory exists and is writeable
cache.debug             = "enable"
cache.refresh-pattern   = (
    "^/mapserver" => "0"      # zero means cache forever
)

Check the lighttpd error log and the mapserver log file in order to verify that your setup works. When sending identical query strings, MapServer should be called only for the first request. All other requests with the same query string should be handled by modcache.

P.S.: The upcoming MapServer 6.2 release includes functionality for tile caching via MapCache. You might want to look into this once it is released.

Fixing wrong direction of segments in pgRouting output

In my last post I covered the steps needed to set up a PostGIS database on Mac OS X together with pgRouting. Now, pgRouting solves my geospatial routing problems quite well, but I have run into a rather annoying problem: pgRouting returns the edges as they are stored in the database and not as traversed, so on some edges the direction (from-to or to-from) needs to be flipped around.

The following QGIS screenshot illustrates the problem (click for bigger version):

Some segments with wrong direction

You can see clearly that the line arrows do not all point in the same direction.

But don’t despair! My SQL skills are admittedly a bit rusty, but I managed to write a PL/pgSQL script that fixes the problem.

-- Wrapper for shortest_path() function to find the shortest route between two points
CREATE OR REPLACE FUNCTION dijkstra(
    geom_table varchar, source int4, target int4) 
    RETURNS SETOF GEOMS AS
$$
DECLARE 
    r record;
    path_result record;
    v_id integer;
    e_id integer;
    geom geoms;
    id integer;
BEGIN    
    id := 0;    
    FOR path_result IN EXECUTE 'SELECT gid,the_geom FROM ' ||
          'shortest_path(''SELECT gid AS id, start_id::integer AS source, end_id::integer AS target, ' || 
          'length(the_geom)::double precision AS cost FROM ' ||
      quote_ident(geom_table) || ''', ' || quote_literal(source) || 
          ' , ' || quote_literal(target) || ' , false, false), ' || 
          quote_ident(geom_table) || ' WHERE edge_id = gid ' 
        LOOP
            geom.gid      := path_result.gid;
            geom.the_geom := path_result.the_geom;
            id            := id + 1;
            geom.id       := id;                 
            RETURN NEXT geom;
        END LOOP;
    RETURN;
END;
$$ LANGUAGE 'plpgsql' VOLATILE STRICT; 

-- Returns records for the route from start to end with correct directions
CREATE OR REPLACE FUNCTION calc_route(
    geom_table varchar, start_id int, end_id int)
    RETURNS SETOF record AS
$$
DECLARE 
    r record;
    id int;
    prev int;
    i int;
BEGIN
    prev := 0;
    id   := start_id;
    FOR r IN EXECUTE 'SELECT start_id, end_id, route.* ' ||
      'FROM ' || quote_ident(geom_table) || ' JOIN ' ||
      '(SELECT * FROM dijkstra(' || quote_literal(geom_table) || ',' || start_id || ', ' || end_id || ')) AS route ' ||
      'ON ' || quote_ident(geom_table) || '.gid = route.gid ORDER BY route.id; '
    LOOP
        IF (r.start_id = id AND r.end_id <> prev) THEN
            RETURN NEXT r;
        ELSIF (r.end_id = id AND r.start_id <> prev) THEN
            i           := r.end_id;
            r.end_id    := r.start_id;
            r.start_id  := i;
            r.the_geom  := ST_Reverse(r.the_geom);
            RETURN NEXT r;
        ELSE
            RAISE NOTICE 'error: record % % %', r.start_id, r.end_id, r.id;
            RETURN;
        END IF;
        prev := r.start_id;
        id   := r.end_id;
    END LOOP;
    RETURN;
END;
$$ LANGUAGE 'plpgsql' VOLATILE STRICT; 

Sample usage:

DROP TABLE IF EXISTS result;
SELECT start_id, end_id, gid, the_geom INTO result FROM calc_route('network', 22353, 21587)
AS (start_id int, end_id int, id int, gid int, the_geom geometry);

As a result, here is a QGIS screenshot with all segments in the correct direction (again, click for bigger version)

All segments with correct direction
All segments with correct direction

A note about line 59: this script assumes that pgRouting returns the segments in the correct order, which means that either the condition from line 50 or from line 52 evaluates to true. So, line 59 should never be executed. However, some people have reported on the pgrouting-users mailing list that the segments returned from pgRouting are not in the right order – in which case you would see the message from line 59.

pgRouting on Mac OS X 10.6.8

Some notes about how to build pgRouting on Mac OS X, in case I have to do this again in the future. 🙂 I’m writing this from my memory and I hope I did not forget something, but it should work similar to what I’ve described.

pgRouting adds geospatial routing functionality to a PostGIS / PostgreSQL geospatial database. So, as a preliminary requirement, we have to start with installing PostgreSQL and PostGIS.

I’ve done this using MacPorts, so

$ sudo port -d install postgresql90 postgresql90-server postgis

If you like, you can also install pgAdmin3 via MacPorts, an administration program to PostgreSQL

To create a database instance, after install do

$ sudo mkdir -p /opt/local/var/db/postgresql90/defaultdb
$ sudo chown postgres:postgres /opt/local/var/db/postgresql90/defaultdb
$ sudo su postgres -c '/opt/local/lib/postgresql90/bin/initdb -D /opt/local/var/db/postgresql90/defaultdb'

When executing the last command from above, I remember I had some trouble with the postgres user that was created during installation. Check that the postgres user account can be executed using the su command, e.g.

$ sudo su - postgres

If not, make sure that the postgres account has a valid shell

$ sudo dscl . -create /Users/postgres UserShell /bin/sh

Note that the installation has created a startup item for starting postgresql90-server with launchd. It is disabled by default. Execute the following command to start the server, and to cause the server to launch at startup:

$ sudo port load postgresql90-server

As an alternative to launchd, you can also start the database server using:

$ /opt/local/lib/postgresql90/bin/postgres -D /opt/local/var/db/postgresql90/defaultdb

or

$ /opt/local/lib/postgresql90/bin/pg_ctl -D /opt/local/var/db/postgresql90/defaultdb -l logfile start

After starting the server, you might want to look at /opt/local/var/log/postgresql90/postgres.log to check if the server is running OK.

Now, let’s try to install pgRouting.

The pgRouting package is also available from MacPorts, but only an older version 1.03. So, I just went to www.pgrouting.org and downloaded the latest source package pgrouting-1.05.tar.gz from 2010/11/22.

However, you won’t be available to successfully build from this archive (at least not on Mac OS X…)! Instead, you will end up with a couple of “undefined symbols” when linking.

For solving this, go to the pgRouting GitHub repository instead. This issue explains the cause why the linking fails on Mac OS X. Fortunately, this issue has already been fixed, so download a current snapshot from the repository and try to build this one instead.

First, you need to edit cmake/FindPostgreSQL.cmake. Add the include path /opt/local/include/postgresql90/server
and the library path /opt/local/lib/postgresql90 to the file so that cmake is able to locate your PostgreSQL installation.

Then build pgRouting with

$ cmake .
$ make
$ sudo make install

You should end up with something like this

[100%] Built target routing
Install the project...
-- Install configuration: ""
-- Installing: /opt/local/lib/postgresql90/librouting.so
-- Installing: /usr/share/pgrouting/routing_core.sql
-- Installing: /usr/share/pgrouting/routing_core_wrappers.sql
-- Installing: /usr/share/pgrouting/routing_topology.sql
-- Installing: /usr/share/pgrouting/matching.sql

Finally, I hope my memory still serves me right:

$ createdb -U postgres -E utf8 template_postgis
$ psql -U postgres -d template_postgis -f /opt/local/share/postgresql90/contrib/postgis-1.5/postgis.sql 
$ psql -U postgres -d template_postgis -f /opt/local/share/postgresql90/contrib/postgis-1.5/spatial_ref_sys.sql 
$ psql -U postgres -d template_postgis -f /usr/share/pgrouting/routing_core.sql 
$ psql -U postgres -d template_postgis -f /usr/share/pgrouting/routing_core_wrappers.sql 
$ psql -U postgres -d template_postgis -f /usr/share/pgrouting/routing_topology.sql 

Create a pgRouting enabled database using the new template_postgis template

$ createdb -U postgres -T template_postgis routing

That’s it! If you are new to pgRouting, you can find a good example about route calculation using pgRouting here.

Best way to merge color relief with shaded relief map

When creating a raster map featuring a shaded relief combined with elevation colors, then chances are that you are using the gdaldem tool to create both the hillshading (gdaldem hillshade) and the color relief map (gdaldem color-relief). But as soon as you have created both images, you might ask yourself how to combine both images into one. First, I want this to be reproducible in a script, so GIMP or Photoshop is not a viable option.

Let’s take a closer look at the results of some different methods I’ve tried out. The color relief and the hillshade relief images were created with the following commands:

$ gdalwarp -te 29 39 31 41 SRTM30.tif clip.tif
$ gdaldem color-relief clip.tif palette.cpt color.tif
$ gdaldem hillshade clip.tif hills.tif -z 3 -s 111120
Color relief
Hillshade


One option for combining these two images is the hsv_merge.py script from Frank Warmerdam. This script merges the two images in a way where the Hue and Saturation come from the color bands, and the Value comes from the panchromatic band.

But as stated here, “[this approach] works ok, but usually turns vegetation from a deep forest green to a minty pistachio green. This is because vegetation is very reflective in the IR band, which the pan band covers.” Let’s take a look at the result of hsv_merge.py:

$ ./hsv_merge.py color.tif hills.tif merged.tif
Result of hsv_merge.py


As you can see very clearly, the resulting image doesn’t match well with the original color relief in the lower green regions. I didn’t find this result satisfactory, so I looked into whether ImageMagick can be of help here.

A very simple operation would be

$ composite -blend 60 color.tif hills.tif output.tif

but the composite image looks, well, a bit gray (no surprise here):

Simple blend operation


A much better approach consists of two steps: adjusting the gamma level of the hillshade image so that the resulting gamma-corrected hillshade has the RGB value (128, 128, 128) in the plain areas, followed by a lighting composition. I’ve used the “overlay” method below, but you can also experiment with the other lighting methods available.

$ convert -gamma .5 hills.tif hills_gamma.tif
$ convert color.tif hills_gamma.tif -compose Overlay -composite output.tif
Gamma-corrected hillshade
Result


Looks good doesn’t it? Note that the colors of the composite image match very well with the original color relief image.

Finally, note that applying any ImageMagick conversion leads to loosing the georeferencing in the TIFF image. To keep the georeferencing information, you can save the GeoTIFF metadata to a file, and then copy the metadata into a new image:

$ listgeo clip.tif > meta.txt
$ geotifcp -g meta.txt output.tif final.tif 

After this, final.tif should be georeferenced (you can check this with gdalinfo).

C++ metaprogramming

A completely useless but fun-to-write example for C++ template metaprogramming. Did you ever wanted to calculate the control digit of a UIC wagon number at compile-time? This serves just as an example that any computer algorithm can be implemented in some way in a template metaprogram (template metaprogramming is generally Turing-complete).

#include <iostream>

template <long long u, int m>
struct S
{
    static long long const val = S<u/10, m%2+1>::val + ((u%10)*m)/10 + ((u%10)*m)%10;
};

template <int m>
struct S<0, m>
{
    static long long const val = 0;
};

template <long long u>
struct UIC
{
    static long long const val = u*10+(10-S<u,2>::val%10)%10;
};

int main()
{
    std::cout << UIC<21812471217LLU>::val << std::endl;
    return 0;
}

As one would expect, the output is 218124712173 🙂

C++ singleton with policies

The traditional Singleton pattern only supports a default constructor (one that takes no arguments), but recently I needed to create a singleton from a class which required an integer argument in its constructor.

This can be solved quite elegantly with a policy-based design. Take a look at the following (non thread-safe!) implementation:

#ifndef SINGLETON_HPP
#define SINGLETON_HPP

// Creation policies
template<class T> struct CreateWithIntParameter
{
    static T* Create()
    {
        return new T(m_intParameter);
    }
    static const int m_intParameter;
};

// Non-thread-safe singleton
template <class T, template<class> class CreationPolicy>
class Singleton
{
public:
    static T& instance()
    {
        if (!m_instance)
        {
            m_instance = CreationPolicy<T>::Create();
        }
        return *m_instance;
    }
private:
    Singleton();
    Singleton(const Singleton&);
    Singleton& operator=(const Singleton&);

    static T* m_instance;
};

template <class T, template<class> class C>
T* Singleton<T, C>::m_instance = 0;

#endif // SINGLETON_HPP

Here is a small example how to use it:

#include "singleton.hpp"
#include <iostream>

using namespace std;

class Foo
{
public:
    Foo(int a) : m_a(a)
    {
        cout << "Foo constructed" << endl;
    }
    void bar()
    {
        cout << "a = " << m_a << endl;
    }
private:
    Foo();
    int m_a;
};

template<> const int CreateWithIntParameter<Foo>::m_intParameter = 5;

typedef Singleton<Foo, CreateWithIntParameter> FooSingleton;

int main()
{
    FooSingleton::instance().bar();
    FooSingleton::instance().bar();
    return 0;
}

The output looks like this:

Foo constructed
a = 5
a = 5

Independent software consultant and contractor