samedi 15 novembre 2014

Streets of Paris Colored by Orientation

Recently, I read an article by datapointed which presented maps of streets of different cities colored by orientation.
The author gave some details about the method, which I tried to reproduce. In this post, I present the different steps from the calculation in my favorite spatial R ToolBox to the rendering in QGIS using a specific blending mode.

The code is given at the end of the article.
Data used in this post is from OpenStreetMap contributors :
Get it on OSM Metro Extract or GeoFabrik

First, polylines are split into lines.





In what we want, the orientation is not directed and orthogonal lines are rendered the same way.


This way, all our segments cover just a quarter of a circle and their sine values range from 0 to 1.

We can affect each segment an angle. I classified my angle values in two classes, with punchy base colors.

Actually, values near 0 and 1 correspond to almost equivalent orientations. They correspond to angles near 0 / 90°. The difference of colors is maximal when two streets have a 45° difference of angles. I chose a kind of diverging palette with dark saturated colors near 45° and bright desaturated ones near 0 / 90°. This way, similar angles are given similar brightness values.

Here is the result for our roads :

Rendering with the inner orientation affected to the segments

The second part is mostly cosmetic. We'll use interpolation and smoothing to get a nice halo effect.

First, we generate the mid point of each line and affect the orientation to it.

Next, we process the interpolation by inverse distance weight based on these points and their orientation values.

interpolation

Finally, we smooth the continuous surface to get something more pleasant.

Interpolated raster & the same, smoothed


We export the smoothed raster for a rendering within QGIS.

Here is how I superimposed my layers :
  1. The raster is put on top with the palette described above and a specific "Overlay" blending mode.
  2. The roads are put below with a grey value. Note that the grey value affects the final color because of the blending mode.
  3. Finally, a near black background lets the halo appear slightly.



How the underlying color affects the final rendering in an overlay blending mode

Here is the result for Paris :

Smoothed result with a kind of halo  (click to see in large scale)

I decided to apply the process to a sample of NYC :

NYC Streets oriented (click to see in large scale)

Technical details and code

The two key functions are :
spatstat::angles.psp which compute the angles of segments
spatstat::idw which interpolates the data.

7 commentaires:

  1. Hi and thanks for this post. I'm trying to reproduce it with this data :

    https://geo.pays-de-brest.fr/donnees/Documents/Public/BMO_Regles_Circulation.zip

    For now, my map looks like this :)

    http://tiles.kupaia.fr/IMG/jpg/oriented_colors_v1.jpg

    The R script is running well for the shapefile, but i get an error on raster generation, does it sound familiar to you ?

    Message d'avis :
    In readOGR("BMO_Regles_Circulation/arr_circ_stat_l.shp", "arr_circ_stat_l") :
    Dropping null geometries: 1895, 2090, 4894
    Erreur : is.matrix(w) n'est pas TRUE
    Exécution arrêtée

    I'm not sure about the palette i'm using in QGIS, here are my settings :

    http://min.us/i/bfSIQmtRreUxc

    Would you share your QGIS project file to see the settings used for styling the vector data ?

    Thanks again.

    RépondreSupprimer
  2. mn <- focal(interp, w=3, fun=mean) # we apply a focal window to the data that will smooth the result.
    Error: is.matrix(w) is not TRUE

    RépondreSupprimer
    Réponses
    1. You're right, the code is
      mn <- focal(interp, w=matrix(1,3,3), fun=mean)

      If you want a bigger matrix, then try mn <- focal(interp, w=matrix(1,5,5), fun=mean)

      The number of pixels for the focal window must be uneven. w=matrix(1,4,4) won't work.

      Supprimer
  3. Wow, super beau travail !! ;-)

    Philippe Garvie

    RépondreSupprimer
  4. I'm having trouble running the code because apparently the "rdal" library is not available for R version 3.1.2. Any ideas on what to do? Thanks!

    RépondreSupprimer
  5. Always nice to bump into a nice code in R :) Thank you!

    RépondreSupprimer