10 days ago Joachim Topf commited the new flex backend into osm2pgsql. This new backend gives you full control on what data goes where, and the opportunity to massage it a little bit before committing to the db. This is useful not only because we could move processing from the rendering time to the import time (which, let's be honest, we could already do some with the --tag-transform-script option), but also because we can move away from the 4 table setting of the original osm2pgsql mode (line, point, polygons, roads) and also ignore data we don't want in the db or even create new one.

This new backend works exactly like osm2pgsql. There are two stages; the first one goes through all the nodes, ways and relations (in that order), and the second one only through ways and relations. For each one of these types you can define a function, process_[type](). process_node() is only going to be called in stage 1, but process_way() and process_relation() are going to be called in both stages. The functions, of course, can figure out in which stage they are; my advise is to split the functions into process_node_stage1() and process_node_stage2() and call them from process_node() to make that more clear. An object is processed in stage 2 only if it was marked in stage 1. For more details, please read the docs, and definitely learn lua.

My first task for this is bus routes. Last time I put some of my few spare time on my fork I managed to display bus routes but it was ugly: overlapping dash based lines that were difficult to follow.

Compare that to a real map of a (quite small) bus network, and you get the idea that my take wasn't going to be useful.

Let's try to think how this can be achieved. What we need is a way to take each way and count how many bus lines are going through it. Bus lines are represented as relations, of which we can take the members. So the first step must be to find all the members of a relation. This can be done in process_relation_stage1(). We mark all the ways in a bus relation, and we also somehow associate the bus to the way.

In stage 2, we process all the marked ways (ways with buses), we count how many buses go through it, we create a new way for rendering each bus line, and we displace this new line based on which 'slot' it belongs to, so the lines are parallel.

I have a first version of such algorithm. It's based on one of the examples. Let's take it apart:

local tables = {}

tables.routes = osm2pgsql.define_way_table('routes', {
    { column = 'tags',   type = 'hstore' },
    { column = 'offset', type = 'float' },
    { column = 'ref',    type = 'text' },
    { column = 'colour', type = 'text' },
    { column = 'way',    type = 'linestring' },

We declare a table that will store ways, called routes. It will store the original tags of the way in hstore format, plus an offset column, that will tell the renderer how much to displace the line, plus a ref and a colour; and of course, the way itself.

pt_ways = {}

function osm2pgsql.process_relation(relation)
    -- Only interested in relations with type=route, route=bus and a ref
    -- TODO: other pt routes too
    if relation.tags.type == 'route' and relation.tags.route == 'bus' and relation.tags.ref then
        for _, member in ipairs(relation.members) do
            if member.type == 'w' then
                -- print(relation.id, member.ref, relation.tags.ref)
                if not pt_ways[member.ref] then
                    pt_ways[member.ref] = {}

                -- insert() is the new append()
                table.insert(pt_ways[member.ref], relation)

We process the relations. Like I said before, we take all the members that are ways, we store this bus line in an array indexed by the original way id, and we mark the way for processing in stage 2. Notice the the bus line's ref is in the tags.ref indexes (more on this later), while it's id is in the ref index. This last part was confusing to me.

function sort_by_ref(a, b)
    return a.tags.ref < b.tags.ref

function osm2pgsql.process_way(way)
    -- Do nothing for ways in stage1, it'll be in relations where the magic starts
    if osm2pgsql.stage == 1 then

    -- We are now in stage2
    local routes = pt_ways[way.id]
    table.sort(routes, sort_by_ref)

We do nothing for ways in stage 1, and in stage 2 we sort the routes to give them a consistent ordering when rendering consecutive ways.

    local line_width = 2.5
    local offset = 0
    local side = 1
    local base_offset
    local offset
    local index
    local ref
    local slot
    local shift

    if #routes % 2 == 0 then
        base_offset = line_width / 2
        shift = 1
        base_offset = 0
        shift = 0

    for index, route in ipairs(routes) do
        -- index is 1 based!
        slot = math.floor((index - shift) / 2)
        offset = (base_offset + slot * line_width) * side

        if side == 1 then
            side = -1
            side = 1

This is the part of the algo that calculates the offset. It was refined after a couple of iterations and it seems to work fine with odd and even amount of bus lines. line_width will be moved to the style later, so I can apply widths depending on ZL. In short, we're assigning slots from the center to the outside, alternating sides.

        -- TODO: group by colour
        -- TODO: generic line if no colour
        row = {
            tags = way.tags,
            ref = route.tags.ref,
            colour = route.tags.colour,
            offset = offset,
            geom = { create = 'line' }

        tables.routes.add_row(tables.routes, row)

And this is the end. We set the row values and add it to our table. It's that simple :) Now we run osm2pgsql with the flex backend:

$ osm2pgsql --cache 1024 --number-processes 4 --verbose --create --database mc --output=flex --style bus-routes.lua --slim --flat-nodes nodes.cache --hstore --multi-geometry --drop
osm2pgsql version 1.2.0 (1.2.0-248-gba17b0c) (64 bit id space)
Reading in file: monaco-latest.osm.pbf
Processing: Node(46k 0.4k/s) Way(3k 3.78k/s) Relation(10 10.00/s)  parse time: 125s
Node stats: total(46745), max(7199027992) in 124s
Way stats: total(3777), max(770935077) in 1s
Relation stats: total(215), max(10691624) in 0s
Entering stage 2...
Creating id indexes...
Creating id index on table 'routes'...
Creating id indexes took 0 seconds
Lua program uses 0 MBytes
Entering stage 2 processing of 458 ways...
Entering stage 2 processing of 0 relations...
Clustering table 'routes' by geometry...
Using native order for clustering
Creating geometry index on table 'routes'...
Analyzing table 'routes'...
All postprocessing on table 'routes' done in 0s.
Osm2pgsql took 129s overall

Now, for the render part, it's twofold; a layer definition...:

    COALESCE(colour, 'purple') as color,
    ref as route,
FROM routes

... and the style itself:

#routes {
  [zoom >= 14] {
    line-width: 1;
    line-color: @transportation-icon;
    line-join: round;
    line-offset: [offset];
    [zoom >= 17] {
      line-color: [color];
      line-width: 2;

Quite simple too! That's because all the data is already prepared for rendering. All the magic happens at import time. Here's the same region as before:

Now, two caveats: This last thing means that if you want to change the style, you most probably will need to reimport the data. It must have taken me some 20 iterations until I got the data in a way I could use for rendering, that's why I tested with an extract of Monaco :) I also used a separate db from the main render database, but maybe just another table would be enough.

Second, you have to specify all the data you want to save, there is no compatibility with the current rendering database, so you will also need to base your code on the compatible example. In my tests, I just imported the data the usual way:

$ osm2pgsql --cache 1024 --number-processes 4 --verbose --create --database mc --output=flex --style bus-routes.lua --slim --flat-nodes nodes.cache --hstore --multi-geometry --drop
osm2pgsql version 1.2.0 (1.2.0-248-gba17b0c) (64 bit id space)
Mid: loading persistent node cache from nodes.cache
Setting up table: planet_osm_nodes
Setting up table: planet_osm_ways
Setting up table: planet_osm_rels
Using lua based tag processing pipeline with script openstreetmap-carto.lua
Setting up table: planet_osm_point
Setting up table: planet_osm_line
Setting up table: planet_osm_polygon
Setting up table: planet_osm_roads
Reading in file: monaco-latest.osm.pbf
Processing: Node(46k 0.8k/s) Way(3k 3.78k/s) Relation(210 105.00/s)  parse time: 61s
Node stats: total(46745), max(7199027992) in 58s
Way stats: total(3777), max(770935077) in 1s
Relation stats: total(215), max(10691624) in 2s
Stopping table: planet_osm_nodes
Stopped table: planet_osm_nodes in 0s
Stopping table: planet_osm_ways
Stopped table: planet_osm_ways in 0s
Stopping table: planet_osm_rels
Stopped table: planet_osm_rels in 0s
Sorting data and creating indexes for planet_osm_point
Sorting data and creating indexes for planet_osm_roads
Sorting data and creating indexes for planet_osm_polygon
Sorting data and creating indexes for planet_osm_line
Using native order for clustering
Using native order for clustering
Using native order for clustering
Using native order for clustering
Copying planet_osm_point to cluster by geometry finished
Creating geometry index on planet_osm_point
Creating indexes on planet_osm_point finished
Copying planet_osm_roads to cluster by geometry finished
Creating geometry index on planet_osm_roads
Creating indexes on planet_osm_roads finished
All indexes on planet_osm_point created in 0s
Completed planet_osm_point
All indexes on planet_osm_roads created in 0s
Completed planet_osm_roads
Copying planet_osm_polygon to cluster by geometry finished
Creating geometry index on planet_osm_polygon
Creating indexes on planet_osm_polygon finished
All indexes on planet_osm_polygon created in 0s
Completed planet_osm_polygon
Copying planet_osm_line to cluster by geometry finished
Creating geometry index on planet_osm_line
Creating indexes on planet_osm_line finished
All indexes on planet_osm_line created in 0s
Completed planet_osm_line
Osm2pgsql took 64s overall

One thing to notice is that it took half of the time of the flex backend.

openstreetmap elevation

Posted sáb 15 feb 2020 10:59:36 CET Tags: openstreetmap

Today I finally sat down and spun off generate_tiles.py to its own repository, so people can follow its development without having to clone my own osm-carto fork. This happened just after I finished making the storage thread optional, because I usualy have as many rendering threads as cores I have, so that extra thread was not only compiting with them, but also it spent some time de/marshaling the metatiles between processes, so I'm not sure it's worth it. Maybe if I had more cores?

The new repo is here. There are some other tools there (hence the tools in the repo's name), but they're not so polished or documented. You're free to look and ask :)

What if I do a release, I hear (it must be the voice in my head)? Why not! I even decided to go bold and tagged it as v1.0.

openstreetmap pyhton

Posted vie 24 may 2019 12:46:20 CEST Tags: openstreetmap

Another aspect I've been looking into with respect to optimizing the rendering speed is data sources that are not in the target projection. Not being in the target projection forces Mapnik to reproject them on the fly, and this for each (meta)tile, whereas it would make more sense to have them already reprojected.

In my case I have 3 datasets, big ones, in EPSG:4258. The three datasets are elevation, hill and slope shade based on EEA's DEM files. The source tiles amount to almost 29GiB, and the elevation layer, being RGB, takes more than that. So I set off to try to reproject the things.

My first, more obvious approach was to reproject every 5x5°, 18000x18000px tile, then derive the data I need, but I started to get gaps between tiles. Notice that the original data is cleanly cut (5° x 3600"/° x 1px/" == 18000px), without any overlapping.

The next approach was to merge them all in a .vrt file, then reproject chunks of it with gdalwarp. What I wanted as output was the same 5x5° tiles, reprojected, but with an extra pixel, so they overlap. This last requirement was the problematic one. See, the final projection makes any square in the original projection a tall rectangle, stretching more and more towards the poles. The closest I could get was to use the -ts option, but that meant that I didn't get any control about how many extra pixels I got in the vertical/latitude direction. My OCD started thrashing :) In fact what happened was that I was not sure how GDAL would handle the possible partial pixel, whether rounding down (meaning excluding it), up (finishing it), or simply leaving the pixel with partial data and impacting the final rendering.

Even Rouault pointed to me that gdalwarp can do something fantastic: it can generate a .vrt file too with all the parameters needed for the reprojection, so reading from there was automatically reprojecting the original data. The resulting dataset is 288,000x325,220px (the original is 288,000x180,000px), so I'm definitely going to cut it down in small tiles. After consulting with a eight-ball, I decided to discard the idea of tiles with boundaries based on coordinates, which might not even make sense anymore, but settle for pixel based sizes, still with an extra pixel. The chosen size is 2**14+1 a.k.a. 16385. For this gdal_translate is perfect.

The final algorithm is like this:

gdalwarp -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 \
    +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over" \
    -r lanczos -tr 30.92208077590933 -30.92208077590933 \
    -of VRT EU-DEM.vrt EU-DEM-corrected.vrt

The values for the -tr option is the pixel size in meters, which is the unit declared in the SRS. Notice that as Mercator stretches towards the poles, this is the size at the origin of the projection; in this case, at 0° lat, 0° lon.

Then the reprojection (by reading from the reprojecting dataset) and cut, in a couple of loops:

tile_size=$((2**14)); \
for i in $(seq 0 17); do
    for j in $(seq 0 5); do
        for k in $(seq 0 3); do
            gdal_translate -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZMA \
                -co LZMA_PRESET=9 \
                -srcwin $(($tile_size*$i)) $(($tile_size*$l)) \
                    $(($tile_size+1)) $(($tile_size+1)) \
                -of GTiff EU-DEM-corrected.vrt \
                $(printf "%03dx%03d-corrected.tif" $i $l) &

There's an extra loop to be able to launch 4 workers at the same time, because I have 4 cores. This doesn't occupy the 4 cores a 100% of the time (cores that already finished stay idle until the other finished), but it was getting akward to express in a Makefile, and this is run only once.

Before deriving the rest of the data there's an extra step: removing those generated tiles that actually have no data. I do a similar thing with empty sea tiles in the rendering process. Notice also that the original data is not tile complete for the covered region (79 tiles instead of the 160 they should be).

openstreetmap gis gdal

Posted mar 23 ene 2018 14:08:42 CET Tags: openstreetmap

Since I started playing with rendering maps locally I've been modifying and extending the original generate_tiles.py script from mapnik-stilesheets. I added option parsing and lots of features; here's the current usage:

$ ./generate_tiles.py --help
usage: generate_tiles.py [-h] [-b BBOX] [-B BBOX_NAME] [-n MIN_ZOOM]
                        [-x MAX_ZOOM] [--tiles Z,X,Y [Z,X,Y ...]]
                        [-i MAPFILE] [-f FORMAT] [-o TILE_DIR]
                        [-m METATILE_SIZE] [-t THREADS]
                        [-p {threads,fork,single}] [-X] [-N DAYS]
                        [-E {skip,link,render}] [-d] [--dry-run]

optional arguments:
-h, --help            show this help message and exit
-b BBOX, --bbox BBOX
-B BBOX_NAME, --bbox-name BBOX_NAME
-n MIN_ZOOM, --min-zoom MIN_ZOOM
-x MAX_ZOOM, --max-zoom MAX_ZOOM
--tiles Z,X,Y [Z,X,Y ...]
-i MAPFILE, --input-file MAPFILE
-f FORMAT, --format FORMAT
-o TILE_DIR, --output-dir TILE_DIR
-t THREADS, --threads THREADS
-p {threads,fork,single}, --parallel-method {threads,fork,single}
-X, --skip-existing
-N DAYS, --skip-newer DAYS
-E {skip,link,render}, --empty {skip,link,render}
-d, --debug

BBoxes are stored in a file called bboxes.ini, so I can say -B Europe instead of remembering the coords. The idea of --format is that I should be supporting slippy maps .png file structure or mbtiles, but the latter support is a little lagging behind because I don't have a use for them yet. You can choose to whether use threads (broken because mapnik cannot handle the situation; I can't find a reference to the problem now), child processes (probably the only one working correctly) or a single main process (so no parallelism). It handles resuming a stopped render by not rendering if the tile exists or it's too new. It also can skip writing empty seas tiles.

I use it to rerender my style everytime I make a modification (or just update to the latest openstreetmap-carto, of which is a fork). I usually bulk render a great part of Europe up to ZL 11 or 14, and them some regions down to ZL 18 or 19 as needed for trips or other interests.

For Europe, it can take a long while, so I've been thinking on ways to optimize the rendering. Besides tuning the database, I first found that rendering big metatiles (8x8, for instance) gave a big boost in rendering time. The next idea is to reuse disk cache. When you render a (meta)tile in ZL N, the same data used for rendering it is going to be used for the 4 sub(meta)tiles of ZL N+1 (except when you remove features, which is rare but exists; city labels come to mind). I don't think something could be done at mapnik level, but one can think of the tiles as a tree: a node in ZL N has 4 subtiles in level N+1 and the leafs are the last ZL rendered. The original algorithm did a breath first traveling of this tree, but if you do a depth first algorithm, it could reuse the kernel's page/disk cache for the data collected by mapnik from the database or other files. Also, I can check whether the subtiles are render worthy: if they're only sea, I don't need to render it or its subtiles; I can cut down whole tile trees. The only point at which this could no longer be true is when we start rendering more things on sea, which currently ammounts to ferry routes at ZL 7.

I finished implementing all these ideas, but I don't have any numbers to prove it works. Definitely not rendering sea tiles should be a great improvement, but I don't really know whether the caching idea works. At least it was fun to implement.

So the rendering batch will be cut in 3: ZLs 0-6 in one run, then 7 and 8 with less threads (these ZLs of my style use so much memory the machine starts thrashing), then 9-18/19 with full threads.

elevation openstreetmap

Posted vie 10 mar 2017 20:44:23 CET Tags: openstreetmap

I just uploaded my first semi-automated change. This change was generated with my hack for generating centerlines for riverbank polygons. This time I expanded it to include a JOSM plugin which will take all the closed polygons from the selection and run the algorithm on them, creating new lines. It still needs some polishing, like making sure they're riverbanks and copying useful tags to the new line, and probably running a simplifying algo at some point. Also, even simple looking polygons might generate complex lines (in plural, and some of these lines could be spurious), so some manual editing might be required afterwards, specially connecting the new line to existing centerlines. Still, I think it's useful.

Like I mentioned last time, its setup is quite complex: The JOSM plugin calls a Python script that needs the Python module installed. That module, for lack of proper bindings for SFCGAL, depends on PostgreSQL+PostGIS (but we all have one, right? :-[ ), and connection strings are currently hardcoded. All that to say: it's quite hacky, not even alpha quality from the installation point of view.

Lastly, as imagico mentioned in the first post about this hack, the algorithms are not fast, and I already made my computer start thrashing the disk swapping like Hell because pg hit a big polygon and started using lots of RAM to calculate its centerline. At least this time I can see how complex the polygons are before handing them to the code. As an initial benchmark, the original data for that changeset (I later simplified it with JOSM's tool) took 0.063927s in pg+gis and 0.004737s in the Python code. More test will come later.

Okey, one last thing: Java is hard for a Pythonista. At some point it took me 2h40 to write 60 lines of code, ~2m40 per line!

openstreetmap gis python

Posted lun 29 ago 2016 20:06:39 CEST Tags: openstreetmap

My latest Europe import was quite eventful. First, I run out of space several times during the import itself, at indexing time. The good thing is that, if you manage to reclaim some space, and reading a little of source code[1], you can replay the missing queries by hand and stop cursing. To be fair, osm2pgsql currently uses a lot of space in slim+flat-nodes mode: three tables, planet_osm_node, planet_osm_way and planet_osm_relation; and one file, the flat nodes one. Those are not deleted until the whole process has finished, but they're actually not needed after the processing phase. I started working on fixing that.

But that was not the most difficult part. The most difficult part was that I forgot, somehow, to add a column to the import.style file. Elevation, my own style, renders different icons for different types of castles (and forts too), just like the Historic Place map of the Hiking and Bridle map[2]. So today I sat down and tried to figure out how to reparse the OSM extract I used for the import to add this info.

The first step is to add the column to the tables. But first, which tables should be impacted? Well, the line I should have added to the import style is this:

node,way   castle_type  text         polygon

That says that this applies to nodes and ways. If the element is a way, polygon will try to convert it to a polygon and put it in the planet_osm_polygon table; if it's a node, it ends in the planet_osm_point table. So we just add the column to those tables:

ALTER TABLE planet_osm_point   ADD COLUMN castle_type text;
ALTER TABLE planet_osm_polygon ADD COLUMN castle_type text;

Now how to process the extract? Enter pyosmium. It's a Python binding for the osmium library with a stream-like type of processing à la expat for processing XML. The interface is quite simple: one subclasses osmium.SimpleHandler, defines the element type handlers (node(), way() and/or relation()) and that's it! Here's the full code of the simple Python script I did:

#! /usr/bin/python3

import osmium
import psycopg2

conn= psycopg2.connect ('dbname=gis')
cur= conn.cursor ()

class CastleTypes (osmium.SimpleHandler):

    def process (self, thing, table):
        if 'castle_type' in thing.tags:
                name= thing.tags['name']
            # osmium/boost do not raise a KeyError here!
            # SystemError: <Boost.Python.function object at 0x1329cd0> returned a result with an error set
            except (KeyError, SystemError):
                name= ''
            print (table, thing.id, name)

            cur.execute ('''UPDATE '''+table+
                         ''' SET castle_type = %s
                            WHERE osm_id = %s''',
                         (thing.tags['castle_type'], thing.id))

    def node (self, n):
        self.process (n, 'planet_osm_point')

    def way (self, w):
        self.process (w, 'planet_osm_polygon')

    relation= way  # handle them the same way (*honk*)

ct= CastleTypes ()
ct.apply_file ('europe-latest.osm.pbf')

The only strange part of the API is that it doesn't seem to raise a KeyError when the tag does not exist, but a SystemError. I'll try to figure this out later. Also interesting is the big amount of unnamed elements with this tag that exist in the DB.

[1] I would love for GitHub to recognize something like https://github.com/openstreetmap/osm2pgsql/blob/master/table.cpp#table_t::stop and be directed to that method, because #Lxxx gets old pretty quick.

[2] I just noticed how much more complete those maps are. more ideas to use :)

openstreetmap gis python

Posted mar 09 ago 2016 17:45:56 CEST Tags: openstreetmap

In this last two days I've been expanding osm-centerlines. Now it not only supports ways more complex than a simple rectangle, but also ones that lead to 'branches' (unfortunately, most probably because the mapper either imported bad data or mapped it himself). Still, I tested it in very complex polygons and the result is not pretty. There is still lots of room for improvements.

Unluckily, it's not as stand alone as it could be. The problem is that, so far, the algos force you to provide now only the polygon you want to process, but also its skeleton and medial. The code extends the medial using info extracted from the skeleton in such a way that the resulting medial ends on a segment of the polygon, hopefully the one(s) that cross from one riverbank to another at down and upstream. Calculating the skeleton could be performed by CGAL, but the current Python binding doesn't include that function yet. As for the medial, SFCGAL (a C++ wrapper for CGAL) exports a function that calculates an approximative medial, but there seem to be no Python bindings for them yet.

So, a partial solution would be to use PostGIS-2.2's ST_StraightSkeleton() and ST_ApproximateMedialAxis(), so I added a function called skeleton_medial_from_postgis(). The parameters are a psycopg2 connection to a PostgreSQL+PostGIS database and the way you want to calculate, as a shapely.geometry, and it returns the skeleton and the medial ready to be fed into extend_medials(). The result of that should be ready for mapping.

So there's that. I'll be trying to improve it in the next days, and start looking into converting it into a JOSM plugin.

openstreetmap gis python

Posted vie 29 jul 2016 00:39:57 CEST Tags: openstreetmap

For a long time now I've been thinking on a problem: OSM data sometimes contains riverbanks that have no centerline. This means that someone mapped (part of) the coasts of a river (or stream!), but didn't care about adding a line that would mark its centerline.

But this should be computationally solvable, right? Well, it's not that easy. See, for given any riverbank polygon in OSM's database, you have 4 types of segments: those representing the right and left riverbanks (two types) and the flow-in and flow-out segments, which link the banks upstream and downstream. With a little bit of luck there will be only one flow-in and one flow-out segment, but there are no guarantees here.

One method could try and identify these segments, then draw a line starting in the middle of the flow-in segment, calculating the middle by traversing both banks at the same time, and finally connect to the middle for the flow-out segment. Identifying the segments by itself is hard, but it is also possible that the result is not optimal, leading to a jagged line. I didn't try anything on those lines, but I could try some examples by hand...

Enter topology, the section of maths that deals with this kind of problems. The skeleton of a polygon is a group of lines that are equidistant to the borders of the polygon. One of the properties this set of lines provides is direction, which can be exploited to find the banks and try to apply the previous algorithm. But a skeleton has a lot of 'branches' that might confuse the algo. Going a little further, there's the medial axis, which in most cases can be considered a simplified skeleton, without most of the skeleton branches.

Enter free software :) CGAL is a library that can compute a lot of topological properties. PostGIS is clever enough to leverage those algorithms and present, among others, the functions ST_StraightSkeleton() and ST_ApproximateMedialAxis(). With these two and the original polygon I plan to derive the centerline. But first an image that will help explaining it:

The green 'rectangle' is the original riverbank polygon. The thin black line is the skeleton for it; the medium red line is the medial. Notice how the medial and the center of the skeleton coincide. Then we have the 4 branches forming a V shape with its vertex at each end of the medial and its other two ends coincide with the ends of the flow in and flow out segments!

So the algorithm is simple: start with the medial; from its ends, find the branches in the skeleton that form that V; using the other two ends of those Vs, calculate the point right between them, and extend the medial to those points. This only calculates a centerline. The next step would be to give it a direction. For that I will need to see if there are any nearby lines that could be part of the river (that's what the centerline is for, to possibly extend existing rivers/centerlines), and use its direction to give it to the new centerline.

For the moment the algorithm only solves this simple case. A slightly more complex case is not that trivial, as skeletons and medials are returned as a MultiLineString with a line for each segment, so I will have to rebuild them into LineStrings before processing.

I put all the code online, of course :) Besides a preloaded PostgreSQL+PostGIS database with OSM data, you'll need python3-sqlalchemy, geoalchemy, python3-fiona and python3-shapely. The first two allows me to fetch the data from the db. Ah! by the way, you will need a couple of views:

CREATE VIEW planet_osm_riverbank_skel   AS SELECT osm_id, way, ST_StraightSkeleton (way)      AS skel   FROM planet_osm_polygon WHERE waterway = 'riverbank';
CREATE VIEW planet_osm_riverbank_medial AS SELECT osm_id, way, ST_ApproximateMedialAxis (way) AS medial FROM planet_osm_polygon WHERE waterway = 'riverbank';

Shapely allows me to manipulate the polygonal data, and fiona is used to save the results to a shapefile. This is the first time I ever use all of them (except SQLAlchemy), and it's nice that it's so easy to do all this in Python.

openstreetmap gis python

Posted mar 26 jul 2016 18:55:14 CEST Tags: openstreetmap

A few weeks ago an interesting PR for osm-carto landed in the project's GitHub page. It adds rendering for several natural relief features, adding ridges, valleys, aretes, dales, coulouirs and others to cliffs, peaks and mountain passes, which were already being rendered. I decided to try it in Elevation (offline for the moment).

I sync'ed the style first with the latest release, applied the patch and... not much. My current database is quite old (re-importing takes ages and I don't have space for updates), so I don't have much features like that in the region I'm interested in. In fact, I went checking and the closest mountain range around here was not in the database, so I added it.

By the way, the range is mostly concurrent with a part of an administrative boundary, but SomeoneElse and SK53 suggested to make a new line. Even when other features are nearby (there's a path close to the crest and it's also more or less the limit between a forest and a bare rock section), which already makes the region a little bit crowded with lines, it makes sense: boundaries, paths, forest borders and ridges change at different time scales, so having them as separate lines makes an update to any of those independent of the rest.

Now I wanted to export this feature and import it in my rendering database, so I can actually see the new part of the style. This is not an straightforward process, only because when I imported my data I used osm2pgsql --drop, which removes the much needed intermediate tables for when one wants to update with osm2pgsql --append. Here's a roundabout way to go.

First you download the full feature (thanks RichardF!). In this case:


This not only exports the line (which is a sequence of references to nodes) with its tags, but the nodes too (which are the ones storing the coords). The next step is to convert it to something more malleable, for instance, GeoJSON. For that I used ogr2ogr like this:

ogr2ogr -f GeoJSON 430573542.GeoJSON 430573542.xml lines

The last parameter is needed because, quoting Even Rouault (a.k.a. José GDAL): «you will always get "points", "lines", "multilinestrings", "multipolygons" and "other_relations" layers when reading a osm file, even if some are empty», and the GeoJSON driver refuses to create layers for you:

ERROR 1: Layer lines not found, and <span class="createlink">CreateLayer</span> not supported by driver.

But guess what, that not the easiest way :) At least we learned something. In fact postgis already has a tool called shp2pgsql that imports ESRIShapeFiles, and ogr2ogr produces by default this kind of file. It creates a .shp file for each layer as discussed before, but again, we're only interested in the line one. So:

ogr2ogr 430573542 430573542.xml lines
shp2pgsql -a -s 900913 -S 430573542/lines.shp > 430573542.sql

We can't use this SQL file directly, as it has a couple of problems. First, you can't tell shp2pgsql the names of the table where you want to insert the data or the geometry column. Second, it only recognizes some attributes (see below), and the rest it tries to add them as hstore tags. So we have to manually edit the file to go from:

INSERT INTO "lines" ("osm_id","name","highway","waterway","aerialway","barrier","man_made","z_order","other_tags",geom)
    VALUES ('430573542','Montagne Sainte-Victoire',NULL,NULL,NULL,NULL,NULL,'0','"natural"=>"ridge"','010500002031BF0D[...]');


INSERT INTO "planet_osm_line" ("osm_id","name","z_order","natural",way)
    VALUES ('430573542','Montagne Sainte-Victoire','0','ridge','010500002031BF0D[...]');

See? s/lines/planet_osm_line/, s/other_tags/"natural"/ (with double quotes, because natural is a keyword in SQL, as in natural join), s/geom/way/ and s/'"natural"=>"ridge"'/'ridge'/ (in single quotes, so it's a string; double quotes are for columns). And I also removed the superfluous values and the ANALIZE line, as I don't care that much. Easy peasy.

A comment on the options for shp2pgsql. -s 900913 declares the SRID of the database. I got that when I tried without and:

ERROR:  Geometry SRID (0) does not match column SRID (900913)

-S is needed because shp2pgsql by default generated MultiLineStrings, but that table in particular has a LineString way column. This is how I figure it out:

ERROR:  Geometry type (MultiLineString) does not match column type (LineString)

Incredibly, after this data massacre, it loads in the db:

$ psql gis < 430573542.sql


openstreetmap gdal postgis

Posted mar 12 jul 2016 17:37:22 CEST Tags: openstreetmap

For a long time I've been searching for a program that would allow me to plan (car) trips with my friends. Yes, I know of the existence of Google Maps, but the service has several characteristics that doesn't make it appealing to me, and lacks a couple of features I expect. This is more or less the list of things I want:

  1. Define the list of points I want to go to. No-brainer.
  2. Define the specific route I want to take. This is normally implemented by adding more control points, but normally they're of the same category as the waypoins of the places you want to visit. I think they shouldn't.
  3. Define stages; for instance, one stage per day.
  4. Get the distance and time of each stage; this is important when visiting several cities, for having an idea of how much time during the day you'll spend going to the next one.
  5. Define alternative routes, just in case you don't really have/make the time to visit some points.
  6. Store the trips in cookies, share them via a URL or central site, but that anybody can easily install in their own server.
  7. Manage several trips at the same time.

So I sat down to try and create such a thing. Currently is just a mashup of several things GIS: my own OSM data rendering, my own waypoints-in-cookies idea (in fact, this is the expansion of what fired that post) and OSRM for the routing. As for the backend, I decided to try flask and flask-restful for creating a small REST API for storing all this. So far some basics work (points #1 and #6, partially), and I had some fun during the last week learning RESTful, some more Javascript (including LeafLet and some jQuery) and putting all this together. Here are some interesting things I found out:

  • RESTful is properly defined, but not for all URL/method pairs. In particular, given that I decide that trip ids are their name, I defined a POST to trips/ as the UPSERT for that name. I hope SQLAlchemy implements it soon.
  • Most of the magic of RESTful APIs happen in the model of your service.
  • Creating APIs with flask-restful could not be more obvious.
  • I still have to get my head around Javascript's prototypes.
  • Mouse/finger events are a nightmare in browsers. In particular, with current leafLet, you get clicked events on double clicks, unless you use the appropriate singleclick plugin from here.
  • Given XSS attacks, same-origin policy is enforced for AJAX requests. If you control the web service, the easiest way to go around it is CORS.
  • The only way to do such calls with jQuery is using the low level function $.ajax().
  • jQuery provides a function to parse JSON but not to serialize to it; use window.JSON.stringify().
  • Javascript's default parameters were not recognized by my browser :(.
  • OSRM's viaroute returns the coordinates multiplied by 10 for precision reasons, so you have to scale it down.
  • Nominatim and OSRM rock!

I still have lots of things to learn and finish, so stay tunned for updates. Currently the code resides in Elevation's code, but I'll split it in the future.


I have it running here. You can add waypoints by clicking in the map, delete them by doublecliking them, save to cookies or the server (for the moment it overwrites what's there, as you can't name the trips or manage several yet) and ask for the routing.

trip-planner elevation openstreetmap osrm python flask leaflet javascript jquery

Posted lun 25 ene 2016 15:52:06 CET Tags: openstreetmap