Manipulation vecteurs

Objectifs du premier volet

Le but de cette partie est de commencer à manipuler des données de type vecteurs et d’y faire des opérations de base mais revenant souvent. Nous allons travailler sur les communes de France métropolitaine et particulièrement sur les populations de ces communes. À la fin de cette partie, vous serez capables de faire les manipulations suivantes en R :

importer une couche vecteur
exporter une couche vecteur
vérifier le système de coordonnées de référence (SCR) d’une couche géographique
effectuer une jointure attributaire
effectuer une requête attributaire sur du texte ou des nombres
regrouper des entités selon un attribut commun
effectuer une requête spatiale
regrouper des entités selon un champ commun en calculant des statistiques

Données fournies

Pour cet exercice, nous utiliserons les données fournies suivantes :

la couche des communes de France en GeoPackage (communes_France.gpkg)
le tableur des populations communales (communes_population.xlsx)
la couche des principaux cours d’eau de France en GeoPackage (hydro_France.gpkg)

Mise en route

Pour le traitement des données spatialisées, nous utiliserons la librairie terra. Il existe différentes librairies pour la géomatique mais celle-ci présente l’avantage de travailler avec les deux types d’objets (vecteur et raster). Vous trouverez la documentation complète en ligne, mais nous nous contenterons ici d’explorer seulement quelques fonctionnalités. Installez cette librairie comme expliqué en introduction. Une fois terra installée, il suffit de l’appeler en ajoutant la ligne suivante en haut de votre script.

# Un script pour effectuer des manipulations fondamentales en SIG
# chargement de la librairie terra
library(terra)

Nous pouvons maintenant charger notre couche vecteur via la librairie terra. Nous la chargeons dans un objet de type spatVector que nous nommons comm. Notons que la librairie terra gère de très nombreux formats de fichiers. Nous aurions pu lui fournir du shapefile si notre couche avait été un .shp.

comm <- vect('./Data/communes_France.gpkg')

Formats

Faites attention au chemin d’accès à vos fichiers. Le chemin peut être relatif ou absolu. De plus, sous Windows, il est nécessaire de doubler le séparateur de chemin comme ceci // ou de le renverser **.

Pour vérifier l’intégrité de notre couche nous pouvons la tracer rapidement et afficher ses propriétés et un extrait de sa table attributaire.

plot(comm)

comm

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 35798, 6  (geometries, attributes)
 extent      : 99217.1, 1242417, 6049646, 7110480  (xmin, xmax, ymin, ymax)
 source      : communes_France.gpkg (communes_ign_geofla_2019)
 coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
 names       : INSEE_COM          NOM_COM CODE_DEPT       NOM_DEPT CODE_REG
 type        :     <chr>            <chr>     <chr>          <chr>    <chr>
 values      :     32216 LOURTIES-MONBRUN        32           GERS       76
                   47033 BOUDY-DE-BEAURE~        47 LOT-ET-GARONNE       75
                   32009    ARMOUS-ET-CAU        32           GERS       76
          NOM_REG
            <chr>
 LANGUEDOC-ROUSS~
 AQUITAINE-LIMOU~
 LANGUEDOC-ROUSS~

Que signifient les différentes propriétés ?

Réponse

L’objet est de classe SpatVector, les géométries stockées sont des polygones et il y en 35798 avec 6 attributs dans la table attributaire. La couche a une certaine étendue géographique (extent) dans le système de coordonnées RGF93 v1 / Lambert-93 (EPSG:2154). Enfin, nous avons le nom des différents attributs ainsi que leurs types et un extrait des 3 premières lignes de la table attributaire.

Nous constatons que nous avons plusieurs champs attributaires. Combien en avons-nous ? Avons-nous une information sur la population de chaque commune ?

Sélections attributaires

Les tables attributaires sont assimilables à des dataframes et le même genre d’opération peut s’y pratiquer. Par exemple, il est possible d’extraire une commune selon son nom et de l’isoler dans un nouvel objet.

# extraction d'une commune
ma_commune <- comm[comm$NOM_COM == 'MALAKOFF', ]
ma_commune

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 1, 6  (geometries, attributes)
 extent      : 646723.3, 649643.8, 6856858, 6858568  (xmin, xmax, ymin, ymax)
 coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
 names       : INSEE_COM  NOM_COM CODE_DEPT       NOM_DEPT CODE_REG
 type        :     <chr>    <chr>     <chr>          <chr>    <chr>
 values      :     92046 MALAKOFF        92 HAUTS-DE-SEINE       11
       NOM_REG
         <chr>
 ILE-DE-FRANCE

plot(ma_commune)

Dans la sélection précédente, nous avons sélectionné la commune de Malakoff et avons stocké son polygone et les attributs associés dans une nouvelle variable nommée ma_commune.

Avertissement

Par défaut, R est sensible à la casse du texte ainsi qu’aux accents, soyez vigilants.

Nous pouvons faire de même pour extraire toutes les communes d’un département donné. Extrayons par exemple les communes de la Nièvre (58). Le code du département est un caractère et pas un nombre.

comm_nievre <- comm[comm$CODE_DEPT == '58', ]
comm_nievre

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 310, 6  (geometries, attributes)
 extent      : 688355, 793165.5, 6616955, 6720868  (xmin, xmax, ymin, ymax)
 coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
 names       : INSEE_COM NOM_COM CODE_DEPT NOM_DEPT CODE_REG          NOM_REG
 type        :     <chr>   <chr>     <chr>    <chr>    <chr>            <chr>
 values      :     58011   ARMES        58   NIEVRE       27 BOURGOGNE-FRANC~
                   58111  FACHIN        58   NIEVRE       27 BOURGOGNE-FRANC~
                   58076 CHOUGNY        58   NIEVRE       27 BOURGOGNE-FRANC~

plot(comm_nievre)

À vous de jouer

Extrayez une commune de votre choix, puis toutes les communes d’un département et enfin toutes les communes d’une région.

Comme nous pouvons le voir, notre couche des communes ne contient pas d’informations sur la population. Mais nous avons cette information dans le tableur fourni communes_population.xlsx. Nous allons donc pouvoir effectuer une jointure attributaire. Nous commençons par charger le tableur.

library(readxl)
tab_comm <- read_excel('./Data/communes_population.xlsx')

Nous créons une nouvelle variable, de type dataframe, contenant les populations communales.

Jointures attributaires

Pour effectuer une jointure attributaire il faut que les deux tables possèdent un champ commun de même type. Quel est le champ commun aux deux tables ?

Réponse

Le champ commun est le champ INSEE_COM, de type texte.

Avant d’effectuer la jointure nous allons créer un dataframe intermédiaire qui ne contiendra que les attributs intéressants pour la jointure. Ceci nous évitera d’avoir des champs en double dans le résultat de la jointure. Depuis le dataframe tab_comm nous ne conservons que les champs INSEE_COM et POPULATION.

tab_comm_xtr <- tab_comm[c('INSEE_COM', 'POPULATION')]

Nous effectuons ensuite la jointure entre notre vecteur comm et ce dataframe expurgé selon le champ commun. La fonction de terra à utiliser est merge.

comm_pop <- terra::merge(comm, tab_comm_xtr, by='INSEE_COM')

Nous disposons maintenant d’une couche vecteur des communes avec une information sur la population. Nous pouvons ainsi interroger les communes selon leur population.

# communes inférieures à 100 habitants
comm_inf_100 <- comm_pop[comm_pop$POPULATION < 100, ]

Pareil avec les communes de plus de 50 00 habitants.

# communes supérieures à 50000 habitants
comm_sup_50000 <- comm_pop[comm_pop$POPULATION > 50000, ]

À vous de jouer

Sélectionnez les communes dont la population est comprise entre 100 000 et 500 000 habitants. Puis les communes de plus de 50 000 habitants des Bouches-du-Rhône.

Réponse

Les deux lignes suivantes répondent respectivement aux deux questions.

comm_xtr <- comm_pop[comm_pop$POPULATION > 100000 & comm_pop$POPULATION < 500000, ]
comm_xtr2 <- comm_pop[comm_pop$CODE_DEPT == '13' & comm_pop$POPULATION > 50000, ]

Calculs de champs attributaires

Il est très fréquent de devoir ajouter des nouveaux champs dans une table attributaire. Lorsque nous travaillons avec des entités administratives, ou des polygones en général, le calcul de la superficie de chaque entité est une manipulation très courante. Nous allons ajouter un attribut nommé sup_km2 dans lequel nous stockerons la superficie en kilomètres carrés de chaque commune. La fonction à utiliser dans la librairie terra se nomme expanse. Cette fonction a la particularité de calculer la superficie de chaque entité directement sur l’ellipsoïde utilisé dans le système de coordonnées de référence. Cette fonction, par défaut, ne prend pas en compte la projection cartographique. Il n’est donc pas nécessaire de se poser la question sur la nature de la projection (conforme ou équivalente). Par contre, l’unité de base est le mètres carrés, il faut donc convertir la superficie en kilomètres carrés.

# ajout d'un champ de superficie
comm_pop$sup_km2 <- expanse(comm_pop) / 1000000

À vous de jouer

Quel est le nom de la commune la plus étendue ?

Réponse

Avec la première ligne nous sélectionnons les attributs de la commune la plus étendue. Avec la seconde ligne, nous ne sélectionnons que son nom.

comm_pop[comm_pop$sup_km2 == max(comm_pop$sup_km2)]

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 1, 8  (geometries, attributes)
 extent      : 815129.1, 851818.9, 6249584, 6297071  (xmin, xmax, ymin, ymax)
 coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
 names       : INSEE_COM NOM_COM CODE_DEPT         NOM_DEPT CODE_REG
 type        :     <chr>   <chr>     <chr>            <chr>    <chr>
 values      :     13004   ARLES        13 BOUCHES-DU-RHONE       93
         NOM_REG POPULATION sup_km2
           <chr>      <num>   <num>
 PROVENCE-ALPES~  5.257e+04     757

comm_pop[comm_pop$sup_km2 == max(comm_pop$sup_km2)]$NOM_COM

[1] "ARLES"

Disposant maintenant de la superficie de chaque commune ainsi que les populations communales, nous mourrons d’envie de calculer la densité d’habitants pour chaque commune. Calculez ces densités dans un champ nommé densite_hab. Quelle est la densité communale moyenne en France ?

Réponse

La densité est simplement le rapport entre le nombre d’habitants et la superficie en kilomètres carrés de chaque commune.

comm_pop$densite_hab <- comm_pop$POPULATION / comm_pop$sup_km2
# densité moyenne
mean(comm_pop$densite_hab)

[1] 176.8681

Nous pouvons exporter la couche finale des communes jointe aux populations.

# export de la couche des communes avec les populations
writeVector(comm_pop, './Data_processed/communes_France_pop.gpkg', overwrite=TRUE)

Regroupements d’entités selon un champ

Le regroupement d’entités selon des valeurs communes dans un champ de la table attributaire est une manipulation courante en SIG. Ici, par exemple, nous disposons d’une couche de polygones des communes de France, pour lesquelles nous connaissons leur département d’appartenance d’après les champs NOM_DEPT et CODE_DEPT. Nous allons créer un nouvel objet SpatVector, nommé dep, qui contiendra les polygones des départements de France métropolitaine. Nous allons ainsi regrouper les polygones des communes selon le champ CODE_DEPT (il est toujours préférable de travailler avec des codes). Il faut simplement faire attention aux champs attributaires que nous souhaitons conserver. Tous les champs textuels liés aux communes n’auront plus de sens à l’échelle des départements, nous pourrons les enlever. Par contre nous pourrons regrouper les champs quantitatifs, comme la population, la superficie et la densité selon une statistique de regroupement pertinente. Selon vous, quelles statistiques allons-nous utiliser pour ces trois champs ?

Réponse

Il faut bien voir que certaines quantités peuvent être moyennées et d’autres plutôt sommées. Par exemple, la population peut être sommée, ainsi nous obtiendrons la population totale de chaque département. Il en va de même pour la superficie. Il ne serait pas faux de calculer les moyennes de ces deux champs, mais nous obtiendrions la population communale moyenne du département et la superficie communale moyenne du département, ce qui n’est pas forcément intéressant. Le champ contenant les densités, quant à lui, ne peut pas être sommé. Une somme de densité n’a pas de sens mathématique. Nous pourrions le moyenner afin d’obtenir la densité communale moyenne pour le département, mais ce n’est pas très intéressant. Par conséquent nous recalculerons la densité départementale à posteriori.

Le regroupement peut se faire via la fonction aggregate de terra.

dep <- aggregate(comm_pop, by='CODE_DEPT', dissolve=TRUE, fun="sum")
dep

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 96, 10  (geometries, attributes)
 extent      : 99217.1, 1242417, 6049646, 7110480  (xmin, xmax, ymin, ymax)
 coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
 names       : CODE_DEPT sum_POPULATION sum_sup_km2 sum_densite_hab INSEE_COM
 type        :     <chr>          <num>       <num>           <num> <logical>
 values      :        01      6.195e+05        5785       5.757e+04      <NA>
                      02      5.401e+05        7420        5.31e+04      <NA>
                      03      3.434e+05        7378       2.193e+04      <NA>
   NOM_COM NOM_DEPT CODE_REG         NOM_REG agg_n
 <logical>    <chr>    <chr>           <chr> <int>
      <NA>      AIN       84 AUVERGNE-RHONE~   410
      <NA>    AISNE       32 NORD-PAS-DE-CA~   805
      <NA>   ALLIER       84 AUVERGNE-RHONE~   318

Les champs relatifs aux communes sum_densite_hab, INSEE_COM, CODE_COM et NOM_COM ne sont plus utiles, nous pouvons les supprimer. Nous pouvons également renommer certains champs avec des noms plus explicites. Nous en profitons pour recalculer la densité de population à l’échelle des départements.

# nous supprimons des champs selon leurs noms
dep$sum_densite_hab <- NULL
dep$INSEE_COM <- NULL
dep$CODE_COM <- NULL
dep$NOM_COM <- NULL
# nous renommons des champs
names(dep)[names(dep) == 'sum_POPULATION'] <- 'Population_dep'
names(dep)[names(dep) == 'sum_sup_km2'] <- 'Superficie_dep_km2'
names(dep)[names(dep) == 'agg_n'] <- 'nb_communes_dep'
# calcul de la densité départementale
dep$densite_hab <- dep$Population_dep / dep$Superficie_dep_km2

Nous avons maintenant une couche des départements avec une table attributaire pertinente.

À vous de jouer

Faites une manipulation similaire pour créer une couche des régions cette fois-ci.

Réponse

Comme pour les départements, nous effectuons la manipulation en deux temps, d’abord le regroupement par région, puis le ménage dans la table attributaire.

# regroupement
regions <- aggregate(comm_pop, by='CODE_REG', dissolve=TRUE, fun="sum")
# ménage
regions$sum_densite_hab <- NULL
regions$INSEE_COM <- NULL
regions$CODE_COM <- NULL
regions$NOM_COM <- NULL
names(regions)[names(regions) == 'sum_POPULATION'] <- 'Population_region'
names(regions)[names(regions) == 'sum_sup_km2'] <- 'Superficie_region_km2'
names(regions)[names(regions) == 'agg_n'] <- 'nb_communes_region'

Une fois que nous disposons de ces deux nouvelles couches, nous pouvons les exporter en GeoPackage afin de pouvoir les ouvrir dans un logiciel de SIG ou autre outil.

writeVector(dep, './Data_processed/dep_France.gpkg', overwrite=TRUE)
writeVector(regions, './Data_processed/regions_France.gpkg', overwrite=TRUE)

Sélections spatiales

Nous allons maintenant sélectionner des communes selon un critère spatial vis-à-vis d’un département donné. Pour l’exemple, nous allons extraire le département de la Nièvre de notre variable dep contenant tous les départements puis nous allons sélectionner toutes les communes des départements voisins qui touchent la Nièvre.

# extraction de la Nièvre
dep_58 <- dep[dep$CODE_DEPT == '58', ]
# export de la couche du département (décommentez la ligne suivante)
#writeVector(dep_58, './Data_processed/dep_58.gpkg', overwrite=TRUE)

La sélection spatiale se fait ensuite en deux temps. Tout d’abord, pour chaque commune nous regardons si oui ou non (True ou False) elle touche notre département. Puis parmi toutes les communes nous sélectionnons celles qui sont à True.

# les communes qui touchent ou non
relation_touche <- relate(comm_pop, dep_58, relation="touches")
# sélection des communes qui répondent True au fait de toucher
comm_adjacentes <- comm_pop[relation_touche[,1], ]

À vous de jouer

Combien de communes touchent la région Île-de-France ?

Réponse

Comme précédemment cette sélection se fait en trois temps : création d’une variable SpatVector de l’Île-de-France, puis attribution True ou False selon qu’une commune touche ou non ce polygone, et enfin extraction des polygones à True.

# le polygone de l'Île-de-France
reg_idf <- regions[regions$NOM_REG == 'ILE-DE-FRANCE', ]
# les communes qui touchent ou non
relation_touche <- relate(comm_pop, reg_idf, relation="touches")
# sélection des communes qui répondent True au fait de toucher
comm_adjacentes_idf <- comm_pop[relation_touche[,1], ]

Croisement avec l’hydrographie

Nous allons maintenant croiser les données communales avec les principaux cours d’eau. Dans un premier temps, chargeons le réseau hydrographique fourni dans le fichier hydro_France.gpkg. Cette couche répertorie tous les cours d’eau dont l’ordre de Strahler est supérieur à 3, donc les rivières d’une certaine importance.

hydro <- vect('./Data/hydro_France.gpkg')

Quel est le type de géométrie de cette couche ?

Réponse

Pour le savoir, il suffit d’entrer le nom de la variable dans la console et les informations basiques s’affichent.

hydro

 class       : SpatVector 
 geometry    : lines 
 dimensions  : 8849, 14  (geometries, attributes)
 extent      : -4.764583, 9.547917, 41.56875, 51.04792  (xmin, xmax, ymin, ymax)
 source      : hydro_France.gpkg (hydro)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       : HYRIV_ID NEXT_DOWN MAIN_RIV LENGTH_KM DIST_DN_KM DIST_UP_KM
 type        :    <int>     <int>    <int>     <num>      <num>      <num>
 values      : 20347376         0 20347376      1.85          0       39.2
               20347507  20347376 20347376      1.13        2.2       37.1
               20348943         0 20348943      4.59          0       92.5
 CATCH_SKM UPLAND_SKM ENDORHEIC DIS_AV_CMS (and 4 more)
     <num>      <num>     <int>      <num>             
      2.43        406         0      6.332             
      1.22      379.3         0      5.899             
      12.3      933.6         0      17.58

La couche contient des lignes (geometry : lines).

Comment se nomme le champ contenant l’ordre de Strahler ?

Réponse

Pour afficher de façon exhaustive le nom de tous les attributs de la couche, nous utilisons la fonction names.

names(hydro)

 [1] "HYRIV_ID"   "NEXT_DOWN"  "MAIN_RIV"   "LENGTH_KM"  "DIST_DN_KM"
 [6] "DIST_UP_KM" "CATCH_SKM"  "UPLAND_SKM" "ENDORHEIC"  "DIS_AV_CMS"
[11] "ORD_STRA"   "ORD_CLAS"   "ORD_FLOW"   "HYBAS_L12"

Le champ contenant l’ordre de Strahler de chaque cours d’eau se nomme ORD_STRA.

Quel est le système de coordonnées de référence de cette couche d’hydrographie ? Que préconisez-vous ?

Réponse

Nous pouvons afficher le SCR d’une couche via la fonction crs de terra et l’option describe=TRUE.

crs(hydro, describe=TRUE)

    name authority code  area             extent
1 WGS 84      EPSG 4326 World -180, 180, -90, 90

Notre couche est non projetée et utilise le système mondial WGS84. Il est donc nécessaire de la reprojeter dans le système français Lambert93 (ESPG:2154) avant de la croiser avec nos autres couches. Cela se fait simplement avec la fonction project de terra.

# nous créons une nouvelle couche hydro en Lambert 93
hydro_L93 <- project(hydro, 'EPSG:2154')

Nous disposons maintenant de la couche des principaux cours d’eau français dans le bon SCR. Nous pouvons donc la comparer à nos autres couches de communes ou départements. Il peut être intéressant de connaître la population totale des communes traversées par un cours d’eau important. Extrayez les communes traversées par un grand cours d’eau, quelle fonction SIG allez-vous utiliser ?

Réponse

Ce type de traitements fait appel à une sélection spatiale. Nous allons croiser la couche des communes et celle des cours d’eau selon le même principe que nous avons suivi précédemment, mais avec le bon prédicat géométrique.

# identification des communes traversées par un cours d'eau
relation_hydro <- relate(comm_pop, hydro_L93, relation="intersects")
# extraction des communes qui répondent à cette relation
relation_hydro <- cbind(relation_hydro, apply(relation_hydro, 1, any))
comm_hydro <- comm_pop[relation_hydro[,ncol(relation_hydro)]]

Vous remarquerez qu’il y a une ligne en plus que précédemment. Elle est nécessaire sinon seule la première commune traversée par un cours d’eau est sélectionnée. Acceptez la syntaxe telle qu’elle est.

Nous pouvons vérifier le résultat en le visualisant et exporter la couche créée.

plot(comm_hydro)

writeVector(comm_hydro, './Data_processed/communes_hydro.gpkg', overwrite = TRUE)

La population totale des communes traversées par des cours d’eau peut être calculée à partir de cette couche.

pop_hydro <- sum(comm_hydro$POPULATION)
pop_hydro

[1] 33876693

Au total, 3.3876693^{7} habitants habitent dans une commune traversée par un cours d’eau d’ordre de Strahler supérieur à 3.

Après ces quelques manipulations fondamentales sur des couches vecteurs, nous allons voir dans le deuxième volet comment effectuer de nouveaux géotraitements et interagir avec de la donnée raster.