Shaded relief maps using mapnik

From OpenStreetMap Wiki
(Redirected from HikingBikingMaps)
Jump to navigation Jump to search
Shaded relief maps using mapnik - Other languages

broom

This article or section may contain out-of-date information. The information may no longer be correct, or may no longer have relevance.
If you know about the current state of affairs, please help keep everyone informed by updating this information. (Discussion)

This page documents set of steps to create a relief map with hillshading, combined with OpenStreetMap data (particularly hiking and biking routes) using the Mapnik renderer. Hillshading is a nice cartographic technique for showing relief, but if you prefer contour lines see Contour relief maps using mapnik.

The page was initially created by numenor in 2009, documenting his development of a hiking biking map. As with all wiki pages, please feel free to help keep it updated as a working set of instructions.

Tour at Lake Lipno

Required features

These features are implemented in the final maps. Some of these are already documented in other parts of the wiki (but summarized here as well, to provide a comprehensive view), while some had to be implemented newly during work on the maps.

  • Topographical map (elevation contours), see Contours
  • Hill shading (new)
  • Rendering of cycle routes and hiking routes. Cycle routes are already implemented in the cycle map, but the styles are not published, and hiking routes are not implemented.
  • Ways with highway=path should be rendered (new)

The image on the right shows the current status of rendering using the information on this page. The transparent red line is part of national cycle route 33 in the Czech Republic. The red-and-blue dotted line is a mixed foot and bicycle path.

The process is mapnik based, and is currently only tested on Linux (Ubuntu 08.04). Anybody having experiences about problems, success, or workarounds for other platforms is welcome to contribute!

Steps

Software setup

Please follow the information provided in the linked pages, to install the necessary software:

  • Mapnik (Mapnik, Mapnik Homepage)
  • PostgreSQL/PostGIS (Mapnik/PostGIS)
  • osm2pgsql (Osm2pgsql)
  • (not needed with recent GDAL) srtm_generate_hdr.sh (from [1])
  • GDAL library/tools (GDAL Homepage)
  • PerryGeo GDAL-based DEM utilities (Homepage)
  • Optional: osmosis (Osmosis), for sorting *.osm files
  • Optional: Quantum GIS (QGIS Homepage), for testing and working on raster data and possibly overlaying OSM data, without using mapnik rendering (e.g. to check whether projection system fits to OSM data)

Note: GDAL 1.7 and newer can read SRTM data files (.hgt) directly and also includes the PerryGeo DEM utilities.

Download and import of OSM data into postgis

The first step is to download a part of the planet as a *.osm-file. One possibility is to use the GeoFabrik server to get a country (at least for Europe). You do not need to decompress the *.bz2 format, as osm2pgsql can decompress it on the fly. Another possibility is to download and save a file using Josm (for smaller areas). Josm can also be useful to edit the data before rendering, to add markers, etc., which can get their own special styles in the mapnik style file.

If the data is not sorted by node/way/relation id (the Geofabrik data seems to be, but not the data saved by Josm), this has to be done, e.g. using Osmosis (bzip2 -d <osm.bz2> before, if necessary): 

java -Xmx1048m -jar <path to osmosis>/osmosis.jar --read-xml enableDateParsing=no file=<input.osm> --write-xml file=<sorted.osm>

For some styles used later on in mapnik style files, some tags are needed, which are not included by default using osm2pgsql. To have them available, the following lines have to be added to the default.style file of osm2pgsql: 

way       network      text     linear
way       surface      text     linear

Then, to import the sorted *.osm file you do 

osm2pgsql --slim -d <postgis db name> <sorted.osm>

The "--slim" is important (I think), as route rendering will not work without it. If you want to add several OSM files to the database, you can provide the flag "--append" to the second and each following call. It is also possible, to import only part of a file by providing the "--bbox <west,south,east,north>" flag.

If you want to visualize OSM data using Quantum GIS, the following lines are important, because a primary key is needed: 

psql -d <postgis db name> -c 'ALTER TABLE planet_osm_point ADD COLUMN pid serial;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_point ADD CONSTRAINT pid_pkey_point PRIMARY KEY (pid);'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_line ADD COLUMN pid serial;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_line ADD CONSTRAINT pid_pkey_line PRIMARY KEY (pid);'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_polygon ADD COLUMN pid serial;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_polygon ADD CONSTRAINT pid_pkey_polygon PRIMARY KEY (pid);'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_roads ADD COLUMN pid serial;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_roads ADD CONSTRAINT pid_pkey_roads PRIMARY KEY (pid);'

It is best, to put these into a script, as they have to be called after each new import. Actually, wait with executing these commands, until all OSM files are imported. Otherwise, before a new import, you have to undo these changes, or osm2pgsql will fail: 

psql -d <postgis db name> -c 'ALTER TABLE planet_osm_point DROP CONSTRAINT pid_pkey_point;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_point DROP COLUMN pid;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_line DROP CONSTRAINT pid_pkey_line;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_line DROP COLUMN pid;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_polygon DROP CONSTRAINT pid_pkey_polygon;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_polygon DROP COLUMN pid;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_roads DROP CONSTRAINT pid_pkey_roads;'
psql -d <postgis db name> -c 'ALTER TABLE planet_osm_roads DROP COLUMN pid;'

It might also be useful to add the SRS (spatial reference system) used by the OpenStreetMap slippy map (and by default by the mapnik rendering performed here) to the list of SRS known to PostGIS: 

psql -d <postgis db name> -c "INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) VALUES (900913, 'EPSG', 900913, 'PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6378137,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]', '+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m');"

Download and import/preparation of elevation data into postgis/for hill shading

The contour lines need data which is not provided in the OSM files, named DEM (digital elevation model). There are some sources for free (some only for non-commercial purposes) data, usually based on SRTM (Shuttle Radar Topography Mission) data by the NASA. This is, in part of the world, available in 1 arc-second resolution (30m at the equator), and in 3 arc-second resolution (90m at the equator) for the whole world (below 60° latitude). For obtaining such data, see Contours. A problem of this data in some areas is that they have gaps, where the data acquisition could not determine the elevation. This occurs e.g. in some regions with snow, such that these gaps occur relatively often in the Alps. An alternative are the data provided by CIAT/CSI, which used other available DEM sources or sophisticated interpolation algorithms to fill these gaps. Unfortunately, this data is available only for non-commercial purposes, and thus contradicts the OSM CC by sa license. It might be okay to use it for private purposes, but the resulting images/maps cannot be published.

First, one or more DEM tiles have to be downloaded for the region in question. Then two different processing steps are necessary, for importing contour levels into the PostGIS database, and for preparing the DEM for the hill shading step.

Contour Lines

This is mainly a summary of Contours, so please also have a look on that page. The available SRTM data files (*.hgt) have to be preprocessed using srtm_generate_hdr.sh. Then region in question has to be extracted from them, and saved to a new GeoTIFF file. This usually is necessary, because all available SRTM data is too large to be imported in one step. A region of the size of one Austrian state (Upper Austria for example) works well with 2Gb of memory. Additionally, the region in question might stretch over more than one DEM image, thus a merging is necessary: 

gdal_merge.py -v -o <srtm.tif> -ul_lr <west> <north> <east> <south> <tif-directory>/*.tif

This extracts the region between <west>-<east> and <north>-<south>, specified in degrees, from the GeoTIFF files in the <tif-directory> into the <srtm.tif>, and provides some processing information ("-v).

Note:
  • My version of gdal_merge (from ms4w) could not read files from a tif-directory, but required all files on the command line
  • If you repeat the process, make sure to delete your <srtm.tif> before running gdal_merge. If the file exists, the data is copied into it, but the GDAL information remains so other tools will use the old information.

Next, the contour lines are extracted into a Shape file: 

gdal_contour -i <step> -snodata 32767 -a height '<srtm.tif>' '<srtm.shp>'

For <step>, specify the elevation difference between adjacent contour lines, in the units used by the DEM. In this case, this is meters. The -snodata option specifies the value used in the image tiles to denote missing elevation values (which is -32768 for current CIAT/CSI data). Now, the contour lines can be imported into the PostGIS database: 

shp2pgsql <-c/-a/-d> -I -g way '<srtm-file without shp extension>' contours | psql -q <postgis db name>

Depending on wheter it is the first import (-c), the contours should be added to existing contour data (-a), or existing contour data should be replaced (-d), different flags have to be specified. This creates a new table (or adds to / replaces the existing table) contours.

In some cases, I noted, that there seems to be either some inaccuracy in the location of the SRTM data, or a "misunderstanding" between the provider of the SRTM data (I noticed such for the CIAT/CSI data) and gdal_contour, because the data imported this way seems to have a certain slight offset (e.g. when comparing known peak or river locations with the contours). A reason for the latter case could be, that the SRTM files could specify the elevation either as average of the pixel area (so, the reference should be the center of the pixel), or at certain grid locations. gdal_contour tool seems to assume the reference location of the elevation to be the center of each pixel.

Whatever the reason, if there is an offset between OSM data and DEM data, the DEM image data can be shifted slightly: 

psql -d <postgis db name> -c 'UPDATE contours SET way = ST_Translate(way, 0.000416666666666, -0.000416666666666, 0);'

In this example the given translation distances are half the width/height of a pixel. In current versions of CIAT/CSI data (i.e. 4.1) they seem to have changed this "misunderstanding". When working with a GeoTIFF file created by gdal_merge.py from the CIAT GeoTIFF files, 

psql -d <postgis db name> -c 'UPDATE contours SET way = ST_Translate(way, -0.000416666666666, -0.000416666666666, 0);'

fixes the contour lines. Another possibility is changing the GeoTIFF metadata: After extracting the data with 

listgeo srtm.tif > srtm.hdr

you can change srtm.hdr: In the ModelTiepointTag section, subtract 0.000416666666666 from the x and y coordinate of the corner (one line). In the "Corner Coordinates" section, shift every point south by 0d00'01.50" and west by 0d00'01.50". After writing the metadata to a new GeoTIFF file using 

geotifcp -g srtm.hdr srtm.tif srtm_shift.tif

and verifying the change using listgeo, this problem should disappear.

Hill Shading

Hill shading works completely different. It uses mapnik's ability to underlay raster images below the rendered lines, points and text. One important point to observe is, that mapnik's support for raster images so far is limited to images which already have the SRS of the final image to generate (it does not reproject or warp raster images). Another limitation, for which there currently seems to be no workaround, is that raster images cannot be rendered transparently, as opposed to polygon or line data. This restricts the hill shading images to be the lowest layer, and the only raster layer possible to show (except if it should be overlaid/hidden by smaller other images). But more on that later, in the mapnik styles section.

First, if necessary, the needed dem *.tif files should again be merged, using gdal_merge.py (see above). Then, the data has to be warped to the "spherical Mercator" projection used by OpenStreetMap rendering by default. Sometimes (e.g. for the CIAT/CSI data), I had problems directly warping the files using gdalwarp, because the SRS in the DEM data seems not to be specified, so it might be necessary to first add it to the file which should be converted into a hill shading grey scale image: 

gdal_translate -of GTiff -co "TILED=YES" -a_srs "+proj=latlong" <input.tif> <input_adapted.tif>

If desired, this step can also be used to extract only a subregion of the original image (-projwin option). The next step is to actually warp the image into the spherical Mercator projection: 

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 <input_adapted.tif> <warped.tif>

This specifies the best interpolation method for warping (-rcs -order 3, this can take a while if the region is large) and a resolution of the target image having 30m per pixel horizontally and vertically. This is higher than the original resolution, but because of the good interpolation, it looks good (at least much better than with less resolution), even for relatively high zoom levels.

You can make it look even better (again at the expense of computing time and resulting file size) by something like these options: -tr 15 15 -wt Float32 -ot Float32 -wo SAMPLE_STEPS=100

Now the hill shading can take place, using the PerryGeo tool: 

hillshade <warped.tif> <hillshade.tif> -z 2

The "-z 2" slightly exaggerates the slopes, leading to an image of better contrast. This looks good after rendering, because often there are partially transparent layers which are put on top of the hill shading image, thus reducing its contrast again. There are also a few other options available for the hillshade tool, which allow you to adapt for different scales between x-/y-coordinates and elevation units, and to adapt direction and altitude of the simulated light source. Alternatively, the color-relief tool from the same toolkit could be used to produce an image colored for the elevation of each pixel (like the background in the cycle map).

If you see artefacts in your hillshaded image that look like plateaus or similar to contour lines, try adding the parameters "-wt Float32 -ot Float32" to your gdalwarp call. This did remove the artefacts for me.

Style adaptations for Mapnik rendering

The following suggestions assume that you start from a fresh personalized copy of osm.xml, see "Rendering with mapnik" at Mapnik. Where you encounter occurrences of svnmapnik in the following code, replace it by the absolute path to a copy of the mapnik subversion directory (https://svn.openstreetmap.org/applications/rendering/mapnik).

Contours

Possible mapnik style adaptations for rendering contours are already described on the aforementioned page Contours. The suggestions given there work very well. Although visualization of 10m contours for a larger area can be by far the most time consuming part of rendering, areas of about 20x20km should complete in more or less a minute rendering time (on a 2-year old laptop with 2GB of RAM).

Hill Shading and other Raster Images

When the process is known, and the necessary preprocessing is completed, the actual rendering of raster images with is quite easy. Some points should be kept in mind, though:

  • The SRS of the image has to match the projection of the map produced by mapnik, which is spherical Mercator by default (i.e. the standard stylesheet in the mapnik directory of the SVN repository).
  • Currently, there seems to be no possibility to render raster images transparently, so usually they will be the lowest layer, or they will be used only in small areas, where they do not cover other important information.
  • The resolution of the raster images should be appropriate for the image resolution and extent. The gdalwarp step mentioned above can be adapted to the final resolution required. At least for hill shading, even a well interpolated (e.g. cubic spline sampling, -rcs option) image constructed from low resolution looks a lot better than the original low resolution image.

The changes for actually including the raster image are given at HikingBikingMaps/HillShading.

Unfortunately the result of all these changes is not yet perfect. If you have a very close look at the image at the top of this page, you will notice a soft yellowish horizontal bar in the lower part of the image. I am not sure, where this originates from, but I suspect one or more of the background coloring layers/styles (world, world-1, coast-poly). These effects occur at regular horizontal and vertical distances. Maybe someone can point out, what to change to get rid of these.

Those bars are from the shoreline shapefiles. The polygons from those shapefiles overlap, so when they're partially transparent the overlapped regions get darker. I solved the problem by putting the raster images on top of the world and coast-poly layers. --Asciiphil 16:00, 13 May 2010 (UTC)

Routes

The design of the following styles for rendering hiking and bicycle routes is inspired by the cycle map. There are separate styles for transparently rendering the ways of the route, for placing the name of the route, and for placing "shields", i.e. the ref tag. The way rendering is repeated for two different zoom level ranges (scale denominators). And all of that is repeated six time, for routes with the network=ncn, network=rcn, network=lcn (national, regional and local cycle networks), and for hiking routes with network=nwn, network=rwn, and network=lwn (national, regional, local walking networks). So far this scheme does not support individual coloring by colors specified as tags (does someone know how to do so?). But it is at least possible, to include a new set of such rendering rules for each route which should get its individual color, and use a filter to apply the rules just to that route.

In the mapnik svn directory there is a tool (perl script) mkshield.pl, which can be used to produce an individually colored set of shields for the routes (it has to be edited, but existing code can be used as a template). Such shields with colors matching the way coloring are used in the following.

The changes for actually rendering routes are given at HikingBikingMaps/RouteRendering.

Path rendering

Technically, path rendering is not very difficult (compared to the expamples already shown above). The challenge is in finding a good way to show the different variations of paths. There are highway=paths with foot=designated or foot=yes, with bicycle=designated or bicycle=yes, and horses can also be allowed. There are existing styles for highway=footway and highway=bicycle, which should be copied by ways with equivalent path tagging. And in future I plan to incorporate specific rendering of paths which have a sac_scale tag.

The following is my suggestion, but by no means yet perfect (if that could ever be the case). It resembles the rendering of the traditional footway and cycleway, when the path has the additional tag surface=paved. The reasoning here is, that path is used for all kinds of paths, most of which probably are not paved. Accordingly, paths without further hints about their surface quality should be rendered slightly less prominent than the existing paths (as thin lines). Paved paths, which have both foot and bicycle tagged with yes or designated are renderd as a combination with alternating salmon and blue dots. Other paths are colored blue (bicycle), salmon (foot), or green (horse), or dashed alternatingly when a combination of these "transportation modes" is allowed. So far, rendering of tunnels and bridges for paths is not implemented, but "left as an exercise to the reader" ...

The changes for actually rendering paths are given at HikingBikingMaps/PathRendering, in addition to a sample image showing how paths are rendered using these styles.

Related

Retrieved from "https://wiki.openstreetmap.org/w/index.php?title=Shaded_relief_maps_using_mapnik&oldid=2238217"