Archive

Cartography

Using OSM data in QGIS is a hot topic but so far, no best practices for downloading, preprocessing and styling the data have been established. There are many potential solutions with all their advantages and disadvantages. To give you a place to start, I thought I’d share a workflow which works for me to create maps like the following one from nothing but OSM:

osm_google_100k

Getting the data

Raw OSM files can be quite huge. That’s why it’s definitely preferable to download the compressed binary .pbf format instead of the XML .osm format.

As a download source, I’d recommend Geofabrik. The area in the example used in this post is part of the region Pays de la Loire, France.

Preparing the data for QGIS

In the preprocessing step, we will extract our area of interest and convert the .pbf into a spatialite database which can be used directly in QGIS.

This can be done in one step using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat 2.59 46.58 -1.44 47.07 noirmoutier.db noirmoutier.pbf

where the -spat option controls the area of interest to be extracted.

When I first published this post, I suggested a two step approach. You can find it here for future reference:

For the first step: extracting the area of interest, we need Osmosis. (For Windows, you can get osmosis from openstreetmap.org. Unpack to use. Requires Java.)

When you have Osmosis ready, we can extract the area of interest to the .osm format:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>..\bin\osmosis.bat --read-pbf pays-de-la-loire-latest.osm.pbf --bounding-box left=-2.59 bottom=46.58 right=-1.44 top=47.07 --write-xml noirmoutier.osm

While QGIS can also load .osm files, I found that performance and access to attributes is much improved if the .osm file is converted to spatialite. Luckily, that’s easy using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES noirmoutier.db noirmoutier.osm

Finishing preprocessing in QGIS

In QGIS, we’ll want to load the points, lines, and multipolygons using Add SpatiaLite Layer:

Screenshot 2014-05-31 11.39.40

When we load the spatialite tables, there are a lot of features and some issues:

  • There is no land polygon. Instead, there are “coastline” line features.
  • Most river polygons are missing. Instead there are “riverbank” line features.

Screenshot 2014-05-31 11.59.58

Luckily, creating the missing river polygons is not a big deal:

  1. First, we need to select all the lines where waterway=riverbank.
    Screenshot 2014-05-31 13.14.00
  2. Then, we can use the Polygonize tool from the processing toolbox to automatically create polygons from the areas enclosed by the selected riverbank lines. (Note that Processing by default operates only on the selected features but this setting can be changed in the Processing settings.)
    Screenshot 2014-05-31 13.40.16

Creating the land polygon (or sea polygon if you prefer that for some reason) is a little more involved since most of the time the coastline will not be closed for the simple reason that we are often cutting a piece of land out of the main continent. Therefore, before we can use the Polygonize tools, we have to close the area. To do that, I suggest to first select the coastline using "other_tags" LIKE '%"natural"=>"coastline"%' and create a new layer from this selection (save selection as …) and edit it (don’t forget to enable snapping!) to add lines to close the area. Then polygonize.

Screenshot 2014-05-31 14.38.48

Styling the data

Now that all preprocessing is done, we can focus on the styling.

You can get the styles used in the map from my Github QGIS-resources repository:

  • osm_spatialite_googlemaps_multipolygon.qml … rule-based renderer incl. rules for: water, natural, residential areas and airports
  • osm_spatialite_googlemaps_lines.qml … rule-based renderer incl. rules for roads, rails, and rivers, as well as rules for labels
  • osm_spatialite_googlemaps_roadshields.qml … special label style for road shields
  • osm_spatialite_googlemaps_places.qml … label style for populated places such as cities and towns

qgis_osm_google_100k

Advertisement

Since I finally managed to download the elevation model of the city of Vienna, I thought I’d share some eye candy with you: The map uses layer blending to combine hillshade and elevation raster, and the elevation raster’s color ramp is a modified “garish14” from QGIS’ cpt-city color ramp collection.

wien_elevation by underdarkGIS
wien_elevation, a photo by underdarkGIS on Flickr.

Update

Here is how you get access to the “garish14” color ramp:

Start by selecting the "new color ramp" option in the raster's style window.

Start by selecting the “new color ramp” option in the raster’s style window.

Chose the "cpt-city" color ramp type.

Chose the “cpt-city” color ramp type.

In the "cpt-city color ramp" window, you will find lots of different premade color ramps. "garish14" is part of the "Topography" collection.

In the “cpt-city color ramp” window, you will find lots of different premade color ramps. “garish14” is part of the “Topography” collection.

Yesterday, I described my process to generate a basic population density map from the city of Vienna’s open government data. In the end of that post, I described some ideas for further improvement. Today, I want to follow-up on those ideas using what is known as dasymetric mapping. GIS Dictionary defines it well (much better than Wikipedia):

Dasymetric mapping is a technique in which attribute data that is organized by a large or arbitrary area unit is more accurately distributed within that unit by the overlay of geographic boundaries that exclude, restrict, or confine the attribute in question.
For example, a population attribute organized by census tract might be more accurately distributed by the overlay of water bodies, vacant land, and other land-use boundaries within which it is reasonable to infer that people do not live.

That’s exactly what I want to do: Based on subdistricts with population density values and auxiliary data – Corine Land Cover to be exact – I want to create an improved representation of population density within the city.

This is the population density map I start out with:

… and this is the Corine Land Cover dataset for the same area:

It shows built-up areas (red), parks and natural areas (green) as well as water-covered regions (blue). For further analysis, I follow the assumption that people only live in areas with Corine code 111 “Continuous urban fabric” and 112 “Discontinuous urban fabric”. Therefore, I use the Intersection tool to clip only these residential areas from the subdistrict polygons. The subdistrict population can now be distributed over these new, smaller areas (use Field Calculator) to create a more realistic visualization of population density:

For easier comparison, I put the original density and the dasymetric map into a looping animation. Some subdistricts change their population density values quite drastically, especially in regions where big parts covered by water or rail infrastructure were removed:

Corine Land Cover is not too detailed but I think it still usable on this scale. One thing to note is that I used data from 2006 with population data from 2012 so some areas in the outer districts will have been turned residential in the meantime. But I hope this doesn’t distort the overall picture too much.

The city of Vienna provides both subdistrict geometries and population statistics. Mapping the city’s population density should be straightforward, right? Let’s see …

We should be able to join on ZBEZ and SUB_DISTRICT_CODE, check! But what about the actual population counts? Unfortunately, there is no file which simply lists population per subdistrict. The file I found contains four lines for each subdistrict: females 2011, males 2011, females 2012 and males 2012. That’s not the perfect format for mapping general population density.

A quick way to prepare our input data is applying pivot tables, eg. in Open Office: The goal is to have one row per subdistrict and columns for population in 2011 and 2012:

Export as CSV, add CSVT and load into QGIS. Finally, we can join geometries and CSV table:

A quick look at the joined data confirms that each subdistrict now has a population value. But visualizing absolute values results in misleading maps. Big subdistricts with only average density will overrule smaller but much denser subdistricts:

That’s why we need to calculate population density. This is easy to do using Field Calculator. The subdistrict file already contains area values but even if they were missing, we could calculate it using the $area operator: "pop2012" / ($area / 10000). The resulting population density in population per ha finally shows which subdistricts are the most densely populated:

One could argue that this is still no accurate representation of population density: Big parts of some subdistricts are actually covered by water – especially along the Danube – and therefore uninhabited. There are also big parks which could be excluded from the subdistrict area. But that’s going to be the topic of another post.

If you want to use my results so far, you can download the GeoJSON file from Github.

Today, I’ve finished my submission for the Hubway Data Visualization Challenge. All parts of the resulting dataviz were created using open source tools. My toolbox for this work contains: QGIS, Spatialite, Inkscape, Gimp and Open Office Calc. To see the complete submission and read more about it, check the project page.

QGIS Cloud is a great cloud hosting service for QGIS Server by Sourcepole. After online registration and installation of an (experimental) plugin, QGIS projects can be uploaded to the cloud quite comfortably.

For a quick test, I tried to recreate the map from “Open Data for Physical Maps”. Right now, one of the limitations is that raster layers cannot be uploaded. Instead of the SRTM data, I therefore chose OCM landscape from OpenLayers plugin to provide some hillshade. The process of uploading data and publishing the project went smoothly and I didn’t encounter any problems.

You can explore the map online at qgiscloud.com/anitagraser/corine_austria.

Considering the complexity of the Corine dataset, rendering is quite fast – certainly much better than on my notebook.

Corine Land Cover is a European program to create a land cover inventory of Europe. The data is freely available and a valuable input for many analyses. In this post, we’ll be using it to create a physical map.

For the background, reused the hillshade presented in “Mapping Open Data With Open GIS”

Instead of the standard grayscale, I defined a sand-colored colormap that looks warmer and more natural:

On top of this hillshade, I put the Corine land cover layer. Instead of the official, rather bright colors I selected a more neutral color palette and varying transparency values: Water areas are drawn with no transparency while forests are set to 50 % and built-up areas to up to 80 % transparency. I also skipped classes such as “bare rock” by setting them to be totally transparent:

On top of the land cover, I added a river dataset and styled it with the same color used for water surfaces in the Corine layer. Obviously, this is an optional step. Big rivers are visible within the land cover data too.

After adding a mask and labels, the map is ready to add the finishing touches in Print Composer: Title, explanatory text and a scale bar. I decided against adding a legend to this particular map since I hope that the color choices are intuitive enough.

If you want to create your own physical map, you can use Corine Land Cover for European regions or the National Vegetation Classification in the U.S.

This post explores some cartographic features of QGIS while mapping the river network of Tirol, Austria. All data used is freely available.

For the background, I downloaded NASA SRTM elevation data from CGIAR-CSI and created a hillshade using Terrain Analysis tools in QGIS 1.8.

To emphasize both state borders and the fact that Tirol consists of two separate areas, I created a mask using the Difference tool and styled it a transparent white.

The river network is too dense to label all rivers on an A4 map. Expression-based labels make it possible to only label selected features. For this dataset, the expression limits labeling to features with certain values of GEW_WRRL attribute:

CASE WHEN (GEW_WRRL = '10.000 km2 Fluss' OR "GEW_WRRL"= '4.000 km2 Fluss' OR "GEW_WRRL"='1.000 km2 Fluss') AND length("GEW_NAME_A") < 10 THEN "GEW_NAME_A" END

Labels of neighboring areas together with map title, descriptions and scale bar were added in Print Composer.

Working with Print Composer, it is useful to know that you can use Copy&Paste to duplicate map components and right-click to lock them from being moved. Also, every new component by default comes with a black outline and white background which can (and should) be disabled/changed in “General options”.

This is the final QGIS Print Composer output – without any further post-processing in Inkscape or Gimp:

This is a follow up post on “Guide to Advanced Labeling for OSM Roads”. This post covers how to create a map that looks similar to the classic Google Maps map based on OSM data.

Styling OSM road data is a bit tricky due to the many different possible options in OSM “type” tag. It takes some fiddling around to get results similar to Google. The easiest way to create such a style is using rule-based renderer with “Symbol levels” enabled.

I’ve excluded a number of road types from being rendered and labeled. This is done by NOT specifying a rule that fits those classes. That’s the reason why I don’t have a “no filter” rule. Instead, I’ve specified

type NOT IN ('footway','footpath','steps','cycleway','pedestrian','track','bridleway')

to cover all the roads I don’t want specifically highlighted but displayed using default style. This way I can make sure that footways and similar are neither drawn nor labeled.

Google-style rules for OSM with multiple zoom levels

Now, we can have a look at how to create those road shields Google uses:

Sample from Google Maps

Currently, it’s not possible to solve this using the labeling engine. And even after labels with “shields” will be implemented, there is still the problem that we want both road names and road numbers labeled – preferably without having to duplicate the layer.

For now, one possible solution for this is to create the shields using the “Marker lines” feature of new symbology. As the name suggests, you can put markers on lines. In this case, the marker will be our road shield. It basically has two parts: the colored shield plus the text. The shield can be created by putting two squares besides each other. (Adjust the offsets to move them besides each other.) Then we can put the text on top, one letter at a time.

Creating the road shield marker

Place the marker on the line once and deactivate “Rotate marker”.

A "road with shield" style

These styles can be assigned to every road you want to decorate with a shield. If you save the style, you just need to change the text next time.

And this is how the result can look like.

Zoomed-out result

Maybe the shields turned out a little too big compared with the original.

Zoomed in, more labels will be displayed and an additional layer with metro stations becomes visible. The metro symbol was created using the same technique described for the road shields.

Zoomed-in result

One disadvantage of creating road shields with this technique – besides the fact that it’s rather tedious – is that there is no collision detection between labels and shields. Nonetheless, it’s a viable solution that allows you to create high-resolution maps that look very similar to Google Maps using free OSM data.

%d bloggers like this: