Geoinformatics

How to filter OpenStreetmap ways by tag from .osm.pbf file using Python & Osmium

In our previous post Minimal example how to read .osm.pbf file using Python & osmium we investigated how to read a .osm.pbf file and count all nodes, ways and relations.

Today, we’ll investigate how to filter for specific ways by tag using osmium working on .osm.pbf files. In this example, we’ll filter by highwaytrunk

#!/usr/bin/env python3
import osmium as osm
from collections import namedtuple
import itertools

Way = namedtuple("Way", ["id", "nodes", "tags"])

class FindHighwayTrunks(osm.SimpleHandler):
    """
    Find ways with "highway: trunk" tag
    """

    def __init__(self):
        osm.SimpleHandler.__init__(self)
        self.ways = []

    def way(self, way):
        # If this way is a highway: trunk, ...
        if way.tags.get("highway") == "trunk":
            # ... add it to self.ways
            # NOTE: We may not keep a reference to the way object,
            # so we have to create a new Way object
            nodes = [node.ref for node in way.nodes]
            self.ways.append(Way(way.id, nodes, dict(way.tags)))

# Find ways with the given tag
way_finder = FindHighwayTrunks()
way_finder.apply_file("kenya-latest.osm.pbf")

print(f"Found {len(way_finder.ways)} ways")

This example uses kenya-latest.osm.pbf downloaded from Geofabrik, but you can use any .osm.pbf file. Parsing said file takes about 10 seconds on my Desktop.

Note that we can’t just self.ways.append(way) because way is a temporary/internal Protobuf object and may not be used directly since. Hence, we create our own object that just wraps the extracted contents of way: Its id, its list of nodes (node IDs) and a dictionary containing the tags. Otherwise, you’ll see this error message:

Traceback (most recent call last):
  File "run.py", line 32, in <module>
    trunk_handler.apply_file("kenya-latest.osm.pbf")
RuntimeError: Way callback keeps reference to OSM object. This is not allowed.
Posted by Uli Köhler in Geography, Geoinformatics, OpenStreetMap, Python

How to install cartopy on Ubuntu

First, you need to install libproj (which is a dependency of cartopy) using

sudo apt -y install libproj-dev libgeos-dev

After that, install cartopy using

sudo pip3 install cartopy

 

Posted by Uli Köhler in Cartopy, Geoinformatics, Python

How fix Cartopy pip install error ‘Proj 4.9.0 must be installed.’

Problem:

When trying to install cartopy using e.g. pip3 install cartopy you see this error message:

Collecting cartopy
  Downloading Cartopy-0.19.0.post1.tar.gz (12.1 MB)
     |████████████████████████████████| 12.1 MB 13.2 MB/s 
  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  ERROR: Command errored out with exit status 1:
   command: /usr/bin/python3 /tmp/tmp2q7tvpo8 get_requires_for_build_wheel /tmp/tmpfwd5htse
       cwd: /tmp/pip-install-dg2i313t/cartopy
  Complete output (1 lines):
  Proj 4.9.0 must be installed.
  ----------------------------------------
ERROR: Command errored out with exit status 1: /usr/bin/python3 /tmp/tmp2q7tvpo8 get_requires_for_build_wheel /tmp/tmpfwd5htse Check the logs for full command output.

Solution:

Install proj using e.g.

sudo apt -y install libproj-dev

and then install cartopy again using e.g.

sudo pip3 install cartopy

 

Posted by Uli Köhler in Cartopy, Geoinformatics, Python

Computing bounding box for a list of coordinates in Python

Problem:

You have a list of X/Y coordinates, for example:

coords = [(6.74219, -53.57835),
          (6.74952, -53.57241),
          (6.75652, -53.56289),
          (6.74756, -53.56598),
          (6.73462, -53.57518)]

For these coordinates you want to compute the minimum bounding box.

Continue reading →

Posted by Uli Köhler in Geoinformatics, Python

Reading a Shapefile directly from a ZIP using pyshp

Problem

You have a ZIP file containing a single ESRI shapefile database (i.e. three files), for example the Natural Earth Countries database. Without unzipping the ZIP you want to use pyshp in order to read the data contain  in the shapefile.

Continue reading →

Posted by Uli Köhler in Geoinformatics