Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
layout title
examples
Geospatial Analysis

An osm2pgsql database and PostGIS are well-suited for geospatial analysis using OpenStreetMap data. PostGIS provides a large number of geometry functions{:.extlink} and a full description of how to perform analysis with them is beyond the scope of this document, but a simple example of finding the total road lengths by classification for a municipality shows the concepts.

To start with, we'll download the data for the region as an extract from Geofabrik{:.extlink}. We could import all of that data, but that's going to take a while. It is better to filter out only what we need using Osmium{:.extlink}. In this case these are the highway geometries and the administrative boundaries.

osmium tags-filter -v -o data.osm.pbf british-columbia-latest.osm.pbf w/highway r/boundary=administrative

We can now import the result with osm2pgsql into a database called gis that we created beforehand:

osm2pgsql -d gis -O flex -S highways.lua data.osm.pbf

Here is the highways.lua config file we are using:

{% include download.html file="highways.lua" %}

{%- include_relative highways.lua -%}

Loading should take only a few seconds. Once this is done we'll open a PostgreSQL terminal with psql -d gis, although a GUI like pgadmin or any standard tool could be used instead.

We'll first find the ID of the polygon we want

SELECT area_id FROM boundaries
    WHERE tags->'boundary' = 'administrative'
      AND tags->'admin_level' = '8'
      AND tags->'name' = 'New Westminster';

The result is this:

  osm_id
----------
 -1377803

The negative sign tells us that the geometry is from a relation, and checking on the OpenStreetMap site{:.extlink} confirms that it is.

We want to find all the roads in the city and get the length of the portion in the city, sorted by road classification. Roads are in the highways table, the administrative areas in the boundaries table:

SELECT
    round(SUM(
      ST_Length(
        ST_Intersection(geom, (SELECT geom FROM boundaries WHERE area_id=-1377803))::geography
    ))) AS "length (meters)", type
  FROM highways
    WHERE ST_Intersects(geom, (SELECT geom FROM boundaries WHERE area_id=-1377803))
    GROUP BY type
    ORDER BY "length (meters)" DESC
    LIMIT 10;

The result is this table:

 length (meters) | type
-----------------+--------------
          118166 | residential
          106993 | service
           69912 | footway
           25540 | tertiary
           23180 | unclassified
           17407 | primary
           17075 | cycleway
           10502 | secondary
            7119 | path
            5306 | motorway

The cast ...::geography is the easiest way to get the result in meters.

More complicated analysises can be done, but this simple example shows how to use the tables and put conditions on the columns.