[GEO1] Manipuler les vecteurs avec R et le package sf
GEO UNIV’R Tunisie 2024
Ce support est inspiré en de nombreux points du manuel Géomatique avec R (Giraud and Pecout (2024)). Les exemples proposés sont tantôt adaptés aux contexte tunisien ou au contexte africain. On les en remercie chaleureusement 🙏
Packages utilisés
readxl: importer des fichiers Excel.sf: importer, manipuler et exporter des données géographiques vectorielles.mapsf: cartographie thématiqueunits: Manipuler les unités des objets Rdplyr: Manipuler les donnéesrmapshaper: simplifier des géométriesterra: importer, manipuler et exporter des fichiers raster.
Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.
Les formats vectoriels - un retour rapide
Les données spatiales vectorielles correspondent à une géométrie qui peut être ponctuelle, linéaire ou zonale. Plusieurs formats de données géospatiales existent : shapefile (shp), GeoPackage (GPKG), GeoJSON (geojson), sqlite (sqlite), geodatabase d’ESRI (FileGDB), etc…
Nous nous attardons rapidement ici aux deux formats utilisés dans cette école d’été, le shapefile (shp) et le GéoPackage (gpkg).
Le shapefile
Le shapefile est un format de données spatiales vectorielles propriétaire édité par l’entreprise ESRI, il se présente comme une collection de fichiers portant le même nom avec les extensions et les usages suivants :
.shp: stocke le format de la forme, sa géométrie.dbf: données attributaires relatives aux géométries du.shpcontenus dans le shapefile ;.shx: index de la géométrie..prj: système de coordonnées au format WKT (well-known text)
Chaque fichier doit porter le même nom et est indispensable au bon fonctionnement du shapefile.
D’autres fichiers annexes peuvent aussi être fournis ou générés (indexes de formes, d’attributs, métadonnées au format xml…).
Le GeoPackage
Le GeoPackage est un format de données spatiales vectorielles et raster non propriétaire et défini selon les standards de l’Open Geospatial Consortium (OGC). Il se comprend comme un fichier base de données qui peut contenir des objets vectoriels, les tuiles d’images, des données raster, des schémas de métadonnées comme ceux des normes INSPIRE.
Hormis le format ouvert et l’interopérabilité, l’avantage du géopackage est qu’il peut se requêter comme une base de données avec du SQL dans QGIs, R ou ailleurs, et qu’il peut stocker dans un même fichier plusieurs couches raster ou vecteur. Cela facilite grandement l’organisation et la gestion de projets en géomatique.
On peut par exemple organiser son travail autour du GeoPackage de différente façon :
- un géopagkage par projet
- un géopackage par echelle ou zone de travail (monde, tunisie)
- un géopagkage par type d’utilisation (géotraitements et versionnage ou cartographie)
La figure suivante illustre la différence d’organisation ou d’arborescence fichiers dans le cas d’utilisation de fichiers au format shapefile ou GéoPackage
Le package sf
Le package sf (Pebesma and Bivand 2023), publié en 2016 par Edzer Pebesma permet l’import, l’export, la manipulation et l’affichage de données spatiales vectorielles. Pour cela sf s’appuie sur une série de bibliothèques spatiales : GDAL (GDAL/OGR contributors, 2022), PROJ (PROJ contributors, 2021) pour les opérations d’import, d’export et de projection, et GEOS (GEOS contributors, 2021) pour les opérations de géotraitement (buffer, intersection…). Ce package propose des objets simples (suivant le standard simple feature) dont la manipulation est assez aisée. Une attention particulière a été portée à la compatibilité du package avec la syntaxe pipe (|> ou %>%) et les opérateurs du tidyverse (Wickham et al., 2019).
Anatomie d’un objet sf
Un objet sf se présente comme un tableau de données data.frame, soit comme une table attributaire, à laquelle on ajoute une colonne geom ou geometry spécifique, de classe sfc (simple feature column) contenant les géométries (simple feature geometry). Chaque ligne, chaque individu est appelé simple feature.
Dans l’écosystème R, le data.frame est conçu comme une structure de vecteurs. Chaque colonne correspond à un vecteur et peut donc prendre tous les types que peuvent prendre des vecteurs (character, integer, numeric, boolean etc…). On peut donc effectuer sur les objets sf les même manipulations que l’on fait sur les vecteurs et tableaux dans R.
La colonne contenant les géométries peut quant à elle prendre les types de géométrie traditionnels pris en charge par les GeoJSON :
| type | description |
|---|---|
POINT |
géométrie à zéro dimension contenant un seul point |
LINESTRING |
séquence de points reliés par des morceaux de lignes droites qui ne se coupent pas entre elles ; géométrie unidimensionnelle |
POLYGON |
Géométrie à aire positive (bidimensionnelle) ; une séquence de points forme un anneau fermé, non auto-intersecté ; le premier anneau désigne l’anneau extérieur, zéro ou plusieurs anneaux suivants désignent des trous dans cet anneau extérieur |
MULTIPOINT |
ensemble de points ; un MULTIPOINT est simple si aucun des points du MULTIPOINT n’est égal. |
MULTILINESTRING |
ensemble de LINESTRING |
MULTIPOLYGON |
ensemble de POLYGON |
GEOMETRYCOLLECTION |
ensemble de géométries de tout type sauf GEOMETRYCOLLECTION |
Cette colonne se prête aux géotraitements typiques d’un SIG comme le buffer, l’intersection ou l’agrégation.
Importer et Exporter des données spatiales vectorielles avec sf
Importer des données
Chargement de la librairie
La fonction st_read() de sf permet d’importer des formats de données géographiques variés (shapefile .shp, GeoPackage .GPKG, GeoJSON .geojson …). Avant l’importation d’un fichier au format Geopackage il est préconnisé de consulter le contenu du fichier.
Pour cela on utilise la fonction st_layers()
Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 delegation Multi Polygon 266 11 ETRS89-extended / LAEA Europe
2 gouvernorat Multi Polygon 24 4 ETRS89-extended / LAEA Europe
3 region Multi Polygon 6 2 ETRS89-extended / LAEA EuropeCette fonction nous donne quelques résumés sur chacune des couches du GeoPackage. Le nom, le type de géométrie, le nombre d’individus, le nombre de champs, et le CRS, la projection de chaque couche.
Nous pouvons à présent importer la couche “delegation” dans l’objet “del”
Reading layer `delegation' from data source
`/Users/claudegrasland1/worldregio/geounivr2024/data/tun/geom/tun_admin.gpkg'
using driver `GPKG'
Simple feature collection with 266 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4089658 ymin: 807127.5 xmax: 4473286 ymax: 1584946
Projected CRS: ETRS89-extended / LAEA EuropeExporter des données
sf peut lire plusieurs formats géospatiaux et peut donc tout autant écrire ces formats à l’aide de la fonction st_write().
Cette fonction prend trois paramètres indispensables en entrée, l’objet à écrire obj =, son chemin d’écriture dsn =, la couche à écrirelayer =.
Pour l’exercice nous allons construire un GeoPackage de projet nommé mar1Vector.gpkg, dans lequel nous allons enregistrer notre couche del.
Remarquez qu’ici le chemin d’accès prend le nom complet du fichier et son extension. Ici le nom du GeoPackage et l’extension .gpkg. Si nous avions voulu exporter un GeoJSON, nous aurions écrit dsn = delegation.geojson
Quelques paramètres additionnels sont utiles à connaître. Par exemple pour “écraser” ou remplacer une couche déja existante on utilise l’argument delete_layer = TRUE, lorsqu’il s’agit de remplacer le fichier delete_dsn = TRUE.
D’autres paramétrages sont possibles pour les autres types d’exports, pour cela référez-vous à la vignette du package
Premieres explorations
L’objet et la table attributaire (le data.frame)
Puisque l’objet sf est un data.frame, on peut faire des opérations typiques telles que :
head() pour en visualiser un extrait :
Simple feature collection with 6 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4184398 ymin: 1282883 xmax: 4380989 ymax: 1539573
Projected CRS: ETRS89-extended / LAEA Europe
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code
1 TN.SF.AG Agareb TS2146 TN347 عقارب TN.SF
2 TN.JE.AD Ain Drahem TS1224 TN225 عين دراهم TN.JE
3 TN.SS.AK Akouda TS2115 TN316 اكودة TN.SS
4 TN.KR.AL Alaa TS2216 TN417 العلاء TN.KR
5 TN.BJ.AM Amdoun TS1212 TN213 عمدون TN.BJ
6 TN.AN.AR Ariana Ville TS1120 TN121 أريانة المدينة TN.AN
gou_nom gou_cap gou_cap_dist reg_code reg_nom
1 Sfax 0 30.1 CE Centre-est
2 Jendouba 0 38.5 NO Nord-ouest
3 Sousse 0 14.3 CE Centre-est
4 Kairouan 0 48.6 CO Centre-ouest
5 Beja 0 18.9 NO Nord-ouest
6 Ariana 1 1.0 NE Nord-est
geom
1 MULTIPOLYGON (((4375308 130...
2 MULTIPOLYGON (((4224769 152...
3 MULTIPOLYGON (((4373533 142...
4 MULTIPOLYGON (((4306097 141...
5 MULTIPOLYGON (((4243838 153...
6 MULTIPOLYGON (((4338370 153...str() pour connaitre le type ou la classe de chaque champ :
Classes 'sf' and 'data.frame': 266 obs. of 12 variables:
$ del_code : chr "TN.SF.AG" "TN.JE.AD" "TN.SS.AK" "TN.KR.AL" ...
$ del_nom_fr : chr "Agareb" "Ain Drahem" "Akouda" "Alaa" ...
$ del_code_riate: chr "TS2146" "TS1224" "TS2115" "TS2216" ...
$ del_code_ins : chr "TN347" "TN225" "TN316" "TN417" ...
$ del_nom_ar : chr "عقارب" "عين دراهم" "اكودة" "العلاء" ...
$ gou_code : chr "TN.SF" "TN.JE" "TN.SS" "TN.KR" ...
$ gou_nom : chr "Sfax" "Jendouba" "Sousse" "Kairouan" ...
$ gou_cap : int 0 0 0 0 0 1 0 0 0 0 ...
$ gou_cap_dist : num 30.1 38.5 14.3 48.6 18.9 1 36.6 15.5 17.1 28.9 ...
$ reg_code : chr "CE" "NO" "CE" "CO" ...
$ reg_nom : chr "Centre-est" "Nord-ouest" "Centre-est" "Centre-ouest" ...
$ geom :sfc_MULTIPOLYGON of length 266; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:204, 1:2] 4375308 4375308 4379546 4379546 4379546 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "geom"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
..- attr(*, "names")= chr [1:11] "del_code" "del_nom_fr" "del_code_riate" "del_code_ins" ...colnames() pour connaitre le nom chaque champ :
[1] "del_code" "del_nom_fr" "del_code_riate" "del_code_ins"
[5] "del_nom_ar" "gou_code" "gou_nom" "gou_cap"
[9] "gou_cap_dist" "reg_code" "reg_nom" "geom" ou summary() pour avoir un résumé statistique de chaque champ :
del_code del_nom_fr del_code_riate del_code_ins
Length:266 Length:266 Length:266 Length:266
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
del_nom_ar gou_code gou_nom gou_cap
Length:266 Length:266 Length:266 Min. :0.00000
Class :character Class :character Class :character 1st Qu.:0.00000
Mode :character Mode :character Mode :character Median :0.00000
Mean :0.09023
3rd Qu.:0.00000
Max. :1.00000
gou_cap_dist reg_code reg_nom geom
Min. : 1.00 Length:266 Length:266 MULTIPOLYGON :266
1st Qu.: 14.32 Class :character Class :character epsg:3035 : 0
Median : 24.40 Mode :character Mode :character +proj=laea...: 0
Mean : 28.66
3rd Qu.: 41.38
Max. :147.40 Les géométries
Pour connaitre le type de géométries on utilise la fonction st_geometry_type(). On peut rajouter l’argument by_geometry = FALSE pour avoir le type de l’ensemble.
[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEOn peut aussi visualiser les géométries avec la fonction plot() native de R.
Dans ce premier cas on obtient une carte par champs de la table attributaire
Il existe deux façons de ne veut visualiser que les géométries
la fonction st_geometry() qui renvoi uniquement aux géométrie de l’objet sf
plot(st_geometry(del))
ou on sélectionne le champs contenant les géométries
plot(del$geom)
Pour plus de détail sur la cartographie avec le package sf voir le notebook de N. Lambert conçu pour l’école thématique : Faire des cartes thématiques avec R
Les systèmes de coordonnées
Pour connaître le détail du système de coordonnées de référence d’un objet sf on utilise la fonction st_crs()
Coordinate Reference System:
User input: ETRS89-extended / LAEA Europe
wkt:
PROJCRS["ETRS89-extended / LAEA Europe",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
MEMBER["European Terrestrial Reference Frame 2020"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["Europe Equal Area 2001",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (Y)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (X)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Statistical analysis."],
AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
BBOX[24.6,-35.58,84.73,44.83]],
ID["EPSG",3035]]Ceci nous renvoi la définition complète du système de coordonnées selon les standards en WKT.
Mais on peut aussi récupérer des informations spécifiques relatives à la projection. Pour cela on utilise l’argument parameters = TRUE.
Il faut alors soit enregistrer le résultat de cette fonction dans un objet puis acceder aux paramètres via l’objet…
crsDel
$SemiMajor
6378137 [m]
$SemiMinor
6356752 [m]
$InvFlattening
[1] 298.2572
$IsGeographic
[1] FALSE
$units_gdal
[1] "metre"
$IsVertical
[1] FALSE
$WktPretty
[1] "PROJCS[\"ETRS89-extended / LAEA Europe\",\n GEOGCS[\"ETRS89\",\n DATUM[\"European_Terrestrial_Reference_System_1989\",\n SPHEROID[\"GRS 1980\",6378137,298.257222101]],\n PRIMEM[\"Greenwich\",0],\n UNIT[\"degree\",0.0174532925199433,\n AUTHORITY[\"EPSG\",\"9122\"]],\n AUTHORITY[\"EPSG\",\"4258\"]],\n PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],\n PARAMETER[\"latitude_of_center\",52],\n PARAMETER[\"longitude_of_center\",10],\n PARAMETER[\"false_easting\",4321000],\n PARAMETER[\"false_northing\",3210000],\n UNIT[\"metre\",1],\n AXIS[\"Northing\",NORTH],\n AXIS[\"Easting\",EAST],\n AUTHORITY[\"EPSG\",\"3035\"]]"
$Wkt
[1] "PROJCS[\"ETRS89-extended / LAEA Europe\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",52],PARAMETER[\"longitude_of_center\",10],PARAMETER[\"false_easting\",4321000],PARAMETER[\"false_northing\",3210000],UNIT[\"metre\",1],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"3035\"]]"
$Name
[1] "ETRS89-extended / LAEA Europe"
$proj4string
[1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
$epsg
[1] 3035
$yx
[1] TRUE
$ProjJson
[1] "{\n \"$schema\": \"https://proj.org/schemas/v0.7/projjson.schema.json\",\n \"type\": \"ProjectedCRS\",\n \"name\": \"ETRS89-extended / LAEA Europe\",\n \"base_crs\": {\n \"name\": \"ETRS89\",\n \"datum_ensemble\": {\n \"name\": \"European Terrestrial Reference System 1989 ensemble\",\n \"members\": [\n {\n \"name\": \"European Terrestrial Reference Frame 1989\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1990\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1991\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1992\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1993\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1994\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1996\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 1997\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 2000\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 2005\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 2014\"\n },\n {\n \"name\": \"European Terrestrial Reference Frame 2020\"\n }\n ],\n \"ellipsoid\": {\n \"name\": \"GRS 1980\",\n \"semi_major_axis\": 6378137,\n \"inverse_flattening\": 298.257222101\n },\n \"accuracy\": \"0.1\"\n },\n \"coordinate_system\": {\n \"subtype\": \"ellipsoidal\",\n \"axis\": [\n {\n \"name\": \"Geodetic latitude\",\n \"abbreviation\": \"Lat\",\n \"direction\": \"north\",\n \"unit\": \"degree\"\n },\n {\n \"name\": \"Geodetic longitude\",\n \"abbreviation\": \"Lon\",\n \"direction\": \"east\",\n \"unit\": \"degree\"\n }\n ]\n },\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 4258\n }\n },\n \"conversion\": {\n \"name\": \"Europe Equal Area 2001\",\n \"method\": {\n \"name\": \"Lambert Azimuthal Equal Area\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 9820\n }\n },\n \"parameters\": [\n {\n \"name\": \"Latitude of natural origin\",\n \"value\": 52,\n \"unit\": \"degree\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 8801\n }\n },\n {\n \"name\": \"Longitude of natural origin\",\n \"value\": 10,\n \"unit\": \"degree\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 8802\n }\n },\n {\n \"name\": \"False easting\",\n \"value\": 4321000,\n \"unit\": \"metre\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 8806\n }\n },\n {\n \"name\": \"False northing\",\n \"value\": 3210000,\n \"unit\": \"metre\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 8807\n }\n }\n ]\n },\n \"coordinate_system\": {\n \"subtype\": \"Cartesian\",\n \"axis\": [\n {\n \"name\": \"Northing\",\n \"abbreviation\": \"Y\",\n \"direction\": \"north\",\n \"unit\": \"metre\"\n },\n {\n \"name\": \"Easting\",\n \"abbreviation\": \"X\",\n \"direction\": \"east\",\n \"unit\": \"metre\"\n }\n ]\n },\n \"scope\": \"Statistical analysis.\",\n \"area\": \"Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State.\",\n \"bbox\": {\n \"south_latitude\": 24.6,\n \"west_longitude\": -35.58,\n \"north_latitude\": 84.73,\n \"east_longitude\": 44.83\n },\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 3035\n }\n}"
$WKT1_ESRI
[1] "PROJCS[\"ETRS_1989_LAEA\",\n GEOGCS[\"GCS_ETRS_1989\",\n DATUM[\"D_ETRS_1989\",\n SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],\n PRIMEM[\"Greenwich\",0.0],\n UNIT[\"Degree\",0.0174532925199433]],\n PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],\n PARAMETER[\"False_Easting\",4321000.0],\n PARAMETER[\"False_Northing\",3210000.0],\n PARAMETER[\"Central_Meridian\",10.0],\n PARAMETER[\"Latitude_Of_Origin\",52.0],\n UNIT[\"Meter\",1.0]]"
$srid
[1] "EPSG:3035"
$axes
name orientation
1 Northing 1
2 Easting 3
$ud_unit
1 [m]
attr(,"class")
[1] "crs_parameters"… soit accéder aux paramètres directement en utilisant le $
[1] "ETRS89-extended / LAEA Europe"Changer de projection
La fonction st_transform() permet de changer la projection d’un objet sf. Il est possible de passer d’un système de coordonnées géodésique à un système de coordonnées projetées et inversement.
Par exemple on peut reprojeter les délégations en dans le système de référence tunisien, Carthage - EPSG:2028.
[1] "Carthage / TM 11 NE"Ou les reprojeter en WGS84 - World Geodetic System 1984
[1] "WGS 84"Code
# Initalisation de la fenêtre graphique
par(mfrow = c(1,3),
mar = c(0, 2, 5, 2),
xaxs='i', yaxs='i')
# Afficher les délégations avec la projection initiale
# On dessine les géométries dans leur projection d'origine
plot(st_geometry(del), border = "lightblue", lwd = 2, col = NA,
graticule = TRUE)
# On ajoute en titre le nom de la projection et les unités de mesure
title(paste0(crsDel$Name, "\n", crsDel$units_gdal))
# Carthage - EPSG:2088
plot(st_geometry(del2088), border = "lightblue", lwd = 2, col = NA,
graticule = TRUE)
title(paste0(st_crs(del2088, parameters = TRUE)$Name, "\n",
st_crs(del2088, parameters = TRUE)$units_gdal)
)
# WGS84 EPSG:4326
plot(st_geometry(del4326), border = "lightblue", lwd = 2, col = NA,
graticule = TRUE)
title(paste0(st_crs(del4326, parameters = TRUE)$Name, "\n",
st_crs(del4326, parameters = TRUE)$units_gdal)
)
Pour conserver la projection adaptée à la tunisie on reprojete del en EPSG:2088
Sélections et jointure attributaire
Sélectionner par attributs
Comme les data.frame on peut sélectionner des lignes et des colonnes des objets sf en utilisant les crochets :
monSF[lignes , colonnes]
Selectionner des lignes
On peut donc afficher les 5 premières lignes de l’objet del grace à leur index
Simple feature collection with 5 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 275294.7 ymin: 3827226 xmax: 467583.4 ymax: 4087746
Projected CRS: Carthage / TM 11 NE
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code gou_nom
1 TN.SF.AG Agareb TS2146 TN347 عقارب TN.SF Sfax
2 TN.JE.AD Ain Drahem TS1224 TN225 عين دراهم TN.JE Jendouba
3 TN.SS.AK Akouda TS2115 TN316 اكودة TN.SS Sousse
4 TN.KR.AL Alaa TS2216 TN417 العلاء TN.KR Kairouan
5 TN.BJ.AM Amdoun TS1212 TN213 عمدون TN.BJ Beja
gou_cap gou_cap_dist reg_code reg_nom geom
1 0 30.1 CE Centre-est MULTIPOLYGON (((462201.4 38...
2 0 38.5 NO Nord-ouest MULTIPOLYGON (((315372.8 40...
3 0 14.3 CE Centre-est MULTIPOLYGON (((461740 3974...
4 0 48.6 CO Centre-ouest MULTIPOLYGON (((394811.9 39...
5 0 18.9 NO Nord-ouest MULTIPOLYGON (((334323 4081...On peut aussi les sélectionner en fonction d’une valeur de champ. Par exemple :
Simple feature collection with 1 feature and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 463418.8 ymin: 3961466 xmax: 469492.5 ymax: 3969566
Projected CRS: Carthage / TM 11 NE
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code
242 TN.SS.SM Sousse Medina TS2110 TN311 سوسة المدينة TN.SS
gou_nom gou_cap gou_cap_dist reg_code reg_nom
242 Sousse 1 1 CE Centre-est
geom
242 MULTIPOLYGON (((464400.5 39...Selectionner des colonnes
On peut sélectionner les colonnes par leur index, par exemple pour les dernières colonnes :
Simple feature collection with 266 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 178843.6 ymin: 3345646 xmax: 556188.8 ymax: 4132739
Projected CRS: Carthage / TM 11 NE
First 10 features:
gou_nom gou_code del_nom_ar del_code_ins del_code_riate
1 Sfax TN.SF عقارب TN347 TS2146
2 Jendouba TN.JE عين دراهم TN225 TS1224
3 Sousse TN.SS اكودة TN316 TS2115
4 Kairouan TN.KR العلاء TN417 TS2216
5 Beja TN.BJ عمدون TN213 TS1212
6 Ariana TN.AN أريانة المدينة TN121 TS1120
7 Kasserine TN.KS العيون TN428 TS2227
8 Tunis TN.TU باب بحر TN113 TS1112
9 Tunis TN.TU باب سويقة TN114 TS1113
10 Jendouba TN.JE بلطة بوعوان TN229 TS1228
del_nom_fr del_code geom
1 Agareb TN.SF.AG MULTIPOLYGON (((462201.4 38...
2 Ain Drahem TN.JE.AD MULTIPOLYGON (((315372.8 40...
3 Akouda TN.SS.AK MULTIPOLYGON (((461740 3974...
4 Alaa TN.KR.AL MULTIPOLYGON (((394811.9 39...
5 Amdoun TN.BJ.AM MULTIPOLYGON (((334323 4081...
6 Ariana Ville TN.AN.AR MULTIPOLYGON (((428017.4 40...
7 Ayoun TN.KS.AY MULTIPOLYGON (((319476.8 39...
8 Bab Bhar TN.TU.BB MULTIPOLYGON (((427125.6 40...
9 Bab Souika TN.TU.BS MULTIPOLYGON (((425857.3 40...
10 Balta Bou Aouane TN.JE.BB MULTIPOLYGON (((323552 4060...On peut aussi les selectionner par leur nom
Simple feature collection with 266 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 178843.6 ymin: 3345646 xmax: 556188.8 ymax: 4132739
Projected CRS: Carthage / TM 11 NE
First 10 features:
del_nom_fr gou_nom reg_nom geom
1 Agareb Sfax Centre-est MULTIPOLYGON (((462201.4 38...
2 Ain Drahem Jendouba Nord-ouest MULTIPOLYGON (((315372.8 40...
3 Akouda Sousse Centre-est MULTIPOLYGON (((461740 3974...
4 Alaa Kairouan Centre-ouest MULTIPOLYGON (((394811.9 39...
5 Amdoun Beja Nord-ouest MULTIPOLYGON (((334323 4081...
6 Ariana Ville Ariana Nord-est MULTIPOLYGON (((428017.4 40...
7 Ayoun Kasserine Centre-ouest MULTIPOLYGON (((319476.8 39...
8 Bab Bhar Tunis Nord-est MULTIPOLYGON (((427125.6 40...
9 Bab Souika Tunis Nord-est MULTIPOLYGON (((425857.3 40...
10 Balta Bou Aouane Jendouba Nord-ouest MULTIPOLYGON (((323552 4060...Combiner les selections
Enfin on peut combiner les selections
Simple feature collection with 16 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 422897.2 ymin: 3933116 xmax: 469492.5 ymax: 4028665
Projected CRS: Carthage / TM 11 NE
First 10 features:
del_nom_fr gou_nom reg_nom geom
3 Akouda Sousse Centre-est MULTIPOLYGON (((461740 3974...
36 Bouficha Sousse Centre-est MULTIPOLYGON (((458390.3 40...
67 Enfidha Sousse Centre-est MULTIPOLYGON (((452021.2 40...
100 Hammam Sousse Sousse Centre-est MULTIPOLYGON (((461739.9 39...
107 Hergla Sousse Centre-est MULTIPOLYGON (((452737.5 39...
124 Kalaa Kebira Sousse Centre-est MULTIPOLYGON (((454231.3 39...
125 Kalaâ Seghira Sousse Centre-est MULTIPOLYGON (((459976.1 39...
138 Kondar Sousse Centre-est MULTIPOLYGON (((439303 3970...
181 Msaken Sousse Centre-est MULTIPOLYGON (((457987.4 39...
221 Sid Bou Ali Sousse Centre-est MULTIPOLYGON (((455032.2 39...Et on peut afficher cette sélection
Jointure attributaire
Avec la fonction merge() on peut joindre les données d’un autre data.frame à un objet sf et inversement via un champ de jointure commun.
On importe d’autres données du projet
On identifie le champ de jointure
Simple feature collection with 2 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 275294.7 ymin: 3827226 xmax: 467583.4 ymax: 4087746
Projected CRS: Carthage / TM 11 NE
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code gou_nom
1 TN.SF.AG Agareb TS2146 TN347 عقارب TN.SF Sfax
2 TN.JE.AD Ain Drahem TS1224 TN225 عين دراهم TN.JE Jendouba
gou_cap gou_cap_dist reg_code reg_nom geom
1 0 30.1 CE Centre-est MULTIPOLYGON (((462201.4 38...
2 0 38.5 NO Nord-ouest MULTIPOLYGON (((315372.8 40... del_code del_nom_fr del_nom_ar gou_code gou_nom gou_cap gou_cap_dist
1 TN.AN.AR Ariana Ville أريانة المدينة AN Ariana 1 1.0
2 TN.AN.ET Cité Ettathamen حي التضامن AN Ariana 0 8.8
reg_code reg_nom popto_2004 popco_2004 immig_2004 emigr_2004 mobil_2004
1 NE Nord-est 97687 97687 16961 15426 32387
2 NE Nord-est 78311 78311 5651 5245 10896
menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014 immig_2014
1 27468 9751 22524 18596 114486 114486 15637
2 11950 430 4505 3824 84312 84312 5028
emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014 popto_2010
1 20448 36085 32498 25474 32308 19942 109500
2 6752 11780 22087 6836 21715 3496 82922
surfa_2010 idr_2011
1 18.612 0.638
2 3.376 0.386Ici les deux colonnes identiques ont aussi le même nom. On choisi d’utiliser les codes pour la jointure :
Attention, le sens de la jointure est important. Ici l’objet “x” est l’objet auquel on joint le second. L’objet final prend le type de l’objet x. Ici nous avons créé un nouvel objet delMerge qui résulte de la jointure de del_def à del. del_mergeprend donc le type de del. C’est un objetsf.
La ligne all.x signifie que l’on conserve tous les individus du tableau “x” meme si la correspondance est manquante dans le tableau “y”.
Simple feature collection with 3 features and 40 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 411876.4 ymin: 4075789 xmax: 429940.1 ymax: 4109337
Projected CRS: Carthage / TM 11 NE
del_code del_nom_fr.x del_code_riate del_code_ins del_nom_ar.x
1 TN.AN.AR Ariana Ville TS1120 TN121 أريانة المدينة
2 TN.AN.ET Cité Ettathamen TS1125 TN126 حي التضامن
3 TN.AN.KA kalaât El Andalous TS1123 TN124 قلعة الاندلس
gou_code.x gou_nom.x gou_cap.x gou_cap_dist.x reg_code.x reg_nom.x
1 TN.AN Ariana 1 1.0 NE Nord-est
2 TN.AN Ariana 0 8.8 NE Nord-est
3 TN.AN Ariana 0 19.8 NE Nord-est
del_nom_fr.y del_nom_ar.y gou_code.y gou_nom.y gou_cap.y
1 Ariana Ville أريانة المدينة AN Ariana 1
2 Cité Ettathamen حي التضامن AN Ariana 0
3 kalaât El Andalous قلعة الاندلس AN Ariana 0
gou_cap_dist.y reg_code.y reg_nom.y popto_2004 popco_2004 immig_2004
1 1.0 NE Nord-est 97687 97687 16961
2 8.8 NE Nord-est 78311 78311 5651
3 19.8 NE Nord-est 23045 15313 728
emigr_2004 mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014
1 15426 32387 27468 9751 22524 18596 114486
2 5245 10896 11950 430 4505 3824 84312
3 528 1256 4709 188 1865 829 26796
popco_2014 immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014
1 114486 15637 20448 36085 32498 25474 32308
2 84312 5028 6752 11780 22087 6836 21715
3 18211 1104 752 1856 6554 1701 6305
telef_2014 popto_2010 surfa_2010 idr_2011 geom
1 19942 109500 18.612 0.638 MULTIPOLYGON (((428017.4 40...
2 3496 82922 3.376 0.386 MULTIPOLYGON (((420275.4 40...
3 428 24367 188.206 0.383 MULTIPOLYGON (((423864.5 41...Sélections et jointure spatiale
Sélections spatiales
Les sélections spatiales s’exécutent avec la fonction st_filter() et selon les prédicats géométriques suivants :
st_intersects()st_disjoint()st_touches()st_crosses()st_within()st_contains()st_contains_properly()st_overlaps()st_equals()st_covers()st_covered_by()st_equals_exact()st_is_within_distance()
Voir la vignette sf
On importe les géométries extraites de OSM dans le module sur la Manipulation des données
Pour plus de détail sur l’extraction de données OSM voir le notebook de R. Ysebaert conçu pour l’école thématique : Acquisition de données géographiques et visualisations de base
Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 deleg Multi Polygon 266 7 WGS 84
2 sect Multi Polygon 2104 7 WGS 84
3 poi Point 11571 4 WGS 84Reading layer `poi' from data source
`/Users/claudegrasland1/worldregio/geounivr2024/data/tun/geom/tun_osm.gpkg'
using driver `GPKG'
Simple feature collection with 11571 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 5.364165 ymin: 31.00284 xmax: 14.74145 ymax: 44.40715
Geodetic CRS: WGS 84Simple feature collection with 6 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 5.364165 ymin: 33.23918 xmax: 10.86563 ymax: 43.31226
Geodetic CRS: WGS 84
osm_id amenity name_ar name_fr
1 255880618 ferry_terminal محطة آرنك البحرية <NA>
2 262728007 parking <NA> <NA>
3 283583078 pharmacy الصيدلية المركزية Pharmacie Centrale
4 283658474 ferry_terminal <NA> <NA>
5 290891538 bank بنك قطر الوطني تونس Qatar National Bank Tunisia
6 297765095 place_of_worship الغريبة Synagogue de la Ghriba
geom
1 POINT (5.364165 43.31226)
2 POINT (10.86563 33.23918)
3 POINT (10.18092 36.83346)
4 POINT (10.30467 36.80743)
5 POINT (10.5971 35.83918)
6 POINT (10.85939 33.81392)Pour que la sélection spatiale fonctionne il faut s’assurer que les deux objets possèdent le meme système de coordonnées de référence
[1] "EPSG:2088"[1] "EPSG:4326"Ici ils ne concordent pas donc :
# Re projection en
poi <- st_transform(poi, crs = "EPSG:2088")
# Verification du SCR
st_crs(poi, parameters = TRUE)$srid[1] "EPSG:2088"On peut maintenant réaliser notre sélection spatiale. Ici on va sélectionner les points qui se trouvent dans l’objet delSousse que l’on a construit plus haut.
On visualise le résultat
Code
# Parametre de l'affichage
par(mar = c(0, 0, 4, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
# Initialisation de la carte à l'emprise de Sousse
plot(st_geometry(delSousse), col = NA, border = NA)
# Fond de carte des délégations
plot(st_geometry(del), col = "gray80", border = "white", lwd = 1, add = TRUE)
# Délégations de Soussz
plot(st_geometry(delSousse), col = "#5B89A3", border = "white", lwd = 2, add = TRUE)
# Points remarquables de Souss
plot(st_geometry(poiSousse), col = "red", border = "white", pch = 19, cex = .3, add = TRUE)
# Titre
title("Points remarquables \ndu Gouvernorat de Sousse")
Opérations sur les géométries
Extraction de centroides
La fonction st_centroid() permet d’extraire les centroides des polygones.
Ici on extrait les centroides des délégations de Sousse :
Agrégation de polygones
Agrégation spatiale
La fonction st_union() permet d’agréger des polygones entre eux. Par exemple pour reconstituer le gouvernorat de Sousse
Code
# Parametre de l'affichage
par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
# Initialisation de la carte à l'emprise de Sousse
plot(st_geometry(delSousse), col = NA, border = NA)
# Fond de carte des délégations
plot(st_geometry(del), col = "gray80", border = "white", lwd = 1, add = TRUE)
# Délégations de Sousse
plot(st_geometry(delSousse), col = "#5B89A3", border = "white", lwd = 1, add = TRUE)
# Gouvernorat de Sousse
plot(st_geometry(gouSousse), border = "darkblue", lwd = 3, add = TRUE)
Agrégation spatiale et attributaire
On peut aussi agréger les polygones et demander un résumé statistique pour un ou plusieurs champ décrivant ces polygones. Plusieurs méthodes permettent de réaliser cet objectif.
La première méthode mobilise la fonction aggregate() de sf. Cette fonction permet d’agréger les polygones et de demander le même résumé statistique pour plusieurs champs.
Dans cet exemple nous repartons de l’objet delMerge issu de la fusion entre les géométries des délégations et le tableau additionnel. L’objectif est de construire un objet gou représentant les gouvernorats et leur population. Pour y arriver, nous fusionnons toutes les délégations via le champ gou_nom.x et calculons la somme de la population pour chacun.
Code
Les fonctions group_by() et summarise() du package dplyr permettent d’agréger les polygones en demandant des résumés statistiques différents selon les champs.
Simple feature collection with 3 features and 3 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 313571.7 ymin: 4026799 xmax: 445046.3 ymax: 4113741
Projected CRS: Carthage / TM 11 NE
# A tibble: 3 × 4
gou_nom.x pop immig_mean geom
<chr> <int> <dbl> <POLYGON [m]>
1 Ariana 576088 11297. ((419782.2 4077626, 419782.2 4077626, 419782.2 40…
2 Beja 303032 1049. ((350635.1 4027850, 350635.1 4027850, 350635.1 40…
3 Ben Arous 631842 6088 ((442640 4056428, 442640 4056428, 442640 4056428,…Construction d’une zone tampon (buffer)
Pour construire la zone tampon il est préférable de connaître le SCR et l’unité de mesure de l’objet sf, par exemple avec cette fonction.
Les mesures sont exprimées en metres, on peut à présent utiliser la fonction st_buffer() pour construire la zone tampon.
Ici on peut construire une zone tampon autour du centroide de la délégation de Sidi Bou Ali
Et on construit la zone tampon de 3000 m soit 5km
Code
# Parametre de l'affichage
par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
# Initialisation de la carte à l'emprise de Sousse
plot(st_geometry(delSousse), col = NA, border = NA)
# Fond de carte des délégations
plot(st_geometry(del), col = "gray80", border = "white", lwd = 1, add = TRUE)
# Délégations de Sousse
plot(st_geometry(delSousse), col = "#5B89A3", border = "white", lwd = 1, add = TRUE)
# Gouvernorat de Sousse
plot(st_geometry(gouSousse), border = "darkblue", lwd = 3, add = TRUE)
# Zone tampon de 5km autour du centroide de Sidi Bou Ali
plot(st_geometry(sidiBou_t), border = "pink", col = "#fac0cb50", lwd = 2, add = TRUE)
# Centroide de Sidi Bou Ali
plot(st_geometry(sidiBou_c), col = "pink", pch = 20, cex = 2, add = TRUE)
Intersection
La fonction st_intersection() permet de découper une couche par une autre.
On peut ici par exemple, découper la couche poi des points remarquables, par le centroide de la délégation de Sidi Bou Ali.
Code
# Parametre de l'affichage
par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
# Initialisation de la carte à l'emprise de Sousse
plot(st_geometry(delSousse), col = NA, border = NA)
# Fond de carte des délégations
plot(st_geometry(del), col = "gray80", border = "white", lwd = 1, add = TRUE)
# Délégations de Sousse
plot(st_geometry(delSousse), col = "#5B89A3", border = "white", lwd = 1, add = TRUE)
# Gouvernorat de Sousse
plot(st_geometry(gouSousse), border = "darkblue", lwd = 3, add = TRUE)
# Zone tampon de 5km autour du centroide de Sidi Bou Ali
plot(st_geometry(sidiBou_t), border = "pink", col = "#fac0cb50", lwd = 2, add = TRUE)
# Centroide de Sidi Bou Ali
plot(st_geometry(sidiBou_c), col = "pink", pch = 20, cex = 2, add = TRUE)
plot(st_geometry(poi_sidiBou), col = "red", border = "white", pch = 19, cex = .5, add = TRUE)
Compter les points
La fonction st_intersects() permet d’intersecter deux couches sans les découper, et de compter les éléments d’une couche (y) contenue dans une autre (x).
L’argument sparse = TRUE nous permet de lister pour chaque élément de x les objets de y.
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ...Le 1: signifie qu’il y a un objet, une 1 liste. Les chiffres suivants sont les identifiants de chacun des points.
Pour vérifier que le nombre d’intersection corresponde bien au nombre d’objets intersectés de x, ici il s’agit de sidiBou_t on fait :
Pour compter le nombre de points intersectés par sidiBou_t on utilise la fonction lengths() qui renvoie la longueurs de chaque liste de l’intersection.
On peut ajouter cette information à la table attributaire de sidiBou_t en créant un nouveau champ nb_poi
Sur la base de cet exemple on peut compter les points dans les gouvernorats créé dans la Section 7.2.2
Sparse geometry binary predicate list of length 24, where the predicate
was `intersects'
first 10 elements:
1: 68, 125, 149, 160, 173, 174, 177, 186, 193, 195, ...
2: 24, 260, 1206, 1322, 1323, 1708, 1710, 1711, 1712, 1713, ...
3: 69, 70, 71, 72, 73, 74, 75, 96, 97, 262, ...
4: 276, 277, 278, 289, 329, 330, 1297, 1597, 1728, 1729, ...
5: 94, 746, 761, 762, 763, 784, 785, 787, 809, 810, ...
6: 695, 1215, 1216, 1217, 1218, 1249, 1746, 1747, 1748, 1780, ...
7: 23, 549, 550, 1544, 1545, 1546, 1547, 1548, 1550, 1551, ...
8: 19, 261, 364, 505, 506, 764, 765, 766, 767, 768, ...
9: 1739, 2175, 2985, 2986, 2987, 2988, 3001, 3002, 3003, 3004, ...
10: 29, 145, 1208, 1219, 1220, 1221, 1222, 1223, 1224, 1225, ...[1] TRUE [1] 292 117 705 183 276 225 143 239 277 165 85 277 56 361 508
[16] 1406 1336 262 100 767 316 90 2505 91Et on peut cartographier ce résultat avec le package mapSf
library(mapsf)
# intitialisation du fond de carte
mf_map(x = gou, border = "white", lwd = 0.5)
# cartographie du nombre de points en cercles proportionnels
mf_map(x = gou,
var = "nb_poi",
type = "prop",
border = "white",
col = "#FF000080",
leg_title = "Nombre de points remarquables",
inches = 0.4, leg_pos = "topright")
# Habillage
mf_layout(title = "Equipements dans les gouvernorats", arrow = TRUE, scale = TRUE, credits = "GeoUnivR 2024 - Tunisie")
Pour plus de détail sur la cartographie avec le package mapsf voir le notebook de N. Lambert conçu pour l’école thématique : Faire des cartes thématiques avec R
Changer de type de géométrie
Il est possible de convertir une géométrie en un autre type, par exemple convertir des géométries de type POLYGON à LINESTRING, avec la fonction st_cast() de sf.
Pour connaitre le type de géométries
[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEIci on converti ces MULTIPOLYGON en MULTILINESTRING
Unités de mesures et calculs
Les mesures sont possibles lorsque l’objet sf à des coordonnées projetées. La première étape est donc de vérifier le SCR de l’objet et son unité de mesure.
Unités de mesures
On peut connaitre l’unité de mesure de la projection avec la fonction st_crs()
st_crs(del)
Coordinate Reference System:
User input: EPSG:2088
wkt:
PROJCRS["Carthage / TM 11 NE",
BASEGEOGCRS["Carthage",
DATUM["Carthage",
ELLIPSOID["Clarke 1880 (IGN)",6378249.2,293.466021293627,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4223]],
CONVERSION["TM 11 NE",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",11,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Oil and gas exploration and production."],
AREA["Tunisia - offshore."],
BBOX[33.22,7.81,38.41,13.67]],
ID["EPSG",2088]]Une méthode plus directe consiste à aller chercher précisément cette information dans ce que renvoie cette fonction :
Calculs
Calcul de superficie
Pour calculer la superficie on utilise la fonction st_area(). Elle renvoie une valeur dans l’unité de mesure de l’objet sf. Ici en metres carrés.
Units: [m^2]
[1] 715991519 462714091 44463342 805121622 238700382On peut enregistrer cette information dans l’objet
Calcul de longueur
Le calcul de longueur ne s’applique qu’aux types LINESTRING et MULTILINESTRING.
Calcul de distance
On peut calculer la distance entre deux points avec la fonction st_distance().
Dans le cas d’objets géographiques de type POLYGONS, le calcul de distance s’effectue automatiquement entre leurs centroides.
Le résultat du calcul est une matrice de distance entre tous les points.
Units: [m]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 244010.5 109399.42 91764.93 229873.3
[2,] 244010.51 0.0 167738.93 116685.25 0.0
[3,] 109399.42 167738.9 0.00 62508.93 150025.5
[4,] 91764.93 116685.2 62508.93 0.00 106360.3
[5,] 229873.27 0.0 150025.46 106360.32 0.0Ici l’unité de mesure de la distance est le metre. On peut modifier cette unité grace au package units et de la fonction set_units(). Il ne s’agit pas ici de modifier l’unité de tout l’objet sf mais seulement des objets créés lors des calculs.
Par exemple en reprenant notre calcul de distances en metres dans un nouvel objet :
On peut les convertir en kilometres (km)
# Chargement du package
library(units)
# Modification de l'unité
set_units(x = distances, value = km)Units: [km]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00000 244.0105 109.39942 91.76493 229.8733
[2,] 244.01051 0.0000 167.73893 116.68525 0.0000
[3,] 109.39942 167.7389 0.00000 62.50893 150.0255
[4,] 91.76493 116.6852 62.50893 0.00000 106.3603
[5,] 229.87327 0.0000 150.02546 106.36032 0.0000Pour que le résultat soit conservé :
Aller plus loin
Simplifier les géométries
La fonction st_simplify() de sfpermet de généraliser des géométries.
Cett fonction n’est pas la plus éfficiace, on peut préférer utiliser la fonction ms_simplify() du package rmapshaper permet de généraliser ou simplifier les géométries en préservant la topologie.
On peut choisir la proportion de sommets à garder avec l’argument keep = ..., et forcer la conservation des formes avec keep_shapes = TRUE
Voici une comparaison de généralisation avec différents paramètres et avec la fonction st_simplify() de sf
par(mfrow = c(1,4),
mar = c(0, 1, 3, 1),
xaxs='i', yaxs='i',
bg = "#F1F3F5")
plot(del$geom, col = "#5B89A3", border = "white")
title("Géométries \ninitiales")
plot(del_simp_sf$geom, col = "#5B89A3", border = "white")
title("Simplification avec sf")
plot(del_simp_rmap$geom, col = "#5B89A3", border = "white")
title("Simplification avec \nrMapshaper")
plot(del_simp_rmap2$geom, col = "#5B89A3", border = "white")
title("Forte simplification \navec rMapshaper")
Agréger des polygones en fonction d’une variable
Digitalisation
La digitalisation est une étape utile de la manipulation de données spatiales mais n’est pas reproductible.
Certains packages de R permettent de réaliser ces opérations mais ne sont pas les plus adaptés car certains de Digitalisation : proposer des choses et préciser que ce n’est pas reproductible et quil peut y avoir des problèmes de topologie et suilhy des outils plus adaptés à ça. Qgis
Construction d’une grille régulière
La fonction st_male_grid() permet la création d’une grille régulière sur l’emprise d’un objet géographique donné.
Cette fonction renvoi un objet de type sfc constitué de listes de cellules. Pour le manipuler facilement on le converti en objet sf avec la fonction st_sf() et en ajoutant un champ d’identifiants.
# Création de la grille
grid <- st_make_grid(gou, cellsize = 35000)
# Ajout d'un identifiant unique et passage en sf
grid <- st_sf(ID = 1:length(grid), geom = grid)
head(grid)Simple feature collection with 6 features and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 178843.6 ymin: 3345646 xmax: 388843.6 ymax: 3380646
Projected CRS: Carthage / TM 11 NE
ID geom
1 1 POLYGON ((178843.6 3345646,...
2 2 POLYGON ((213843.6 3345646,...
3 3 POLYGON ((248843.6 3345646,...
4 4 POLYGON ((283843.6 3345646,...
5 5 POLYGON ((318843.6 3345646,...
6 6 POLYGON ((353843.6 3345646,...par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
plot(st_geometry(gou), col = "#5B89A3", border = "white", lwd = 1)
plot(st_geometry(grid), col = NA, border = "black", lwd = 1, add = TRUE)
Il est possible de créer des grilles hexagonales avec l’argument square = FALSE
grid_hex <- st_make_grid(gou, cellsize = 35000, square = FALSE)
# Ajout d'un identifiant unique et passage en sf
grid_hex <- st_sf(ID = 1:length(grid_hex), geom = grid_hex)
# Cartographie
par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
plot(st_geometry(gou), col = "#5B89A3", border = "white", lwd = 1)
plot(st_geometry(grid_hex), col = NA, border = "black", lwd = 1, add = TRUE)
Ou de récuperer le centroide de ces polygones avec l’argument what = centers ou les angles avec what = corners
par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i', bg = "#F1F3F5")
plot(st_geometry(gou), col = "#5B89A3", border = "white", lwd = 1)
# Les centres
plot(st_make_grid(gou, cellsize = 35000, what = "centers"), col = "red", pch = 20, add = TRUE)
# Les angles
plot(st_make_grid(gou, cellsize = 35000, what = "corners"), col = "pink", pch = 3, add = TRUE)
Intersecter la grille avec les points
Comme présenté dans la Section 7.4 on peut intersecter des points dans des polygones et les compter.
# Intersection
inter <- st_intersects(grid, poi, sparse = TRUE)
# vérifier l'intersection
length(inter) == nrow(grid)[1] TRUESimple feature collection with 6 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 178843.6 ymin: 3345646 xmax: 388843.6 ymax: 3380646
Projected CRS: Carthage / TM 11 NE
ID geom nb_poi
1 1 POLYGON ((178843.6 3345646,... 0
2 2 POLYGON ((213843.6 3345646,... 0
3 3 POLYGON ((248843.6 3345646,... 0
4 4 POLYGON ((283843.6 3345646,... 0
5 5 POLYGON ((318843.6 3345646,... 0
6 6 POLYGON ((353843.6 3345646,... 0On peut affiner cette grille en ne sélectionnant que les carreaux qui intersectent le fond de carte…
…et cartographier le résultat avec mapsf
# intitialisation du fond de carte
mf_map(x = grid_f, border = "white", lwd = 0.5)
# cartographie du nombre de points en cercles proportionnels
mf_map(x = grid_f,
var = "nb_poi",
type = "prop",
border = "white",
col = "#FF000080",
leg_title = "Nombre de points remarquables",
inches = 0.4, leg_pos = "topright")
Conversion vecteur –> raster
Le package terra permet la manipulation de données raster mais aussi de données vecteur pour certains traitements.
On peut convertir un objet vectoriel sf vers un objet vectoriel terra de format spatVector avec la fonction vect() de terra
[1] "SpatVector"
attr(,"package")
[1] "terra"On peut aussi convertir un objet vectoriel sf vers un objet raster terra de format spatRast pour cela voir le module de raster
Ce support a été créé pour la semaine de formation franco-tunisienne GEO UNIV’R Tunisie 2024 - “Enseigner la statistique, la cartographie et l’analyse spatiale avec R qui se tient à Sousse en mai 2024.
Références
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS 15.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.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: Europe/Paris
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] terra_1.8-29 rmapshaper_0.5.0 units_0.8-5 mapsf_0.12.0
[5] dplyr_1.1.4 sf_1.0-19
loaded via a namespace (and not attached):
[1] jsonlite_1.9.1 compiler_4.4.1 tidyselect_1.2.1 Rcpp_1.0.13-1
[5] yaml_2.3.10 fastmap_1.2.0 lattice_0.22-6 R6_2.6.1
[9] generics_0.1.3 curl_6.2.1 classInt_0.4-10 s2_1.1.7
[13] knitr_1.49 htmlwidgets_1.6.4 maplegend_0.1.0 tibble_3.2.1
[17] DBI_1.2.3 pillar_1.10.1 rlang_1.1.5 sp_2.1-4
[21] utf8_1.2.4 V8_6.0.0 xfun_0.49 cli_3.6.4
[25] magrittr_2.0.3 class_7.3-22 wk_0.9.4 digest_0.6.37
[29] grid_4.4.1 rstudioapi_0.17.1 geojsonsf_2.0.3 lifecycle_1.0.4
[33] vctrs_0.6.5 KernSmooth_2.23-24 proxy_0.4-27 evaluate_1.0.1
[37] glue_1.8.0 codetools_0.2-20 e1071_1.7-16 rmarkdown_2.29
[41] tools_4.4.1 pkgconfig_2.0.3 htmltools_0.5.8.1