Formation au carroyage et lissage spatial sur R

Kim Antunez - Solène Colin

2024-03-01

Plan :

• 1. Introduction
  • 1.1 Objectifs du TP
  • 1.2 Avertissements
• 2. Configurations
  • 2.1 Chargement des librairies
  • 2.2 Chargement de la base
  • 2.3 Chargement et préparation des données cartographiques
  • 2.4 Carroyage de données
• 3. Lissage
  • 3.1 Calcul de densité
  • 3.2 Calcul de moyenne
  • 3.3 Calcul de taux
• 4. Indicateurs sur zonage à façon
• 5. Conclusion

1. Introduction

1.1 Objectifs du TP

• En 2018, création du package btb (PSAR analyse urbaine, Arlindo Dos Santo 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é,
  • de moyennes,
  • des lissages de taux,

• À calculer un indicateur sur une zone à façon

Liens utiles

• Code de la formation
• Site web des supports de formation

1.2 Avertissements

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

Qualité des cartes

Pour simplifier : on prend des libertés avec la sémiologie cartographique

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

Système de projection

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

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

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("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)
    }
  }
)

Sur le SSPCloud :

• Utiliser le package aws.s3 pour charger les données stockées dans Minio ;
• Plus rapide qu’un chargement classique.

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

Base « Demandes de Valeurs Foncières »,

• Produite par la Direction générale des finances publiques (actes notariés).

Elle recense :

• l’ensemble des ventes de biens fonciers,
• Au cours des 5 dernières années,
• Hors Mayotte et en Alsace-Moselle.
• Biens bâtis ou terrains

Constitution de la base à partir de cette source

• Uniquement sur le périmètre de la petite couronne parisienne
• Pour l’année 2021.

8 variables utilisées

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

Remarques

• Code de production de la base.

• Base non représentative de la réalité : usage purement pédagogique.

En dehors du SSPCloud

Importation de la base ventesImmo_couronneParis.RDS, stockée sous Minio

Téléchargement via URL :

# 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"

url_file <- url(paste0(url_bucket, bucket, "/", object))
dfBase <- readRDS(url_file)

Dans le SSPcloud

dans le “bucket public” : s3/projet-formation/r-lissage-spatial/.

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

Manipulation de notre base chargée en mémoire

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 et préparation des données cartographiques

Fond de carte du territoire à étudier [Chargement]

• Chargement de la couche vectorielle des départements de la petite couronne parisienne.
• Type de fichiers : .shp ou .gpkg (format recommmandé, voir ici)

# Récupération du fond de carte grâce à sf::st_read
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

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

Fond de carte du territoire à étudier [Visualisation]

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)

Contours de Paris

# Sélection de Paris
paris_sf <- depCouronne_sf[depCouronne_sf$code == "75", ]
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)

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

Projections

# 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]]

• Si besoin, on reprojette les deux bases dans le même système de projection.
• Ici, en Lambert 93 (epsg 2154)

depCouronne_sf <- sf::st_transform(depCouronne_sf, crs = 2154)
paris_sf <- sf::st_transform(paris_sf, crs = 2154)

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.

Territoires englobants

☠️ Eviter les effets de bord

➡️ Toujours électionner des données à lisser au-delà de la 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.

Méthode 1 : sélection non-géométrique [1/4]

• filtrer les x et y compris dans un grand rectangle englobant Paris intramuros

Avantages / Inconvénients

• ✔️ Très efficace computationnellement

• ❌ Mais requiert de construire au préalable un rectangle adapté…

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

Méthode 1 : sélection non-géométrique [2/4]

# 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 

• … puis buffer autour de celle-ci (marge : 2000m)

marge <- 2000
bufferBbox <- bbox
bufferBbox[["xmin"]] <- bufferBbox[["xmin"]] - marge
bufferBbox[["xmax"]] <- bufferBbox[["xmax"]] + marge
bufferBbox[["ymin"]] <- bufferBbox[["ymin"]] - marge
bufferBbox[["ymax"]] <- bufferBbox[["ymax"]] + marge

Méthode 1 : sélection non-géométrique [3/4]

# 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
m <- 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)

m

• 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 = 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.

Filtrage des coordonnées sont comprises à l’intérieur du rectangle jaune.

• très efficace computationnellement (filtre de valeurs numériques)
• Pas de traitement vectorielle (plus chronophages).

Méthode 1 : sélection non-géométrique [4/4]

Filtrer (via les x,y) les logements du rectangle jaune

# 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

Méthode 2 : sélection géométrique [1/4]

• Utiliser nos données individuelles comme un ensemble de points géolocalisés
• et procéder à des intersections géographiques.

Avantages / Inconvénients

• ✔️ plus court à coder et plus logique

• ❌ Potentiellement lourd d’un point de vue calculatoire.

➡️ Avec des données volumineuses, faire au minimum un premier filtrage non géométrique.

Méthode 2 : sélection géométrique [2/4]

1. On transforme nos observations en points vectoriels ;

sfBase <- sf::st_as_sf(dfBase, coords = c("x", "y"), crs = 2154)

2. 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, prendre une marge légèrement plus grande que le rayon de lissage envisagé.

Méthode 2 : sélection géométrique [3/4]

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

sfBase_filtre <- st_join(sfBase, buffer_sf, left = FALSE)
head(sfBase_filtre, 2)

Simple feature collection with 2 features and 7 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 646929.9 ymin: 6864730 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
  valeur_fonciere surface_reelle_bati        nom                 geometry
3          384000                  41 territoire POINT (648001.8 6866153)
5          407200                  24 territoire POINT (646929.9 6864730)

Méthode 2 : sélection géométrique [4/4]

Zone, buffer et 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, #col.regions = "black",alpha.regions=0.5,
          alpha = 0, cex = 2)

2.4 Carroyage de données

Les étapes du carroyage [1/7]

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

➡️ Le carroyage nécessite plusieurs étapes

Les étapes du carroyage [2/7]

1. Associer chaque point (= vente géolocalisée) au centroïde du carreau auquel il appartient (btb_add_centroids).

➡️ Le territoire est découpé en carreaux de 200 mètres à partir de l’origine du référentiel.

iCellSize = 200 # Carreaux de 200m
# On repart de la base filtrée selon la première méthode
points_carroyage <- btb::btb_add_centroids(pts = dfBase_filtre,
                                        iCellSize = iCellSize) 

head(points_carroyage,1)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y x_centro y_centro
1          480000                  64 647357.3 6868635   647300  6868700

Les étapes du carroyage [3/7]

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

points_centroides <- points_carroyage %>%
  group_by(x_centro, y_centro) %>% 
  count(name = "nbVentes")
head(points_centroides, 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

Les étapes du carroyage [4/7]

3. Passer d’une table de centroïdes à une table de carreaux vectoriels via btb::btb_ptsToGrid.

Paramètres obligatoires :

• df : un tableau avec les colonnes x et y représentant les coordonnées des centroïdes de la grille ;
• sEPSG : une chaîne de caractères indiquant 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_centroides,
                          sEPSG = "2154", iCellSize = iCellSize)

Les étapes du carroyage [5/7]

4. Se restreindre au champ des carreaux intersectant Paris

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

Les étapes du carroyage [6/7]

On obtient le carroyage des ventes dans Paris intramuros

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

Les étapes du carroyage [7/7]

Remarque :

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.

3. Lissage

3.1 Calcul de densité

Lissage le plus simple à réaliser

• Densité = Quantité par unité de surface (un carreau)
• Ici : densité de ventes de logements dans la ville de Paris au cours de l’année 2021

☠️ Effets de bord

Toujours pour éviter les effets de bord :

• Lissage sur une zone plus large que la zone d’intérêt

➡️ Utilisation de dfBase_filtre (construite précédemment).

Création une variable nbObsLisse

• On crée la variable nbObsLisse = 1 pour chaque observation

# Variable à lisser = nombre d'observations = nombre de ventes
dfBase_filtre$nbObsLisse <- 1
head(dfBase_filtre, 2)

  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
  valeur_fonciere surface_reelle_bati        x       y nbObsLisse
1          480000                  64 647357.3 6868635          1
2          345000                  43 644483.5 6867695          1

➡️ On lissera ensuite la variable nbObsLisse avec la fonction btb::btb_smooth

Lissage et enregistrement dans sfCarrLiss

Paramètres de btb_smooth

• 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 : chaîne 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.

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

Visualisation du lissage [1/2]

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

Visualisation du lissage [2/2]

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

Remarques

• La grille lissée épouse le périmètre des points en entrée.
• Elle va même un peu au-delà (voir infra).

‼️ 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.

Le lissage par densité est conservatif

# 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 sur le secret statistique

‼️ Toujours vérifier le respect du secret statistique avant publication de cartes !

Par exemple, avec données fiscales, pas de carreaux comportant moins de 11 observations, même lissées…

sfCarrLiss <- sfCarrLiss[sfCarrLiss$nbObsLisse >= 11, ]

Représentation cartographique de la densité lissée

• Carte choroplèthe ;
• Discrétisation quantile de la densité lissée ;
• Avec 5 classes de carreaux (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.

Récapitulatif : lissage de densité

On souhaite obtenir le nombre de ventes lissé, en 2021 dans Paris intramuros :

• On sélectionne les ventes au-delà de Paris intramuros (effets de bord)
• On lisse grâce à btb_smooth
• On obtient une grille carroyée bien plus vaste que Paris intramuros
• On filtre uniquement les carreaux dans Paris
• On cartographie la variable lissée (nbObsLissee)

Faire varier le rayon de lissage

Taille optimale du rayon de lissage ?

• Rayon grand (1km) : aspects structurels des données
• Rayon petit (600m) : spécificités locales des données.

➡️ Compromis à faire par le statisticien géographe entre précision et généralisation (connaissance des données et objectifs recherchés)

Faire varier le rayon de lissage (600m)

dataLissage <- dfBase_filtre[, c("nbObsLisse", "x", "y")]
sfCarrLiss <- btb::btb_smooth(pts = dataLissage, 
                              sEPSG = "2154",
                              iCellSize = 200, 
                              iBandwidth = 600)
# Filtrage des carreaux lissés dans 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 600m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Faire varier le rayon de lissage (1000m)

dataLissage <- dfBase_filtre[, c("nbObsLisse", "x", "y")]
sfCarrLiss <- btb::btb_smooth(pts = dataLissage, 
                              sEPSG = "2154",
                              iCellSize = 200, 
                              iBandwidth = 1000)
# Filtrage des carreaux lissés dans 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 1000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Faire varier le rayon de lissage (1000m)

Sans grille apparente !

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

Faire varier la taille des carreaux (50m)

• ❌ Trop gros carreaux = effet granuleux : arêtes des carreaux visibles
• ❌ Trop petits carreaux = temps de calcul important

dataLissage <- dfBase_filtre[, c("nbObsLisse", "x", "y")]
sfCarrLiss <- btb::btb_smooth(pts = dataLissage, 
                              sEPSG = "2154",
                              iCellSize = 50, 
                              iBandwidth = 1000)
# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[, "geometry"], left = FALSE)

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

Illustration des effets de bord

Sans effet de bord

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)

mf_init(x=sfCarLiss_paris)
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = c(0, 3, 5, 7, 9, 17),
       # nbreaks = 5,
       border = NA,
       leg_val_rnd = 1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Sans effets de bord, avec R=2000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Avec effets de bord

sfBase_intramuros <- sfBase %>% st_join(paris_sf, left = FALSE)
dfBase_intramuros <- sfBase_intramuros %>% 
  mutate(x = st_coordinates(geometry)[, 1],
         y = st_coordinates(geometry)[, 2],
         nbObsLisse =1) %>% 
  st_drop_geometry() %>% 
  select(nbObsLisse, x, y)

sfCarLiss_intramuros <- btb::btb_smooth(pts = dfBase_intramuros, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 2000)
sfCarLiss_intramuros_paris <- sfCarLiss_intramuros %>%
  st_join(paris_sf, left = FALSE)

mf_init(x = sfCarLiss_intramuros_paris)
mf_map(x = sfCarLiss_intramuros_paris, 
       type = "choro",
       var = "nbObsLisse",
       breaks = c(0, 3, 5, 7, 9, 17),
       # nbreaks = 5,
       border = NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd = 4,
       col = "black", add = TRUE)
mf_layout(title = "Avec effets de bord, avec R=2000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

➡️ Périphérie de Paris : densités artificiellement plus faibles avec une marge

Les grilles automatiques de btb [1/2]

La grille par défaut produite par btb::btb_smooth dépasse les limites de la zone d’étude choisie (ici Paris intramuros) et contient :

• l’ensemble des carreaux contenant au moins une vente…
• + 4 carreaux limitrophes (règle automatique)

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

Les grilles automatiques de btb [2/2]

Modification possible de cette règle avec le paramètre iNeighbor :

• iNeighbor = 2 ➡️ 2 carreaux limitrophes aux carreaux contenant au moins une vente.
• iNeighbor = 0 ➡️ grille uniquement constituée des carreaux contenant au moins une vente (attention à l’effet « gruyère »)

sfCarLiss_intramuros <- btb::btb_smooth(
  pts = dfBase_intramuros, 
  sEPSG = "2154",
  iCellSize = 200, 
  iBandwidth = 2000,
  iNeighbor = 0)

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

Spécificités locales : ville cotière…

On peut souhaiter utiliser une grille de lissage particulière et adaptée à un territoire donné.

Par exemple, sur une ville cotière, on ne veut :

• ni de carreau dans la mer
• ni l’effet « gruyère » du paramétrage iNeighbor = 0.

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.

3.2 Calcul de moyenne

Calcul de moyenne [1/2]

Nous allons nous intéresser désormais au prix moyen des logements vendus à Paris en 2021.

Il faut :

1. Lisser les prix de ventes des biens vendus
2. Lisser le nombre de biens vendus (exactement comme précédemment)
3. Faire le ratio des deux précédents lissages pour chaque carreau de la grille produite.

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

Calcul de moyenne [2/2]

Les prix moyens sont particulièrement élevés dans le centre et l’Ouest de la capitale

Remarque

le lissage en moyenne est potentiellement sensible aux valeurs atypiques (contrairement au lissage quantile).

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

3.3 Calcul de taux

Taux de grands logements

Proportion de ventes comportant plus de 4 pièces principales ?

• indicatrice quatrePieces
• lissage numérateur et dénominateur
• calcul taux
• intersection géographique avec Paris

# Creation indicatrice quatrePieces
dfBase_filtre$quatrePieces <- ifelse(
  dfBase_filtre$nombre_pieces_principales >= 4, 1, 0)

# Lissage des variables `quatrePieces` et `nbObsLisse`
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)

# Carto
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")

➡️ Carte similaire à la carte des prix moyens

➡️ Lissage du prix au mètre-carré

Prix au m²

Il faut lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur).

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

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. Indicateurs sur zonage à façon

Agréger des données individuelles sur un zonage

Calculons le nombre et le prix moyen des transactions immobilières en 2021 dans le triangle d’or ¹.

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

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

1. Ce triangle a été créé manuellement sur le site du Géoportail.

Intersection géographique

2. Intersection entre des points (les transactions) et un polygone (le quartier).

# Attention au système de projection.
st_crs(triangleOr)[["input"]]

[1] "WGS 84"

# Transformation en Lambert93
triangleOr <- st_transform(triangleOr, crs = 2154) 

# Intersection géographique
transac_triangle <- sfBase %>%
  st_join(triangleOr[, "geometry"], left = FALSE)

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

Calcul des statistiques recherchées

3. 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. Conclusion

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

Autres logiciels

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

Références

• Manuel d’analyse spatiale de l’Insee
• Introduction à la géomatique pour le statisticien, Document de travail, Insee
• Github du package BTB