class: center, middle, inverse, title-slide # Le package sp ### CEFE-CNRS ### 20 mars 2019 --- # 1. Présentation du package sp --- ## 1.1. sp: "Classes et Méthodes pour données spatiales" Le package **sp** propose des structures de données pour gérer des données spatiales dans R, un peu à la manière de ce que l'on trouve dans les SIG - des données ponctuelles, linéaires, surfaciques (= données vectorielles) - des données matricielles (= données raster) C'est le principal package pour la gestion de données vectorielles ! C'est pourquoi beaucoup d'autres packages dépendent de sp (exemples: rgdal, adehabitat, biomod2, maptools). Pour les données raster, on utilise plus fréquemment le package ... raster (dont on parlera à la séance 2). --- ## 1.2. Ce que contient et ne contient pas sp Les classes proposées par sp permettent de gérer la **partie géométrique** et la **partie attributaire** des données spatiales, ainsi que des métadonnées essentielles telles que le **système de coordonnées** employé dans les données. Il est possible de créer de toutes pièces des entités géographiques avec du code R ! En revanche, sp propose peu de fonctions pour **analyser** les données spatiales : ce sont d'autres packages dépendant de sp qui offrent cela. Les entrées/sorties depuis ou vers des fichiers SIG (exemple: .shp, .kml, .gpx, .tif, .asc) ne sont gérées par sp. On utilise d'autres packages tel que rgdal ou maptools. --- ## 1.3. Interactions sp <--> raster Le package raster est compatible avec le package sp, notamment lorsqu'il faut croiser des données vectorielles et raster. Exemples : - extraire la valeur d'un raster pour N points (`raster::extract`) - calculer la moyenne, l'écart type, etc. des pixels recouverts par un polygone (`raster::zonal`) --- ## 1.4. sp, rgeos, sf Le package **rgeos** a longtemps été, en complément de sp, le seul moyen de faire certains calculs tels que : _zones tampons, intersection, calcul de longueur et de surface_. La logique de rgeos est difficile à appréhender, son efficacité n'est pas optimale. Heureusement le package **sf** est arrivé ! Il est beaucoup plus rapide et simple d'utilisation. Mais il introduit un nouveau modèle de données différent de sf. --- ## 1.5. Conclusion - Pour manipuler des données raster, les commandes de base sont dans le package **raster** - Pour manipuler des données vectorielles, le package **sf** est simple et efficace - Le package **sp** reste utilisé en 2019 car ses classes de données spatiales se retrouvent dans d'autres packages tierces ... --- ## 1.6. Livres, sites - Bivand R., Pebesma E., Gomez-Rubio V. Applied Spatial Data Analysis 2nd ed. (bibliothèque du CEFE, et voir aussi ASDAR book (<http://www.asdar-book.org/>) - Lovelace R., Nowosad J., Muenchow J. Geocomputation with R (<https://geocompr.robinlovelace.net/>) --- # 2. TP : les classes de données spatiales dans sp Chargement des packages pour le TP 2 ```r library(sp) ``` --- ## Répertoire de travail Nous allons créer un répertoire **data1** dans lequel seront téléchargées les données nécessaires au TP ```r # afficher le répertoire de destination getwd() ``` ``` ## [1] "D:/GitRepo/RSpatial/formation_2018" ``` ```r # créer un sous répertoire path_data <- "data1" dir.create(path_data) ``` --- ## 2.1. Les classes Spatial* - `SpatialPoints` : pour les données spatiales ponctuelles - `SpatialLines` : pour les données spatiales linéaires - `SpatialPolygons` : pour les données spatiales surfaciques - `SpatialPixels` et `SpatialGrid` : pour les données spatiales matricielles ou raster (sur une grille) --- ## 2.2. Créer un objet SpatialPoints à partir d'une matrice de coordonnées GPS Définissons un `data.frame` avec 5 lignes et 4 colonnes. Les vecteurs `lon` and `lat` donnent la position GPS des 5 lieux, en degrés décimaux (système de coordonnées WGS84) ```r fid <- c(1, 2, 3, 4, 5, 6, 7, 8) name <- c("Montpellier Saint-Roch", "Saint-Aunès", "Villeneuve-lès-Maguelone", "Vic - Mireval", "Baillargues", "Occitanie", "Sabines", "Mosson") lon <- c(3.88067, 3.9629, 3.84991, 3.79945, 4.00716, 3.84785, 3.86032, 3.81933) lat <- c(43.60474, 43.63544, 43.54382, 43.50077, 43.65332, 43.63523, 43.58374, 43.61668) categ <- c("SNCF", "SNCF", "SNCF", "SNCF", "SNCF", "Hérault Transport", "Hérault Transport", "Hérault Transport") coul <- c("red", "red", "red", "red", "red", "orange", "orange", "orange") df <- data.frame(fid, name, lon, lat, coul, categ) ``` --- ## 2.3. La classe SpatialPoints La classe `SpatialPoints` est une structure de données pour stocker des points : seulement la partie "spatiale", pas la partie "attributs". Pour construire un objet de type `SpatialPoints`, nous avons besoin de : - une matrice à 2 colonnes (avec des coordonnées X Y, ou "longitude latitude") - si possible, un objet CRS généré avec la **définition proj 4** du système de coordonnées. Ici la définition vient du site epsg.io : <http://epsg.io/4326>. ```r matcoords <- as.matrix(df[,c("lon","lat")]) sp_pts <- SpatialPoints(matcoords, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs")) # la definition proj4string avec le code EPSG ID est equivalente à la definition complète sp_pts <- SpatialPoints(matcoords, proj4string = CRS("+init=EPSG:4326")) slotNames(sp_pts) ``` ``` ## [1] "coords" "bbox" "proj4string" ``` --- ## 2.4. Les classes Spatial*DataFrame - `SpatialPointsDataFrame` : pour les données spatiales ponctuelles et leurs attributs - `SpatialLinesDataFrame` : pour les données spatiales linéaires et leurs attributs - `SpatialPolygonsDataFrame` : pour les données spatiales surfaciques et leurs attributs - `SpatialPixelsDataFrame` et `SpatialGridDataFrame` : pour des données spatiales matricielles ou raster accompagnées de leur(s) valeur(s) --- ## 2.5. Créer un SpatialPointDataFrame à partir d'un data.frame Créer un SpatialPointDataFrame à partir d'un data.frame **muni de coordonnées**: utiliser la fonction `coordinates` pour créer l'objet, en désignant les 2 colonnes qui contiennent les coordonnées géographiques. ```r # copier le data.frame sp_pts_df <- df # transformer un data.frame en SpatialPointsDataFrame coordinates(sp_pts_df) <- c("lon","lat") # ceci fonctionne aussi : sp_pts_df <- df coordinates(sp_pts_df) <- ~lon+lat # definir le CRS (optionnel) proj4string(sp_pts_df) <- CRS("+init=EPSG:4326") ``` --- --- ### 2.5.1. Exercice: étudiez la structure d'1 objet SpatialPointDataFrame Exécutez les lignes suivantes et répondez aux questions : - Quels sont les __slots__ (propriétés) de l'objet `sp_pts_df` ? - Quel slot contient la définition du système de coordonnées ? - Quel slot contient les coordonnées des points, et sous quelle forme ? ```r slotNames(sp_pts_df) str(sp_pts_df) ``` --- ### 2.5.2. Exercice: .csv -> SpatialPointDataFrame Créer un objet SpatialPointDataFrame à partir du fichier MMM_MMM_VeloParc.csv (cf. URL dans le code ci-dessous) ```r URL <- "http://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_VeloParc.csv" download.file(URL, "data1/MMM_MMM_VeloParc.csv") sp_exercice <- read.csv("data1/MMM_MMM_VeloParc.csv", encoding = "UTF-8") coordinates(sp_exercice) <- ______________ proj4string(sp_exercice) <- ______________ ``` --- ## 2.6. Relation d'héritage entre les classes Spatial\*, SpatialPoints\* et SpatialPointsDataFrame .pull-left[ Le schema suivant tiré du "ASDAR Book" (<http://www.asdar-book.org/>), p.35 nous montre la composition de la class `SpatialPoints`. La classe S4 `Spatial` est la plus générique. Elle a 2 "slots" : bbox (matrix) et proj4string (CRS) La classe S4 `SpatialPoints` hérite de la classe `Spatial` et l'étend ] .pull-right[  ] --- ```r showClass("SpatialPoints") ``` ``` ## Class "SpatialPoints" [package "sp"] ## ## Slots: ## ## Name: coords bbox proj4string ## Class: matrix matrix CRS ## ## Extends: "Spatial" ## ## Known Subclasses: ## Class "SpatialPointsDataFrame", directly ## Class "SpatialPixels", directly ## Class "SpatialPixelsDataFrame", by class "SpatialPixels", distance 2 ``` --- ```r showClass("Spatial") ``` ``` ## Class "Spatial" [package "sp"] ## ## Slots: ## ## Name: bbox proj4string ## Class: matrix CRS ## ## Known Subclasses: ## Class "SpatialPoints", directly ## Class "SpatialMultiPoints", directly ## Class "SpatialGrid", directly ## Class "SpatialLines", directly ## Class "SpatialPolygons", directly ## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2 ## Class "SpatialPixels", by class "SpatialPoints", distance 2 ## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2 ## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2 ## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2 ## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3 ## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2 ``` --- ## 2.7. SpatialLines et SpatialPolygons .pull-left[ Le schema suivant tiré du "ASDAR Book" (<http://www.asdar-book.org/>), p.40 nous montre la composition de la class `SpatialPolygons` et `SpatialLines` Un objet **`SpatialLines`** est composé d'une **liste de `Lines`**. Un objet **`Lines`** est une **liste of `Line`**. Un objet **`Line`** est une matrice **matrice de coordonnées** à 2 colonnes, qui représente un ensemble ordonné de points. **Polygons** = équivalent d'un *MultiPolygone* dans un Shapefile = *MULTIPOLYGON* en notation WKT. **Lines** = équivalent d'une *MultiPolyligne* dans un Shapefile = *MULTILINESTRING* en notation WKT : <https://en.wikipedia.org/wiki/Well-known_text#Geometric_objects>) ] .pull-right[  ] --- ## 2.8. Créer des objets SpatialLines et SpatialLinesDataFrame avec R. ### 2.8.1. Exemple (1/3) : créer 2 polylignes (ensembles de 7 et 5 points) ```r # cree 2 objets Lines avec ID slot = L1 and L2 x1 <- c(3.8, 3.8, 3.9, 3.9) y1 <- c(43.7, 43.6, 43.6, 43.7) x2 <- c(3.84, 3.84, 3.85, 3.85, 3.86) y2 <- c(43.65, 43.64, 43.64, 43.63, 43.63) mat_1 <- cbind(x1, y1) mat_2 <- cbind(x2, y2) line_1 <- Line(mat_1) line_2 <- Line(mat_2) lines_1 <- Lines(list(line_1), "L1") # L1 est l'identifiant de la polyligne (obligatoire) lines_2 <- Lines(list(line_2), "L2") # L2 est l'identifiant ``` --- ### 2.8.2. Exemple (2/3) : SpatialLines ```r sp_lines <- SpatialLines(list(lines_1, lines_2)) plot(sp_lines, axes=T) ``` <!-- --> ```r #str(sp_lines) ``` --- ### 2.8.3. Exemple (3/3) : SpatialLinesDataFrame Un objet **`SpatialLines`** résulte de la combinaison entre un **`SpatialLines`** et un **`data.frame`**. Utiliser le slot **`ID`** du `SpatialLines` et le nom des lignes (**`row.names`**) du `data.frame` pour les mettre en correspondance. ```r # créons la table attributaire (data.frame, 2 colonnes) NAME=c("RANDOM1", "RANDOM2") LENGTH_M = SpatialLinesLengths(sp_lines, longlat=T) * 1000 df_demo <- data.frame(NAME, LENGTH_M) row.names(df_demo) <- c("L1","L2") sp_lines_df <- SpatialLinesDataFrame(sp_lines, df_demo) proj4string(sp_lines_df) <- CRS("+init=EPSG:4326") # take a look at data sp_lines_df@data ``` ``` ## NAME LENGTH_M ## L1 RANDOM1 30295.598 ## L2 RANDOM2 3836.049 ``` --- # 3. TP: Lecture / écriture shapefile Chargement des packages pour le TP 3 ```r library(rgdal) ``` --- ## 3.2. Téléchargement des données utilisées dans ce TP * Communes de Montpellier Métropole Méditerranée (Multipolygones) * Lignes de tramway (Multilignes) * Arrêts de tramway (Points) ```r # changer le répertoire de destination URL1 <- "http://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_Limites.zip" URL2 <- "http://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_LigneTram.zip" URL3 <- "http://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_ArretsTram.zip" download.file(URL1, "data1/MMM_MMM_Limites.zip") download.file(URL2, "data1/MMM_MMM_LigneTram.zip") download.file(URL3, "data1/MMM_MMM_ArretsTram.zip") unzip("data1/MMM_MMM_Limites.zip", exdir="data1") unzip("data1/MMM_MMM_LigneTram.zip", exdir="data1") unzip("data1/MMM_MMM_ArretsTram.zip", exdir="data1") ``` --- ## 3.3. Lire / écrire des shapefiles 2 packages sont utilisables * **rgdal** avec les fonctions `readOGR` et `writeOGR` * **maptools** avec les fonctions `readShapePoints`, `readShapeLines`, `readShapePoly`, `writeShapePoints`, `writeShapeLines`, `writeShapePoly` (obsolètes) **rgdal** peut aussi lire les formats ESRI File GDB, PostGIS, MapInfo, GRASS ... Remarque : le package **sf** est beaucoup plus performant, avec les fonctions `st_read` et `st_write` --- ### 3.3.1. Paramètres de readOGR / writeOGR pour les shapefiles : - `dsn` = répertoire du .shp - `layer` = nom du .shp **sans l'extension** - `driver` (pour writeOGR) = "ESRI Shapefile" - `p4s` (optionnel) = définition système coordonnées - `overwrite_layer` = FALSE / TRUE - `encoding` ("latin1" or "UTF-8") Remarque : pour la lecture / écriture de données volumineuses, le package sf est nettement plus performant ! --- ## 3.4. Lire un shapefile avec rgdal::readOGR .pull-left[ ```r shp_communes <- readOGR(dsn=path_data, layer="MMM_MMM_Limites") ``` ``` ## OGR data source with driver: ESRI Shapefile ## Source: "D:\GitRepo\RSpatial\formation_2018\data1", layer: "MMM_MMM_Limites" ## with 31 features ## It has 2 fields ``` ```r # shp_communes : données géométriques plot(shp_communes, axes=T, main="Montpellier Méditerranée Métropole", col="beige") ``` ] .pull-right[  ] --- ## 3.4. Lire un shapefile avec rgdal::readOGR ```r # shp_communes : données attributaires head(shp_communes@data) ``` ``` ## codcomm nom ## 0 340090 LE CRES ## 1 340120 JACOU ## 2 340027 BEAULIEU ## 3 340057 CASTELNAU LE LEZ ## 4 340169 MONTFERRIER-SUR-LEZ ## 5 340295 SAUSSAN ``` --- ## 3.4. Lire un shapefile avec rgdal::readOGR ```r shp_tramlignes <- readOGR(dsn=path_data, layer="MMM_MMM_LigneTram") ``` ``` ## OGR data source with driver: ESRI Shapefile ## Source: "D:\GitRepo\RSpatial\formation_2018\data1", layer: "MMM_MMM_LigneTram" ## with 5 features ## It has 2 fields ``` ```r shp_tramarrets <- readOGR(dsn=path_data, layer="MMM_MMM_ArretsTram", encoding = "UTF-8") ``` ``` ## OGR data source with driver: ESRI Shapefile ## Source: "D:\GitRepo\RSpatial\formation_2018\data1", layer: "MMM_MMM_ArretsTram" ## with 180 features ## It has 1 fields ``` --- ### 3.4.1. Exercices Téléchargez le shapefile "Contour des départements français" dans le répertoire data Remarque : pour connaître l'encodage des données, vous pouvez procéder par essai/erreur, ou ouvrir le fichier .cpg ```r # cf https://www.data.gouv.fr/fr/datasets/contours-des-departements-francais-issus-d-openstreetmap/ URL <- "http://osm13.openstreetmap.fr/~cquest/openfla/export/departements-20140306-100m-shp.zip" download.file(URL, "data1/departements-20140306-100m-shp.zip") unzip("data1/departements-20140306-100m-shp.zip", exdir="data1") shp_departements <- readOGR(dsn=__________, layer=__________, encoding=__________) names(shp_departements) ``` --- ## 3.5. Ecrire un shapefile avec rgdal::writeOGR L'objet à écrire doit être un Spatial*DataFrame. ```r writeOGR(sp_pts_df, dsn=path_data, layer="test_points", driver="ESRI Shapefile", overwrite_layer=TRUE) ``` --- ## 3.6. Comment convertir des données spatiales d'un système de coordonnées vers un autre ? Transformer les coordonnées d'un système vers un autre requiert `rgdal`. Le package `rgdal` ne fournit pas seulement les drivers nécessaires pour lire et écrire les formats de fichier raster et vecteur. **Le package `rgdal` est aussi nécessaire pour convertir les coordonnées de données spatiales d'un système vers un autre.** La fonction `rgdal::spTransform` s'applique sur tous les objets de classe `Spatial*` et `Spatial*DataFrame`. --- ### 3.6.1. La fonction spTransform Le système de coordonnées de l'objet **origine** doit bien sûr être défini (slot `proj4string`). Le système de coordonnées de l'objet **destination** est passé en paramètre de `spTransform` ```r # input SpatialPointsDataFrame : check CRS proj4string(sp_pts_df) ``` ``` ## [1] "+init=EPSG:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ``` ```r # transformation vers RGF93 / Lambert93 sp_pts_df_l93 <- spTransform(sp_pts_df, CRS("+init=EPSG:2154")) sp_pts_df_l93@coords ``` ``` ## lon lat ## 1 771123.6 6278806 ## 2 777722.6 6282295 ## 3 768712.4 6272009 ## 4 764681.4 6267181 ## 5 781269.6 6284327 ## 6 768436.7 6282166 ## 7 769505.6 6276454 ## 8 766156.1 6280080 ``` --- ### 3.6.2. Exercice : colonnes X et Y à partir d'un `SpatialPointsDataFrame` Nous venons de convertir les points en Lambert 93 dans *sp_pts_df_l93*. Complétez le code ci-dessous pour ajouter au jeu de données 2 nouvelles colonnes *x_l93* et *y_l93* avec les coordonnées des points en Lambert 93. ```r matrix_coords <- _________ sp_pts_df_l93$x_l93 <- _________ sp_pts_df_l93$y_l93 <- _________ ``` --- ### 3.6.3. Exercice avec spTransform et writeOGR Dans cet exercice, vous convertirez l'objet *sp_lines_df* (`SpatialLinesDataFrame`) en Lambert 93 dans un nouvel objet *sp_lines_df_l93*. Puis vous enregistrerez le résultat dans un nouveau fichier shapefile *test_lines_l93.shp*. ```r # transformation en RGF93 / Lambert93 sp_lines_df_l93 <- spTransform(sp_lines_df, CRS(__________)) writeOGR(__________, dsn=__________, layer=__________, driver=__________, overwrite_layer=TRUE) ``` --- # 4. TP: Visualiser les données sur une carte Chargement des packages pour le TP 4 ```r library(leaflet) library(tmap) ``` --- ## 4.1. plot .pull-left[ Avec `plot`, les systèmes de coordonnées doivent être **homogènes**. Ici tout est en Lambert 93 ! ```r plot(shp_communes, axes=T, main="Montpellier Méditerranée Métropole", col="beige") #plot(sp_lines_df_l93, col="red",add=TRUE) plot(sp_pts_df_l93, col="blue",add=TRUE) ``` Pour en savoir plus sur sp + plot : Bivand, ASDAR Book, chap. 3 ] .pull-right[  ] --- ## 4.2. leaflet Avec le package `leaflet`, visualisez facilement vos données sur des cartes web interactives, avec un fond de carte OpenStreetMap. Leaflet est utilisable avec **Shiny** et **RPubs**. --- ## 4.2.1. leaflet ```r # cliquer sur les marqueurs pour afficher les popups ! leaflet() %>% addTiles() %>% addMarkers(data=sp_pts_df, popup=name) ``` <div id="htmlwidget-a7f11f6914624fc603d4" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-a7f11f6914624fc603d4">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addMarkers","args":[[43.60474,43.63544,43.54382,43.50077,43.65332,43.63523,43.58374,43.61668],[3.88067,3.9629,3.84991,3.79945,4.00716,3.84785,3.86032,3.81933],null,null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},["Montpellier Saint-Roch","Saint-Aunès","Villeneuve-lès-Maguelone","Vic - Mireval","Baillargues","Occitanie","Sabines","Mosson"],null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]}],"limits":{"lat":[43.50077,43.65332],"lng":[3.79945,4.00716]}},"evals":[],"jsHooks":[]}</script> RMarkdown users ! installer webshot et phantomJS pour imprimer les cartes leaflet en PDF https://bookdown.org/yihui/bookdown/html-widgets.html --- ### 4.2.2. leaflet, exemple 2 ```r leaflet() %>% addTiles() %>% addPolylines(data=sp_lines_df, label=paste(NAME,"=",LENGTH_M,"m.")) %>% addCircleMarkers(data=sp_pts_df, color=coul) ``` <div id="htmlwidget-466e25d57ba062f527b7" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-466e25d57ba062f527b7">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addPolylines","args":[[[[{"lng":[3.8,3.8,3.9,3.9],"lat":[43.7,43.6,43.6,43.7]}]],[[{"lng":[3.84,3.84,3.85,3.85,3.86],"lat":[43.65,43.64,43.64,43.63,43.63]}]]],null,null,{"interactive":true,"className":"","stroke":true,"color":"#03F","weight":5,"opacity":0.5,"fill":false,"fillColor":"#03F","fillOpacity":0.2,"smoothFactor":1,"noClip":false},null,null,["RANDOM1 = 30295.5982157777 m.","RANDOM2 = 3836.04939585801 m."],{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]},{"method":"addCircleMarkers","args":[[43.60474,43.63544,43.54382,43.50077,43.65332,43.63523,43.58374,43.61668],[3.88067,3.9629,3.84991,3.79945,4.00716,3.84785,3.86032,3.81933],10,null,null,{"interactive":true,"className":"","stroke":true,"color":["red","red","red","red","red","orange","orange","orange"],"weight":5,"opacity":0.5,"fill":true,"fillColor":["red","red","red","red","red","orange","orange","orange"],"fillOpacity":0.2},null,null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]}],"limits":{"lat":[43.50077,43.7],"lng":[3.79945,4.00716]}},"evals":[],"jsHooks":[]}</script> --- ## 4.3. tmap .pull-left[ ```r tm_shape(sp_pts_df) + tm_symbols(col="categ") ``` ``` ## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3 ``` ] .pull-right[  ] --- ## 4.3. tmap .pull-left[ ```r tm_shape(shp_communes) + tm_polygons() + tm_shape(sp_pts_df) + tm_symbols(size=0.4, col="categ", palette=c("orange", "red"), shape=17) ``` Graph parameter reminder https://www.r-graph-gallery.com/6-graph-parameters-reminder/ http://www.r-graph-gallery.com/42-colors-names/ ] .pull-right[  ] --- ## 4.4. étendue, barre d'échelle, flèche du nord, labels dans tmap ```r tm_shape(shp_communes) + tm_polygons(col="beige") + tm_shape(sp_pts_df) + tm_symbols(size=0.4, col="categ", palette=c("orange", "red"), title.col="Types de gare", shape=17) + tm_text("fid", fontface="bold", ymod=1) + tm_style(style="white") + tm_layout(main.title="Gares routières et SNCF dans la métropole de Montpellier", main.title.position = c("center","top"), main.title.size = 1, legend.title.size = 1, legend.title.fontface = "bold", inner.margins = c(0.1, 0.1, 0.1, 0.1)) + tm_compass(type = "arrow", size=1.5, fontsize=0.6, position = c("right","top")) + tm_scale_bar(size=0.4, breaks=c(0,5,10)) + tm_credits("Données : Montpellier 3M, SNCF", position=c("left","bottom")) ``` --- ## 4.4. étendue, barre d'échelle, flèche du nord, labels dans tmap  --- ## 4.5. Exercice: lignes de tramways Dans cet exercice, vous allez partir du code ci-dessous avec tmap, pour l'améliorer progressivement. 1. Laisser les communes de la métropole en fond de carte, mais zoomer sur l'étendue des 4 lignes. 2. Avec tmap, afficher chacune des 4 lignes de tramway avec une couleur particulière. ```r tm_shape(shp_communes) + tm_polygons() + tm_shape(shp_tramlignes) + tm_lines() + tm_shape(shp_tramarrets) + tm_symbols(size=0.2) + tm_layout(main.title="Lignes de tramway dans la métropole de Montpellier", main.title.position = c("center","top"), main.title.size = 1.) ``` --- # 5. "COMMENT FAIRE POUR ... ?" Chargement des packages pour le TP 5 ```r library(maptools) ``` --- ## 5.1. Comment lire des coordonnées en degrés minutes secondes ? `char2dms` convertit une chaîne en objet DMS. La chaîne doit se terminer par N,S,W ou E. Précisez le caractère utilisé pour les degrés, les minutes, les secondes (par défaut : d,',"). `as.numeric` convertit l'objet DMS en degrés décimaux. ```r villes <- c("Paris", "Madrid", "Buenos-Aires", "New-York") char_lat <- c("48°51' 24\" N", "40°26'00\"N", "34°36'29\"S", "40°42' 52\"N") char_lon <- c("2°21'07\"E", "3°41' W", "58°22' 13\" W", "74°0'22\"W") dd_lat <- as.numeric(char2dms(char_lat, "°", "'", "\"")) dd_lon <- as.numeric(char2dms(char_lon, "°", "'", "\"")) ``` --- ## 5.1. Comment lire des coordonnées en degrés minutes secondes ? ```r leaflet() %>% addTiles() %>% addMarkers(lng=dd_lon, lat=dd_lat) ``` <div id="htmlwidget-57a568007aec20936f10" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-57a568007aec20936f10">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addMarkers","args":[[48.8566666666667,40.4333333333333,-34.6080555555556,40.7144444444444],[2.35194444444444,-3.68333333333333,-58.3702777777778,-74.0061111111111],null,null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},null,null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]}],"limits":{"lat":[-34.6080555555556,48.8566666666667],"lng":[-74.0061111111111,2.35194444444444]}},"evals":[],"jsHooks":[]}</script> --- ## 5.2. Comment calculer une matrice de distance entre des points ? `spDists` calcule les distances entre 2 ensembles de points. Si les coordonnées des points sont dans un système métrique (exemple UTM), alors utiliser longlat=FALSE pour obtenir la distance euclidienne en mètres. Si les coordonnées des points sont en degrés (exemple WGS84), alors utiliser longlat=TRUE. On obtient une distance en km. ```r # distance entre Madrid, Buenos Aires, Paris, New-York (cf exemple précédent) sp_pts_villes <- data.frame(villes, dd_lat, dd_lon, char_lat, char_lon) coordinates(sp_pts_villes) <- c("dd_lon","dd_lat") proj4string(sp_pts_villes) <- CRS("+init=EPSG:4326") ``` --- ```r # matrice de distance : toutes les combinaisons m <- spDists(sp_pts_villes, longlat=TRUE) dimnames(m) <- list (c("PAR","MAD","BAI","NYC"), c("PAR","MAD","BAI","NYC")) m ``` ``` ## PAR MAD BAI NYC ## PAR 0.000 1050.571 11024.886 5852.808 ## MAD 1050.571 0.000 10022.522 5783.334 ## BAI 11024.886 10022.522 0.000 8492.632 ## NYC 5852.808 5783.334 8492.632 0.000 ``` ```r # distances consécutives v <- spDists(sp_pts_villes, longlat=TRUE, segments=TRUE) names(v) <- c("PAR-MAD","MAD-BAI", "BAI-NYC") v ``` ``` ## PAR-MAD MAD-BAI BAI-NYC ## 1050.571 10022.522 8492.632 ``` --- ## 5.3. Comment filtrer des données par attribut ? Pour filtrer les données d'un objet Spatial*DataFrame d'après les valeurs d'une colonne, nous procédons de la même manière d'avec un data.frame. Exemple : à partir de shp_communes (cf 3.4) créez un nouvel objet shp_communes2 (toujours un SpatialPolygonsDataFrame) avec les 4 communes suivantes : COURNONSEC, COURNONTERRAL, LAVERUNE, PIGNAN. ```r communes_match <- shp_communes$nom %in% c("COURNONSEC", "COURNONTERRAL", "LAVERUNE", "PIGNAN") which(communes_match) ``` ``` ## [1] 11 14 17 31 ``` --- ```r shp_communes2 <- shp_communes[communes_match,] plot(shp_communes, col="grey") plot(shp_communes2, col="red", add=T) ``` <!-- --> --- ## 5.4. Comment joindre à un objet Spatial*DataFrame les données d'une autre table ? Dans l'exemple suivant nous joignons aux communes les données d'un tableau avec la population. Avec `match`, nous redéfinissons l'ordre des données du tableau population ; puis nous ajoutons la nouvelle colonne *poptot* à *shp_communes*. Voir aussi la fonction `maptools::spCbind` pour fusionner des données attributaires. *Remarque :* avec `sf` et `dplyr`, les jointures seront plus faciles !! ```r # la population des communes (recensement 2015) dans un data.frame codcom <- c("340022", "340027", "340057", "340058", "340077", "340087", "340088", "340090", "340095", "340116", "340120", "340123", "340129", "340134", "340164", "340169", "340172", "340179", "340198", "340202", "340217", "340227", "340244", "340249", "340256", "340259", "340270", "340295", "340307", "340327", "340337") poptot <- c(7202, 1778, 19886, 6162, 5527, 3262, 5844, 9343, 6860, 8027, 6881, 10783, 16567, 3211, 996, 3833, 282143, 1928, 9177, 6761, 5300, 1754, 2925, 2365, 1884, 5498, 9233, 1518, 2768, 6247, 9744) dfpop <- data.frame(codcom, poptot) ``` --- ## 5.4. Comment joindre à un objet Spatial*DataFrame les données d'une autre table ? ```r # index des lignes dans dfpop qui correspondent à shp_communes, en comparant codcom row_indices <- match(shp_communes$codcomm, dfpop$codcom) # nouvelle colonne poptot dans shp_communes shp_communes$poptot <- dfpop[row_indices,"poptot"] head(shp_communes@data) ``` ``` ## codcomm nom poptot ## 0 340090 LE CRES 9343 ## 1 340120 JACOU 6881 ## 2 340027 BEAULIEU 1778 ## 3 340057 CASTELNAU LE LEZ 19886 ## 4 340169 MONTFERRIER-SUR-LEZ 3833 ## 5 340295 SAUSSAN 1518 ``` --- ## 5.5. Comment trouver dans quel polygone est situé un point ? Dans `sp`, la fonction `over` permet facilement de faire certaines jointures spatiales. Exemple : pour chaque arrêt de tramway, trouver dans quel commune il est situé. Combien y a-t-il d'arrêts dans chaque commune ? ```r commune_tramarrets <- over(shp_tramarrets, shp_communes) nb_stations <- as.data.frame(table(commune_tramarrets$nom)) head(nb_stations) ``` ``` ## Var1 Freq ## 1 BAILLARGUES 0 ## 2 BEAULIEU 0 ## 3 CASTELNAU LE LEZ 18 ## 4 CASTRIES 0 ## 5 CLAPIERS 0 ## 6 COURNONSEC 0 ``` --- ## 5.5.1. Exercice : jointure spatiale + sélection Utilisez le data.frame *commune_tramarrets* obtenu dans l'exercice précédent pour déterminer et afficher les noms des arrêts de tramway situés sur la commune "PEROLS". ```r # indices_perols <- which(__________) __________$nom[indices_perols] ``` --- ## 5.6. Comment enregistrer vos données en KML ? Enregistrer des points sous forme de fichier KML est un autre moyen intéressant pour visualiser des données, notamment avec Google Earth . Le moyen le plus direct est d'utiliser la fonction `maptools::kmlPoints`. Cette fonction permet d'associer à chaque point un nom et l'URL d'une icone. Il est possible d'associer des icones de différentes couleurs suivant une variable qualitative. Exemple : créez un fichier KML avec les gares routières et SNCF (objet sp_pts_df, cf. 2.5). Nous nous baserons sur la colonne *coul* pour les couleurs (red = SNCF, orange = Hérault Transport) URL des icones : https://sites.google.com/site/gmapsdevelopment/ --- ## 5.6. Comment enregistrer vos données en KML ? ```r # url for red & orange markers url_color_markers <- paste0("http://maps.google.com/mapfiles/ms/micons/", sp_pts_df$coul, ".png") kmlPoints(sp_pts_df, kmlfile="data1/exemple_formation.kml", name=sp_pts_df$name, icon=url_color_markers) ``` --- # 6. Révision --- ## 6.1. Exercice récapitulatif Dans cet exercice, nous allons réviser tous les chapitres vus précédemment. Nous utiliserons une fonction qui, à partir d'un point de départ, réalise une marche avec des pas de 50 m dans une direction aléatoire. - Vous partirez d'un point P1 dont les coordonnées, en Lambert 93, sont : *X=766350, Y=6313470*. - A partir de ce point de départ, réalisez 5 marches aléatoires de 100 pas de 50 m. Nous appellerons ces 5 marches : *RW1, RW2, RW3, RW4, RW5*. - Pour la marche *RW1*, mesurez la distance entre le point de départ et les points suivants. Vous utiliserez pour cela la fonction `sp::spDistsN1`. - Créez un objet `SpatialLinesDataFrame` nommé *sp_rw* avec 5 entités Lines, à partir des 5 marches aléatoires. Le système de coordonnées de *sp_rw* sera Lambert 93 (EPSG:2154). Chaque entité aura pour attribut : un identifiant unique et la longueur renvoyée par `SpatialLinesLengths`. - Transformer *sp_rw* dans le système WGS84 ; vous créerez pour cela un nouvel objet `SpatialLinesDataFrame` nommé *sp_rw_wgs84* - Avec leaflet, afficher *sp_rw_wgs84* avec une couleur différente pour chaque trajet, sur un fond de carte OpenStreetMap --- --- ### 6.1.1. A compléter ```r # random walk function library(sp) library(leaflet) f_point_in_circle <- function(x_from, y_from, step_dist=50, step_number=100) { angle <- runif(step_number + 1, 0, 360) x_diff <- step_dist * cos(angle) y_diff <- step_dist * sin(angle) x_pos <- x_from + cumsum(x_diff) y_pos <- y_from + cumsum(y_diff) cbind(x=x_pos, y=y_pos) } # point de depart p1 <- c(x=766350, y=6313470) # appeler 5 fois la fonction f_point_in_circle rw_matlist <- replicate(5, f_point_in_circle(p1["x"], p1["y"]), simplify = F) rw_names <- c("RW1", "RW2", "RW3", "RW4", "RW5") # distance depart - points suivants f_distance_suivants <- function(matcoords) { spDistsN1(matcoords, matcoords[1,], longlat=F) } f_distance_suivants(rw_matlist[[5]]) # creer SpatialLines = list de Lines list_rw_lines <- list() for (i in 1:5) { rw_coords <- rw_matlist[[i]] rw_name <- rw_names[i] rw_line <- Line(__________) rw_lines <- Lines(__________, rw_name) list_rw_lines[[i]] <- rw_lines } sp_l <- SpatialLines(__________, proj4string = __________) df_attr <- data.frame(rw_names, long_m = SpatialLinesLengths(sp_l, longlat=F), row.names=__________) sp_rw <- SpatialLinesDataFrame(__________, __________) # transformation Lambert 93 -> WGS 84 sp_rw_wgs84 <- spTransform(sp_rw, __________) coul5 <- c("cyan", "deeppink" ,"brown", "darkorange", "chartreuse") # carto leaflet avec une couleur differente pour chaque ligne leaflet(data=sp_rw_wgs84) %>% addTiles() %>% addPolylines(_____ = coul5) ``` --- # Conclusion de la séance 1 --- ## Ce que nous avons vu aujourd'hui - les structures de données spatiales introduites dans sp : SpatialPoints, SpatialLines, SpatialPolygons, SpatialPointsDataFrame, SpatialLinesDataFrame, SpatialPolygonsDataFrame - comment lire et écrire des fichiers SIG vectoriel avec rgdal - gérer le système de coordonnées des données - présenter les données sur une carte simple --- ## Ce que nous n'avons pas abordé aujourd'hui ### ... mais que nous verrons avec le package sf - calcul de surface - jointures attributaires et spatiales - zones tampons - intersecter des données vectorielles (lignes et polygones, polygones et polygones) - requêtes spatiales complexes - cartes choroplèthes