User:NicolasDumoulin/PlanStatique/DICRIM

I've been solicitated by my municipality council for building a map that shows major risks on the administrative entity for illustrating a legal document. I've used mapnik with a python script and a postgis database with OSM data and some subjectives data.

How-to

Software and data requirements

• Install a qemu virtual debian squeeze that will contain the postgis database (french tutorial)
• Install osm2pgsql and mapnik
• Download data of the region from geofabrik and import data into the database (named "osm" in my scripts)

I've developped a little python module with some helpfull functions for having a nicer script.

Subjective datas

The map need some subjective datas for representing the risk areas. These additionnal data had been stored in a separated database, named "dicrim" in the script.

Legend

I've used inkscape for the legend, but it should not be so hard to generate from a selection of layers using cairo.

The script

#!/usr/bin/python
# -*- coding: utf-8 -*-

import mapnik
import sys, os
# personnal module with some helper functions
from easymap.easymap import *

if __name__ == "__main__":
    ll = (3.098,45.739,3.1543,45.764)
    srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over"
    prj = mapnik.Projection(srs)
    c0 = prj.forward(mapnik.Coord(3.098, 45.739))
    c1 = prj.forward(mapnik.Coord(3.1543, 45.764))
    z = 8
    imgx = 542 * z
    imgy = 360 * z
    m = mapnik.Map(imgx,imgy)
    m.background = mapnik.Color('white')
    m.srs = srs

    # SQL for ways inside the administrative boundary
    sqlWays="(select l.way,l.highway,l.name,l.waterway,l.tunnel from planet_osm_polygon as c join planet_osm_line as l on (ST_Contains(c.way, l.way) OR ST_Intersects(c.way, l.way)) where c.name='Aubière')"
    # SQL for polygons inside the administrative boundary
    sqlPolygons="(select p.way,p.building from planet_osm_polygon as c join planet_osm_polygon as p on ST_Contains(c.way, p.way) where c.name='Aubière')"
    # buildings
    addPGLayer(m, 'buildings', "(select way from "+sqlPolygons+" as polygons where building='yes') as roads", [line('#333', 2)])
    addPGLayer(m, 'buildings-fill', "(select way from "+sqlPolygons+" as polygons where building='yes') as roads", [polygon('#444')])
    # river
    addLinesLayer(m, 'river-tunnel', "(select way from "+sqlWays+" as ways where waterway='river' and tunnel in ('yes','true','1')) as roads", '#D5EDFF',17)
    addLinesLayer(m, 'river', "(select way from "+sqlWays+" as ways where waterway='river' and (tunnel is null or tunnel not in ('yes','true','1'))) as roads", '#57BFFF',17)
    # highway
    addLinesLayer(m, 'service-ways', "(select way from "+sqlWays+" as ways where highway in ('service')) as roads", '#AAA',3)
    addLinesLayer(m, 'minor-ways', "(select way from planet_osm_line as ways where highway in ('footway','path','track')) as roads", '#DDD',4)
    addLinesLayer(m, 'classic-ways', "(select way from planet_osm_line as ways where highway in ('residential','unclassified','road')) as roads", '#666',3)
    addLinesLayerWithContour(m, 'major-ways', "(select way from planet_osm_line as ways where highway in ('primary','secondary','tertiary','trunk', 'trunk_link', 'motorway', 'motorway_link')) as roads", '#FFAB35',25,'#674515',3)
    # schrub area
    addPGLayer(m, 'puy', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='scrub') as roads",[line('#7BA673', 5)],"dicrim")
    addPGLayer(m, 'puy-fill', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='scrub') as roads",[polygon('#49CD2B', 0.3)],"dicrim")
    addPGLayer(m, 'puy-text', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='scrub') as roads",[text(42,'black','white',4)],"dicrim")
    # collasing risk area
    addPGLayer(m, 'cave', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='cave') as roads",[line('#672F00', 5)],"dicrim")
    addPGLayer(m, 'cave-fill', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='cave') as roads",[polygon('#C05800', 0.4)],"dicrim")
    addPGLayer(m, 'cave-text', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='cave' AND name='Les Grandes Caves') as roads",[text(42,'black','white',8)],"dicrim")
    ts = text(42,'black','white',8)
    ts.displacement(20,-110) # shift the text because the area is tiny and is invisible if covered by the label
    addPGLayer(m, 'cave-text-2', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='cave' AND name<>'Les Grandes Caves') as roads",[ts],"dicrim")
    # administrative boundary
    addLinesLayer(m, 'admin', "(select way from planet_osm_polygon where name='Aubière') as roads", '#BB2200',14)
    # highway labels
    addPGLayer(m, 'major-ways-text', "(select way,name from planet_osm_line as ways where highway in ('primary','secondary','tertiary') and name<>'Rue de Romagnat') as roads", [text(36, 'black', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)])
    # two street name are shifted because the street are covered by flooding area
    ts = text(36, 'black', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)
    ts.displacement(20,-20)
    addPGLayer(m, 'rue-romagnat-text', "(select way,name from planet_osm_line as ways where name='Rue de Romagnat') as roads", [ts])
    ts.displacement(0,-20)
    addPGLayer(m, 'rue-gergovie-text', "(select way,name from planet_osm_line as ways where name='Rue de Gergovie') as roads", [ts], 'dicrim')
    # river label
    addPGLayer(m, 'river-text', "(select way, name from planet_osm_line where waterway='river') as roads", [text(36, 'blue', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)])
    # administrative boundary clip
    bb_geometry="ST_GeomFromText('POLYGON(("
    for c in [[c0.x,c0.y],[c0.x,c1.y],[c1.x,c1.y],[c1.x,c0.y],[c0.x,c0.y]]:
      bb_geometry = bb_geometry + str(c[0]) + " "+str(c[1]) + ","
    bb_geometry = bb_geometry[:-1] + "))', 900913)"
    addPGLayer(m, 'admin-clip', "(select ST_BUFFER(ST_DIFFERENCE(ST_BUFFER("+bb_geometry+",100),way),-12) as way from planet_osm_polygon where name='Aubière') as roads",[polygon('#FFFFFF',0.7)])
    # motorway labels
    addPGLayer(m, 'motorways-text', "(select way,ref as name from planet_osm_line where highway='motorway') as roads", [text(36, 'black', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)])
    addPGLayer(m, 'motorways-text-dicrim', "(select way,name from planet_osm_line where highway in ('motorway', 'trunk')) as roads", [text(36, 'black', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)],"dicrim")
    addPGLayer(m, 'motorways-text-dicrim-ref', "(select way,ref as name from planet_osm_line where highway in ('motorway', 'trunk')) as roads", [text(36, 'black', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)],"dicrim")
    # an approximative positioning for a buried water stream
    addPGLayer(m, 'river-text-dicrim', "(select way, name from planet_osm_line where waterway='river') as roads", [text(36, 'blue', 'white', 4, mapnik.label_placement.LINE_PLACEMENT, True)], 'dicrim')
    # a neighbour town
    addPGLayer(m, 'romagnat-text', "(select way,name from planet_osm_polygon where planet_osm_polygon.natural='village') as roads",[text(52,'black','white',4)],"dicrim")

    m.zoom_to_box(mapnik.Envelope(c0.x,c0.y,c1.x,c1.y))
    # flooding area from shapefiles
    m_zi = mapnik.Map(imgx,imgy)
    m_zi.srs = srs
    addLayer(m_zi,"zone-innondable",[polygon('#60CFFF', 0.5)], mapnik.Shapefile(file="zone_innondable/zone_innondable"),projection="+proj=lcc +lat_1=45.898918964419 +lat_2=47.696014502038 +lat_0=46.8 +lon_0=0 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +towgs84=-168,-60,320,-0,-0,-0,0 +pm=2.337229166667 +units=m +no_defs")
    m_zi.zoom_to_box(mapnik.Envelope(c0.x,c0.y,c1.x,c1.y))
    
    cairoDraw([m, m_zi], "aubiere_dicrim.png", imgx, imgy)