Category: Tutorials

Creating OpenStreetMap History Visualisations

I created a couple of OSM visualisations for my talk at the OSGeo Ireland conference.

See: History of OpenStreetMap in Ireland

These are pretty easy to make, but take a fair bit of time. I did mine for Ireland, but should work with any part of the world.

Required software:

  • PostgreSQL with PostGIS
  • Python
  • QGIS
  • osmium-tools

This is the trickiest part, installing osmium-tools: here.

Data:

An OSM full history export. The best source for these is GEOFABRIK.

For Ireland:

http://download.geofabrik.de/europe/ireland-and-northern-ireland.html

Due to GDPR, you will have to log in with an OSM id to download the full history extracts. User ID’s are personal data.

Process:

The workflow is pretty simple. Osmium-tools provides pretty easy API access to the history files, where you can provide a data, and it will extract what OSM was like at that date. We simply need to loop through the desired dates we want to extract, and pipe the results into a workflow that loads the data into PostgreSQL. The final step is simply rendering in QGIS using the time manager plugin.

Python Script:

Github GIST:

https://gist.github.com/HeikkiVesanto/f01ea54cca499a6a144d18cf8909c940

The tables in the database will be:

  • lines
  • multilinestrings
  • multipolygons
  • other_relations
  • points

Each feature will be tagged with the date it is associated with.

Visualisation:

To visualise the data in QGIS we use simply use the excellent time manager plugin, filtering on the load_date field and with a monthly interval.

Result:

Storm Harvey QGIS Geometry Generator

Storm Harvey produced some extremely high levels of rainfall. Some areas of Texas received over 50 inches of rain over 9 days.

The National Oceanic and Atmospheric Administration (NOAA) provided some really great real time datasets to map the progress of the storm.

Among these were:
Hourly Precipitation
and
Hurricane Path

From these we can produce a GIF of hourly precipitation:

Hourly precipitation.

And total precipitation:
Hurricane Harvey Total Precipitation

Particularly the hurricane path was possible to create in QGIS using the Atlas Generator, and the excellent new:ish geometry generator. This can be found as an option for any layers symbology, as one of the renderers.

For my map I had a non spatial table that drove my atlas. This was a log table of all of the hours of precipitation I had loaded into my database. So I looped through each entry and showed the corresponding points of hourly precipitation for the corresponding hour. I also had hurricane path data as points for every 6 hours. So I could use the geometry generator to interpolate points in between known points.

While the query ended up being pretty long it is pretty straightforward.

It only needs to be run when the hour being generated does not end with a 00, 06, 12, or 18, because those are the positions I already know.

For the rest I need to generate two points. One for the previously known point, and one for the next known point.

Then I would create a line between those two points, measure the line, and place a point on the line x times one sixth of the way for the start of the line depending on the hour from the last known point.

Overall I am very impressed and happy with the result. With a bit of data defined rotation the storm progress looks great.

line_interpolate_point(
Make_line(
geometry(
case when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12', '18') then
    get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') )
else
    get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - (attribute(  @atlas_feature , 'id') % 100 % 6  ))
end)
,
geometry(
case
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 6 )
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('18') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 100 - 18 )
    when to_int(right(to_string(attribute(@atlas_feature , 'id')),2)) > 18 then
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 100 - 18)
    else
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 6)
end)
),
length(Make_line(geometry(
case when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12', '18') then get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') )
else get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - (attribute(  @atlas_feature , 'id') % 100 % 6  ))
end
)
,
geometry(
case
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 6 )
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('18') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 100 - 18 )
    when to_int(right(to_string(attribute(@atlas_feature , 'id')),2)) > 18 then
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 100 - 18)
    else
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 6)
end)))
*
((attribute(  @atlas_feature , 'id') % 100 % 6) * 0.16666666666666666))