Archives de catégorie : urbanism

[Sextante Plugin] An R script for Neighborhood Detection

Des bâtiments regroupés en bloc à Paris

Artice in English

[Plugin Sextante] Un script R pour détecter les voisinages : Cet article explique comment réaliser un petit script R via le plugin Sextante pour QGIS afin d’établir des relations de voisinage entre objets en fonction de la distance qui les sépare. Il comporte une petite introduction sur l’écosystème QGIS, sur R aussi. Et puis bien sûr, des détails sur l’écriture d’un script R pour le plugin Sextante GIS Plugin.


Content:

  • An opensource spatial constellation
    • R
    • QGIS as a solar planet
    • Sextante GIS plugin as a glue
    • R and QGIS
  • Performing a neighborhood analysis with Sextante plugin
    • The R example scripts
    • Neighborhood by distance: my R script
    • Writing an R script
    • The complete Sextante R script
  • Usage in QGIS
    •  Quiet places where to sleep/live in Paris
    • Styling and visualizing
    • Statistics
    • Convex hulls
Download the R script directly here.

An opensource constellation

http://commons.wikimedia.org/wiki/File:Solar_sys.jpg

R

R is a statistical language from which a huge number of libraries, representing a multicoloured palette of domains, have proliferated, from finance to ecology, from spatial to morphometry.

Its popularity is increasing very fast. Even the New York Times wrote about it. Also, an article describes how it is used at google.

More and more companies search for data scientists that master R, besides SAS, the proprietary equivalent tool.

In the spatial domain, there are some functions that you’ll never find in any opensource GIS software, be it GRASS or Sextante.

In this spatial R page, you might discover some tasks you’ve never heard of in the geospatial world, like spatial modeling, point patterns and small area estimations.

QGIS as a solar planet

More and more, QGIS tends to integrate famous tools. The first one was GRASS, then there was R, Sextante GIS, recently there was the Orfeo Toolbox for Imagery and Remote Sensing.

All of this constitutes an opensource spatial galaxy/ecosystem in which every planet (here, software) develops itself in contact with other planets, also stars (the technologies that gain popularity like Mapnik, Leaflet).

Sextante GIS in QGIS as a glue

Initially, Sextante GIS is a geospatial library.
Then, it became a plugin for QGIS, developed by Victor Olaya. It not only included Sextante functionalities but also functions from Orfeo Toolbox, R, GDAL, GRASS,…

The goal of this plugin is to indifferently execute functions from different libraries in a graphical environment since not everybody is familiar with command lines.

Also, the Sextante QGIS Plugin integrates a modeler in which you can chain functions and make QGIS a kind of automated factory for spatial tasks. Remember GRASS also has a modeler.

I think that with the improvement of his Sextante modeler, QGIS could become a very powerful spatial ETL and could compete with Spatial Data Integrator/Talend Open Studio, GeoKettle or FME (when it deals with spatial treatments).

R and QGIS

There are many plugins that use the R library: manageR by Carson Farmer, spqr by Barry Rowlingson, SDA4PP (Spatial Data analysis For Point Patterns, Home Range Estimation

manageR is an R editor inside QGIS. Knowing R is a prerequisite for using this plugin.

spqr produces graphics for stat visualizations .

– SDA4PP is a suite of functions for the analysis of point patterns: are the points clustered, seggregated? 

I think that the most interesting R plugins for QGIs are those that guide the user through a graphical interface and perform complex R tasks behind.

R in Sextante QGIS plugin

When you install the Sextante GIS plugin, you must activate the R scripts to get them available in the Sextante dock

– Also, you have to mention where R is installed.

– I preferred to move all the R scripts that initially were in C:Documents and Settings[username].qgispythonpluginssextanterscripts to C:/R/QGIS_SEXTANTE to make them more accessible. I changed the R scripts folder in the configuration dialog as you can see in the image above.

The Sextante R Scripts

If you don’t know anything about R and if you want to use an R script but you don’t have the approriate package to use it, all you have to know is how to install this package. It is very simple. I recommend you RStudio for this (if a library has not been installed in your standard R, it won’t be available in Sextante QGIS R scripts)

Using an R script doesn’t require any knowledge in R programming but it’s preferrable to understand how it is made, and to see which function it uses. You can refer yourself to the documentation to get some details, by using R-Seek. Be studious: some functions sometimes require some theoretical knowledge “that could lead either to years of therapy, or to your Ph.D”.

Performing a neighborhood analysis with Sextante plugin

The R example scripts

There a already 9 R scripts in the package as you can see in the image above. The system of these scripts are input-output systems where the inputs are always spatial data and where outputs can either be spatial data, plots (images), or console outputs.

– With the Ripley Rason script, you can create some envelopes that are a bit like convex hulls but stick more to the group of elements.

– You can create create either regular, either random points based on a polygon (generally, we use a grid)

– You can analyze the clustering of points using the F, G, K Functions scripts

– You can check whether you values are random – that means following a normal distribution – with the Kolgomorov-Smirnov normality test script

– You can make kind of same analysis but with spatial quadrats: consider a study zone with points. You divide this zone into a certain number of squares (quadrats). You count the points inside each square. Finally, you make a test to see if your points are randomly located (it’s called the null hypothesis of Complete Spatial Randomness).

– There is an R script that makes a raster histogram. It illustrates that you can have rasters as inputs, not only vectors.

Neighborhood by distance: my R script

As you can see, the R scripts mainly focus on point patterns. I wanted to implement an R script that I had already written, but in an R environment, that focused on areas.

The goal of this R script is to specify a distance below which some polygons (for instance, buildings) are considered neighbors. Once it finds the neighborhood relationships, it counts the neighbors and affects each building to the group of neighbors it belongs to.

The output of this script is the same layer as the base one but enriched with two more attributes: neighbors count and neighbor group.

This script can be useful to determine the structure of the habitat and detect the buildings that are isolated or highly grouped based on a mimimum distance threshold.

Here, you’ll find a complete description of the spdep package.

Here is the script, in R language:

Writing an R script

The goal was to transform this R script into a Sextante-QGIS R scripts.

If you click on “create a new R script”, an edition window opens and you can have an assistant through the “edit script help” button. Unfortunately, I didn’t understand how this script help worked. So, I inspired myself from the example scripts.

To create a script, create a new file, in my case, neighborhood_by_distance.rsx and edit it.

An R script is structured in two parts:
– the configuration part where you set the input, output, options.
– the code part where you use the variables set in the configuration part inside a standard R script.

The configuration part

Group

Each script belongs to a group of scripts. It’s the group key, for instance,

##[datagistips]=group

means that the script belongs to the “datagistips” group. If the group doesn’t exist, it will automatically be created.

Input

The inputs are either vector, either raster. I don’t think you can mention you’d like a certain type of features: points, lines or polygons (maybe an improvement to be done as some analysises can only be done on a certain type of features?)

For vectors, you write:

##polygons=vector

polygons will be the label in the graphical window

For rasters:

##layer = raster

Options

In options, the user can indicate a numeric value or a string value.

For the numeric value

##distance=number 100

Here, the label will be distance and the default value is set to 100.
When the user will specify the value, he will type it manually or choose it from the different layer statistics

If the value is character, the user will have to type it.
 ##title=string France

Here, the label will be title and the default value is set to France

Fields

One very interesting thing is that you can use a field in a script. In the case of the Kolgomorov-Smirnov example script, you choose the numeric field from which to test the normality.

##field=field layer

Then, in a script

>lillie.test(layer[[field]])     

will integrate the field variable into the script

Outputs

For spatial data outputs, simply write
##output=output vector

Graphical outputs and console outputs will both be in HTML format

For graphical plots outputs, you must mention the tag

##showplots

(you could imagine producing some lovely plots with ggplot2)

For console output, you must precede the command by >

>lillie.test(layer[[field]])    
The complete Sextante R script

Below is the final R script. You can also download it here.

As you can see, you can instantly port your R code to a graphical environment so that every QGIS user, even a beginner, can use your process. Also, it’s easier than developing a python QGIS plugin.

The plugin then appears in the Sextante Toolbox

Usage example in QGIS

http://commons.wikimedia.org/wiki/File:Paris_sunset.JPG

Quiet places where to sleep/live of Paris

Here, I executed the script on an extract of buildings in the 1st district of Paris that I downloaded from the Geofabrik OpenStreetMap data repository

You could use it to determine which buildings are the most isolated, where you would probably sleep one night or live, far from the busy, crowdy life of the City of Lights.
The distance that I chose was 15 meters. In such a big city like Paris, that’s quite a lot.

The R script in action


Styling and visualizing

Once I had my neighbors layer, I could style it differently depending on the neighbors counts or the group it belonged to..

Clustered and isolated buildings on top of OSM data (with OpenLayers plugin)

Statistics

With the standard statistics tools incorporated in QGIS, you can analyze the spatial structure.

The results mention that there 24 groups of neighbors. The biggest group contains 23 buildings. Half of the neighbors contain between 0 and 7 neighbors.
You can also produce histograms of these values with the Statist plugin and see the distribution of counts for the neighbor groups

By using the request functionnality, I could locate buildings that where relatively isolated.

Buildings with less than 3 neighbors (distance of 15m) on top of Bing Aerial Imagery

Convex hulls

Finally, why not creating convex hulls to visualize the neighborhood groups and put some OSM data behind to see how it looks?

The output data on top of OpenStreetMap data (15 meters threshold)

When on top of OSM data, you can easily visualize the buildings groups. Remember that it all depends on the distance thereshold you used. With a 5m threshold, you get this:

The same but with a 5 meters threshold (nicer than with a 15m threshold, isn’t it?)

These spatial groups could probably reflect the different residences that exist in this place. They could be useful to plan a more detailed OSM investigation and mapping on the ground (to map residences in OpenstreetMap, use landuse=residential, name=[residence name], or a relation of type site)

Also, you could use these blocks as a base layer for some aggregated statistics: number of buildings, area, perimeter, lacunarity, fractal dimension, and so forth..

Le carroyage avec QGIS et le plugin QMarxan, c’est facile!

Article en français 

Gridding Paris with OpenStreetMap data
“Gridding for the dummies” 😉 with Quantum GIS and QMarxan: In this article, I explain the purpose of gridding, why it is such an efficient way to analyse data within a large community of actors. Actually, gridding can be considered as a transcription of the raster format to the vector one where one tile <=> one cell. Thanks to QMarxan, gridding with QGIS is really simple:  you don’t need any knowledge in PostGIS, GRASS and R to do it. That’s what we’ll see: in the end of the post, I detail how to constitute some gridded data through a case study aiming to process a building density analysis on the city of Paris. The article is in French, but the printscreens will be helpful for the understanding. 

All the material used in this article is free: Quantum GIS, OpenStreetMap for the buildings and INSEE (equivalent of US Census) for the population data.

Maillage / Carroyage

Guide du CERTU sur les traitements en carroyage

Le CERTU (Centre d’Etudes sur les Réseaux, les Transports, l’Urbanisme et les Constructions) a publié l’année dernière un guide en téléchargement gratuit appelé Traitements géomatiques par carreaux pour l’observation des territoires qui décrit, par le biais de cas d’étude, tout le potentiel que permet la constitution de données carroyées.

Le carroyage est un mode de représentation des données vecteur par carreau. Il s’effectue en croisant d’un côté une grille (soit rectangulaire, soit hexagonale en nid d’abeilles) avec des données métier. La précision avec laquelle s’effectuera l’analyse sera alors relative au pas choisi pour la grille. Ce pas diffère selon l’échelle d’analyse.

Densité d’habitants selon l’INSEE

On peut d’ores et déjà disposer de données carroyées sur Internet. C’est le cas des données INSEE sur le nombre d’habitants disponible en téléchargement avec un pas de 1 kilomètre sur la France entière et 200 mètres sur chaque région. Le Muséum National d’Histoire Naturelle propose, lui, une grille-type pour effectuer un calcul, par carreau de 5 ou 10km, du nombre d’espèces d’un certain type de plantes ou d’animaux. Inspire propose aussi des grilles-types.

En quoi ce mode de visualisation est-il si intéressant?

Potentiel du carroyage

– Le format vecteur est généralement choisi comme base de traitements car les utilisateurs du SIG sont plus familiers de l’exploitation de ce format. Pourtant, parfois, les calculs sur raster sont plus appropriés. Le carroyage permet de transcrire au format vectoriel des méthodes rasters. L’augmentation de la taille d’un carreau s’apparente à une augmentation de la résolution (rééchantillonnage). Il est tout à fait possible d’effectuer des dénombrements, de calculer des moyennes, des écarts-types au niveau de chaque carreau, ce qui s’apparente aux calculs focaux utilisés sur les rasters.

– La transformation de données fines, précises sous la forme de carreaux permet de vulgariser la donnée, de la simplifier, en sorte. Les mailles obtenues peuvent être confrontées à la fois à un niveau spatial: un carreau par rapport à ses voisins (par ex. étude de clusters), de même que dans le temps: état d’un carreau dans un état antérieur et dans un état actuel (par ex. étude de l’évolution de la population). C’est un mode de visualisation efficace qui “élimine le bruit” et découpe le territoire en unités égales et comparables.



– Une troisième raison qui en fait un support intéressant est la possibilité d’associer plusieurs dimensions aux carreaux, comme autant de colonnes et d’attributs. On pourra alors adjoindre à chaque carreau des valeurs issues de différents domaines: environnement, transport, habitat, et de différents types: quantitative, qualitative, ordinale. On pourra faire apparaître différentes thématiques sur les choroplèthes au gré des besoins. Cette multi-dimensionnalité est difficile à atteindre avec les rasters dans la mesure où les pixels ne peuvent être associés qu’à un nombre limité de valeurs, de plus, positives. Sous R, des objets appelés SpatialGridDataFrame et SpatialPixelsDataFrame permettent d’associer à chaque pixel de nombreux attributs, comme les formats vecteur, mais R est un outil complexe qui nécessite un certain temps de prise en main.

– Cette multi-dimensionnalité montre bien que le carroyage peut être un support d’analyse partagé entre des thématiciens de différents métiers. Des urbanistes, des botanistes peuvent se servir d’une même grille comme support pour leurs calculs puis confronter leurs analyses.

En définitive, l’accessibilité, la simplification, la multi-dimensionnalité font du carroyage un support intéressant pour des analyses et des représentations de données partagées.

Carroyage sous QGIS

Sous Quantum GIS, le carroyage est à la portée de tout géomaticien. On peut en réaliser sans utiliser ni PostgreSQL, GRASS, ni R même si je ne saurais que vous conseiller de vous investir dans l’apprentissage de ces derniers.

Nativement, on peut l’exécuter par le biais de la fonction “Création de Grille” puis “Points dans un Polygone” pour des cartes de densité.

Cependant, l’année dernière, un plugin est apparu sur le dépôt d’extension appelé QMarxan qui permet d’aller plus loin. Il s’agit d’un plugin QGIS d’interfaçage avec le logiciel Marxan. 


Voici la définition que l’on trouve sur le site des parcs du Canada quant au logiciel Marxan:
“Marxan a été mis au point par des chercheurs de l’Université de Queensland pour faciliter le rezonage du parc marin de la Grande Barrière et a par la suite été utilisé un peu partout dans le monde dans le cadre de projets de planification d’aires marines. Au moyen d’un algorithme d’optimisation, ce programme sélectionne des séries de secteurs qui permettent d’atteindre des objectifs de conservation tout en réduisant au minimum les incidences socio-économiques sur les utilisateurs des milieux marins et des collectivités côtières. Les données qui y sont entrées peuvent comporter un large éventail d’éléments/de valeurs liés aux aspects écologiques, socio-économiques, culturels et à l’expérience du visiteur, y compris les espèces, les habitats, les activités récréatives et les sites culturels.”

Dans cet article, nous n’exploiterons pas tout le potentiel du plugin. Nous verrons sa fonctionnalité principale, celle permettant de réaliser une grille et d’effectuer des calculs par carreau. Nous ne couvrirons pas l’exploitation des résultats sous Marxan.

Le bâti en Île-de-France

Le cas d’étude sera celui de l’analyse du bâti à Paris avec la constitution d’un carroyage décrivant la densité du bâti en terme de nombre et de surface de bâtiments. Puis, nous exploiterons le carroyage de l’INSEE afin de coupler les données concernant le bâti avec la population.

Téléchargement des données

Dans un premier temps, téléchargeons les limites administratives de la ville de Paris:
http://export.openstreetmap.fr/contours-administratifs/export-communes/

Puis les données du bâti à partir du site GeoFabrik:
http://download.geofabrik.de/osm/europe/france/ile-de-france.shp.zip

Effectuons une extraction des donnés sur la ville de Paris par le biais de l’extension “Requête Spatiale”. Le temps de traitement est très long!..

Par la suite, les données bâtiments buildings seront transformées en Lambert 2 Etendu afin de les coupler par la suite (fin du tutoriel) avec les données INSEE qui sont dans ce système de projection. La couche transformée se nomme buildings_LA2E.

Dans un premier temps, il faut installer l’extension QMarxan depuis les dépôts-tiers.

Réalisation d’une grille en nid d’abeilles

Une fois le plugin installé, lançons le module de génération de grille

J’aime bien les grilles en nid d’abeilles qui sont assez esthétiques. Réalisons-en une. Nous choisirons un pas de 200 mètres comme pour l’INSEE.

Création d’une grille hexagonale sur la base du fichier de bâti

Voici le résultat du carroyage. Notre grille est pour l’instant vide de données statistiques sur le bâti.

Bâti + grille en nid d’abeilles

Calcul de la surface, du nombre de bâtiments et de la surface moyenne par carreau

Effectuons le calcul de la surface de bâti par carreau.

Le résultat est intégré à la grille:

Afin d’avoir une bonne visualisation de l’ensemble, nous avons extrait les nids d’abeille comportant une surface supérieure à 0:



Surface de bâti en 10 classes selon la méthode de Jenks

Pour avoir le nombre de bâtiments par carreau, on doit tout d’abord générer une couche de points. Nous l’appellerons buildings_ct

Centroïdes des bâtiments superposés sur la surface bâtie


Le nombre de bâtiments par carreau peut être indicatif d’un type de bâti: collectif, résidentiel, etc…

Nous pouvons calculer la surface moyenne de bâti par carreau: surface/(nombre de bâtiments)

Utilisation de la calculatrice de champ

Voici différentes choroplèthes obtenues en superposant le résultat du carroyage avec les données OpenStreetMap grâce au plugin OpenLayers de QGIS.

Nombre de bâtiments, surface de bâti et surface moyenne sur couches OpenStreetMap Paris

Couplage données nombre d’habitants et bâti sur la base d’une grille pré-existante

On peut aussi lancer les calculs en se basant sur une grille autre, par exemple celle de l’INSEE, et ensuite coupler les données habitants et bâti.

Les zones d’activité pourraient ainsi apparaître comme des zones de bâti comprenant peu d’habitants.

Téléchargeons le carroyage de l’INSEE qui se trouve en Lambert 2 Etendu. Il contient les estimations du nombre d’habitants. Un préalable, comme au départ, est de ne prendre que les carreaux situés sur Paris.

Associons ces données du bâti avec ces données de population en utilisant le carroyage existant:

 Le résultat est alors intégré dans la grille de l’INSEE:

ajout de la colonne area_built au carroyage de l’INSEE

Ouvrons le fichier dbf issu du carroyage sous OpenOffice, enregitrons-le au format ods et réalisons un graphique bivarié avec en abcisse la surface de bâtiments et en ordonnée le nombre d’habitants.

fichier ods obtenu à partir du dbf

Nous pouvons voir quelle relation existe entre le nombre d’habitants et la surface de bâti, ainsi qu’identifier des carreaux qui échappent à la tendance principale.

Identification d’un point échappant à la “tendance”

Remarque: A vrai dire, il est difficile de dégager une quelconque tendance à partir de ce graphique car dans le cas du bâti, la surface moyenne ne suffit pas à prédire le nombre d’habitants. Il faut aussi intégrer la hauteur du bâtiment ou son nombre d’étages. Cela dit, le but de cet exemple est de montrer comment le carroyage permet de confronter deux (voire plus) variables de manière efficace par le biais d’une analyse bivariée (ou multi-variée).



 Nota Bene: dans la première version de QMarxan, l’utilisation de grilles rectangulaires génère des erreurs. Aussi la dernière version a-t-elle été utilisée pour l’intégration des informations du bâti dans le carroyage de l’INSEE.

Pour ceux qui souhaitent aller plus loin dans le carroyage, je vous conseille d’utiliser R avec la librairie sp. La fonction over( ) permet de croiser dans n’importe quel ordre des données vectorielles et rasters ainsi que d’effectuer une infinité de calculs à l’intérieur de chaque carreau, ceci grâce à une fonction personnelle fn: over(x, y, returnList = FALSE, fn = NULL, …). On peut alors envisager de calculer l’écart-type ou même l’indice écologique de diversité de Shannon. Absolument tout est possible!
Plus d’infos sur les méthodes d’agrégation spatiale avec R ici: http://cran.r-project.org/web/packages/sp/vignettes/over.pdf