Archives de catégorie : R

Le bâti de Paris façon réseau ou comment faire des gribouillis avec R comme avec Processing

Paris sous forme de réseau. Non, on ne parle pas de rues, mais de bâti. Dans cet article, nous proposerons le passage d’une représentation à plat des bâtiments en une figuration sous forme de lignes, aux multiples évocations. On verra qu’avec R, on peut faire certaines choses comme Processing, l’outil de référence en design génératif.

La carte


Cliquer sur l’image pour la voir en grand sur flickr
Sur cette carte est représenté le bâti de Paris sous la forme d’une toile. Chaque point de couleur représente un bâtiment et la couleur de chacun indique la période de construction du bâtiment. 



La rampe de couleurs est celle utilisée par Etienne Côme sur sa carte interactive BATIPARIS.

Sur ma carte, un arc relie deux bâtiments et il est d’autant plus clair qu’il relie deux bâtiments de périodes de construction différentes.

Ce mode de visualisation peut être intéressant pour explorer les différences de valeur dans une variable dans un voisinage, comme dans le cas de cette carte de la densité de population en France :

Différence de densité de population au niveau départemental (RGP INSEE 2013)

Sur ma carte, en suivant le réseau, on trouve toujours un chemin pour aller d’un bâtiment à un autre, même si l’on part de la place de l’Étoile pour se diriger vers Bercy; ceci a été obtenu grâce à une technique appelée arbre couvrant de poids minimal. Cette représentation a pour objet d’illustrer la densité d’une ville où tout serait interconnecté.

Paris dresse sa toile


La représentation du bâti est le plus souvent faite avec l’emprise des bâtiments, ce qui peut donner de beaux résultats comme celui de BATIPARIS.
Afin d’illustrer l’évolution champignonnesque d’une ville et se départir de cette représentation classique des bâtiments, on pourrait être tentés de représenter cela sous forme de réseau.

toile d’araignée, mycélium et réseau neuronal

Sur ma carte, les ramifications qui relient chaque bâtiment peuvent évoquer de multiples choses : les synapses de neurones rappelant les smart cities, les filaments d’un mycélium, voire même, si on regarde tout ça en noir et blanc, une sorte d’écriture qui évoquerait le tissage, au fil du temps, d’un tissu urbain.

Les données


Les données utilisées sont celles de l’Atelier Parisien d’URbanisme (APUR). Les données de bâti de l’APUR sont très riches. Outre la période de construction, on y trouve l’année de réhabilitation, les caractéristiques de la toiture, la surface d’ensoleillement,…
Une variable est particulièrement intéressante, que nous exploiterons peut-être dans un nouveau post : la hauteur des bâtiments.

La colonne que nous avons utilisée est la même que celle utilisée par Etienne Côme dans son application : la colonne C_PERCONST.

La transformation


L’illustration ci-dessus, créée par Benedikt Gross, montre le processus en cours lors d’un design de données. Je lui ai juste ajoutée les outils utilisés.
Pour obtenir notre carte, nous avons créé un algorithme exploitant la théorie des graphes ainsi que des techniques de design génératif. Pour ce dernier domaine, Processing constitue l’outil de référence. Il est adopté par nombre de data artists, d’artistes multimédia. Nous verrons grâce à cet article que R permet aussi de faire du design génératif. Néanmoins, pour peindre nos données, nous nous sommes tournés vers un second outil : QGIS, plutôt que R, en raison de capacités graphiques assez particulières.

Dans un premier temps, les centroïdes des bâtiments ont été générés

L’île de la Cité

Ensuite, nous avons relié tous ces points grâce à la fonction de Delaunay

De l’eau naît

J’ai trouvé que le nombre de segments était bien trop grand pour le rendu final. Afin de l’épurer, j’ai utilisé une technique de la théorie des graphes appelée Arbre couvrant de poids minimal, ou minimum spanning tree en anglais. Cet arbre est tel qu’on peut trouver un chemin entre chaque sommet avec un nombre de branches minimal

MST (Minimum Spanning Tree)

Les liaisons entre bâtiments étaient selon moi trop rectilignes. Afin d’obtenir quelque chose de plus “smooth”, j’ai utilisé les courbes de Bézier. C’est une forme visuelle très usitée dans les infographies de réseaux.

Avec Bézier, tout devient plus doux

Le style


  • J’ai adopté un fonds noir pour mettre en contraste les données. 
  • J’ai posé le réseau de segments avec un petit effet flou, avec le mode de fusion “addition” et une valeur de luminosité différente selon la différence de période reliant deux bâtiments. Grâce au mode addition ajouté entre les objets, les endroits où les segments se superposent semblent s’éclairer.
  • Par-dessus, j’ai posé les centroïdes de bâtiments coloriés selon la période de construction, avec un mode de fusion “revêtement”. Pour rappel, la rampe de couleurs est celle utilisée dans BATIPARIS. J’avais utilisé le mode de fusion “revêtement” pour la carte de l’orientation des rues de Paris. Avec ce mode, c’est comme si on teintait une couche avec les couleurs d’une autre.
  • Enfin, en dessous, j’ai disposé la Seine afin de donner un contexte. Celle-ci se trouve dans une couleur hors du spectre utilisé pour le bâti, et avec un peu de flou pour donner un côté organique et en atténuer la présence. 

J’aurais pu rajouter des données contextuelles supplémentaires mais j’ai préféré garder un côté minimaliste pour cette carte.

Voici un tableau récapitulant l’agencement des différentes couches dans QGIS avec leur style et leur ordre de superposition :

Variations


J’ai réalisé des petites variations autour de cette carte.

Tout d’abord, en générant des lots de données par période, on peut voir dans quelles zones Paris a dressé sa Toile au cours du temps. C’est comme si on voyait la signature de chaque époque.

Jusqu’en 1800, on construit surtout au centre
Georges Eugène Hausmann est dans la place. Ca construit à tout va, y compris des places (comme celle de l’Étoile)
Le rythme de construction semble plus modéré dans les 80’s, à part peut-être dans le NE de la capitale


Le code


Décortiquage

La librairie spdep a permis de réaliser la triangulation de Delaunay :

coords = coordinates(ile_cite)
nb = tri2nb(coords)

ainsi que de calculer l’arbre couvrant de poids minimal :

nb = tri2nb(coords)
nbw = nb2listw(nb, zero.policy = FALSE)
mst.bh <- mstree(nbw,5)

J’ai créé une fonction appelée MakeLineFromCoords qui permet de tracer une ligne à partir d’une série de coordonnées. D’ailleurs, je me suis dit qu’elle pouvait être utile pour superposer tout type de graphique en ligne sur des cartes : par exemple, évolution du chômage par département au cours du temps :

makeLineFromCoords <- function(coords, i) {

  Sl1 = Line(coords)
S1 = Lines(list(Sl1), ID=as.character(i))
Sl = SpatialLines(list(S1))
return(Sl)
}

R inclut une librairie appelée bezier qui, comme son nom l’indique, permet de dessiner…des courbes de Bézier.

b = bezier(t=t <- seq(0, 1, length=100), p=rbind(coords[from, ], midpoint, coords[to, ]))

midpoint désigne le point d’ancrage situé entre les deux points à connecter. La position de ce dernier détermine la courbure finale de l’arc. J’ai injecté un peu de hasard dans la position de ce dernier.
Dans un premier temps, j’ai calculé la distance séparant les deux points :

d = sqrt((coords[from, 1] – coords[to, 1])^2 + (coords[from, 2] – coords[to, 2])^2)

Ensuite, j’ai bougé sa position selon certains seuils. Tout d’abord, j’ai calculé un facteur de transformation qui déterminera dans quelle direction se fait le déplacement du point :

ifelse(runif(1,-1,1) > 0, 1, -1)


Ensuite, j’ai appliqué un ratio tel que le point bougera dans un rayon allant du milieu situé entre les deux points à une distance égale au tiers de la distance séparant les deux points.

Voici le code correspondant au déplacement en X :

(coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) > 0, 1, -1

Et celui en Y est similaire :

(coords[from, 2] + coords[to, 2])/2 + ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3)

Au final, le point intermédiaire d’ancrage de la courbe de Bézier se situera ici :

midpoint = c((coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3),(coords[from, 2] + coords[to, 2])/2 + ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3))

Le code dans son intégralité



Corine Disco Cover

This post is about the making of an experimental video that dives into the land cover of Paris & its surroundings with a hint of RGB magic. 

You’ll see the purpose of gridding, how to cross a vector and a raster layer in R. Finally we’ll cover one of the latest functionnality of QGIS 2.12 : data-defined layer combinations in atlases.


If the world were square


do you prefer ice creams or sandwiches ?
Matrix is the first movie approaching grids

Gridding is a way to simplify the representation of data on a large scale. It can be a useful way to “resample” the vision of a territory if you work on bigger units. These units, usually squares, are comparable, on the contrary to administrative units. If you can omit the administrative units, gridding is a way to get an homogeneous view of your territory. The larger your area of interest will be, the bigger the squares from your grid you might configure.

hexagons look ecologic, triangles ecclesiastic

Also, you could prefer hexagonal units to square units for more harmony. Some prefer triangles. Anyway, the number of possible shapes are limited.

Land Cover of Paris & surroundings

I wanted to visualize the land cover of Ile-De-France width a grid. Ile-de-France is a region that includes Paris at its center. Precisely, I decided to study the repartition of urban , agricultural and natural (forests, gardens, etc…) areas, that is 3 variables over the whole area.

the best diagram I’ve ever made, even if pie charts are the worst

The question was : how to represent 3 variables at a time on each square unit ? One answer could have been “diagrams”, that is pie charts, histograms on the map
But I found it cooler to use an RGB representation.

A Freemason way of seeing data

Let’s suppose R means urban, G Natural and B Agricultural. If the color of a square approaches red, then it means urban; if it is purple, it means urban & agricultural areas dominate but not natural; if it is desaturated, it means the 3 types of cover are equally mixed.

The recipe

I remember a mushroom im my parents’garden had the same shape

For the land cover of Ile-de-France, we have a nice 2012 Corine Land Cover GeoTiff.

Laurent Manadou exploring data

I wondered if one could adopt a dynamic view that goes from the global view to the zonal one in a sort of “data diving”.

beware of the logarithmic launching pad

I wondered what series of sizes to adopt. I decided to generate 100 images. 1km would be the minimum unit and we would reach the BIG square covering all the region, the width of Ile-De-France being precisely 154 451 km according to GÉOFLA. The distribution, as illustrated, will be logarithmic.

The same sensation appears when you dive head first in a municipal pool

Above is a GIF presenting the different calculation canvases.

Here is the video with land cover percentages encoded as colors, with reversed order, from the biggest to the smallest unit. The changing colors make it kind of animated dancefloor.

Dataekwondo

This illustration, that I had made for a GéoTribu post, illustrates the usual workflow I adopt, from production to design. I generated the distinct grids associated with the land cover percentage values with R. I designed all the layers with QGIS. With the same tool, I generated an atlas. With this atlas, I created a video with ffmpeg.

Technical details

The QGIS part may interest the QGIS fans as it uses a new QGIS 2.12 functionnality : data defined control over map layers and style presets

Generation : R

I created a function in R that generates a grid. The hack consists in polygonizing a raster layer, of which the resolution determinates the square size :

creerGrille = function(pol, size=50000) {
r = raster(extend(extent(pol), size)); res(r) = size; r[]=1
grille = rasterToPolygons(r)
grille = grille[which(gIntersects(grille, pol, byid=T)), ]
grille = grille[-which(gTouches(grille, reg, byid=T)), ]
return(grille)
}

The extract function crosses the vector grid and the 2012 CLC raster. The table function gives the number of pixels of each category :

  counts = extract(clc2, grille, 
method="simple",
fun=function(x,...)table(x),
na.rm=T)

The prop.table calculates the proportion of urban, natural and agricultural areas :

RGB = round(prop.table(counts[, 1:3], 1)*255)

Design : QGIS

color RGB code from columns, named “R”, “G”, and “B”

An expression defines the color of each square depending on its column values :

MultiQML plugin

I saved the QML style file, then used the multiQML plugin to apply it to the 99 other layers.

Data defined control over map layers and style presets

To get a static coverage layer, I duplicated my area of interest as many times as the number of pages (or layers) and assigned the square size value to the layer names. In the latest version of QGIS, you can generate an atlas with a list of layers on each page defined by an expression.

Also, I used some generative design tricks to get the animated halo effect around the IDF region.

Shapeburst fill and data-defined color

The color follows the chromatic color circle, covering all the rainbow colors :
color_hsv(scale_linear(@atlas_featurenumber,1,100,0,360),40,27)

riding the wave of generative design



sin wave shade distance


while a sin function defines the shade distance. This way, we get a “flame” effect with increasing then decreasing shade width. scale_linear(sin(scale_linear(@atlas_featurenumber,1,100,0,pi()*2)),0,1,5,10)

Video making: ffmpeg

Finally, this command line generates the video from images :
ffmpeg -framerate 25 -r 10 -i output_%d.jpg -vcodec mpeg4 -b 50000k -qscale:v 1 dancingVideo.avi

The code

The functions

The processing code

Credits

Images from The Noun Project : diver by Claire Jones, honey by Angela Horton, sandwich by Gareth Servant, biker by Les vieux garçons, ice cream by Rafaël Massé.

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.