Openptmap/Installation

From OpenStreetMap Wiki
Jump to: navigation, search

This page describes how to install the software you need for creating a map like openptmap.org (see also openptmap Wiki page).

Prerequisites

Hardware

A weak CPU or a virtual Internet server, will suffice if you limit the geographical region (e.g. France, Germany, Australia). For rendering a larger area, at least 4 or 8 GB RAM are recommended.

Operating System

We assume that you use Ubuntu >=16.04 as operating system. All of the following steps have been tested with version 16.04.

Prepare your System

This chapter describes the steps which are to be done once. After having run through this work, you will only need to update the OSM data for your map from time to time, or you choose to implement an automatic update procedure which will be explained in the next chapter.

Create a User

For all openptmap purposes you can use your usual user account if you like. However, it is recommended to create a separate user for that. In this example, we decide to create a user account with the name pt. Please log in with a user with administrator rights and enter the following commands:

sudo adduser pt
sudo usermod -aG sudo pt

You will have noticed that we gave sudo rights to user pt. These rights should be removed after the whole installation process has been completed: sudo deluser pt sudo Now logoff, then login with the new user name pt.

Software Packages

You must install all packages necessary to run your own web server. Additionally, you will need some other packages. Ensure the installation of all necessary packages by performing these commands:

sudo apt-get update
sudo apt-get upgrade
sudo apt-get install nano unzip gcc zlib1g-dev
sudo apt-get install postgresql-9.5-postgis-2.2 postgresql-contrib-9.5
sudo apt-get install osm2pgsql
sudo apt-get install python-mapnik mapnik-utils
sudo apt-get install apache2 php php-pgsql libapache2-mod-php libjs-openlayers

At this point, it does not hurt to reboot the system to ensure all services have been started in proper order.

PostgreSQL Database

The GIS database must be intitialized. Please follow these steps:

sudo -u postgres -i -H
createuser -SdR ptuser
createdb -E UTF8 -O ptuser ptgis
createlang plpgsql ptgis
psql -d ptgis -f /usr/share/postgresql/9.5/contrib/postgis-2.2/postgis.sql
psql -d ptgis -f /usr/share/postgresql/9.5/contrib/postgis-2.2/spatial_ref_sys.sql
psql ptgis -c "ALTER TABLE geography_columns OWNER TO ptuser"
psql ptgis -c "ALTER TABLE geometry_columns OWNER TO ptuser"
psql ptgis -c "ALTER TABLE spatial_ref_sys OWNER TO ptuser"
exit

To enable easy database login by user ptuser we need to edit one of the database configuration files. In case you are running Ubuntu with a graphical interface, you could use a more comfortable editor, e.g. gedit instead of nano.

sudo nano /etc/postgresql/9.5/main/pg_hba.conf

Near to the bottom of the file you will find these lines:

local   all         all                               peer
host    all         all         127.0.0.1/32          md5
host    all         all         ::1/128               md5

In these lines change the word "peer" and the words "md5" each to "trust" and close the editor (for nano: Ctrl-O, Enter, Ctrl-X). Now reload the database configuration:

sudo service postgresql reload

For a short test, login to the database by using the previously created database user ptuser:

psql ptgis ptuser

Type \d to see a list with the two tables which we have created some steps above (geography_columns, geometry_columns, and spatial_ref_sys). Then logout with \q

If you encounter any problems, you may find a solution here: PostGIS.

Choose a Project Directory

Our example directory for downloading OSM data and generating the data base contents will be the subdirectory pt of the user pt, "/home/pt/pt", which can be abbreviated as "~/pt" if being logged in with this user. Let's create this directory:

mkdir /home/pt/pt

If you want to use a different directory, please create it and grant all access rights to the user pt. You also will have to replace all cd /home/pt/pt commands in the following examples with cd and the full path to your alternative directory.

Tools osmconvert, osmfilter and osmupdate

These tools are mainly used to accelerate the PostgreSQL import process. If we filter out all information we want to use and drop everything else, osm2pgsql will run faster. On top of this, database queries will be faster too, so the rendering process will be accelerated as well. You can install these programs from Ubuntu repository. They are supplied by package osmctools. However, it is recommended to download and compile the source code because the binary may be out of date. To do this – downloading and compiling – from the command line, enter these commands:

cd /home/pt/pt
wget -O - http://m.m.i24.cc/osmconvert.c |cc -x c - -lz -O3 -o osmconvert
wget -O - http://m.m.i24.cc/osmfilter.c |cc -x c - -O3 -o osmfilter
wget -O - http://m.m.i24.cc/osmupdate.c |cc -x c - -o osmupdate

For further details on both tools refer to osmconvert, osmfilter, and osmupdate.

Get Coastlines

Mapnik needs external shape-file data for coastlines at lower zoom levels; let's get them now:

cd /home/pt/pt
wget http://svn.openstreetmap.org/applications/rendering/mapnik/get-coastlines.sh
chmod +x get-coastlines.sh
./get-coastlines.sh
rm *.tar.bz2 *.tgz *.zip

Since more than 500 MB will be downloaded and decompressed in this step, it may take a while.

Get Icons

Mapnik renderer will need certain icons to represent the objects we want to display on the map. You can create these icons by yourself or download sets of icons from various Internet sources. In this example we take the icon set from openptmap.org:

cd /home/pt/pt
wget -r -l 1 -nd -A \*.png -P symbols http://openptmap.org/f/symbols/

Mapnik Renderer Initialization

For the rendering tool Mapnik, some initializations are to be done. Afterwards, a Mapnik style file needs to be created. You can create the file mapnik_pt.xml on your own or use one of the examples from the Internet. You also may download the public-transport style file from here: openptmap.org/f/mapnik_pt.xml

cd /home/pt/pt
wget -r -l 1 -nd -R index.html -P inc http://svn.openstreetmap.org/applications/rendering/mapnik/inc/
wget http://svn.openstreetmap.org/applications/rendering/mapnik/generate_xml.py
chmod +x generate_xml.py
./generate_xml.py --host localhost --user ptuser --dbname ptgis --symbols ./symbols/ --world_boundaries ./world_boundaries/ --port '' --password ''
wget http://openptmap.org/f/mapnik_pt.xml

Create the Map

This chapter will show you how to create or update the tiles for your map. Tiles are small bitmaps which will be assembled to a whole map by the map browser. We assume that all tasks of the previous chapter have been completed successfully.

The process of map creation involves two steps – filling the database and rendering the map tiles.

Fill the Database

All information we need to create the map tiles must be written into our database. To do this, we need to develop a script which performs several tasks, step by step. We will do this step by step, entering every single command at the command line terminal. That makes it much easier to find errors.

Get the Planet File

It is not recommended to use the file of the entire planet. Please choose the file of an area you are interested in, in this example a part of Germany. The first task will be to download this file and to store it using .o5m file format.

cd /home/pt/pt
wget http://download.geofabrik.de/europe/germany/bayern/mittelfranken-latest.osm.pbf -O - |./osmconvert - -o=a.o5m

Filter the Planet File

We need to do a hierarchical filtering because there will be nodes in ways, ways in relations and relations in other relations. For performance reasons a pre-filtered file will be used for interrelational filtering. In this example we decided to filter public-transport specific information because we want to create a public transport map. The filtering criteria need to be specified in a file. This can be done via commandline like this:

cd /home/pt/pt
cat <<eof >toolchain_pt.filter
--keep=
route=ferry =train =light_rail =subway =tram =monorail =trolleybus =bus =aerialway =funicular
line=ferry =rail =train =light_rail =subway =tram =trolleybus =bus =funicular
public_transport=stop_area =platform
railway=station =halt =tram_stop
highway=bus_stop
amenity=bus_station
aerialway=station

--drop=
railway=platform

--keep-tags=all
railway=station =halt =tram_stop
highway=bus_stop
amenity=bus_station
disused=yes
public_transport=stop_area =platform
aerialway=station =yes
ferry=yes train=yes subway=yes
monorail=yes tram=yes trolleybus=yes bus=yes
ref= uic_ref= ref_name= ref:ibnr= ref:IBNR= name= website= wheelchair=

--keep-way-relation-tags=all
route= line= type= state=
eof

Then the OSM raw data can be filtered:

./osmfilter a.o5m --parameter-file=toolchain_pt.filter -o=gis.o5m

Some seconds or minutes later – depending on the size of the chosen area – we will get the file gis.o5m, containing only that information we need.

Note: We are using .o5m file format, because this will accelerate the filter process significantly. For further information see osmconvert and .o5m.

Transfer the Data to Postgres Database

Now we can transfer the OSM data from the file gis.o5m into the Postgres database. The program osm2pgsql will do this job. To fit the needs of our specialised map, we need to create our own osm2pgsql style file first. This can be done with an editor or by executing these commands:

cd /home/pt/pt
cat <<eof >osm2pgsql_pt.style
node,way aerialway text linear
node,way amenity text polygon
node,way area text polygon
node,way bus text linear
node,way disused text linear
node,way ferry text linear
node,way highway text linear
node,way monorail text linear
node,way name text linear
node,way public_transport text polygon
node,way railway text linear
node,way ref text linear
node,way ref_name text linear
node,way ref:ibnr text linear
node,way ref:IBNR text linear
node,way route text linear
node,way train text linear
node,way tram text linear
node,way trolleybus text linear
node,way uic_ref text linear
node,way z_order int4 linear
node,way website text linear
node,way wheelchair text linear
way way_area real linear
node,way line text linear
node,way state text linear
eof

Now you can start the data import:

osm2pgsql -s -C 1500 -d ptgis -U ptuser -S osm2pgsql_pt.style gis.o5m

Because the .o5m file had been filtered in the previous step, this transfer should take only a few minutes. If you have more main memory to spend, you can increase the number of available MB increasing the parameter -C 1500.

If you encounter any problems, you may find a solution here: osm2pgsql.

Rendering Test

At first, we should test if the renderer works as expected. Thereby, the nik4 program will help us. Let's pick an area for which we have already imported GIS data – in this example: Nürnberg, Germany – generate a test image and examine it with the viewer (here Eye of Gnome).

cd /home/pt/pt
wget https://raw.githubusercontent.com/Zverik/Nik4/master/nik4.py
chmod +x nik4.py
./nik4.py -c 11 49.5 -z 12 -d 800 400 -f png256 mapnik_pt.xml image.png
eog image.png &

Render the Tiles

If the test-run has been successful, the next step should be generating map tiles. For this purpose, we use the script generate_tiles.py. First, the bounding box must be adjusted to our map area. Open the script with nano or gedit and delete all the lines below the line # Start with an overview.

cd /home/pt/pt
wget http://svn.openstreetmap.org/applications/rendering/mapnik/generate_tiles.py
chmod +x generate_tiles.py
nano generate_tiles.py

Replace the deleted lines with these (do not change the 4 spaces indent):

    x1= float(os.environ['X1'])
    y1= float(os.environ['Y1'])
    x2= float(os.environ['X2'])
    y2= float(os.environ['Y2'])
    minZoom= int(os.environ['Z1'])
    maxZoom= int(os.environ['Z2'])
    bbox= (x1,y1,x2,y2)
    render_tiles(bbox,mapfile,tile_dir,minZoom,maxZoom,"pt")

The next step would be to render map tiles and store them into our webserver's public directory /var/www/html/tiles. To do so it is necessary to have sufficient access rights to that directory. Let's care about these rights:

sudo groupadd www
sudo adduser pt www
sudo chgrp -R www /var/www/html
sudo chmod -R g+w /var/www/html

To get these changes work, you will need to relogin. If you had chosen another user account, not pt, than please change the last command accordingly. We need to create a directory where all the rendered map tiles will be stored in:

mkdir /var/www/html/tiles

Before generating a huge bulk of tiles, it's best to try the rendering with a single tile and see if it works. The following commands will create a map overlay tile for the city of Nürnberg (Germany) at zoom level 10.

cd /home/pt/pt
X1=11 Y1=49.5 X2=11 Y2=49.5 Z1=10 Z2=10 MAPNIK_MAP_FILE="/home/pt/pt/mapnik_pt.xml" MAPNIK_TILE_DIR="/var/www/html/tiles/" ./generate_tiles.py

Now open the created tile with an image viewer, e.g. eog:

eog /var/www/html/tiles/10/543/349.png &

If the rendering of this test tile has been successful, we can risk to start the rendering of a larger map, the whole city of Nürnberg, for example. In dependence of the hardware, it may take 15 minutes up to a few hours to render all the tiles. In case you want to abort the rendering, you will have to use a system monitor program like top or the command kill to kill the rendering process because Ctrl-C does not work here.

cd /home/pt/pt
F=1 X1=10.9 Y1=49.3 X2=11.3 Y2=49.6 Z1=4 Z2=17 MAPNIK_MAP_FILE="/home/pt/pt/mapnik_pt.xml" MAPNIK_TILE_DIR="/var/www/html/tiles/" ./generate_tiles.py

Build the Website

The web server Apache will expect our web content at "/var/www/html/". In one of the previous steps we have granted us the necessary access rights to this directory. To display, move and zoom the map on the screen a specialized framework will be needed: openlayers. For this, some additional files must be copied resp. downloaded:

cd /var/www/html
cp /usr/share/javascript/openlayers/OpenLayers.js .
cp /usr/share/openlayers/theme/default/style.css .
mkdir img
cp /usr/share/openlayers/img/* img/
wget http://www.openstreetmap.org/openlayers/OpenStreetMap.js

Now open the existing example index file and replace its contents by the following text. You can use gedit or nano editor for this task: nano /var/www/html/index.html

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html lang="de">
<head>
  <meta http-equiv="Content-Type" content="text/html;charset=utf-8" >
  <title>Public Transport Lines - &Ouml;V-Linien - openptmap.org</title>
  <link rel="shortcut icon" href="favicon_pt.png">
  <link rel="stylesheet" href="style.css" type="text/css">
  <style>  /* avoid empty tiles */ .olImageLoadError {display: none;} </style>
  <script src="OpenLayers.js" type="text/javascript"></script>
  <script src="OpenStreetMap.js" type="text/javascript"></script>
  <script type="text/javascript">
    var zoom=4; var lon=11.08; var lat=49.45; var map;

    OpenLayers.Protocol.HTTPex= new OpenLayers.Class(OpenLayers.Protocol.HTTP, {
      read: function(options) {
        OpenLayers.Protocol.prototype.read.apply(this,arguments);
        options= OpenLayers.Util.applyDefaults(options,this.options);
        options.params= OpenLayers.Util.applyDefaults(
        options.params,this.options.params);
        options.params.resolution= map.getResolution();
        options.params.zoom= map.getZoom();
        if(options.filter) {
          options.params= this.filterToParams(
            options.filter, options.params);
          }
        var readWithPOST= (options.readWithPOST!==undefined)?
          options.readWithPOST: this.readWithPOST;
        var resp= new OpenLayers.Protocol.Response({requestType: "read"});
        if(readWithPOST) {
          resp.priv= OpenLayers.Request.POST({
            url: options.url,
            callback: this.createCallback(this.handleRead,resp,options),
            data: OpenLayers.Util.getParameterString(options.params),
            headers: {"Content-Type": "application/x-www-form-urlencoded"}
            });
        } else {
          resp.priv= OpenLayers.Request.GET({
            url: options.url,
            callback: this.createCallback(this.handleRead,resp,options),
            params: options.params,
            headers: options.headers
            });
          }
        return resp;
        },
      CLASS_NAME: "OpenLayers.Protocol.HTTPex"
      });

    function init() {
      var args= OpenLayers.Util.getParameters();

      map= new OpenLayers.Map("map", {
        controls:[
          new OpenLayers.Control.Navigation(),
          new OpenLayers.Control.PanZoomBar(),
          new OpenLayers.Control.LayerSwitcher(),
          new OpenLayers.Control.Permalink(),
          new OpenLayers.Control.ScaleLine(),
          new OpenLayers.Control.Permalink('permalink'),
          new OpenLayers.Control.MousePosition(),                    
          /*new OpenLayers.Control.Attribution()*/],
        maxExtent: new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34),
        maxResolution: 156543.0399,
        numZoomLevels: 17,
        units: 'm',
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG:4326")
        } );

    map.addLayer(new OpenLayers.Layer.OSM.Mapnik("Map",
      { maxZoomLevel: 17, numZoomLevels: 18 }));
 
    map.addLayer(new OpenLayers.Layer.OSM.Mapnik("Map, pale",
      { maxZoomLevel: 17, numZoomLevels: 18, opacity: 0.5 }));
 
    map.addLayer(new OpenLayers.Layer.OSM.CycleMap("Cycle map",
      { maxZoomLevel: 17, numZoomLevels: 18 }));

    map.addLayer(new OpenLayers.Layer.OSM.CycleMap("Cycle map, pale",
      { maxZoomLevel: 17, numZoomLevels: 18, opacity: 0.5 }));

    map.addLayer(new OpenLayers.Layer.OSM("No Background","blank.gif",
      { maxZoomLevel: 17, numZoomLevels: 18 }));

    map.addLayer(new OpenLayers.Layer.OSM("Public transport lines","tiles/${z}/${x}/${y}.png",
      { maxZoomLevel: 17, numZoomLevels: 18, alpha: true, isBaseLayer: false }));

    if(!map.getCenter()){
      var lonLat= new OpenLayers.LonLat(lon, lat).transform(new OpenLayers.Projection("EPSG:4326"),
        map.getProjectionObject());
      map.setCenter(lonLat,zoom);
      }
    }
  </script>
</head>
<body onload="init();">
<div style="width:100%; height:100%;" id="map"></div><br></body>
</html>

Further information about OpenLayers can be found here.

Test your new Map

Open a web browser and try to access your new website. If your browser is installed on the same computer as the Apache server, type localhost or 127.0.0.1 as URL. If you did not render the complete map yet but have rendered that single tile from the rendering test example above, you will find it here: 127.0.0.1/?zoom=10&lat=49.5&lon=11&layers=B0000T All tiles which have not been rendered are displayed as pink squares.

Please be aware that this map is an overlay map. That means in this case, the background tiles are fetched directly from openstreetmap.org map server, and just the foreground tiles – the public transport layer – comes from your server.

Create and update the Map Automatically

Did all previous steps work without any errors? Then it's time to pack them into a single script and to care for map updates on a regular basis. For the updates we will us a so-called diff-based rendering strategy.

With diff-based rendering only that tiles will be rendered whose OSM data have been changed. With the right call, the program osm2pgsql will supply us with a list of these tiles. This is usually not very helpful for a thematic overlay map because 90% of the changed data will not affect the overlay image, but there is a way to get exactly the information we need: update the PostgreSQL database with OSM-Change files which have been retrieved from intensively filtered data.

Getting Differential OSM Data

To do the map updates there is no need to completely reimport all the OSM data. It will suffice to download and process just that data which have been changed. There are several sources of this kind of data – so-called OSC files (OSM-Change files, .osc) – in the Internet. For updating a local copy of OSM raw data, we will use the program osmupdate. The same job could be done with the well-known multy-purpose program Osmosis als well, of course.

Indirectly filtering OSM-Change File

At first view we can save a lot of effort if we directly filtered the differential data. Unfortunately this does not work as expected. The reason is that OSM-Change file usually contain only that OSM objects which have been changed. If – as in this example – all bus lines shall be rendered and a single node of one of these bus lines had been moved, only this node would be in the OSM-Change file. This node does not have any bus-line information, thus it would be filtered out and the bus-line relevant information would be lost.

To prevent this from happening the filtered data from after an update need to be compared against against the filtered data from before this update. On the base of this comparison a new OSM-Change file is created. This can be done by the --diff function of osmconvert. The resulting OSM-Change file should not bee very large because it is retrieved from filtered data.

Getting a List of all Dirty Tiles

Map tiles which are not up-to-date are referred to as dirty tiles. During map data import the program osm2pgsql generates a list of all dirty tiles. Unfortunately the program does not provide us with the names of all affected tiles because this list is stripped of all redundant information. If there is, for example, the tile 3/0/1 (zoom, x, y) in the list, the also affected tiles 2/0/0 and 1/0/0 will be omitted. Furthermore, a list entry 3/0/1 means that all four tiles of the next zoom level have also been omitted although affected: 4/0/2, 4/0/3, 4/1/2, and 4/1/3.

This is an effective way to reduce the dirty-tiles list length but we still need a list of all effected tiles. For this reason, the dirty tiles list must be expanded accordingly. This can be done with a small C program which needs to be added to our toolchain:

cat <<eof |cc -x c - -O3 -o dtexpand
// dtexpand 2017-03-10 18:00
//
// expands the dirty-tiles list as provided by osm2pgsql
// example: dtexpand 4 17 <dirty_tiles >dirty_tiles_ex
// (c) 2017 Markus Weber, Nuernberg, License LGPLv3
#include <ctype.h>
#include <inttypes.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

int main(int argc,char** argv) {
  int fieldedge[20];
  uint8_t* field[20],*fieldend[20];
  int zoom_min,zoom_max,zoom_sub;
  int z,x,y,zz,xx,yy,xd,yd,q; uint64_t u;
  uint8_t b; uint8_t* bp,*bpa,*bpe;
  char line[40]; char* sp;
  typedef struct {int zs,tx_min,tx_max,ty_min,ty_max;} sub_t;
  sub_t* sub,*subp; int subn;

  if(argc<3 || (argc-3)%5!=0) { fprintf(stderr,
    "Usage: dtexpand ZOOM_MIN ZOOM_MAX <DT_IN >DT_OUT\n"
    "Optional value sets each: ZOOM_SUB TX_MIN TX_MAX TY_MIN TY_MAX\n"
    "         ZOOM_SUB: max zoom level for blindly adding subtiles\n"
    "         TX...TY: concerning tile range at ZOOM_MAX level\n"
      ); return 1; }
  zoom_min= atoi(argv[1]); zoom_max= atoi(argv[2]);
  subn= (argc-3)/5;
  sub= (sub_t*)malloc(sizeof(sub_t)*(subn+1));
  argv+= 2; subp= sub;
  for(z= 0;z<subn;z++) {  // for each subtile value set
    subp->zs= atoi(*++argv);
    subp->tx_min= atoi(*++argv); subp->tx_max= atoi(*++argv);
    subp->ty_min= atoi(*++argv); subp->ty_max= atoi(*++argv);
    if(zoom_min<0 || zoom_min>19 || zoom_max<zoom_min || zoom_max>19 ||
        subp->zs<=zoom_max || subp->zs>19) { fprintf(stderr,
        "Unsupported zoom value(s).\n"); return 2; }
    subp++; }  // for each subtile value set
  subp->zs= 0;
  for(z= zoom_max-1;z>=zoom_min;z--) {  // for each zoom level
    fieldedge[z]= u= UINT32_C(1)<<z; u= u*u+7>>3;
    if((field[z]= (uint8_t*) malloc(u))==NULL) { fprintf(stderr,
      "Not enough main memory.\n"); return 3; }
    memset(field[z],0,u);
    fieldend[z]= field[z]+u;
    }  // for each zoom level
  while(fgets(line,sizeof(line),stdin)!=NULL) {  // for each input line
    sp= line;
    z= 0; while(isdigit(*sp)) z= z*10+*sp++-'0'; if(*sp=='/') sp++;
    x= 0; while(isdigit(*sp)) x= x*10+*sp++-'0'; if(*sp=='/') sp++;
    y= 0; while(isdigit(*sp)) y= y*10+*sp++-'0';
    if(z<zoom_min || z>zoom_max)
  continue;
    zz= z-1; xx= x>>1; yy= y>>1;
    while(zz>=zoom_min) {  // for all lower zoom levels
      u= xx; u*= fieldedge[zz]; u+= yy; bp= &field[zz][u/8];
      if(bp<fieldend[zz]) *bp|= 1<<(u&7);  // (preventing overflow)
      zz--; xx>>= 1; yy>>= 1;
      }  // for all lower zoom levels
    q= 1;
    do {  // for this and all higher zoom levels
      for(xd= 0;xd<q;xd++) for(yd= 0;yd<q;yd++) {  // all these tiles
        xx= x+xd; yy= y+yd;
        if(z<zoom_max) {
          u= xx; u*= fieldedge[z]; u+= yy; bp= &field[z][u/8];
          if(bp<fieldend[z]) *bp|= 1<<(u&7);  // (preventing overflow)
          }
        else {  // at zoom_max level
          printf("%i/%i/%i\n",z,xx,yy);
          // determine relevant subtile zoom level 
          zoom_sub= zoom_max;
          subp= sub;
          while(subp->zs!=0) {  // for each subtile value set
            if(subp->zs>zoom_sub &&
                xx>=subp->tx_min && xx<=subp->tx_max &&
                yy>=subp->ty_min && yy<=subp->ty_max) zoom_sub= subp->zs;
            subp++; }  // for each subtile value set
          // write subtiles
          int xxd,yyd,qq;
          zz= z; qq= 1;
          while(zz<zoom_sub) {  // for each subtile level
            zz++; xx<<= 1; yy<<= 1; qq<<= 1;
            for(xxd= 0;xxd<qq;xxd++) for(yyd= 0;yyd<qq;yyd++)
              printf("%i/%i/%i\n",zz,xx+xxd,yy+yyd);
            }  // for each subtile level
          }  // at zoom_max level
        }  // all these tiles
      z++; x<<= 1; y<<= 1; q<<= 1;
      } while(z<=zoom_max);  // for this and all higher zoom levels
    }  // for each input line
  for(z= zoom_max-1;z>=zoom_min;z--) {  // for each zoom level
    bp= bpa= field[z]; bpe= fieldend[z];
    do {  // for each byte in bitfield
      b= *bp;
      if(b) {  // at least one bit in byte is set
        int be= fieldedge[z];
        for(int bi=0;bi<8;bi++) {  // for each bit in byte
          if(b&1) {  // bit is set
            u= (bp-bpa)*8+bi; x= u/be; y= u&be-1;
            printf("%i/%i/%i\n",z,x,y);
            }  // bit is set
          b>>= 1;
          }  // for each bit in byte
        }  // at least one bit in byte is set
      } while(++bp<bpe);  // for each byte in bitfield
    free(field[z]);
    }  // for each zoom level
  free(sub); return 0;
  }  // main()
eof

Besides its main function this program can also be used to add subtile names for certain regions to the dirty-tiles list. If you want to create a global map up to zoom level 17, for example, and having a small region rendered up to zoom level 18, you may later add this subtile zoom level as well as the region's level-17 tile range to the command line:

dtexpand 4 17 18 69000 70000 44000 45000

If there are two or more of such regions, just enter all of their 5-value sets. For example:

dtexpand 4 17 19 80000 81000 48000 49000 18 69000 70000 44000 45000

Feeding Mapnik with a Tiles List

Mapnik must be enabled to read and process dirty-tiles lists. The following python script will exactly do this. Enter /home/pt/pt directory and create a file with the name mapnik_pt.py and the following contents:

#!/usr/bin/python
# mapnik_pt 2017-02-19 18:30
# reads dirty-tiles file and renders each tile of this list;
# afterwards, deletes the rendered dirty-tiles file;
# call: DTF="dirty_tiles" mapnik_pt.py
# parallel call: DTF="dirty_tiles" PID=123 mapnik_pt.py
# (c) Markus Weber, Nuernberg, License: AGPLv3
from math import pi,cos,sin,log,exp,atan
from subprocess import call
import sys, os
from Queue import Queue
import mapnik
import threading
import shutil

DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi

def minmax (a,b,c):
    a = max(a,b)
    a = min(a,c)
    return a

class GoogleProjection:
    def __init__(self,levels=18):
        self.Bc = []
        self.Cc = []
        self.zc = []
        self.Ac = []
        c = 256
        for d in range(0,levels):
            e = c/2;
            self.Bc.append(c/360.0)
            self.Cc.append(c/(2 * pi))
            self.zc.append((e,e))
            self.Ac.append(c)
            c *= 2

    def fromLLtoPixel(self,ll,zoom):
         d = self.zc[zoom]
         e = round(d[0] + ll[0] * self.Bc[zoom])
         f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
         g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
         return (e,g)

    def fromPixelToLL(self,px,zoom):
         e = self.zc[zoom]
         f = (px[0] - e[0])/self.Bc[zoom]
         g = (px[1] - e[1])/-self.Cc[zoom]
         h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
         return (f,h)

if __name__ == "__main__":

  # read environment variables
  dt_file_name = os.environ['DTF']
  try:
    temp_file = "/dev/shm/tile" + os.environ['PID']
  except:
    temp_file = "/dev/shm/tile"

  # do some global initialization
  print "Rendering " + dt_file_name + ": started."
  mapfile = "mapnik_pt.xml"
  tile_dir = "/var/www/html/tiles/"
  max_zoom = 18
  tile_count = 0
  empty_tile_count = 0

  # create tile directory and all possible zoom directories
  if not os.path.isdir(tile_dir):
    os.mkdir(tile_dir)
  for z in range(0, max_zoom):
    d = tile_dir + "%s" % z + "/"
    if not os.path.isdir(d):
      os.mkdir(d)

  # open the dirty-tiles file
  try:
    dt_file = file(dt_file_name, "r")
  except:
    dt_file = None
  if (dt_file != None):

    # do some Mapnik initialization
    mm = mapnik.Map(256, 256)
    mapnik.load_map(mm, mapfile, True)
    # Obtain <Map> projection
    prj = mapnik.Projection(mm.srs)
    # Projects between tile pixel co-ordinates and LatLong (EPSG:4326)
    tileproj = GoogleProjection(max_zoom + 1)

    # process the dirty-tiles file
    while True:

      # read a line of the dirty-tiles file
      line = dt_file.readline()
      if (line == ""):
        break
      line_part = line.strip().split('/')

      # create x coordinate's directory - if necessary
      d = tile_dir + line_part[0] + "/" + line_part[1] + "/"
      if not os.path.isdir(d):
        os.mkdir(d)

      # get parameters of the line
      z = int(line_part[0])
      x = int(line_part[1])
      y = int(line_part[2])
      tile_name= tile_dir + line[:-1] + ".png"

      # render this tile -- start

      # Calculate pixel positions of bottom-left & top-right
      p0 = (x * 256, (y + 1) * 256)
      p1 = ((x + 1) * 256, y * 256)

      # Convert to LatLong (EPSG:4326)
      l0 = tileproj.fromPixelToLL(p0, z);
      l1 = tileproj.fromPixelToLL(p1, z);

      # Convert to map projection (e.g. mercator co-ords EPSG:900913)
      c0 = prj.forward(mapnik.Coord(l0[0], l0[1]))
      c1 = prj.forward(mapnik.Coord(l1[0], l1[1]))

      # Bounding box for the tile
      if hasattr(mapnik,'mapnik_version') and mapnik.mapnik_version() >= 800:
        bbox = mapnik.Box2d(c0.x,c0.y, c1.x,c1.y)
      else:
      	bbox = mapnik.Envelope(c0.x,c0.y, c1.x,c1.y)
      render_size = 256
      mm.resize(render_size, render_size)
      mm.zoom_to_box(bbox)
      mm.buffer_size = 256  # (buffer size was 128)

      # render image with default Agg renderer
      im = mapnik.Image(render_size, render_size)
      mapnik.render(mm, im)
      im.save(temp_file, 'png256')
      tile_count = tile_count + 1

      # render this tile -- end

      # copy this tile to tile tree
      l= os.stat(temp_file)[6]
      if l>116:
        shutil.copyfile(temp_file,tile_name)
      else:
        try:
          os.unlink(tile_name)
        except:
          l= l  # (file did not exist)
        empty_tile_count = empty_tile_count + 1

    # close and delete the dirty-tiles file
    dt_file.close()
    try:
      os.unlink(dt_file_name)
    except:
      l= l  # (file did not exist)

  # print some statistics
  print "Rendering " + dt_file_name + ":", tile_count, "tiles (" + "%s" % empty_tile_count + " empty)."

This Python file must be made executable:

chmod +x mapnik_pt.py

Omitting Empty Tiles

Empty tiles on a transparent map layer are completely transparent graphics and not of any use for our map, they just waste disk space. Therefore it is better to refrain from saving them in the first place. This will be accomplished by a few additional lines in the Mapnik Python script (see above).

OpenLayers must be advised what to do if an tile image is to be loaded which does not exist on the server. We already did the necessary changes to the file index.html (see above, look for .olImageLoadError {display: none;}).

Schematic Diagram

Here is an overview of the programs and files we are using. This diagram is somewhat simplified; it has been created to show all main steps of the strategy in one picture.

                   hourly.osc
                        |
                        v
                  +-----------+
                  | osmupdate |
      a.o5m ----> | -B=p.poly | ----> b.o5m
        |         +-----------+         |
        v               ^               v
+---------------+       |       +---------------+
|   osmfilter   |     p.poly    |   osmfilter   |
| --keep=route= |               | --keep=route= |
+---------------+               +---------------+
        |                               |
        v         +------------+        v
      fa.o5m ---> | osmconvert | <--- fb.o5m
                  |   --diff   |
                  +------------+
                        |
                        v
                     gis.osc
                        |
                        v
                  +------------+
                  | osm2pgsql  | ---> database
                  | -e 5-12 -a |      update
                  +------------+
                        |
                        v
                   dirty_tiles
                        |
                        v
                 +--------------+
                 | mapnik_pt.py | ---> new tiles
                 +--------------+
Involved Files
  • hourly.osc: hourly OsmChange file, downloaded from planet.openstreetmap.org
  • p.poly: border polygon of our graphical region
  • a.o5m: previous OSM data file (all data of our geographical region)
  • b.o5m: updated OSM data file (all data of our geographical region)
  • fa.o5m: previous filtered OSM data file (only thematic data, e.g. public transport)
  • fb.o5m: filtered new OSM data file (only thematic data, e.g. public transport)
  • gis.osc: the OsmChange file (so-called diff file) you would need to update fa.o5m to fb.o5m
  • dirty_tiles: a list of all tiles which are affected by the gis.osc and therefore have to be rerendered
Used Programs
  • osmupdate: download .osc files and use them to update an .o5m file
  • osmfilter: filter OSM data and discard all the information we do not need
  • osmconvert: compare the difference between two .o5m files and create an .osc diff file
  • osm2pgsql: update a postgreSQL geo database and calculate a dirty-tiles list
  • mapnik_pt.py: read the dirty-tiles file and (re)render all listed tiles

Most of the task also can be done – maybe more reliable – with Osmosis. If you already have installed Osmosis, feel free to use it. Osmosis offers a much wider variety of functions but will be slower in processing the data.

Toolchain Script

The following toolchain script will link all the previously introduced steps. Create a file named toolchain_pt.sh and put it into main directory /home/pt/pt:

#!/bin/bash
# toolchain_pt.sh
# (c) Markus Weber, 2017-02-23 02:50, License: AGPLv3
#
# This script cares about creating an overlay map.
# It loads and updates the map data and renders the tiles.
# The script will run endless in a loop.
# Start it with "nohup ./toolchain_pt.sh &"
# To terminate it, delete the file "toolchain_pt_running.txt".

PLANETURL=http://planet.openstreetmap.org/pbf/planet-latest.osm.pbf
PLANETMINSIZE=10000000000  # minimum size of OSM data file in .o5m format
BORDERS=
# ALTERNATE MAP DATA SOURCE (will overwrite previous settings):
PLANETURL=http://download.geofabrik.de/europe/germany/bayern/mittelfranken-latest.osm.pbf
PLANETMINSIZE=40000000  # minimum size of OSM data file in .o5m format
BORDERS="-B=mittelfranken.poly"
#
MAXPROCESS=2  # maximum number of concurrent processes for rendering
MAXMERGE=5  # number of files to be merged in parallel by osmupdate
OSM2PGSQLPARAM="-s -C 1500 -d ptgis -U ptuser -S osm2pgsql_pt.style"
  # main parameters to be passed to osm2pgsql

# enter working directory and do some initializations
cd /home/pt/pt
echo >>tc.log
echo $(date)"  toolchain started." >>tc.log
PROCN=1000  # rendering-process number (range 1000..1999)
mkdir d_t 2>/dev/null

# publish that this script is now running
rm toolchain_pt_ended.txt 2>/dev/null
echo -e "The toolchain script has been started at "$(date)"\n"\
"and is currently running. To terminate it, delete this file and\n"\
"wait until the file \"toolchain_pt_ended.txt\" has been created.\n"\
"This may take some minutes." \
>toolchain_pt_running.txt

# clean up previous Mapnik processes
killall "mapnik_pt.py" 2>/dev/null
while [ $(ls d_t/at* 2>/dev/null |wc -l) -gt 0 ]; do
    # there is at least one incompleted rendering process
  AT=$(ls -1 d_t/at* 2>/dev/null|head -1)  # tile list
  echo $(date)"  Cleaning up incomplete rendering: "$AT >>tc.log
  DT=${AT:0:4}d${AT:5:99}  # new name of the tile list
  mv $AT $DT  # rename this tile list file;
    # now the tile is is marked as 'to be rendered'
  done

# download and process planet file - if necessary
if [ "0"$(stat --print %s a.o5m 2>/dev/null) -lt $PLANETMINSIZE ]; then
  echo $(date)"  Missing file a.o5m, downloading it." >>tc.log
  rm -f fa.o5m
  wget -nv $PLANETURL -O - 2>>tc.log |./osmconvert - \
    $BORDERS -o=a.o5m 2>>tc.log
  echo $(date)"  Updating the downloaded planet file." >>tc.log
  rm -f b.o5m
  ./osmupdate -v a.o5m b.o5m --day --hour \
    --max-merge=$MAXMERGE $BORDERS --drop-author >>tc.log 2>&1
  if [ "0"$(stat --print %s b.o5m 2>/dev/null) -gt $PLANETMINSIZE ]; then
    echo $(date)"  Update was successful." >>tc.log
  else
    echo $(date)"  No update available. Dropping meta data." >>tc.log
    ./osmconvert -v a.o5m b.o5m --drop-author >>tc.log 2>&1
    fi
  mv -f b.o5m a.o5m
  if [ "0"$(stat --print %s a.o5m 2>/dev/null) -lt $PLANETMINSIZE ]; then
    echo $(date)"  toolchain Error: could not download"\
      $PLANETURL >>tc.log
    exit 1
    fi
  fi

# refill the database - if necessary
if [ "0"$(stat --print %s fa.o5m 2>/dev/null) -lt 5 ]; then
  echo $(date)"  Missing file fa.o5m, creating it." >>tc.log
  rm dirty_tiles d_t/* 2>/dev/null
  echo $(date)"  Filtering the downloaded planet file." >>tc.log
  ./osmfilter a.o5m --parameter-file=toolchain_pt.filter \
    -o=fa.o5m 2>>tc.log  # filter the planet file
  ./osmconvert fa.o5m -o=gis.o5m 2>>tc.log
    # convert to .o5m format
  echo $(date)"  Writing filtered data into the database." >>tc.log
  rm dirty_tiles 2>/dev/null
  osm2pgsql $OSM2PGSQLPARAM -c gis.o5m -e 4-17 >/dev/null 2>&1
    # enter filtered planet data into the database
  echo $(date)"  All tiles need to be rerendered." >>tc.log
  echo $(date)"  If the tile directory is not empty, please remove" >>tc.log
  echo $(date)"  all outdated tile files." >>tc.log
  fi

# main loop
while [ -e "toolchain_pt_running.txt" ]; do
  echo $(date)"  Processing main loop." >>tc.log

  # limit log file size
  if [ "0"$(stat --print %s tc.log 2>/dev/null) -gt 6000000 ]; then
    echo $(date)"  Reducing logfile size." >>tc.log
    mv tc.log tc.log_temp
    tail -c +4000000 tc.log_temp |tail -n +2 >tc.log
    rm tc.log_temp 2>/dev/null
    fi

  #care about entries in dirty-tile list
  while [ $(ls dirty_tiles d_t/?t* 2>/dev/null |wc -l) -gt 0 -a \
      -e "toolchain_pt_running.txt" ]; do
      # while still tiles to render

    # start as much rendering processes as allowed
    while [ $(ls d_t/dt* 2>/dev/null |wc -l) -gt 0 -a \
        $(ls -1 d_t/at* 2>/dev/null |wc -l) -lt $MAXPROCESS ]; do
        # while dirty tiles in list AND process slot(s) available
      DT=$(ls -1 d_t/dt* |head -1)  # tile list
      AT=${DT:0:4}a${DT:5:99}  # new name of the tile list
      touch -c $DT  # remember rendering start time
      mv $DT $AT  # rename this tile list file;
        # this is our way to mark a tile list as 'being rendered now'
      #echo $(date)"  Rendering "$DT >>tc.log
      PROCN=$(($PROCN + 1))
      if [ $PROCN -gt 1999 ]; then
        PROCN=1000;  # (force range to 1000..1999)
        fi
      N=${PROCN:1:99}  # get last 3 digits
      DTF=$AT PID=$N nohup ./mapnik_pt.py >/dev/null 2>&1 &
        # render every tile in list
      echo $(date)"  Now rendering:"\
        $(ls -m d_t/at* 2>/dev/null|tr -d "d_t/a ") \
        >>tc.log
      sleep 2  # wait a bit
      tail -30 tc.log >/var/www/html/status1
      done

    # determine if we have rendered all tiles of all lists
    if (ls d_t/?t* >/dev/null 2>&1); then  # still tiles to render
      sleep 11  # wait some seconds
  continue  # care about rendering again
      fi
    # here: we have rendered all tiles of all lists

    # care about dirty-tiles master list and split it into parts
    if (ls dirty_tiles >/dev/null 2>&1); then
        # there is a dirty-tiles master list
      echo $(date)"  Expanding \"dirty_tiles\" file." >>tc.log
      ./dtexpand 4 17 <dirty_tiles >dirty_tiles_ex 2>/dev/null
      mv -f dirty_tiles_ex dirty_tiles 2>/dev/null
      echo $(date)"  Splitting dirty-tiles list" \
        "("$(cat dirty_tiles |wc -l)" tiles)" >>tc.log
      grep --color=never -e "tiles)" -e "newest hourly timestamp" tc.log >/var/www/html/status
      split -l 1000 -d -a 6 dirty_tiles d_t/dt
      echo "*** "$(date) >>dt.log
      cat dirty_tiles >>dt.log  # add list to dirty-tiles log
      rm dirty_tiles 2>/dev/null
      # limit dirty-tiles log file size
      if [ "0"$(stat --print %s dt.log 2>/dev/null) -gt 750000000 ]; then
        echo $(date)"  Reducing dirty-tiles logfile size." >>tc.log
        mv dt.log dt.log_temp
        tail -c +500000000 dt.log_temp |tail -n +2 >dt.log
        rm dt.log_temp 2>/dev/null
        fi
      fi

    done  # while still tiles to render
  if [ ! -e "toolchain_pt_running.txt" ]; then  # script shall be terminated
continue;  # exit the main loop via while statement
    fi
  # here: all tiles have been rendered

  # update the local planet file
  echo $(date)"  Updating the local planet file." >>tc.log
  rm b.o5m fb.o5m 2>/dev/null
  # for maintenance only:
  #echo $(date)"---> Waiting 24 hours." >>tc.log
  #sleep 86400
  #echo $(date)"---> Resuming work after wait." >>tc.log
  ./osmupdate -v a.o5m b.o5m --day --hour \
    --max-merge=$MAXMERGE $BORDERS --drop-author >>tc.log 2>&1

  if [ "0"$(stat --print %s b.o5m 2>/dev/null) -lt \
      $(expr "0"$(stat --print %s a.o5m 2>/dev/null) \* 9 / 10) ]; then
      # if new osm file smaller than 90% of the old file's length
    # wait a certain time and retry the update
    I=0
    while [ $I -lt 33 -a -e "toolchain_pt_running.txt" ]; do  # (33 min)
      sleep 60
      I=$(( $I + 1 ))
      done
continue
    fi

  # filter the new planet file
  echo $(date)"  Filtering the new planet file." >>tc.log
  ./osmfilter b.o5m --parameter-file=toolchain_pt.filter -o=fb.o5m \
    2>>tc.log

  # calculate difference between old and new filtered file
  echo $(date)"  Calculating diffs between old and new filtered file."\
    >>tc.log
  ./osmconvert fa.o5m fb.o5m --diff-contents --fake-lonlat \
    --out-osc|gzip -1 >gis.osc.gz 2>>tc.log

  # check if the diff file is too small
  if [ "0"$(stat --print %s gis.osc.gz 2>/dev/null) -lt 10 ]; then
    echo $(date)"  Error: diff file is too small." >>tc.log
    exit 2
    fi

  # enter differences into the database
  echo $(date)"  Writing differential data into the database." >>tc.log
  osm2pgsql $OSM2PGSQLPARAM -a gis.osc.gz -e 4-17 >/dev/null 2>&1

  # replace old files by the new ones
  echo $(date)"  Replacing old files by the new ones." >>tc.log
  ls -lL a.o5m fa.o5m b.o5m fb.o5m gis.o5m gis.osc.gz >>tc.log
  mv -f a.o5m a_old.o5m
  mv -f fa.o5m fa_old.o5m
  mv -f b.o5m a.o5m
  mv -f fb.o5m fa.o5m

  # check if there are any tiles affected by this update
  if [ "0"$(stat --print %s dirty_tiles 2>/dev/null) -lt 5 ]; then
    echo $(date)"  There are no tiles affected by this update." >>tc.log
    rm -f dirty_tiles
    fi

  # wait a certain time to give other maps a chance for rendering
  I=0
  while [ $I -lt 30 -a -e "toolchain_pt_running.txt" ]; do  # (30 min)
    sleep 60
    I=$(( $I + 1 ))
    done

  done  # main loop

# wait until every rendering process has terminated
while (ls d_t/at* 2>/dev/null) ; do sleep 30; done

# publish that this script has ended
rm toolchain_pt_running.txt 2>/dev/null
echo "The toolchain script ended at "$(date)"." \
  >toolchain_pt_ended.txt

As usual, the file must be made executable:

cd /home/pt/pt
chmod +x toolchain_pt.sh

There is one border polygon file you will need if you run the script for Mittelfranken region (thats the region where Nürnberg is in). Download that polygon file and then start the script as background process:

cd /home/pt/pt
wget http://download.geofabrik.de/europe/germany/bayern/mittelfranken.poly
nohup ./toolchain_pt.sh &

You can watch the rendering process(es) by invoking the process monitor top:

top -u pt

Other Strategies

Mapnik experts will have noticed that the rendering strategy which has been introduced on this page does not represent the standard solution. Most Slippy Map installations will use the more comfortable Tirex/Modtile tools to define which tiles when to render. Furthermore, of course, you usually do not render regularly sized tiles but so-called meta tiles instead, which are 16 or 64 times larger. On standard maps this increases rendering speed significantly. It would be interesting to know if there are advantages for diff-based rendered thematic maps too because you would have to render a lot of areas which have not changed.

Please feel free to add your experiences here if you have tried different methods and had the opportunity to compare them.

Benchmarks

Be prepared that in the first run every affected tile of your thematic layer will be rendered. This may take between one hour and a few weeks, depending on the available hardware, the size of the chosen geographical region and the density of the thematic data.

The planet-wide public transport map, for example, took a bit more than two days to be initially rendered on a quad core CPU. The statistics in detail:

  • 35 minutes downloading and converting the planet file
  • 15 minutes initially updating the planet file
  • 15 minutes filtering the planet file
  • 15 minutes writing the filtered data into the database
  • 50 hours rendering all tiles of the public transport overlay (18 Mio. tiles)

Updates run faster. An average update takes less than a day. In detail (example):

  • 15 minutes updating the planet file
  • 15 minutes filtering the updated planet file
  • 10 seconds calculating the differences between old and new filtered data
  • 20 minutes writing the differences of the filtered data into the database
  • 12 hours rerendering the modified tiles of the public transport overlay

Troubleshooting

Too many links

If this message occurs during rendering, you most likely are using a file system which supports only a limited number of subdirectories (e.g. ext2 or ext3). Please upgrade to ext4. The migration from ext3 to ext4 should be possible without any loss of data. Please refer to the users guide of your operating system.