Formation au carroyage et lissage spatial sur R

TUTORIEL

Author

Kévin Milin et Solène Colin

Modified

September 2026

1 Introduction

1.1 Objectifs du TP

En 2018, le PSAR Analyse Urbaine (ancêtre de la section Analyse Urbaine à la direction générale de l’Insee), a développé un package R, nommé btb (auteurs : Arlindo Dos Santos et François Sémécurbe).

Sa principale fonction, btb_smooth, permet de réaliser très facilement un carroyage et un lissage sur des données géolocalisées.

À partir de données ponctuelles, nous allons apprendre, en utilisant le langage R :

À carroyer les informations.
À lisser des densités, des moyennes, des taux et des quantiles.
À calculer un indicateur sur une zone à façon à partir des données ponctuelles et de données carroyées de l’Insee.

Liens utiles

Code de la formation : https://github.com/InseeFrLab/formation-r-lissage-spatial
Site web des supports de formation : https://inseefrlab.github.io/formation-r-lissage-spatial

Crédits (Division Statistiques et Analyses Urbaines)

Cette formation s’inspire d’une formation élaborée par Arlindo Dos Santos en 2018-2019
Elle a été refondue par Julien Pramil principalement, et Kim Antunez en co-animatrice, en 2021 en y intégrant, notamment, les données en open-data DVF
Elle est animée depuis 2024 par Solène Colin

1.2 Avertissements

1.2.1 Secret statistique

Avant toute diffusion auprès des partenaires, il faut bien veiller à respecter :

le secret
- primaire
- secondaire
- fiscal
les conventions établies avec les fournisseurs des données

1.2.2 Qualité des cartes

Pour simplifier les programmes présentés dans ce TP, les représentations graphiques ne respectent pas toutes les règles élémentaires de la sémiologie cartographique, résumées par cette infographie :

Auteur : Timothée Giraud, auteur de la librairie mapsf

Il faudra bien veiller à appliquer ces règles de sémiologie sur les cartes que vous réaliserez ultérieurement en diffusion.

1.2.3 Système de projection

Voici les systèmes de projection que vous pouvez régulièrement rencontrer pour la métropole :

Nom	Description	Code EPSG
Lambert93	Système de projection officiel pour la métropole	2154
LAEA	Système de projection européen	3035
WGS84	GPS (utile pour utiliser Leaflet)	4326

2 Configurations

2.1 Chargement des librairies

Cinq librairies principales sont nécessaires pour ce TP.

sf pour manipuler des fichiers spatiaux : importation, projections, géotraitements… (fonctions commençant par st_)
dplyr pour le traitement des données, en particulier l’agrégation géographique
mapsf pour réaliser des cartes dans RStudio (fonctions commençant par mf_)
mapview (reposant sur leaflet) pour réaliser des cartes interactives (fond de carte OpenStreetMap)
btb pour le carroyage et lissage (fonctions commençant par btb_).

Remarque : Le choix de dplyr plutôt que data.table se justifie ici du fait de sa forte compatibilité avec les objets géomatiques.

Code

## Liste des librairies utilisées
packages <-  c("devtools", "sf", "dplyr", "mapsf", "mapview", "leaflet", "btb")

## Vérifier si la librairie est installée, si non l'installer, puis la charger
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE, quiet = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

2.2 Chargement de la base

Pour ce TP, nous utilisons la base « Demandes de Valeurs Foncières » (DVF), qui recense l’ensemble des ventes de biens fonciers réalisées au cours des cinq dernières années (hors Mayotte et Alsace-Moselle). Les biens concernés peuvent être bâtis (appartement et maison) ou non bâtis (parcelles et exploitations). Les données sont produites par la Direction Générale des Finances Publiques (DGFiP). Elles proviennent des actes enregistrés chez les notaires et des informations contenues dans le cadastre.

À partir de cette source, nous avons constitué une base de données qui s’intéresse uniquement au périmètre de la petite couronne parisienne (départements 75, 92, 93 et 94) et au millésime 2021.

Remarque importante : La base a été filtrée de manière à être la plus pédagogique possible pour cette formation. Elle n’est donc ni représentative de la réalité, ni fidèle au champ proposé par le producteur. En somme, mieux vaut donc de pas appuyer vos investissements immobiliers sur les résultats de ce TP !

Pour information, le code ayant permis de constituer la base de données est disponible sur le répertoire github de la formation.

La base contient les 8 variables suivantes :

id_mutation : identifiant unique de la vente
date_mutation : date de la vente
type_local : appartement ou maison
nombre_pieces_principales : nombre de pièces dans le logement
valeur_fonciere : prix de vente
surface_reelle_bati : surface du logement
x : longitude (en projection Lambert 93)
y : latitude (en projection Lambert 93)

La base pédagogique utilisée dans ce tutoriel se nomme ventesImmo_couronneParis.RDS, et est stockée sous Minio, dans le “bucket public” : s3/projet-formation/r-lissage-spatial/. Nous téléchargeons cette base grâce à son URL public :

Code

url_bucket <- "https://minio.lab.sspcloud.fr/projet-formation/r-lissage-spatial/"
object <- "ventesImmo_couronneParis.RDS"

# Charger la source de données (variable nommée dfBase) depuis l'URL
url_file <- url(paste0(url_bucket, object))
dfBase <- readRDS(url_file)

On peut ensuite manipuler la base :

Code

# Visualiser les premières lignes de la base
head(dfBase, 3) ; dim(dfBase)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y
1          480000                  64 647357.3 6868635
2          345000                  43 644483.5 6867695
3          384000                  41 648001.8 6866153

[1] 34489     8

2.3 Chargement des données cartographiques

2.3.1 Fond de carte du territoire à étudier

Nous allons charger le contour géographique des départements de la petite couronne parisienne. En cartographie, ces fichiers sont appelés des “couches vectorielles” (à ne pas confondre avec les couches “raster” qui correspondent, par exemple, aux fonds de carte satellite OpenStreetMap ou Google Maps).

Différents types de fichiers permettent de stocker des couches vectorielles. Les deux principaux sont le Geopackage (.gpkg) et le Shapefile (.shp). Nous recommandons l’utilisation du .gpkg pour les raisons suivantes :

Un .gpkg est un fichier unique alors qu’un fichier .shp tient dans 5 fichiers séparés et interdépendants.
Le format Géopackage est libre et ouvert, contrairement au Shapefile.
Un Shapefile impose des noms de variables de moins de 10 caractères, ce qui peut être pénible en pratique.
Le géopackage est devenu le format standard dans Qgis depuis la version 3 (logiciel libre de cartographie préconisé à l’Insee).

Pour aller plus loin sur la comparaison des formats, voir ici.

Récupérons le contour du territoire étudié avec la fonction st_read du package sf :

Code

# Charger le fond de carte du territoire étudié
chemin_file <- paste0(url_bucket, "depCouronne.gpkg")
depCouronne_sf <- st_read(chemin_file)

Reading layer `depCouronne' from data source 
  `https://minio.lab.sspcloud.fr/projet-formation/r-lissage-spatial/depCouronne.gpkg' 
  using driver `GPKG'
Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 v1 / Lambert-93

Code

# Visualisation de la couche vectorielle
head(depCouronne_sf)

Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 v1 / Lambert-93
  code           libelle reg surf                           geom
1   75             Paris  11  105 MULTIPOLYGON (((660897 6860...
2   92    Hauts-de-Seine  11  176 MULTIPOLYGON (((648796 6847...
3   93 Seine-Saint-Denis  11  237 MULTIPOLYGON (((659428 6861...
4   94      Val-de-Marne  11  245 MULTIPOLYGON (((656908 6846...

Code

plot(depCouronne_sf$geom)

Code

# On renomme la variable "geom" en "geometry" par convention
depCouronne_sf <- depCouronne_sf %>% rename(geometry = geom)

Une fois les départements chargés, nous allons sélectionner le contour de la commune de Paris.

Code

# Sélection de Paris
paris_sf <- depCouronne_sf[depCouronne_sf$code == "75", ]

# Visualisation de la nouvelle couche vectorielle
head(paris_sf)

Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
Projected CRS: RGF93 v1 / Lambert-93
  code libelle reg surf                       geometry
1   75   Paris  11  105 MULTIPOLYGON (((660897 6860...

Code

plot(paris_sf$geometry)

Pour être certain que les données et le territoire soient dans le même système de projection, il est possible de transformer ce dernier à l’aide de la fonction st_transform.

Code

# Connaître le système de correction d'une couche cartographique
st_crs(depCouronne_sf)$epsg ; st_crs(paris_sf)$epsg

[1] 2154

[1] 2154

Code

# Transformer la projection en Lambert 93 (epsg 2154)
depCouronne_sf <- st_transform(depCouronne_sf, crs = 2154)
paris_sf <- st_transform(paris_sf, crs = 2154)

2.3.2 Territoires englobants, sélection des données à lisser

Pour éviter les effets de bord, il faut sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros). Autrement, si on ne sélectionne que les ventes immobilières situées dans Paris intramuros, et qu’on lisse ce nuage de points, les zones situées à proximité du périphérique seront artificiellement peu denses en ventes immobilières. En effet, elles seront situées à proximité de zones vides en ventes immobilières, à savoir les communes limitrophes. Pourtant, dans la réalité, il y a bien des ventes réalisées au delà de Paris, et qui influencent la densité lissée près du périphérique.

2.3.2.1 Première méthode : sélection géométrique des ventes

La première méthode est géométrique/vectorielle. Elle consiste à utiliser nos données comme un ensemble de points géolocalisés et procéder à des intersections géographiques.

Plus précisément :

On transforme nos observations en points vectoriels ;

Code

sfBase <- dfBase %>% mutate(lon = x, lat = y) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 2154)

On crée une zone tampon (buffer) autour du territoire d’intérêt, avec une marge (ici 2 000m), sous la forme d’un objet sf vectoriel ;

Code

buffer_sf <- st_buffer(paris_sf, dist = 2000)

Remarque: Pour la zone tampon, il convient de prendre une marge légèrement plus grande que le rayon de lissage envisagé, ceci afin d’éviter les effets de bord tout en limitant les temps de calcul.

On repère les observations comprises dans cette zone tampon par intersection géographique.

Code

sfBase_filtre <- st_join(sfBase, buffer_sf, left = F)
nrow(sfBase_filtre) ; head(sfBase_filtre, 3)

[1] 22221

Simple feature collection with 3 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 645132.7 ymin: 6864646 xmax: 648001.8 ymax: 6866153
Projected CRS: RGF93 v1 / Lambert-93
  id_mutation date_mutation  type_local nombre_pieces_principales
3 2021-447025    2021-01-08 Appartement                         3
5 2021-447029    2021-01-05 Appartement                         2
6 2021-447030    2021-01-15 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y code libelle reg surf
3          384000                  41 648001.8 6866153   75   Paris  11  105
5          407200                  24 646929.9 6864730   75   Paris  11  105
6         1040000                  90 645132.7 6864646   75   Paris  11  105
                  geometry
3 POINT (648001.8 6866153)
5 POINT (646929.9 6864730)
6 POINT (645132.7 6864646)

Schéma explicatif des zones sélectionnées :

Code

# Échantillon (en dehors et dans le buffer)
sfBase_sample <- sfBase[sample(1:nrow(sfBase), 2000), ] # échantillon 2000 observations
sfBase_filtre_sample <- st_join(sfBase_sample, buffer_sf, left = F) # dans le buffer

# Cartographie pédagogique
mf_map(buffer_sf, col = "yellow")
mf_map(paris_sf, col = "blue", add = T)
mf_map(sfBase_sample, col = "red", add = T)
mf_map(sfBase_filtre_sample, col = "green", add = T)

Ce schéma présente :

en bleu, le territoire d’intérêt (Paris)
en jaune, la zone tampon du territoire d’intérêt, avec une marge de 2 km
en rouge et vert, 2000 points tirés au sort dans la base initiale
en vert, les points (parmis les 2000 tirés) qui sont dans le buffer d’intérêt et qu’on va garder.

Cette méthode filtre efficacement, mais elle peut-être lourde d’un point de vue calculatoire. Avec des données volumineuses, il convient au minimum de faire un premier filtrage (par exemple avec la seconde méthode).

2.3.2.2 Seconde méthode : sélection non-géométrique des ventes

Ici, nous allons sélectionner les ventes immobilières appartenant à un grand rectangle englobant Paris.

L’idée est de pouvoir réduire la quantité de ventes à considérer avant le lissage, uniquement en filtrant les coordonnées x et y comprises dans ce rectangle. Cette méthode est très efficace d’un point de vue calculatoire, mais filtre moins, et requiert de construire au préalable un rectangle (“bbox”) adapté.

Code

# Création d'une bbox autour du territoire
bbox <- st_bbox(paris_sf) ; bbox

   xmin    ymin    xmax    ymax 
 643076 6857499  660897 6867034

Code

# Création d'un buffer de la bbox, avec une marge de 2000 mètres
marge <- 2000
bufferBbox <- bbox
bufferBbox[["xmin"]] <- bufferBbox[["xmin"]] - marge
bufferBbox[["xmax"]] <- bufferBbox[["xmax"]] + marge
bufferBbox[["ymin"]] <- bufferBbox[["ymin"]] - marge
bufferBbox[["ymax"]] <- bufferBbox[["ymax"]] + marge
bufferBbox

   xmin    ymin    xmax    ymax 
 641076 6855499  662897 6869034

Pour sélectionner les données d’intérêt pour le lissage, on filtre les observations dont les coordonnées sont comprises à l’intérieur de la grande bbox. Ce filtre est très efficace computationnellement car il ne dépend que des valeurs numériques prises par les variables de longitude et de latitude et n’utilise pas les propriétés vectorielles des données qui sont des attributs plus chronophages à utiliser.

Code

# Ne garder que les données dans le rectangle englobant, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= bufferBbox["xmin"] &
                          dfBase$x <= bufferBbox["xmax"] &
                          dfBase$y >= bufferBbox["ymin"] &
                          dfBase$y <= bufferBbox["ymax"], ]

# nb de lignes avant/après le filtre
nrow(dfBase_filtre) ; nrow(dfBase)

[1] 24419

[1] 34489

Schéma explicatif des zones sélectionnées :

Code

# Échantillon (en dehors et dans le buffer)
sfBase_sample <- sfBase[sample(1:nrow(sfBase), 2000), ] # échantillon 2000 observations
sfBase_sample_filtre <- sfBase_sample[sfBase_sample$x >= bufferBbox["xmin"] &
                                        sfBase_sample$x <= bufferBbox["xmax"] &
                                        sfBase_sample$y >= bufferBbox["ymin"] &
                                        sfBase_sample$y <= bufferBbox["ymax"], ] # dans le buffer

# Petit rectangle vectoriel
bbox_sf = st_sf(geometry = st_as_sfc(bbox), crs = 2154)

# Grand rectangle vectoriel
bufferBbox_sf = st_sf(geometry = st_as_sfc(bufferBbox), crs = 2154)

# Cartographie pédagogique
mf_map(bufferBbox_sf, col = "yellow")
mf_map(bbox_sf, col = "grey", add = T)
mf_map(paris_sf, col = "blue", add = T)
mf_map(sfBase_sample, col = "red", add = T)
mf_map(sfBase_sample_filtre, col = "green", add = T)

Ce schéma présente :

en bleu, le territoire d’intérêt (Paris)
en gris, la bbox du territoire d’intérêt (à savoir le plus petit rectangle englobant cette zone)
en jaune, la bbox élargie avec une marge de 2km. Autrement dit, la zone tampon autour de la bbox permettant de prendre en compte des observations au-delà de la seule zone d’intérêt, et d’ainsi éviter les effets de bord au moment du lissage
en rouge et vert, 2000 points tirés au sort dans la base initiale
en vert, les points (parmis les 2000 tirés) qui sont dans le buffer d’intérêt et qu’on va garder.

3 Carroyage de données

Avant de lisser les données ponctuelles, on peut représenter ces données sous forme carroyée afin de se les approprier. Il faut pour cela :

Découper le territoire en carreaux de 200 mètres à partir de l’origine du référentiel, puis associer chaque point (1 point = 1 vente géolocalisée) au centroïde de son carreau. La fonction btb_add_centroids réalise ce travail, elle prend pour paramètre :

pts : un tableau avec les colonnes x et y (coordonnées des points) ;
iCellSize : la taille des carreaux (longueur du côté, en mètres).

Code

iCellSize = 200 # carreaux de 200m
points_carroyage <- btb_add_centroids(pts = dfBase_filtre, iCellSize = iCellSize) 
head(points_carroyage, 3)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y x_centro y_centro
1          480000                  64 647357.3 6868635   647300  6868700
2          345000                  43 644483.5 6867695   644500  6867700
3          384000                  41 648001.8 6866153   648100  6866100

Agréger les données sur chaque centroïde de la grille. En d’autres termes, compter le nombre de ventes par carreau :

Code

points_carroyage <- points_carroyage %>% 
  group_by(x_centro, y_centro) %>% count(name = "nbVentes")
head(points_carroyage, 3)

# A tibble: 3 × 3
# Groups:   x_centro, y_centro [3]
  x_centro y_centro nbVentes
     <dbl>    <dbl>    <int>
1   641100  6857700        1
2   641100  6857900        5
3   641100  6858300        1

Passer d’une table de centroïdes à une table de carreaux vectoriels grâce à la fonction btb_ptsToGrid qui attend comme paramètres :

pts : un tableau avec les colonnes x_centro et y_centro représentant les coordonnées des centroïdes de la grille ;
sEPSG : le code epsg du système de projection utilisé ;
iCellSize : la taille des carreaux (longueur du côté, en mètres).

Code

carreaux_large <- btb_ptsToGrid(pts = points_carroyage,
                                sEPSG = 2154, iCellSize = iCellSize)
head(carreaux_large, 3)

Simple feature collection with 3 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 641000 ymin: 6857600 xmax: 641200 ymax: 6858400
Projected CRS: RGF93 v1 / Lambert-93
# A tibble: 3 × 4
# Groups:   x_centro, y_centro [3]
  x_centro y_centro nbVentes                                            geometry
     <dbl>    <dbl>    <int>                                       <POLYGON [m]>
1   641100  6857700        1 ((641000 6857800, 641200 6857800, 641200 6857600, …
2   641100  6857900        5 ((641000 6858000, 641200 6858000, 641200 6857800, …
3   641100  6858300        1 ((641000 6858400, 641200 6858400, 641200 6858200, …

Se restreindre au champ des carreaux intersectant Paris

Code

carreaux <- carreaux_large %>% st_join(paris_sf, left = FALSE)

# nb de carreaux avant/après filtre
nrow(carreaux_large) ; nrow(carreaux)

[1] 4107

[1] 1754

On peut alors carroyer les ventes dans Paris :

Code

mf_map(x = paris_sf, col = NA, lwd = 4)
mf_map(x = carreaux,
       type = "choro",
       var = "nbVentes",
       breaks = "quantile",
       nbreaks = 5,
       leg_val_rnd = 1,
       add = TRUE)
mf_layout(title = "Carroyage du nombre de ventes",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Le carroyage permet d’avoir un premier aperçu des données le plus fidèle aux données initiales

Il peut être également utilisé pour simplifier les données avant lissage, afin de rendre le calcul de lissage moins long en cas de très grand nombre d’observations (cf 4.1.6).

4 Lissage

4.1 Calcul de densité

Un lissage simple à réaliser correspond au calcul d’une densité, à savoir une quantité par unité de surface (un carreau). Dans cette section, nous allons calculer la densité de ventes de logements dans la ville de Paris au cours de l’année 2021.

4.1.1 Grille automatique ; carreaux 200m ; rayon de lissage 400m

Pour ce premier lissage, nous choisissons des carreaux de 200 mètres et un rayon de lissage de 400 mètres.

Toujours pour éviter les effets de bord, le lissage sera effectué sur une zone plus large que la zone d’intérêt (utilisation d’un buffer), via la base de données dfBase_filtre construite précédemment.

Nous créons préalablement une variable nbObsLisse qui vaut 1 pour chaque observation (chaque logement vendu en 2021) et qui sera lissée avec la fonction btb_smooth.

Code

# Variable à lisser = nombre d'observations = nombre de ventes
dfBase_filtre$nbObsLisse <- 1

# Visualiser les premières lignes de la base
head(dfBase_filtre, 3)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y nbObsLisse
1          480000                  64 647357.3 6868635          1
2          345000                  43 644483.5 6867695          1
3          384000                  41 648001.8 6866153          1

La fonction de lissage btb_smooth intègre automatiquement l’étape de génération de la grille de carreaux. Puis, chaque carreau est affecté de la densité lissée en son centroïde.

Cette fonction nécessite l’utilisation de 4 paramètres :

pts : le tableau des données à lisser. Il doit nécessairement contenir les coordonnées (soit les colonnes x et y, ou bien une colonne geometry), et 1 à n colonnes numériques (variables à lisser) ;
sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé ;
iCellSize : un entier indiquant la taille des carreaux ;
iBandwidth : un entier indiquant le rayon de lissage.

Code

sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 400)

# visualiser les premières lignes des données lissées, et le nombre de lignes
head(sfCarrLiss, 3) ; nrow(sfCarrLiss)

Attention, la fonction retourne une erreur en cas de :

présence d’une variable non-numérique ;
valeur(s) absente(s) dans les coordonnées.

Ci-dessous la grille carroyée obtenue (hors données lissées). On remarque qu’elle épouse le périmètre des points en entrée de la fonction. Par défaut, elle va même un peu au-delà (voir partie 4.1.5).

Code

mf_map(x = paris_sf, lwd = 8, col = NA, border = "wheat")
mf_map(x = sfCarrLiss, col = NA, lwd = 1, add = T)
mf_layout(title = "btb_smooth génère une grille carroyée")

Une propriété particulièrement importante de la fonction de lissage btb_smooth est qu’elle est conservative. Cela signifie que nous avons la même somme des variables additives sur le champ géographique concerné avant et après lissage.

Code

# nb de ventes dans la petite couronne vs nb lissé de ventes dans les carreaux
sum(dfBase_filtre$nbObsLisse) ; sum(sfCarrLiss$nbObsLisse)

[1] 24419

[1] 24419

Remarque : veillez toujours à respecter le secret statistique au moment de la publication de vos cartes ! Par exemple, si vous utilisez des données issues des bases fiscales, vous ne pouvez représenter des informations sur des carreaux comportant moins de 11 observations, même si ce nombre est lissé. Pour ce faire, vous pouvez utiliser le filtre ci-dessous :

Code

# Exclure les carreaux ne respectant pas le secret statistique
sfCarrLiss <- sfCarrLiss[sfCarrLiss$nbObsLisse >= 11, ]

(Pour garder une carte complète, nous conservons la table sans filtre dans ce tutoriel).

Affichons maintenant le nombre lissé de ventes de logements à Paris en 2021 suite à ce premier lissage.

Quelques informations préalables :

On cherche à analyser le phénomène sur le périmètre de la ville de Paris ;
Ainsi, on lisse les données de ventes sur un périmètre plus large que la ville de Paris afin d’éviter les effets de bord aux frontières de la commune. En effet, la base dfBase_filtre contient toutes les ventes comprises dans un buffer de 2 km autour de Paris.
Suite au lissage, une grille carroyée est produite comportant les données lissées. Les carreaux vont bien au-delà de la commune de Paris. Néanmoins, pour la représentation, on ne sélectionne que les carreaux intersectant Paris (grâce à la fonction st_join).
Pour cartographier les carreaux :
- On réalise une carte de type choroplèthe ;
- Avec une méthode de discrétisation de la densité lissée (par quantiles par exemple) ;
- Avec 5 classes (c’est-à-dire 5 couleurs).

Code

# Filtrage des carreaux lissés intersectant la ville de Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# Carte lissée
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = "quantile",
       nbreaks = 5)
mf_map(x = paris_sf, lwd = 4, col = NA, border = "black", add = TRUE)
mf_layout(title = "Lissage avec rayon de 400m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque : la variable lissée nbObsLiss correspond au nombre de ventes (observations) par carreau (ici de 200 mètres). Ainsi, il est possible d’obtenir une densité au km² en multipliant la variable nbObsLiss par 25 dans le cas présent.

4.1.2 Grille automatique ; carreaux 200m ; rayons de lissage 600m et 1 km

Il n’est pas possible de connaître a priori la taille optimale du rayon de lissage. Il faut donc en essayer plusieurs pour trouver le meilleur compromis entre précision et généralisation. Nous allons donc faire le même lissage que précédemment avec un rayon de 600m, puis de 1 km, au lieu de 400m.

Code

# reproduction pour 600m
sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 600)

# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = "quantile",
       nbreaks = 5)
mf_map(x = paris_sf, 
       lwd = 4,
       col = NA, border = "black", add = TRUE)
mf_layout(title = "Lissage avec rayon de 600m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Code

sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 1000)

# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# Carte lissée
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = "quantile",
       nbreaks = 5)
mf_map(x = paris_sf, 
       lwd = 4,
       col = NA, border = "black", add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

À mesure que nous augmentons le rayon de lissage, la carte révèle des aspects structurels des données. Mais ceci se fait au détriment des spécificités locales visibles seulement avec un petit rayon de lissage. Le choix revient dont au statisticien-géographe, en fonction de sa connaissance des données et des objectifs recherchés.

Enfin, la grille carroyée est ici visible à des fins pédagogiques, il est possible et même conseillé de la supprimer sur les cartes publiées.

Code

# pour automatiser les cartes, on fait une fonction
creation_carte <- function(data, var, titre){
  mf_map(x = data, 
         type = "choro",
         var = var,
         breaks = "quantile",
         nbreaks = 5, 
         leg_val_rnd = 1, 
         border = NA) # enlever la grille
  mf_map(x = paris_sf, 
         lwd = 4,
         col = NA, border = "black", add = TRUE)
  mf_layout(title = titre,
            credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
}

creation_carte(sfCarrLiss_paris, "nbObsLisse", titre = "Lissage avec rayon de 1000m, sans grille")

4.1.3 Grille automatique ; carreaux 50m ; rayon de lissage 1 km

Pour éviter l’aspect “granuleux”, on pourrait aussi réduire la taille des carreaux (ci-dessous à 50 m). Attention toutefois aux temps de calcul.

Code

sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 50, 
                         iBandwidth = 1000)

# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# carte lissée
creation_carte(sfCarrLiss_paris, "nbObsLisse", 
               titre = "Lissage avec pas de 50m")

4.1.4 Exemple d’effet de bord à éviter

Dans cette section, nous illustrons les effets de bord produits dans le cas où on lisserait uniquement les ventes réalisées dans Paris intramuros.

Dans les deux cartes ci-dessous, le rayon de lissage (2 000 mètres) et les bornes de discrétisation sont communes pour assurer la comparabilité des deux cartes :

Dans le premier cas, on lisse les ventes situées dans Paris et ses alentours (comme précédemment).

Code

sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 2000)
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

creation_carte(sfCarrLiss_paris, var = "nbObsLisse", 
               titre = "Sans effets de bord, avec R = 2000m")

Dans le second cas, on lisse uniquement les ventes réalisées dans Paris intramuros pour observer les effets de bord (attention : les seuils des classes sont légérement différents).

Les effets de bord se matérialisent à la périphérie de Paris (en particulier au sud), où les densités locales sont impactées si on enlève des points.

4.1.5 Petit point sur les grilles automatiques de `btb`

La grille produite automatiquement par btb_smooth dépasse automatiquement les limites de la zone d’étude choisie.

Prenons l’exemple précédent du lissage (rayon 2000 m, pas 200 m) réalisé uniquement sur les ventes à l’intérieur de Paris intramuros.

La grille obtenue dépasse Paris intramuros et contient :

l’ensemble des carreaux contenant au moins une vente…
… mais aussi les 4 carreaux limitrophes de ces derniers (règle automatique¹)

Code

sfCarrLiss_extramuros <- btb_smooth(pts = sfBase_intramuros[, c("nbObsLisse", "x", "y")], 
                                   iCellSize = 200, 
                                   iBandwidth = 2000)

mf_map(paris_sf, lwd = 2, col = NA)
mf_map(sfCarrLiss_extramuros, col = NA, border = "red", add = T)

Il est possible de modifier la règle de création de cette grille automatique en rensignant le paramètre iNeighbor. Pour iNeighbor = 2, la grille du lissage contiendra 2 carreaux limitrophes aux carreaux contenant au moins une vente.

On peut aussi vouloir une grille uniquement constituée des carreaux contenant au moins une vente (iNeighbor = 0), mais attention à l’effet “gruyère” :

Code

sfCarrLiss_intramuros <- btb_smooth(pts = sfBase_intramuros[, c("nbObsLisse", "x", "y")], 
                                   iCellSize = 200, 
                                   iBandwidth = 2000,
                                   iNeighbor = 0)

mf_map(paris_sf, lwd = 2, col = NA)
mf_map(sfCarrLiss_intramuros, col = NA, border = "red", add = T)

De même, on peut souhaiter utiliser une grille de lissage particulière, adaptée à un territoire donné. Dans ce cas, il convient de la confectionner soi-même et de la renseigner dans le paramètre dfCentroids de la fonction btb_smooth. Par exemple, en faisant un lissage sur une ville cotière, on peut vouloir une grille qui ne déborde pas sur la mer tout en ne souhaitant pas de l’effet “gruyère” du paramétrage iNeighbor = 0.

4.1.6 Partir de données ponctuelles ou partir de données carroyées ?

Pour garantir la meilleure qualité, le mieux est de lisser en partant des données ponctuelles. Parfois, on ne dispose pas de ces données, ou bien leur traitement est trop long (souci de volumétrie par exemple). Il est possible de repartir de données agrégées au carreau pour lisser. En effet, chaque centroïde est une donnée ponctuelle en tant que telle.

Les points utilisés :

Code

# données déjà carroyées
points_carroyage <- btb_add_centroids(pts = dfBase_filtre, iCellSize = 200) %>% 
  group_by(x_centro, y_centro) %>% count(name = "nbVentes") %>% 
  st_as_sf(coords = c("x_centro", "y_centro"), crs = 2154) |> 
  st_join(paris_sf |> st_buffer(2000), left = F)

par(mfrow = c(1, 2)) # pour afficher les 2 cartes côte à côté

mf_map(x = paris_sf, col = NA, lwd = 4)
mf_map(x = sfBase_filtre, cex = 0.3, col = "red", add = TRUE)
mf_layout(title = "Données ponctuelles")

mf_map(x = paris_sf, col = NA, lwd = 4) 
mf_map(x = points_carroyage, cex = 0.3, col = "red", add = TRUE)
mf_layout(title = "Données carroyées")

Les données lissées :

Code

# lissage classique
sfBase_filtre$nbObsLisse <- 1
sfCarrLiss <- btb_smooth(pts = sfBase_filtre[, c("nbObsLisse", "x", "y")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 400)
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# lissage en partant des données déjà carroyées
sfCarrLiss_car <- btb_smooth(pts = points_carroyage[, c("nbVentes", "geometry")], 
                         sEPSG = 2154,
                         iCellSize = 200, 
                         iBandwidth = 400)
sfCarrLiss_car_paris <- sfCarrLiss_car %>% st_join(paris_sf, left = FALSE)

par(mfrow = c(1, 2))
creation_carte(data = sfCarrLiss_paris, var = "nbObsLisse", 
               titre = "Lissage des données ponctuelles")

creation_carte(data = sfCarrLiss_car, var = "nbVentes", 
               titre = "Lissage des données carroyées")

L’agrégation préalable par carreau réduit le nombre de points pour le lissage, d’où la perte d’information : le centroïde résume les ventes de son carreau, on perd la diversité des positions dans l’espace, et donc en précision. Mais c’est un gain de simplicité de calcul, en cas de données volumineuses.

Code

# comparaison de points 
nrow(points_carroyage) ; nrow(sfBase_filtre)

[1] 3235

[1] 22221

4.2 Calcul de moyenne

Dans cette section, nous allons nous intéresser au prix moyen des logements vendus à Paris en 2021.

RAPPEL : Une moyenne n’est pas le lissage du rapport (liss(variable / nbObsLisse)) mais le rapport des lissages : liss(variable) / liss(nbObs).

Ainsi, il convient de :

Lisser les prix de ventes des biens vendus
Lisser le nombre de biens vendus
Faire le ratio des deux précédents lissages pour chaque carreau de la grille produite.

Dans l’exemple ci-dessous, on prend un rayon de lissage de 800m et un pas de 100m :

Code

# Lissage des 2 variables
sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "valeur_fonciere", 
                                                "x", "y")] , 
                        sEPSG = 2154,
                        iCellSize = 100, 
                        iBandwidth = 800)
head(sfCarrLiss, 3)

Code

# Calculer le taux
sfCarrLiss$prixMoyen <- sfCarrLiss$valeur_fonciere / sfCarrLiss$nbObsLisse
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

creation_carte(data = sfCarrLiss_paris, var = "prixMoyen", 
               titre = "Lissage des prix de vente avec un rayon de 800m")

Les prix moyens sont particulièrement élevés dans le centre et l’Ouest de la capitale, ce qui semble conforme à l’intuition.

Attention cependant aux éventuelles valeurs atypiques (ventes extrêmement chères, erreurs dans les données) : le lissage en moyenne est potentiellement sensible à ces valeurs atypiques (contrairement au lissage quantile, voir partie 4.4).

4.3 Calcul de taux

Nous nous intéressons ici à la proportion de ventes portant sur des grands logements (comportant plus de 4 pièces principales).

Code

# Créer une indicatrice renseignant si le logement comporte plus de 4 pièces
dfBase_filtre <- dfBase_filtre %>% 
  mutate(quatrePieces = ifelse(nombre_pieces_principales >= 4, 1, 0))

# Lissage des 2 variables
sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("nbObsLisse", "quatrePieces",
                                                "x", "y")], 
                        sEPSG = 2154,
                        iCellSize = 100, 
                        iBandwidth = 800)

# Calcul du taux
sfCarrLiss$txQuatrePieces <- sfCarrLiss$quatrePieces / sfCarrLiss$nbObsLisse

# Intersection géographique avec Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# cartographie
creation_carte(data = sfCarrLiss_paris, var = "txQuatrePieces",
               titre = "Logements avec 4 pièces ou plus (rayon de 800m)")

Logiquement, la carte des prix moyens et la carte des grands logements ont des similitudes importantes.

On peut alors avoir envie de lisser le prix au mètre-carré.

Pour ce faire, n’oubliez pas qu’il convient de lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur) des logements vendus. Autrement, si on lisse le prix au m² de chaque logement vendu, on sur-pondère artificiellement les prix au m² des petits logements (car l’unité statistique devient alors le logement, et non le m² vendu).

Code

# lissage des 2 variables
sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("surface_reelle_bati", 
                                                "valeur_fonciere", "x", "y")],
                        sEPSG = 2154,
                        iCellSize = 100,
                        iBandwidth = 800)

# Calcul du prix au m²
sfCarrLiss$prixM2 <- sfCarrLiss$valeur_fonciere / sfCarrLiss$surface_reelle_bati

# Intersection géographique avec Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# Cartographie
creation_carte(data = sfCarrLiss_paris, var = "prixM2",
               titre = "Prix au m² (rayon de 800m)")

4.4 Calcul de quantiles géographiquement pondérés

Le lissage quantile (médiane, déciles…) permet d’avoir des indicateurs moins sensibles aux valeurs extrêmes et d’enrichir l’analyse avec des indicateurs de dispersion (comme par exemple l’écart interquantile).

Ce lissage se fait toujours avec la fonction btb_smooth mais il faut ajouter un paramètre, vQuantiles, qui contient la liste des quantiles souhaités.

Remarque : Le calcul de quantiles géographiquement pondérés présente certaines différences par rapport au lissage classique vu jusqu’à présent :

il n’est pas conservatif
la fonction crée :
- une colonne nbObs indiquant le nombre réel (non lissé) d’observations ayant contribué au calcul du carreau.
- pour chaque variable, autant de colonnes qu’il y a de quantiles souhaités

Nous allons calculer le 1er décile, la médiane et le 9ème décile du prix de vente des logements à Paris.

Code

# calcul des déciles "locaux"
sfCarrLiss <- btb_smooth(pts = dfBase_filtre[, c("valeur_fonciere", "x", "y")],
                        sEPSG = 2154,
                        iCellSize = 100,
                        iBandwidth = 1500,
                        vQuantiles = c(0.1, 0.5, 0.9))

# Visualiser les premières lignes des données lissées
head(sfCarrLiss, 3)

Simple feature collection with 3 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 651500 ymin: 6854700 xmax: 651800 ymax: 6854800
Projected CRS: RGF93 v1 / Lambert-93
  nbObs valeur_fonciere_01 valeur_fonciere_05 valeur_fonciere_09      x       y
1    57             182000             310000             610000 651550 6854750
2    53             182000             310000             620000 651650 6854750
3    52             182000             310000             620000 651750 6854750
                        geometry
1 POLYGON ((651500 6854800, 6...
2 POLYGON ((651600 6854800, 6...
3 POLYGON ((651700 6854800, 6...

Code

# Intersection géographique avec Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf, left = FALSE)

# Cartographie 1er décile
creation_carte(data = sfCarrLiss_paris, var = "valeur_fonciere_01", 
               titre = "Premier décile des prix de vente (rayon de 1500m)")

Code

# Cartographie médiane
creation_carte(data = sfCarrLiss_paris, var = "valeur_fonciere_05", 
               titre = "Médiane des prix de vente (rayon de 1500m)")

Un intérêt majeur du “lissage quantile” est d’offrir des indicateurs lissés de disparité des distributions. En effet, on peut par exemple calculer le rapport interdéciles lissé afin de représenter l’ampleur des écarts de prix dans les différents quartiers de Paris intramuros.

Code

sfCarrLiss_paris <- sfCarrLiss_paris %>% 
  mutate(rapp_interdecile = valeur_fonciere_09 / valeur_fonciere_01)

# Cartographie du rappot interdécile
creation_carte(data = sfCarrLiss_paris, var = "rapp_interdecile", 
               titre = "Rapport interdéciles des prix de vente (rayon de 1500m)")

(attention, il est ici question des prix de vente, et non des prix au m²).

Enfin, pour que le quantile lissé ait du sens, il faut qu’un nombre conséquent d’observations ait participé à sa construction. Ainsi, il est conseillé d’utiliser un rayon de lissage suffisamment élevé.

Code

summary(sfCarrLiss_paris$nbObs)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     13     725    1040    1106    1440    2917

Code

# Centroïdes dont le quantile lissé a été calculé avec moins de 100 observations
sfCarrLiss_paris %>% st_drop_geometry() %>% group_by(nbObs < 100) %>% count()

# A tibble: 2 × 2
# Groups:   nbObs < 100 [2]
  `nbObs < 100`     n
  <lgl>         <int>
1 FALSE         10200
2 TRUE            229

5 Indicateurs sur zonage à façon

Pour terminer, nous allons produire des indicateurs sur des zones à façon. En d’autres termes, nous importerons une couche cartographique du zonage qui nous intéresse et agrégerons les données sur ce zonage.

5.1 Indicateurs sur zones à façon à partir de données ponctuelles

Dans la continuité du précédent exercice, nous calculons le nombre et le prix moyen des transactions immobilières en 2021 dans le “triangle d’or²”, situé dans le 8e arrondissement de Paris. Ce quartier est connu pour ses immeubles cossus, ses hôtels et ses commerces de luxe.

Nous procédons de la manière suivante :

Import et visualisation de la zone à façon “Triangle d’or”. Ce fichier géographique est au format .kml et se lit avec la fonction st_read.

Code

chemin_file <- paste0(url_bucket, "triangle_or.kml")
triangleOr <- st_read(chemin_file)

Reading layer `Layer #0' from data source 
  `https://minio.lab.sspcloud.fr/projet-formation/r-lissage-spatial/triangle_or.kml' 
  using driver `KML'
Simple feature collection with 1 feature and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2.300717 ymin: 48.86447 xmax: 2.310474 ymax: 48.87203
Geodetic CRS:  WGS 84

Code

mapview(triangleOr, col.regions = "yellow")

Intersection géographique entre des points (les transactions) et un polygone (le quartier). Attention à bien s’assurer que les deux couches géographiques soient bien dans le même système de projection.

Code

# Système de projection de triangleOr
st_crs(triangleOr)$epsg ; st_crs(sfBase)$epsg

[1] 4326

[1] 2154

Code

# Transformation du système de projection de triangleOr 
triangleOr <- st_transform(triangleOr, crs = 2154)

# Intersection géographique avec les transactions
transac_triangle <- sfBase %>% st_join(triangleOr, left = F)

# Visualisation des points retenus
head(transac_triangle, 3)

Simple feature collection with 3 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 648903.5 ymin: 6863145 xmax: 649194 ymax: 6863601
Projected CRS: RGF93 v1 / Lambert-93
      id_mutation date_mutation  type_local nombre_pieces_principales
18629 2021-489256    2021-01-15 Appartement                         5
18676 2021-489336    2021-01-25 Appartement                         3
18677 2021-489337    2021-01-28 Appartement                         4
      valeur_fonciere surface_reelle_bati        x       y Name Description
18629         3300000                 180 648918.9 6863145                 
18676          590000                  53 648903.5 6863601                 
18677         1511345                  65 649194.0 6863377                 
                      geometry
18629 POINT (648918.9 6863145)
18676 POINT (648903.5 6863601)
18677   POINT (649194 6863377)

Code

# Cartographie des points dans le triangle
mapview(transac_triangle) + mapview(triangleOr, col.regions = "yellow")

Calcul du nombre et du prix moyen des transactions immobilières réalisées en 2021 dans ce quartier

Code

# Nombre de transactions en 2021 et prix moyen
nrow(transac_triangle) ; mean(transac_triangle$valeur_fonciere)

[1] 20

[1] 1637256

Remarque : On ne voit que 12 points sur la cartes, alors qu’on a trouvé 20 transactions dans le triangle : ceci s’explique car certaines transactions sont géolocalisées au même endroit (appartements situés dans un même immeuble par exemple).

5.2 Indicateurs sur zones à façon à partir de données carroyées

On cherche ici à connaître la proportion de logements sociaux dans un rayon de 1 km autour de la Direction générale de l’Insee.

Nous mobilisons pour cela les données carroyées à 200 mètres disponibles sur insee.fr.

Remarque : un tutoriel généralisant cet exercice est disponible en ligne, sur insee.fr, en accompagnement des données carroyées.

Le principe de la démarche consiste à :

Charger les données carroyées de Filosofi à 200m, uniquement en région parisienne, en ne conservant qu’un nombre réduit de variables, dont le nombre de logements sociaux.

Code

# fonction permettant de faire le chargement souhaité
st_read_maison <- function(chemin_tab){
  requete <- "SELECT idcar_200m, lcog_geo, i_est_200, ind, men, log_soc, geom
            FROM carreaux_200m_met
            WHERE SUBSTR(lcog_geo, 1, 2) IN ('75','92','93','94')"
  st_read(chemin_tab, query = requete)
}

# Chargement des données
chemin_file <- paste0(url_bucket, "carreaux_200m_met.gpkg")
carreaux <- st_read_maison(chemin_file)
carreaux <- carreaux %>% rename(Log_soc = log_soc, Men = men)
head(carreaux, 3)

Simple feature collection with 3 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 653706.4 ymin: 6862684 xmax: 654721.1 ymax: 6866380
Projected CRS: RGF93 v1 / Lambert-93
                       IdINSPIRE Depcom I_est_cr Men Log_soc
1 CRS3035RES200mN2893400E3763200  75119        0 990     937
2 CRS3035RES200mN2890000E3762200  75111        0 926      80
3 CRS3035RES200mN2893400E3762400  75119        0 508     323
                            geom
1 MULTIPOLYGON (((654521.9 68...
2 MULTIPOLYGON (((653843.4 68...
3 MULTIPOLYGON (((653725.1 68...

Créer un disque de 1 km autour du White

Code

# géolocalistion du White et buffer autour
white <- st_sfc(st_point(c(649218.36, 6857569.14)), crs = 2154)
buffer_white <- st_buffer(white, 1000) %>% st_as_sf()

Déterminer, pour ce disque, les carreaux qui le recouvrent

Remarque : L’ensemble de carreaux étant en général plus large que la zone, les agrégats obtenus seront des estimations des valeurs réelles, plus ou moins précises. Le problème ne se pose pas si vous travaillez directement sur des données disponibles au niveau “individu statistique” (avec des coordonnées x, y) comme en partie 5.1.

Code

carreaux_select <- carreaux %>% st_join(buffer_white, left = F)

Visualisation des carreaux issus de Filosofi carroyé qui intersectent notre zone d’intérêt :

Code

mapview(carreaux_select) +
  mapview(buffer_white, color = "#ffff00", lwd = 10, alpha.regions = 0, legend = F)

À noter que ces carreaux sont “penchés” : ils ont été produits en projection LAEA (standard européen), puis reprojeté en Lambert 93 (standard français).

Calculer l’indicateur sur la zone par somme des données carroyées.

Code

tx_logSoc_white <- sum(carreaux_select$Log_soc) / sum(carreaux_select$Men)
cat("On compte", round(tx_logSoc_white * 100),
    "% de logements sociaux dans un rayon de 1 km autour du White.")

On compte 33 % de logements sociaux dans un rayon de 1 km autour du White.

Vous êtes désormais capables de carroyer et lisser toutes les données que vous aurez à votre disposition !

Reproductibilité

Code

sessionInfo()

R version 4.6.1 (2026-06-24)
Platform: aarch64-apple-darwin23
Running under: macOS Tahoe 26.6.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: UTC
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] btb_0.2.2      leaflet_2.2.3  mapview_2.11.4 mapsf_1.2.2    dplyr_1.2.1   
[6] sf_1.1-2       devtools_2.5.2 usethis_3.2.1 

loaded via a namespace (and not attached):
 [1] xfun_0.60               raster_3.6-32           htmlwidgets_1.6.4      
 [4] lattice_0.22-9          leaflet.providers_3.0.0 vctrs_0.7.3            
 [7] tools_4.6.1             crosstalk_1.2.2         generics_0.1.4         
[10] stats4_4.6.1            tibble_3.3.1            proxy_0.4-29           
[13] pkgconfig_2.0.3         KernSmooth_2.23-26      satellite_1.0.6        
[16] RColorBrewer_1.1-3      uuid_1.2-2              RcppParallel_6.2.1     
[19] lifecycle_1.0.5         compiler_4.6.1          farver_2.1.2           
[22] textshaping_1.0.5       terra_1.9-51            codetools_0.2-20       
[25] maplegend_0.6.3         htmltools_0.5.9         class_7.3-23           
[28] yaml_2.3.12             jquerylib_0.1.4         pillar_1.11.1          
[31] ellipsis_0.3.3          classInt_0.4-11         cachem_1.1.0           
[34] wk_0.9.5                sessioninfo_1.2.4       brew_1.0-10            
[37] tidyselect_1.2.1        digest_0.6.39           purrr_1.2.2            
[40] fastmap_1.2.0           grid_4.6.1              cli_3.6.6              
[43] magrittr_2.0.5          base64enc_0.1-6         pkgbuild_1.4.8         
[46] leafem_0.2.5            e1071_1.7-17            withr_3.0.3            
[49] scales_1.4.0            sp_2.2-3                rmarkdown_2.32         
[52] otel_0.2.0              png_0.1-9               memoise_2.0.1          
[55] evaluate_1.0.5          knitr_1.51              s2_1.1.11              
[58] rlang_1.3.0             Rcpp_1.1.2              leafpop_0.1.0          
[61] glue_1.8.1              DBI_1.3.0               pkgload_1.5.3          
[64] svglite_2.2.2           jsonlite_2.0.0          R6_2.6.1               
[67] systemfonts_1.3.2       fs_2.1.0                units_1.0-1

6 Autre logiciels

Grâce à différentes extensions, le package btb peut être utilisé avec d’autres logiciels que R :

Python : libraryBTBpy
Qgis grâce à un plug-in développé par Lionel Cacheux (DR Insee Grand Est).

Footnotes

cf. séquence “théorie du lissage” de la formation “Comment utiliser les outils de l’analyse urbaine ?”↩︎
Pour information, ce triangle a été créé manuellement sur le site du Géoportail.↩︎

1 Introduction

1.1 Objectifs du TP

1.2 Avertissements

1.2.1 Secret statistique

1.2.2 Qualité des cartes

1.2.3 Système de projection

2 Configurations

2.1 Chargement des librairies

2.2 Chargement de la base

2.3 Chargement des données cartographiques

2.3.1 Fond de carte du territoire à étudier

2.3.2 Territoires englobants, sélection des données à lisser

2.3.2.1 Première méthode : sélection géométrique des ventes

2.3.2.2 Seconde méthode : sélection non-géométrique des ventes

3 Carroyage de données

4 Lissage

4.1 Calcul de densité

4.1.1 Grille automatique ; carreaux 200m ; rayon de lissage 400m

4.1.2 Grille automatique ; carreaux 200m ; rayons de lissage 600m et 1 km

4.1.3 Grille automatique ; carreaux 50m ; rayon de lissage 1 km

4.1.4 Exemple d’effet de bord à éviter

4.1.5 Petit point sur les grilles automatiques de btb

4.1.6 Partir de données ponctuelles ou partir de données carroyées ?

4.2 Calcul de moyenne

4.3 Calcul de taux

4.4 Calcul de quantiles géographiquement pondérés

5 Indicateurs sur zonage à façon

5.1 Indicateurs sur zones à façon à partir de données ponctuelles

5.2 Indicateurs sur zones à façon à partir de données carroyées

6 Autre logiciels

Footnotes

4.1.5 Petit point sur les grilles automatiques de `btb`