L’objectif de ce projet pilote est d’évaluer la faisabilité technique d’une analyse croisée entre la base de données d’occupation des sols OCS-GE et la base de données des permis d’urbanisme de Sitadel. La BDD Sitadel ne disposant pas d’informations géospatiales, un recoupement intermédiaire avec les données cadastrales doit être réalisé. Ce projet pilote est menée à l’échelle du département de la Dordogne (24).
1. Préparation des analyses
1.1 Chargement des packages
# Nettoyage de l'espace de travailrm(list=ls())#Chargement des packageslibrary(readxl)library(rmarkdown)library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attachement du package : 'scales'
L'objet suivant est masqué depuis 'package:purrr':
discard
L'objet suivant est masqué depuis 'package:readr':
col_factor
library(maps)
Attachement du package : 'maps'
L'objet suivant est masqué depuis 'package:purrr':
map
library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(geojsonsf)library(leaflet)library(units)
udunits database from C:/Users/jsiracusaphilippe/AppData/Local/R/win-library/4.4/units/share/udunits/udunits2.xml
library(ggrepel)library(viridis)
Le chargement a nécessité le package : viridisLite
Attachement du package : 'viridis'
L'objet suivant est masqué depuis 'package:maps':
unemp
L'objet suivant est masqué depuis 'package:scales':
viridis_pal
library(units)library(lubridate)library(plotly)
Attachement du package : 'plotly'
L'objet suivant est masqué depuis 'package:ggplot2':
last_plot
L'objet suivant est masqué depuis 'package:stats':
filter
L'objet suivant est masqué depuis 'package:graphics':
layout
1.2 Chargement des données
# Données cadastrale (CA) de la Dordogne (24)cadastre24 <-st_read("Data/Cadastre/parcelles.shp", crs =st_crs(2154))
Reading layer `parcelles' from data source
`C:\Users\jsiracusaphilippe\Documents\Datas\R\OCS-GE x Sitadel\Data\Cadastre\parcelles.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1852595 features and 8 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 460301 ymin: 6387603 xmax: 577760.8 ymax: 6515556
Projected CRS: RGF93 v1 / Lambert-93
# Données d'occupation du sol à grande échelle (OCSGE)ocsge24 <-st_read("Data/OCS-GE/OCS-GE_Dordogne_2021/1_DONNEES_LIVRAISON_2024-07-00021/OCSGE_2-0_SHP_LAMB93_D24-2021/OCCUPATION_SOL.shp", crs =st_crs(2154))
Reading layer `OCCUPATION_SOL' from data source
`C:\Users\jsiracusaphilippe\Documents\Datas\R\OCS-GE x Sitadel\Data\OCS-GE\OCS-GE_Dordogne_2021\1_DONNEES_LIVRAISON_2024-07-00021\OCSGE_2-0_SHP_LAMB93_D24-2021\OCCUPATION_SOL.shp'
using driver `ESRI Shapefile'
Simple feature collection with 537120 features and 8 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 460274.3 ymin: 6387574 xmax: 577772.4 ymax: 6515566
Projected CRS: RGF93 v1 / Lambert-93
# Données de la BDD Sitadel sur les permis d'aménager (PA)PA <-read_delim("Data/Sitadel/Liste-des-permis-damenager.2024-07.csv", delim =";", escape_double =FALSE, trim_ws =TRUE)
Rows: 89633 Columns: 17
── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (12): REG_CODE, REG_LIBELLE, DEP_CODE, DEP_LIBELLE, COMM, NUM_PA, SEC_C...
dbl (4): ETAT_PA, AN_DEPOT, SUPERFICIE_TERRAIN, ZONE_OP
date (1): DATE_REELLE_AUTORISATION
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
1.3 Filtrage et préparation des données
# Filtrage des données Dordogne dans le fichier Sitadel des PAPA24 <- PA %>%filter(DEP_CODE =="24") # on conserve les PA avec une superficie de terrain aménagé > 0 et les projets non-annulésPA24 <- PA24 %>%filter(SUPERFICIE_TERRAIN >"0")%>%rename(commune = COMM) %>%mutate(CA1 =paste(SEC_CADASTRE1, NUM_CADASTRE1, sep="."),CA2 =paste(SEC_CADASTRE2, NUM_CADASTRE2, sep="."),CA3 =paste(SEC_CADASTRE3, NUM_CADASTRE3, sep="."))%>%filter(ETAT_PA %in%c("2", "5", "6"))# On retire les variables inutiles à l'analyse PA24_long <- PA24 %>%pivot_longer(cols =c("CA1", "CA2", "CA3"), names_to ="CA_num", values_to ="CA")%>%select(-c("REG_CODE", "REG_LIBELLE", "AN_DEPOT","SEC_CADASTRE1", "NUM_CADASTRE1","SEC_CADASTRE2", "NUM_CADASTRE2","SEC_CADASTRE3", "NUM_CADASTRE3", "ZONE_OP"))# On supprime les données manquantesPA24_long <- PA24_long %>%filter(CA !="NA.NA")# Création d'une nouvelle variable CA combinant "section" et "numéro" cadastre24 <-mutate(cadastre24, CA =paste(section, numero, sep ="."))
1.4 Combinaison des données CA et PA
# Fusion des dataframe CA et PACAxPA24 <-inner_join(x=PA24_long, y=cadastre24, by=c("CA", "commune"))
Warning in inner_join(x = PA24_long, y = cadastre24, by = c("CA", "commune")): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 36 of `x` matches multiple rows in `y`.
ℹ Row 1204309 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
"many-to-many"` to silence this warning.
# Transformation des données en objet simple feature (sf)CAxPA24_sf <-st_as_sf(CAxPA24, crs =2154)
1.4.1 Export des données créées au format geopackage (.gpkg)
Deleting layer `CAxPA24' using driver `GPKG'
Writing layer `CAxPA24' to data source `Output/CAxPA24.gpkg' using driver `GPKG'
Writing 600 features with 16 fields and geometry type Polygon.
Deleting layer `ocsge24' using driver `GPKG'
Writing layer `ocsge24' to data source `Output/ocsge24.gpkg' using driver `GPKG'
Writing 537120 features with 8 fields and geometry type Polygon.
1.5 Combinaisation avec les données OCSGE
inter <-st_intersection(CAxPA24_sf, ocsge24)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
inter <-st_cast(inter, "MULTIPOLYGON")# Renommer les variablesnames(inter) <-c("dep_code", "dep_nom", "CP", "num_PA", "etat_PA", "date_autorisation","PA_superficie_terrain", "CA_num", "CA_id", "CA_id_full","prefixe", "section", "numero", "CA_contenance", "created", "updated","ocsge_id","code_cs", "code_us", "millesime", "source", "ossature" ,"id_origine", "code_or", "geometry")# On réorganise l'ordre des variablesCAxPAxOCS24_all <- inter %>%relocate(num_PA, .before =everything())%>%select(num_PA, etat_PA, dep_code,dep_nom,CP,date_autorisation,PA_superficie_terrain, CA_contenance,ocsge_id,code_cs,code_us,geometry)
1.5.1 Export des données créées au format geopackage (.gpkg)
Deleting layer `CAxPAxOCS24_all' using driver `GPKG'
Writing layer `CAxPAxOCS24_all' to data source
`Output/CAxPAxOCS24_all.gpkg' using driver `GPKG'
Writing 2038 features with 11 fields and geometry type Multi Polygon.
1.6 On sélectionne les observations situées sur des couvertures de sols non bâtis
Code de couverture des sols OCSGE : “CS1.1.1.1” = “Zones baties”, “CS1.1.1.2” = “Zones non baties”, “CS1.1.2.1” = “Zones à matériaux minéraux”, “CS1.1.2.2” = “Zones à autres matériaux composites”, “CS1.2.1” = “Sols nus”, “CS1.2.2” = “Surfaces d’eau”, “CS1.2.3” = “Névés et glaciers”, “CS2.1.1.1” = “Peuplement de feuillus”, “CS2.1.1.2” = “Peuplement de conifères”, “CS2.1.1.3” = “Peuplement mixtes”, “CS2.1.2” = “Formations arbustives et sous-arbrisseaux”, “CS2.1.3” = “Autres formations ligneuses”, “CS2.2.1” = “Formations herbacées”, “CS2.2.2” = “Autres formations non ligneuses”)
1.6.1 Export des données créées au format geopackage (.gpkg)
st_write(obj=CAxPAxOCS24_NAT, # on exporte au format .gpkgdsn="Output/CAxPAxOCS24_NAT.gpkg",delete_layer =TRUE)
Deleting layer `CAxPAxOCS24_NAT' using driver `GPKG'
Writing layer `CAxPAxOCS24_NAT' to data source
`Output/CAxPAxOCS24_NAT.gpkg' using driver `GPKG'
Writing 1318 features with 11 fields and geometry type Multi Polygon.
2. Analyses descriptives
2.1 Analyses des données OCSGE
2.1.1 Calcul de l’aire totale de la couche OCSGE
# Calcul de l'aire des polygones de chauqes observation en m^2ocsge24$surface <-st_area(ocsge24)# Calcul de l'aire totale de la couche OCS-GEocsge24_all_surface <-sum(ocsge24$surface, na.rm =TRUE)print(ocsge24_all_surface)
9210844764 [m^2]
2.1.2 On créer un subset de données
La puissance de calcul de mon PC ne permet pas de procéder à la suite des analyses sur l’ensemble des données. On selectionne les premières 2000 observations parmi les 537120.
ocsge24_subset <- ocsge24[1:2000, ]
2.1.3 Calcul de l’aire totale par type d’occupation du sol (CODE_CS)
# Calcul de l'aire en m^2ocsge24_group_surface_subset <- ocsge24_subset %>%group_by(CODE_CS) %>%summarize(grouped_surface =sum(surface, na.rm =TRUE))# Calcul de l'aire en hectares (ha)ocsge24_group_surface_subset$grouped_surface_ha <- ocsge24_group_surface_subset$grouped_surface /10000print(ocsge24_group_surface_subset$grouped_surface_ha)
Comment faire pour réaliser ces analyses à l’échelle du territoire français ?
Evolution du nombre de permis d’aménager par région/département au cours du temps
Faudra t’il mener les mêmes analyses pour les permis de construire ?
Périmêtre d’analyse : France métropolitaine ou avec les territoires ultramarins ?
Comment prendre en compte le décalage temporel entre OCSGE (fichier figé à une date donnée, 2021 ou 2017 pour la Dordogne ; il existe également un fichier différentiel 2017-2021, les années utilisées varient selon le département) ?
Ce problème peut toutefois conduire à une analyse intéressante de l’évolution de la couverture/usage du sol en fonction des permis d’urbanisme autorisés
l’utilisation des données Theia d’occupation des sols pourraient permettre de résoudre ce problème ? (car actualisation plus fréquente et mise à disposition des différents millésimes ?)
pas de données sur les territoires ultramarins au format vecteur
