Formation au carroyage et lissage spatial sur R

TUTORIEL

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 avec R.

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

À carroyer les informations.
À réaliser des lissages de densité, des lissages de moyennes, des lissages de taux et des lissages 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 avec l’appui de Kim Antunez

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 très bien 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

Pour information, 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 (importer des .shp, transformer des projections, et réaliser des géotraitements)
mapsf pour réaliser des cartes dans RStudio
mapview (reposant sur leaflet) pour réaliser des cartes interactives (fond de carte OpenStreetMap)
btb pour le carroyage et lissage
dplyr pour le traitement des données, en particulier l’agrégation géographique.

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

# Charger les librairies nécessaires

## Liste des librairies utilisées
packages <-  c("devtools", "dplyr", "sf", "mapsf", "leaflet", "mapview", "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)
    }
  }
)

Pour les personnes utilisant l’environnement de développement du SSPCloud, il est possible d’utiliser le package aws.s3 pour charger les données stockées dans le cloud Minio de l’Insee. Ce package permet un chargement des données plus rapide qu’un chargement classique.

# Chunk permettant d'utiliser les données S3. Trouver une autre solution à termes

if (!require("aws.s3", character.only = TRUE)) {
      install.packages("aws.s3",repos = "https://cloud.R-project.org",
                       dependencies = TRUE, quiet = TRUE)
      library("aws.s3", character.only = TRUE)
    }

2.2 Chargement de la base

Pour ce TP, nous utilisons la base « Demandes de Valeurs Foncières », ou DVF, qui recense l’ensemble des ventes de biens fonciers réalisées au cours des cinq dernières années, en métropole et dans les départements et territoires d’outre-mer — sauf à Mayotte et en 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. 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.

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)

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

Remarque importante : Cette 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.

Le code ci-dessous permet d’importer la première base de données ventesImmo_couronneParis.RDS utilisée dans ce tutoriel. Elle est stockée sous Minio, dans le “bucket public” : s3/projet-formation/r-lissage-spatial/.

# Charger la source de données (variable nommée dfBase)
url_bucket <- "https://minio.lab.sspcloud.fr/"
bucket <- "projet-formation"
object <- "r-lissage-spatial/ventesImmo_couronneParis.RDS"

Nous téléchargeons cette base grâce à son URL public, puis on la charge dans R.

# Charger la source de données (variable nommée dfBase) avec son URL
url_file <- url(paste0(url_bucket, bucket, "/", object))
dfBase <- readRDS(url_file)

L’équivalent en utilisant la library aws.s3 (uniquement depuis une session SSPCloud) :

dfBase <-
  aws.s3::s3read_using(
    FUN = base::readRDS,
    object = object,
    bucket = bucket
    ,
    opts = list("region" = "")
  )

On peut ensuite manipuler notre base de données chargée en mémoire.

# Visualiser les premières lignes de la base
head(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
## 4 2021-447027    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
## 1          480000                  64 647357.3 6868635
## 2          345000                  43 644483.5 6867695
## 3          384000                  41 648001.8 6866153
## 4          261900                  44 647414.8 6868937
## 5          407200                  24 646929.9 6864730
## 6         1040000                  90 645132.7 6864646

dim(dfBase)

## [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 donc le contour du territoire étudié avec la fonction st_read du package sf.

# Charger le fond de carte du territoire étudié
chemin_file <- paste0(url_bucket, bucket, "/r-lissage-spatial/depCouronne.gpkg")
depCouronne_sf <- 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

# Visualisons cette 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...

plot(depCouronne_sf$geom)

# On renomme la variable "geom" en "geometry"
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.

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

# Visualisons cette 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...

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.

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

## Coordinate Reference System:
##   User input: RGF93 v1 / Lambert-93 
##   wkt:
## PROJCRS["RGF93 v1 / Lambert-93",
##     BASEGEOGCRS["RGF93 v1",
##         DATUM["Reseau Geodesique Francais 1993 v1",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4171]],
##     CONVERSION["Lambert-93",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",46.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",44,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",700000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica)."],
##         BBOX[41.15,-9.86,51.56,10.38]],
##     ID["EPSG",2154]]

st_crs(paris_sf)

## Coordinate Reference System:
##   User input: RGF93 v1 / Lambert-93 
##   wkt:
## PROJCRS["RGF93 v1 / Lambert-93",
##     BASEGEOGCRS["RGF93 v1",
##         DATUM["Reseau Geodesique Francais 1993 v1",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4171]],
##     CONVERSION["Lambert-93",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",46.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",44,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",700000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica)."],
##         BBOX[41.15,-9.86,51.56,10.38]],
##     ID["EPSG",2154]]

# Transformer la projection en Lambert 93 (epsg 2154)
depCouronne_sf <- sf::st_transform(depCouronne_sf, crs = 2154)
paris_sf <- 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 également sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros). Autrement dit, 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 sont situées à proximité de zones artificiellement vides en ventes immobilières dans notre nuage de points, à savoir dans les communes limitrophes. Pourtant, dans la réalité, il y a bien des ventes réalisées à Malakoff, Vanves, Montrouge, etc. Et ces ventes influencent la densité lissée aux niveau de la Porte de Vanves.

2.3.2.1 Première 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 intramuros.

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 le rectangle. Cette méthode est très efficace d’un point de vue calculatoire, mais requiert de construire au préalable un rectangle adapté.

# Mise en forme de la couche "territoire" : sélection des variables et renommage
paris_sf$nom <- "territoire"
paris_sf <- paris_sf[, c("nom", "geometry")]

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

##    xmin    ymin    xmax    ymax 
##  643076 6857499  660897 6867034

# 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

Ci-dessous le code permettant de créer un petit schéma explicatif des zones sélectionnées.

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

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

# Options globales pour les cartes avec mapview     
mapviewOptions(
  basemaps = c(
    "OpenStreetMap",
    "CartoDB.Positron",
    "CartoDB.DarkMatter",
    "Esri.WorldImagery"
  )
)
# Cartographie pédagogique avec mapview
mapview(paris_sf, col.regions = "#26cce7")+
  mapview(bufferBbox_sf %>% st_cast("MULTILINESTRING"),
          color = "#FFC300", lwd = 6) +
  mapview(bbox_sf %>% st_cast("MULTILINESTRING"),
          color = "#229954", lwd = 6)

Ce schéma présente :

en bleu, le territoire d’intérêt (Paris)
en vert, 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.

Pour sélectionner les données individuelles d’intérêt pour le lissage, on peut tout d’abord filtrer les observations dont les coordonnées sont comprises à l’intérieur du rectangle jaune. 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 (cf. infra).

# Repérer les coordonnées extrêmes, avec une marge (ici 2km)
bufferBbox

##    xmin    ymin    xmax    ymax 
##  641076 6855499  662897 6869034

xMin <- bufferBbox["xmin"]
xMax <- bufferBbox["xmax"]
yMin <- bufferBbox["ymin"]
yMax <- bufferBbox["ymax"]

# Ne garder que les données dans le rectangle englobant, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= xMin &
                          dfBase$x <= xMax &
                          dfBase$y >= yMin &
                          dfBase$y <= yMax, ]
dim(dfBase_filtre)

## [1] 24419     8

dim(dfBase)

## [1] 34489     8

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

La seconde méthode est géométrique ou vectorielle. Elle consiste à utiliser nos données individuelles 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 ;

sfBase <- sf::st_as_sf(dfBase,
                       coords = c("x", "y"),
                       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 ;

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. Une fois le lissage réalisé, on peut réduire la carte à la zone étudiée.

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

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

## [1] 22221

head(sfBase_filtre)

## Simple feature collection with 6 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 643891.9 ymin: 6864587 xmax: 648239 ymax: 6866235
## 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
## 9  2021-447034    2021-01-18 Appartement                         4
## 10 2021-447035    2021-01-18 Appartement                         4
## 11 2021-447037    2021-01-07 Appartement                         4
##    valeur_fonciere surface_reelle_bati        nom                 geometry
## 3           384000                  41 territoire POINT (648001.8 6866153)
## 5           407200                  24 territoire POINT (646929.9 6864730)
## 6          1040000                  90 territoire POINT (645132.7 6864646)
## 9           625000                  68 territoire POINT (647494.6 6866019)
## 10          721250                  94 territoire POINT (643891.9 6864587)
## 11         1035000                  94 territoire   POINT (648239 6866235)

Cette méthode, bien qu’élégante, peut-être lourde d’un point de vue calculatoire. Avec des données volumineuses, il convient au minimum de faire un premier filtrage avec la précédente méthode.

Ci-dessous le code permettant d’avoir un aperçu de la zone, du buffer et de 2000 points tirés aléatoirement dedans.

# Mise en forme de la couche buffer
buffer_sf$nom <- "buffer"

# Échantillon de 2000 observations dans le buffer
sfBase_sample <- sfBase_filtre[sample(1:nrow(sfBase_filtre), 2000), ]

# Cartographie pédagogique avec mapview
mapview(paris_sf, col.regions = "#26cce7")+
  mapview(buffer_sf %>% st_cast("MULTILINESTRING"),
          color = "#FFC300", lwd = 6) +
  mapview(sfBase_sample, 
          alpha = 0, cex = 2)

3 Carroyage de données

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

Associer chaque point (= vente géolocalisée) au centroïde du carreau auquel il appartient (le territoire est découpé en carreaux de 200 mètres à partir de l’origine du référentiel). La fonction btb_add_centroids est faite pour cela.

iCellSize = 200 # Carreaux de 200m
points_carroyage <- btb::btb_add_centroids(pts = dfBase_filtre,
                                           iCellSize = iCellSize) 
head(points_carroyage)

##   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
## 4 2021-447027    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 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
## 4          261900                  44 647414.8 6868937   647500  6868900
## 5          407200                  24 646929.9 6864730   646900  6864700
## 6         1040000                  90 645132.7 6864646   645100  6864700

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

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

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

pts : un tableau avec les colonnes x et y 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).

carreaux <- btb::btb_ptsToGrid(pts = points_carroyage,
                               sEPSG = 2154, iCellSize = iCellSize)

Se restreindre au champ des carreaux intersectant Paris

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

On obtient alors ce carroyage des ventes dans Paris intramuros :

# Cartographie
contourParis <- st_cast(paris_sf[, c("geometry")], "MULTILINESTRING")

mf_init(x=carreaux)
mf_map(x = carreaux,
       type = "choro",
       var = "nbVentes",
       breaks = "quantile",
       nbreaks = 5,
       lwd = 1,
       leg_val_rnd = 1,
       add = TRUE)
mf_map(x = contourParis,
       lwd = 4,
       col ="black",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 à la réalité.

Il est é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.

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 allons préalablement créer une variable nbObsLisse qui vaudra 1 pour chaque observation (chaque logement vendu en 2021) et qui sera lissée avec la fonction btb_smooth du package btb.

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

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

##   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
## 4 2021-447027    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 nbObsLisse
## 1          480000                  64 647357.3 6868635          1
## 2          345000                  43 644483.5 6867695          1
## 3          384000                  41 648001.8 6866153          1
## 4          261900                  44 647414.8 6868937          1
## 5          407200                  24 646929.9 6864730          1
## 6         1040000                  90 645132.7 6864646          1

La fonction de lissage btb::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 quatre paramètres :

pts : le tableau des données à lisser. Il doit nécessairement contenir une colonne x, une colonne y, 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.

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

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

## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 646000 ymin: 6855400 xmax: 647800 ymax: 6855600
## Projected CRS: RGF93 v1 / Lambert-93
##        x       y nbObsLisse                       geometry
## 1 646100 6855500   1.712851 POLYGON ((646000 6855600, 6...
## 2 646500 6855500   1.261316 POLYGON ((646400 6855600, 6...
## 3 647100 6855500   2.301557 POLYGON ((647000 6855600, 6...
## 4 647300 6855500   2.808277 POLYGON ((647200 6855600, 6...
## 5 647500 6855500   2.395913 POLYGON ((647400 6855600, 6...
## 6 647700 6855500   1.949819 POLYGON ((647600 6855600, 6...

# Nombre de lignes lissées
nrow(sfCarrLiss)

## [1] 4107

Attention, la fonction retourne une erreur en cas de :

présence d’une variable non-numérique ;
valeur(s) absente(s) dans les colonnes x ou y.

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).

mf_init(x = paris_sf)
mf_map(x = st_cast(sfCarrLiss[, c("geometry")], "LINESTRING"), 
       lwd = 1, add = T)
mf_map(x = contourParis, 
       lwd = 8,
       col = "wheat", add = TRUE)
mf_layout(title = "btb_smooth génère une grille carroyée",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Une propriété particulièrement importante de la fonction de lissage classique du package btb 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.

# Nombre de ventes dans la petite couronne
sum(dfBase_filtre$nbObsLisse)

## [1] 24419

# Après lissage, nombre lissé de ventes dans les carreaux
sum(sfCarrLiss$nbObsLisse)

## [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 réaliser le filtre ci-dessous :

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

Voyons 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 sf::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 : ici, par quantiles par exemple ;
- Avec 5 classes (c’est-à-dire 5 couleurs).

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

# Carte lissée
mf_init(x = sfCarrLiss_paris)
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd = 1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "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 1km

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 1km, au lieu de 400m.

À 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 bien-entendu possible et conseillé de la supprimer sur les cartes publiées comme dans cet exemple :

mf_init(x = sfCarrLiss_paris)
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       border = NA, # C'est ici que ça se passe
       # lwd =1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m, sans grille",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

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

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

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).

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

Dans le second cas, on lisse uniquement les ventes réalisées dans Paris intramuros pour observer les effets de bord.

# Lissage sans gestion des effets de bord (intramuros uniquement)
sfBase_intramuros <- sfBase %>% st_join(paris_sf, left = FALSE) %>% 
  mutate(nbObsLisse = 1L) %>% 
  select(nbObsLisse, geometry)

sfCarLiss_intramuros <- btb::btb_smooth(pts = sfBase_intramuros, 
                                    iCellSize = 200, 
                                    iBandwidth = 2000)

sfCarLiss_intramuros_paris <- sfCarLiss_intramuros %>%
  st_join(paris_sf, left = FALSE)

Les effets de bord se matérialisent à la périphérie de Paris, où les densités lissées sont artificiellement plus faibles que quand on lisse en prenant une marge.

4.1.5 Petit point sur les grilles automatiques de `btb`

La grille produite automatiquement par btb::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¹)

mapview(sfCarLiss_intramuros %>%
          st_cast("MULTILINESTRING"), lwd = 2,
        layer.name = "Grille automatique") +
  mapview(contourParis, color = "red", lwd = 10)

Il est possible de modifier la règle de création de cette grille automatique en rensignant le paramètre iNeighbor. Par exemple, si iNeighbor = 2, la grille du lissage contient 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 » :

De même, on peut souhaiter utiliser une grille de lissage particulière et 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.2 Calcul de moyenne

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

Remarque : 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 (exactement comme précédemment)
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 :

# Ajout de la variable "valeur_fonciere"
dataLissage <- dfBase_filtre[, c("nbObsLisse", "valeur_fonciere", "x", "y")] 

# Lissage
sfCarLiss <- btb::btb_smooth(pts = dataLissage, 
                                    sEPSG = 2154,
                                    iCellSize = 100, 
                                    iBandwidth = 800)
sfCarLiss

# Calculer le taux
sfCarLiss$prixMoyen <- sfCarLiss$valeur_fonciere / sfCarLiss$nbObsLisse

# Cartographie
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf, left = FALSE)
mf_init(x=sfCarLiss_paris)
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var = "prixMoyen",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Lissage des prix de vente avec un rayon de 800m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

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).

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

# Lissage
dataLissage <- dfBase_filtre[, c("nbObsLisse", "quatrePieces", "x", "y")]
sfCarLiss <- btb::btb_smooth(pts = dataLissage, 
                                    sEPSG = 2154,
                                    iCellSize = 100, 
                                    iBandwidth = 800)

# Calculer le taux
sfCarLiss$txQuatrePieces <- sfCarLiss$quatrePieces / sfCarLiss$nbObsLisse

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

# Cartographie
mf_init(x = sfCarLiss_paris)
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var = "txQuatrePieces",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Logements avec 4 pièces ou plus (rayon de 800m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Fort 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).

dataLissage <- dfBase_filtre[, c("surface_reelle_bati",
                            "valeur_fonciere", "x", "y")]

sfCarLiss <- btb::btb_smooth(pts = dataLissage,
                                    sEPSG = 2154,
                                    iCellSize = 100,
                                    iBandwidth = 800)

# Calculer le prix au m²
sfCarLiss$prixM2 <- sfCarLiss$valeur_fonciere / sfCarLiss$surface_reelle_bati

# Cartographie
sfCarLiss_paris <- sfCarLiss[unlist(st_intersects(paris_sf, sfCarLiss)),]
mf_init(x = sfCarLiss_paris)
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var = "prixM2",
       breaks = "quantile",
       nbreaks = 5,
       border = NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = st_cast(paris_sf[, c("geometry")], "MULTILINESTRING"),
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Prix au m² (rayon de 800m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

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.

dataLissage <- dfBase_filtre[, c("valeur_fonciere", "x", "y")]

sfCarLiss <- btb::btb_smooth(pts = dataLissage,
                                    sEPSG = 2154,
                                    iCellSize = 100,
                                    iBandwidth = 1500,
                                    vQuantiles = c(0.1, 0.5, 0.9))

## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that

# Visualiser les premières lignes des données lissées
head(sfCarLiss)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651500 ymin: 6854700 xmax: 652100 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
## 4    55             182000             321000             620000 651850 6854750
## 5    48             182000             324000             620000 651950 6854750
## 6    55             182000             324000             620000 652050 6854750
##                         geometry
## 1 POLYGON ((651500 6854800, 6...
## 2 POLYGON ((651600 6854800, 6...
## 3 POLYGON ((651700 6854800, 6...
## 4 POLYGON ((651800 6854800, 6...
## 5 POLYGON ((651900 6854800, 6...
## 6 POLYGON ((652000 6854800, 6...

# Cartographie
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf, left = FALSE)

mf_init(x = sfCarLiss_paris)
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var = "valeur_fonciere_01",
       breaks = "quantile",
       nbreaks = 5,
       border = NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Premier décile des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

# Cartographie
mf_init(x = sfCarLiss_paris)
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var = "valeur_fonciere_05",
       breaks = "quantile",
       nbreaks = 5,
       border = NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Médiane des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

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.

sfCarLiss_paris$rapp_interdecile <- sfCarLiss_paris$valeur_fonciere_09 / sfCarLiss_paris$valeur_fonciere_01

# Cartographie
mf_init(x = sfCarLiss_paris)
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var = "rapp_interdecile",
       breaks = "quantile",
       nbreaks = 5,
       border = NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Rapport interdéciles des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

(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é.

summary(sfCarLiss_paris$nbObs)

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

# Part de centroïdes dont le quantile lissé 
# a été calculé avec moins de 100 observations
sfCarLiss_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 sans problème avec la fonction sf::st_read.

chemin_file <- paste0(url_bucket, bucket, "/r-lissage-spatial/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

mapview(triangleOr, col.regions = "#ffff00")

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.

# Système de projection de triangleOr
st_crs(triangleOr)[["input"]]

## [1] "WGS 84"

# 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[, "geometry"], left = F)

# Visualisation de points retenus
transac_triangle

## Simple feature collection with 20 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 648767.5 ymin: 6863055 xmax: 649194 ymax: 6863601
## Projected CRS: RGF93 v1 / Lambert-93
## First 10 features:
##       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
## 18698 2021-489369    2021-02-05 Appartement                         6
## 18723 2021-489406    2021-01-20 Appartement                         3
## 18725 2021-489410    2021-02-12 Appartement                         3
## 18748 2021-489452    2021-02-24 Appartement                         2
## 18757 2021-489464    2021-02-04 Appartement                         2
## 18768 2021-489483    2021-02-19 Appartement                         2
## 18777 2021-489495    2021-02-19 Appartement                         4
##       valeur_fonciere surface_reelle_bati                 geometry
## 18629         3300000                 180 POINT (648918.9 6863145)
## 18676          590000                  53 POINT (648903.5 6863601)
## 18677         1511345                  65   POINT (649194 6863377)
## 18698         2450000                 133 POINT (648814.5 6863055)
## 18723          997000                  72 POINT (648870.5 6863135)
## 18725         1605000                 100 POINT (648767.5 6863434)
## 18748          821500                  60 POINT (648771.3 6863503)
## 18757         1200000                  74 POINT (648767.5 6863434)
## 18768           74900                  44 POINT (648814.2 6863144)
## 18777         2500000                 100 POINT (648970.8 6863181)

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

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

# Nombre de transactions en 2021 : 
nrow(transac_triangle)

## [1] 20

# Prix moyen
mean(transac_triangle$valeur_fonciere)

## [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 1km 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 : une 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 2019 à 200m, uniquement en région parisienne, en ne conservant que 6 variables dont le nombre de logements sociaux.

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

# Chargement des données
chemin_file <- paste0(url_bucket, bucket, 
                    "/r-lissage-spatial/carreaux_200m_met.gpkg")
carreaux <-  st_read_maison(chemin_file)
head(carreaux)

## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 666880.8 ymin: 6843259 xmax: 667515.6 ymax: 6843714
## Projected CRS: RGF93 v1 / Lambert-93
##                       idcar_200m   lcog_geo i_est_200  men log_soc
## 1 CRS3035RES200mN2869400E3773400 9405691097         1  1.1       0
## 2 CRS3035RES200mN2869400E3773600      94056         0 40.0       0
## 3 CRS3035RES200mN2869400E3773800      94056         0 13.0       0
## 4 CRS3035RES200mN2869600E3773400 9405691097         0 12.0       0
## 5 CRS3035RES200mN2869600E3773600      94056         0 32.0       1
## 6 CRS3035RES200mN2869600E3773800      94056         0 11.0       0
##                             geom
## 1 POLYGON ((666918.1 6843259,...
## 2 POLYGON ((667117.2 6843278,...
## 3 POLYGON ((667316.4 6843297,...
## 4 POLYGON ((666899.5 6843458,...
## 5 POLYGON ((667098.6 6843477,...
## 6 POLYGON ((667297.8 6843496,...

Créer un disque de 1km autour du White

white <- st_sfc(st_point(c(649218.36, 6857569.14)), crs = 2154)
buffer_white <- st_buffer(white, 1000) %>% st_as_sf()

Déterminer, pour cette zone, l’ensemble des carreaux qui la 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.

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

La carte ci-dessous représente les carreaux issus de Filosofi carroyé en 2019 qui intersectent notre zone d’intérêt. À noter que ces carreaux sont « penchés » : ceci est la résultante de leur reprojection en Lambert 93. En effet, ces carreaux ont été produits en projection LAEA (standard européen), ce qui leur donne cet aspect une fois reprojetés dans le standard français.

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

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

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 1km autour du white")

## On compte  33 % de logements sociaux dans un rayon de 1km autour du white

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

Reproductibilité

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.7.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## 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] aws.s3_0.3.21  btb_0.2.0      mapview_2.11.2 leaflet_2.2.1  mapsf_0.9.0   
## [6] sf_1.0-15      dplyr_1.1.4    devtools_2.4.5 usethis_2.2.3 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0        fastmap_1.1.1           promises_1.2.1         
##  [4] digest_0.6.34           mime_0.12               lifecycle_1.0.4        
##  [7] ellipsis_0.3.2          terra_1.7-73            magrittr_2.0.3         
## [10] compiler_4.3.2          rlang_1.1.3             sass_0.4.8             
## [13] tools_4.3.2             leafpop_0.1.0           utf8_1.2.4             
## [16] yaml_2.3.8              knitr_1.45              brew_1.0-10            
## [19] htmlwidgets_1.6.4       sp_2.1-3                pkgbuild_1.4.3         
## [22] classInt_0.4-10         curl_5.2.0              aws.signature_0.6.0    
## [25] RColorBrewer_1.1-3      xml2_1.3.6              pkgload_1.3.4          
## [28] KernSmooth_2.23-22      miniUI_0.1.1.1          withr_3.0.0            
## [31] purrr_1.0.2             grid_4.3.2              stats4_4.3.2           
## [34] fansi_1.0.6             urlchecker_1.0.1        profvis_0.3.8          
## [37] xtable_1.8-4            e1071_1.7-14            leafem_0.2.3           
## [40] colorspace_2.1-0        scales_1.3.0            cli_3.6.2              
## [43] rmarkdown_2.25          generics_0.1.3          remotes_2.4.2.1        
## [46] RcppParallel_5.1.7      httr_1.4.7              sessioninfo_1.2.2      
## [49] DBI_1.2.2               cachem_1.0.8            proxy_0.4-27           
## [52] stringr_1.5.1           s2_1.1.6                base64enc_0.1-3        
## [55] rmdformats_1.0.4        vctrs_0.6.5             jsonlite_1.8.8         
## [58] bookdown_0.37           systemfonts_1.0.5       crosstalk_1.2.1        
## [61] maplegend_0.1.0         jquerylib_0.1.4         units_0.8-5            
## [64] glue_1.7.0              leaflet.providers_2.0.0 codetools_0.2-19       
## [67] stringi_1.8.3           later_1.3.2             raster_3.6-26          
## [70] munsell_0.5.0           tibble_3.2.1            pillar_1.9.0           
## [73] htmltools_0.5.7         satellite_1.0.5         R6_2.5.1               
## [76] wk_0.9.1                evaluate_0.23           shiny_1.8.0            
## [79] lattice_0.21-9          highr_0.10              png_0.1-8              
## [82] memoise_2.0.1           httpuv_1.6.14           bslib_0.6.1            
## [85] class_7.3-22            uuid_1.2-0              Rcpp_1.0.12            
## [88] svglite_2.1.3           xfun_0.42               fs_1.6.3               
## [91] pkgconfig_2.0.3

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).

cf. séquence 5 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.↩︎

Formation au carroyage et lissage spatial sur R TUTORIEL