Contexte

L’objectif de ce script est de produire un lot de données identifiant les différences entre deux millésimes de consommation d’ENAF calculée à partir d’OCS GE donnés pour un département.

En fonction du paramétrage du fichier config.yml, le script peut être exécuté pour un département de la région ou successivement pour tous les départements.

Attention : les temps de traitement sont longs dans ce dernier cas (plusieurs heures).

Configuration

Chargement des librairies

library(datalibaba)
library(DBI)
library(yaml)

Chargement des variables

Attention : il faut vérifier que les paramètres spécifiés dans le fichier config.yml correspondent bien à l’environnement de travail souhaité.

On récupère ensuite les paramètres suivants via le fichier config.yml :

millésimes par département des livrables OCS GE
base de données de travail
schéma de stockage de l’OCS GE
schéma de travail
préfixe du nom de la table à créer
rôle de connexion

# Charger le fichier YAML
config <- yaml::read_yaml("../config.yml")

# Récupérer la liste annee_dept et les paramètres de connexion à la base
annee_dept <- config$annee_dept
database <- config$database
ocs_ge_schema <- config$ocs_ge_schema
work_schema <- config$work_schema
prefix_table <- config$prefix_table
role <- config$role

Paramétrage des scripts

On configure ici :

le préfixe du nom de la table ou des tables de différentiel à générer,
la liste des couples millésimes / département à traiter.

# Préfixe final incluant "diff_"
prefix_table_diff <- paste0(prefix_table, "diff_")
# Exemple : "r_consommation_enaf_ocsge_diff_"

# Construction d'une liste : département → c(year1, year2)
dept_years <- list()

for (year in names(annee_dept)) {
  for (dept in annee_dept[[year]]) {
    dept <- as.character(dept)
    dept_years[[dept]] <- sort(c(dept_years[[dept]], as.integer(year)))
  }
}

Création des tables du différentiel par département

Pour chaque couple millésimes / département :

une table temporaire du recouvrement entre les deux millésimes choisis est créée par analyse spatiale,
la différence entre les valeurs du champs consommation_enaf entre les deux millésimes est ensuite stockée dans une table.

# Boucle : chaque département possède exactement 2 millésimes
for (dept in names(dept_years)) {
  
  first_year <- dept_years[[dept]][1]
  last_year  <- dept_years[[dept]][2]
  
  # Nom final de la table
  table_name <- paste0(work_schema, ".", prefix_table_diff, first_year, "_", last_year, "_s_", dept)

  # Connexion
  connexion <- datalibaba::connect_to_db(db = database, user = role)

  # Requête overlay
  overlay_query <- sprintf(
    "DROP TABLE IF EXISTS overlay;
    CREATE TEMP TABLE overlay AS
    SELECT 
      row_number() OVER () AS uid,
      a.id AS id_%s,
      b.id AS id_%s,
      a.consommation_enaf AS conso_%s,
      b.consommation_enaf AS conso_%s,
      ST_Intersection(a.the_geom, b.the_geom) AS the_geom
    FROM foncier_sol.r_consommation_enaf_ocsge_%s_s_%s a
    JOIN foncier_sol.r_consommation_enaf_ocsge_%s_s_%s b
      ON ST_Intersects(a.the_geom, b.the_geom)
    WHERE NOT ST_IsEmpty(ST_Intersection(a.the_geom, b.the_geom));",
    first_year, last_year,
    first_year, last_year,
    first_year, dept,
    last_year,  dept)

  # Requête création table diff
  create_table_query <- sprintf(
    "CREATE TABLE %s AS
    SELECT *
    FROM overlay
    WHERE conso_%s IS DISTINCT FROM conso_%s;",
    table_name,
    first_year, last_year )

  # Exécution
  DBI::dbExecute(connexion, overlay_query)
  DBI::dbExecute(connexion, create_table_query)

  DBI::dbDisconnect(connexion)

  message(sprintf(
    "Table créée : %s (département %s, années %s → %s)",
    table_name, dept, first_year, last_year))
}

Vérification des géométries (optionnel)

L’opération de recouvrement génère la plupart du temps des scories de géométrie.

Il faut donc vérifier le type de géométrie créée par la requête pour chaque enregistrement :

# Boucle de vérification
for (dept in names(dept_years)) {
  
  first_year <- dept_years[[dept]][1]
  last_year  <- dept_years[[dept]][2]
  
  # Nom complet de la table diff
  table_name <- paste0(
    "zz_traitements.", 
    prefix_table_diff, first_year, "_", last_year, "_s_", dept)

  # Connexion
  connexion <- datalibaba::connect_to_db(db = database, user = role)

  # Requête de vérification
  check_geom_query <- sprintf(
    "SELECT DISTINCT ST_GeometryType(the_geom) AS geom_type
    FROM %s;",
    table_name)

  # Exécution
  geom_types <- DBI::dbGetQuery(connexion, check_geom_query)

  # Déconnexion
  DBI::dbDisconnect(connexion)

  # Message résultat
  message(sprintf(
    "Types de géométrie dans %s : %s",
    table_name,
    paste(unique(geom_types$geom_type), collapse = ", ")
  ))
}

Correction des géométries (optionnel)

Si la vérification indique la création de de points (simples ou multi), de lignes (simples ou multi) ou de collection de géométries, les géométries concernées sont corrigées à l’aide de la requête suivante :

# Boucle de correction des géométries
for (dept in names(dept_years)) {
  
  first_year <- dept_years[[dept]][1]
  last_year  <- dept_years[[dept]][2]
  
  # Nom complet de la table diff
  table_name <- paste0(
    "zz_traitements.", 
    prefix_table_diff, first_year, "_", last_year, "_s_", dept)
  
  # Connexion
  connexion <- datalibaba::connect_to_db(db = database, user = role)

  # Mise à jour des Points / LineStrings -> MultiPolygon via buffer(0)
  update_point_linestring_query <- sprintf(
    "UPDATE %s
    SET the_geom = ST_Multi(ST_Buffer(the_geom, 0))
    WHERE ST_GeometryType(the_geom) IN 
      ('ST_LineString', 'ST_MultiLineString', 'ST_Point', 'ST_MultiPoint');",
    table_name)

  # Mise à jour des GeometryCollection -> extraction polygonale
  update_geom_collection_query <- sprintf(
    "UPDATE %s
    SET the_geom = ST_Multi(ST_CollectionExtract(the_geom, 3))
    WHERE ST_GeometryType(the_geom) = 'ST_GeometryCollection';",
    table_name)

  # Exécution des requêtes
  DBI::dbExecute(connexion, update_point_linestring_query)
  DBI::dbExecute(connexion, update_geom_collection_query)

  # Déconnexion
  DBI::dbDisconnect(connexion)

  # Message résultat
  message(sprintf(
    "Correction géométrique effectuée pour %s",
    table_name))
}

Ajout des contraintes

Une fois les géométries uniquement de type polygone (simple ou multi), on ajoute sur la table du différentiel par département les contraintes géomatiques suivantes :

# Boucle d’ajout des contraintes
for (dept in names(dept_years)) {
  
  first_year <- dept_years[[dept]][1]
  last_year  <- dept_years[[dept]][2]
  
  # Nom complet de la table diff
  table_name <- paste0(
    "zz_traitements.",
    prefix_table_diff, first_year, "_", last_year, "_s_", dept)
  
  # Connexion
  connexion <- datalibaba::connect_to_db(db = database, user = role)
  
  # Requête dynamique d'ajout des contraintes
  add_constraints_query <- sprintf(
    "ALTER TABLE %s
    ADD CONSTRAINT enforce_dims_geom CHECK (st_ndims(the_geom) = 2),
    ADD CONSTRAINT enforce_geotype_geom CHECK (
        geometrytype(the_geom) = 'MULTIPOLYGON'
        OR geometrytype(the_geom) = 'POLYGON'
        OR the_geom IS NULL
    ),
    ADD CONSTRAINT enforce_srid_geom CHECK (st_srid(the_geom) = 2154);",
    table_name
  )

  # Exécution
  DBI::dbExecute(connexion, add_constraints_query)

  # Déconnexion
  DBI::dbDisconnect(connexion)

  # Message d'information
  message(sprintf(
    "Contraintes ajoutées sur la table : %s",
    table_name))
}

Ajout d’un index spatial

On ajoute enfin un index spatial sur la géométrie afin d’optimiser le chargement de la table :

# Boucle de création d'index GIST sur the_geom
for (dept in names(dept_years)) {
  
  first_year <- dept_years[[dept]][1]
  last_year  <- dept_years[[dept]][2]
  
  # Nom complet de la table diff
  table_name <- paste0(
    "zz_traitements.",
    prefix_table_diff, first_year, "_", last_year, "_s_", dept
  )
  
  # Nom de l’index
  index_name <- paste0(
    prefix_table_diff, first_year, "_", last_year, "_s_", dept, "_the_geom_idx"
  )
  
  # Connexion
  connexion <- datalibaba::connect_to_db(db = database, user = role)

  # Requête CREATE INDEX dynamique
  create_idx_query <- sprintf(
    "CREATE INDEX %s
    ON %s
    USING GIST(the_geom);",
    index_name, table_name
  )

  # Exécution
  DBI::dbExecute(connexion, create_idx_query)

  # Déconnexion
  DBI::dbDisconnect(connexion)

  # Message
  message(sprintf(
    "Index créé : %s sur %s",
    index_name, table_name))
}

Création des tables du différentiel entre les millésimes par département

Ronan Vignard

31/03/2026