Hill shading OpenStreet maps

I've been fooling around with OpenStreetmaps again, this time just trying out a couple of experiments that others did and hopefully in the process develop some interesting ideas to use the data. I've blogged about Openstreetmaps before. A couple of things have changed as Karmic came out. Some build steps require you to download more software and some packages may have been removed since then. I've had bigger problems getting mapnik to work, but that was related to the projection of the import of OSM into Postgis. When starting 'osm2pgsql', don't use the "-l" flag, because that indicates conversion to latlong. The standard mapnik projection however is mercator. Most of the times, when nothing shows up in the rendered image, it's a projection problem somewhere. Just verify it all makes sense going backwards from the mapnik process and you should find it quick enough.

Anyway, a very interesting experiment I tried relates to hillshading. This is an effect where elevation data from different sources are merged (nowadays people use 'mashed') with map data in order to shade particular areas of the map to show where the hills are located. Elevation data for most of the world can be downloaded for free from the Shuttle Radar Topography Mission. You should get version 2.1. For the US, there's elevation data available with 30m horizontal accuracy, the rest of the world only gets 90m (or 3 arcseconds). Some tools exist to help you out in managing that data: srtm2osm (I didn't use this) or a script called srtm_generate_hdr.sh (which I did use). Downloaded the data first, then produced the HDR files.

Then run a simple script:

cd /home/gt/srtm/
for X in *.hgt.zip
do
yes | /home/gt/srtm/srtm_generate_hdr.sh $X
done

This basically generates a number of TIF files for each particular region. These TIF files are height maps. When you look into them, you may be able to recognize certain parts of your country, but don't over-interpret the images you see, because the actual height doesn't really show up there.

In the Ubuntu repositories, there is a package called "python-gdal", which is a set of utilities in python to drive gdal utilities. This can easily be used to stitch the generated TIF files together to generate one huge height map for the entire region. I've used this to get a height map for the Benelux for example, resulting in a 55.0MB file.

Problem is that the projection is incorrect, so you need to assign the latlong projection first, then warp the image, reprojecting it to mercator (what mapnik uses) and in the process attempt to take out any invalid data points (set to 32767).

#> gdal_translate -of GTiff -co "TILED=YES" -a_srs "+proj=latlong" combined.tif combined_adapted.tif
#> gdalwarp -of GTiff -co "TILED=YES" -srcnodata 32767 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -rcs -order 3 -tr 30 30 -multi combined_adapted.tif warped.tif

Cool! So we now have a warped.tif in mercator projection that is basically an elevation map of the terrain. This is however not yet an image which uses the elevation data to 'shade' particular areas of the map. Think of putting the elevation data under a spotlight and then moving the spotlight to some particular origin (somewhere NW), located high above the map. This would create shady areas on the parts of the map where there are hills facing away from the lightsource. The other sides would be normally visible.

Perry Geo luckily did this work for me. You can download his utility here. I had to modify the Makefile to set the GDAL include path (-I/usr/include/gdal) and modify the included gdal library (gdal1.5.0). These tools create a 400-500MB hillshade file (18555x24128 pixels) if you use it for the entire Benelux region. That's still manageable, but in the process of getting there you may need 4GB of memory and a fast processor. The following image shows a minute detail of the Dutch coastline (scaled down 33%).

So that's cool. This image can be used as a bottom layer for some mapnik generated tiles. If mapnik then paints other layers above it slightly translucent, the shading will slightly shine through the image, creating a sense of height throughout the map.

A standard mapnik image from a region in Belgium, close to Maastricht which has a number of hills looks like this:

The same region in Google Maps, rendering terrain data, looks like this:

So, my rendered image looks different from Google's, but the styling wasn't optimized for terrain data in the mapnik renderer.

There's also a reason why the zoom level wasn't increased to about 18 or so. Remember that we downloaded the SRTM3 set, which provides 90m accuracy in the horizontal plane. The height map therefore interpolates between subplanes of equal elevation. This doesn't necessarily improve quality :). Google Maps doesn't go beyond zoomlevel 15, which is a reasonable one. Trying to render images at increased zoomlevels creates pixelated images that don't look nice at all. The height map can improve slightly with better algorithms, but no matter what you do, it will never be perfect. Google Maps even has a couple of areas with hills that look a bit blotchy, irregular or pixelated. The idea is just to get a sense of height of a mountain range anyway and drawing the contours around this map certainly helps to find out where the mountains are (especially for hiking and cycling maps, which use these techniques regularly).

Original sources with the tutorials are here and here.

Building a mapserver with a karmic koala

I've updated my Linux system to Karmic Koala over the weekend. It seems to work quite well. For the first time, I decided to kill all the binaries that somehow made it to my machine over a course of 2 years and do a fresh binaries install, keeping my home mount with data. That worked out well and the machine even booted my RAID-0 array with dmraid without any problems this time. Ubuntu 9.10 works like a treat, you should try it!


Getting down to business, if you want to find out how Google / TeleAtlas renders maps, here's a page that gives you an idea how the process works. A mapserver is basically an image server with a HUGE, HUGE, HUGE number of tiles behind it. Each zoomlevel maintains its own set of images basically, so that's why adding a zoomlevel to a finer-grained level will be costly space-wise. The tiles are constructed by adding GIS information of different types together from a rather large database. The tutorial that I found here is very easy to follow and the most comprehensive on the subject.

In this tutorial, they end up building a little world map, which I attempted and worked out, as you can see. The world map on the right was constructed by the information of SHAPE files on my server. The overview is a very generic image, but the information in the SHAPE file is so great, that I can zoom in to great extents and produce the complete coastline of Antarctica or any country in Europe. Thanks to the efforts of the Openstreetmaps project that is and people and organizations that collaborate with them. Notice also how the world seems to appear slightly distorted. I think this is related to the chosen projection method.

Well, once that is working, you can load your own spatial data in postgis and postgres and start drawing specific parts of your interest on detailed images. Instead of writing your own programs to do that, just use the utilities and scripts of the mapnik project. An example of that is here, central Amsterdam:

So this is great. I can now produce very detailed streetmaps of any region in Holland, reason about those places and ways through a database with spatial reasoning predicates, find out extents of regions, and so forth. Mapnik also provides scripts to generate tiles from a given 'planet' database. Tiling a country however can produce a very large amount of data on your PC, so use with care :). The above image was produced using standard styling rules. It is possible to adjust this styling or replace it entirely, such that it becomes more personal. These generated images, together with a bit of JavaScript and the original PostGis database as a backup to the locations of points of interest are at the core of Google Maps.

Well, another interesting application of GIS information is by super-imposing data from different sources over the data in the database, or not rendering specific sets of data from the source in the database so that they have less information, making it easier to focus on the important bits. You can see how Bjørn Sandvik made a thematic mapper for Google Earth by generating KML from thematic data merged with (simplified versions of) world boundaries. Although KML takes some time to render, especially when in 3D (he wrote a nice, detailed paper about the techniques though), you can generate 2D images by loading your thematic data in PostGis first, then relating your data rows with the geographical data. Using a clever query and the pgsql2shp tool, it should be possible to output a file with the attributes you require for rendering. The last step is then to spit out an XML rendering file for mapnik, which basically filters your attributes, assigns colors or other styling measures and then run it through the mapnik renderer.

There's lots of things one can do here. Be reminded that dealing with these tools can be a bit daunting at first. There's generally no need to write complicated mergers/processors, because you can use PostGis as an intermediate data store, which can output .shp files (the most portable format I reckon), which other tools can visualize or process further.

Drawing your country with openstreetmap

In some previous project, I worked with GIS data to show demographic data on Google Maps. The underlying database is PostGreSQL. This post shows how you can use (part of) the database of OpenStreetMap (verify license!) to extract features from the dataset to paint this onto an SVG or PNG image. Using the PostGIS database in combination with the loaded data, you can extract features (these are like gmaps overlays, or GIS layers) and include them in the image. You want to include railroads and regular roads? Not a problem! You only want to include waterways and the general shape of the country? Can do! The image above was constructed by taking water, waterways, train railroads and forest areas of Holland. Then I zoomed in a bit to show the level of detail.

Since I mostly use Ubuntu, I'll explain the steps used to get things done.

Prerequisites:

• Postgres 8.3 + Postgres 8.3 server-dev package
  
• Postgis
• cmake
• qmake
• libqt4-dev
  
• bz2-dev library
• geos library
• gdal libraries (libgdal, gdal binaries *and* libgdal1-dev)
  
• libxml2-dev header files and library
• Benelux or other data: ( http://planet.openstreetmap.nl/ ) (no need to unpack! 170M packed, 1.8G unpacked )
  

Then, you probably need to download the SVN version of a utility called osm2pgsql. This utility allows you to load in the Benelux data into your postgis database. You can get it using:

svn co http://svn.openstreetmap.org/applications/utils/export/osm2pgsql/

Then, make and make install. This probably doesn't copy the default.style across from the svn checkout directory to where it is expected. So:

# mkdir /usr/share/osm2pgsql
# cp default.style /usr/share/osm2pgsql

Now, you're ready to load the Benelux data:

osm2pgsql --slim -c -l -d  planet-benelux-latest.osm.gz

Ok, so let's start visualizing this for a bit. You'll want to get qgis, compile that, install it, run it and then connect to your database:

# svn co https://svn.osgeo.org/qgis/branches/Release-1_1_0 qgis_1.1.0
# cd qgis-1.1.0
# mkdir build
# cd build
# ccmake ..
(verify output, resolve errors). Press 'g' to generate scripts

Then you should be able to run qgis from the command line: qgis. This starts up the application. If you look for a couple of ESRI shape files, you can load them up and play around with them to see how things work. For the open street map data we have downloaded, we can connect to the database using "Add Postgis Layer". This allows you to select the host, database and tables to load in. It takes a while to get data out of the database, but eventually you get there and it shows all of the Benelux in a particularly bad colorful display :).

A better way to get your data out is by using multiple layers, instead of one. There are a set of utilities that can be useful to load ESRI shapes into the database, but also to get them out. Since we can use where statements, disjunctions and conjunctions in SQL, it is simple to pick out what you're looking for, put this into a shape file and load it into qgis for visualization. Note that the polygons are great for showing a bit of volume and color and that lines are more useful for borders:

pgsql2shp -f forests  "select osm_id, landuse, "name", way from planet_osm_polygon where landuse='forest'"
pgsql2shp -f water  "select name, way from planet_osm_polygon where \"natural\"='water'"
pgsql2shp -f waterwegen  "select name, way from planet_osm_polygon where not waterway is null"
pgsql2shp -f railways  "select name, way from planet_osm_line where railway = 'rail'"

These exported shape files can now be imported into QGis. Then adjust the properties per layer, specifying pen for drawing a black line around polygon areas or not (I didn't in this case). The lines only use a pen and don't have fill colors. Match wood / forest areas with green and water with blue, then select particular types for roads and railways and you're set!

Using the Print Composer function in QGis, you can now export what you've been creating. Make sure to use the Add Map button first, then export to either SVG (Inkscape?) or PNG.

Here's a page that explains the features in the database, but note that the features are not consistently used throughout the database, so you should always check your results:

http://wiki.openstreetmap.org/wiki/Map_Features

More detailed zoom-in near Amsterdam area using the same restricted set of layers and a comparison with Google Maps:


Notice how Google maps have painted areas in a more generalized way, whereas the openstreetmap image is still the original ESRI format. It should not be very difficult to start painting images of GIS data in a format like the Google one above, then overlay the original high-resolution data over those images to indicate the position of cafés, cinema's, etc...

Good luck!