This book is in Open Review. I want your feedback to make the book better for you and other readers. To add your annotation, select some text and then click the on the pop-up menu. To see the annotations of others, click the in the upper right hand corner of the page

Cours 5 Spatial data cours

Contient pas mal d’infos : https://r-spatial.org/book/07-Introsf.html

Géomatique avec R : https://rcarto.github.io/geomatique_avec_r/

library(tidyverse)
library(sf)

5.1 Tips

  • Le CRS est souvent une source de bug : absent, pas le bon, contenant des accents

5.2 Mapper des données

Mapper des points dont on a la longitude et la latitude :

ggplot(mapping = aes(x = lon, y = lat)) +
  geom_point()

ggmap : ajoute des couches d’images (fond de carte) téléchargées depuis des services de cartorgaphie en ligne.

library("ggmap")

east_canada <- get_stamenmap(bbox = c(left=-81, right = -59, bottom = 44, top = 51),#boîte de coordonnées délimitant la carte à produire
                             zoom = 6, #zoom = niv de détails (2 suffisant pour une carte du monde)
                             maptype = "terrain") # type de carte

ggmap(east_canada) + 
geom_point(data = weather, mapping = aes(x = lon, y = lat))

## ou 

ggmap(east_canada,
      base_layer = ggplot(weather, aes(x = lon, y = lat))) + 
  geom_point()
##base_layer permet d’effectuer des facettes et d’éviter de spécifier la source des données dans toutes les couches subséquentes.

Maptype : - “terrain” : efficace mais peu esthétique - “toner-lite” : pour l’impression - “watercolor” : pour le web

5.2.1 Types de données spatiales

Les données spatiales sont des données localisées dans un système de coordonnées de référence (CRS)

  • Données vectorielles (en général en 2D, 3D si on prend en compte l’altitude:

    • Données ponctuelles : points
    • Données linéaires : série de points (route, rivière)
    • Données de polygone : aire délimitée par des points (champ, bassin versant)

  • Données raster: grille (image satellite où chaque pixel est associé à un recouvrement foliaire.)

  • Données génériques : shapefiles et geojson

  • Données spécialement conçus pour R : sf

5.2.2 ggplot

geom_polygon : : pour créer des polygones coord_map geom_path() : pour créer des lignes geom_tile() : visualiser une grille -> graphique de type “heatmap” *geom_sf() : pour afficher les objets sf

Cartes intéractives (mode leaflet) : tmap_mode(“view”) Pour revenir en mode statique : tmap_mode(“plot”)

5.2.3 ggspatial

  • ajouter une échelle : annotation_scale(location = “br”) +
  • ajouter une boussole : annotation_north_arrow(pad_y = unit(1, “cm”),style = north_arrow_nautical())

5.2.4 Les rasters

= données (variable(s)) associées à une grille comprenant les combinaisons de longitudes et latitudes.

Raster = * une en-tête : - syst de coordonnées de référence - l’origine : généralement les coordonnées du coin bas-gauche de la matrice (mais le package ‘raster’ prends le coin haut-gauche par défaut) - dimensions : nbr de colonnes, de lignes, la résolution de la taille des cellules * une matrice de cellules équidistantes (~pixels) : la matrice ne stocke qu’1 coordonnée de la cellule : l’origine -> + efficace et rapide que le traitement des données vectorielles.

1 cellule = 1 ou plsrs valeurs, ou NA (numérique ou catégorique) = la valeur moyenne (ou majoritaire) de la zone qu’elle couvre. Cependant, dans certains cas, les valeurs sont en fait des estimations pour le centre de la cellule.

l’information spatiale est implicitement donnée par l’étendue spatiale et le nombre de lignes et de colonnes dans lesquelles la zone est divisée.

Les données raster sont en général continues (altitude, T°, densité de pop), mais peuvent-être aussi catégoriques (type de sol) mais dans ce cas une représentation vectorielle pourrait être plus appropriée.

expand.grid() : créer une grille de paramêtres : associer des coordonées à une variable

grid <- expand.grid(lon = seq(from = -80, to = -60, by = 0.25),
                    lat = seq(from = 45, to = 50, by = 0.25))
grid <- grid %>% 
  mutate(z = 10*sin(lon*lat) - 0.01*lon^2 + 0.05*lat^2) # créer une variable spatialisée

grid %>% head()

ggplot(grid, aes(lon, lat)) +
  geom_tile(aes(fill = z))

Géostatistique = étude statistique des variables spatiales

5.2.4.1 Geometry

Intersection between a raster and a geometry (sf) : mask(raster, geometry)

5.2.5 Les objets spatialisés en R

5.2.5.1 Données vectorielles

Objets géoréférencés R (=variables (data.frame) liés à des coordonées (sfc)): sf Packages : sp, sf (mieux adapté au tidy)

Transformer un dataframe en objet sf : st_as_sf(, coords = c(‘X’,‘Y’,‘Z’)) : - type de géométrie (geometry type: POINT (= sfg objets)), - les limites des objets (bbox: …), - le système de référence (epsg ou proj4string: …) - le tableau descriptif

Revenir à une table non spatialisée : st_drop_geometry()

Les objets sfg peuvent être créés avec : * un vecteurs numérique * une matrice * une liste

plot(st_point(c(5, 2)))           # XY point (2D)
st_point(c(5, 2, 3))              # XYZ point (3D)
st_point(c(5, 2, 1), dim = "XYM") # XYM point
st_point(c(5, 2, 3, 1))           # XYZM point

## the 'rbind' function simplifies the creation of matrices
### MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
plot(st_multipoint(multipoint_matrix))

### LINESTRING
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
plot(st_linestring(linestring_matrix))

### POLYGON
polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
plot(st_polygon(polygon_list))

### POLYGON with a hole
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
plot(st_polygon(polygon_with_hole_list))

### MULTILINESTRING
multilinestring_list = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
plot(st_multilinestring((multilinestring_list)))

### MULTIPOLYGON
multipolygon_list = list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
                         list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
plot(st_multipolygon(multipolygon_list))

### GEOMETRYCOLLECTION
gemetrycollection_list = list(st_multipoint(multipoint_matrix),
                              st_linestring(linestring_matrix))
plot(st_geometrycollection(gemetrycollection_list))
##> GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2),
##>   LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))
5.2.5.1.1 Simple feature columns (sfc)

une sfc est une liste d’objets sfg, de meme type géométrique ou non.

point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2) #st_sfc() crée une liste d'objets sfg
points_sfc
st_geometry_type(points_sfc)
  • sfc to sfg : st_point(as.numeric(unlist(sfc))) (ex ici point)

st_read() : pour charger des données (ex:shapefiles) au format sf

Les tableaux sf sont manipulables sous tydiverse.

plot() : par défaut crée un multi-panel plot : un graph par variable.

Pour ne ploter que les contours (la géométrie):

quebec %>%
  st_geometry() %>%
  plot() # ou bien plot(st_geometry(quebec)), ou bien plot(quebec %>% select(geometry))

  • Jointures spatiales (intersections de polygones, unions, différences) -> 1 nouvelle géométrie

  • Attribuer un crs à un objt sf ou sfc : st_set_crs(objtsanscrs, st_crs(objtaveccrs))

  • st_crs<- : replacing crs does not reproject data; use st_transform for that

  • Exporter un tableau sf en format csv incluant la géométrie : st_write(obj = tableau,dsn = “tableau.csv”, layer_options = “GEOMETRY=AS_XY”) Si la géométrie n’est pas consituée de points, il faudra préalablement transformer les polygones en points avec st_cast()

5.2.5.2 Données raster

Package : raster grilles svt enchasées dans des images tif géoréférencées.

1 image tif = 1 ou plsrs bandes (variables)

raster() : importe des données raster à 1 bande (“names” dans les infos) brick() : importe des données raster à plsrs bandes (“names” dans les infos) stack() : permet de connecter plusieurs objets raster stockés dans différents fichiers ou plusieurs objets en mémoire

Les informations des objets RasterLayer et RasterBrick peuvent être extraites par les fonctions : extent(), ncell(), nlayers() et crs(). La fonction plot() permet d’explorer les données en créant 1 graphique par bande.

Les tif sont trés volumineux. Si l’on n’as pas besoin d’une si haute résolution, on peut simplifier le raster avec raster::aggregate()

  • Attribuer un crs à raster : crs(mnt) <- crs(obj2) ou crs(mnt) <- crs(“EPSG:2972”)

5.2.5.3 Travailler sur du vectoriel & du raster

Le package “raster” ne supporte pas le format sf, il faut donc préalablement convertir en sp : poly_sp <- as(poly, “Spatial”) (poly étant ici un polygone)

5.2.5.4 Opérations sur Raster

Le package terra : https://rcarto.github.io/geomatique_avec_r/les-donn%C3%A9es-raster-le-package-terra.html#

  • Intersection (ce qu’il y a en commun) entre un polygone et un raster : mask()
  • Découper (rectangulaire selon les limites de l’objet) : crop()
  • Découper un raster selon un autre raster, ou mettre à la meme résolution :resample()
  • Extraction : extract()
##Pour effectuer un calcul sur l’intérieur du polygone avec extract()… on spécifie le raster, le polygone et la fonction!

extract(canopy, poly_sp, fun = mean) # "canopy" =  raster, "poly_sp" = polygone
  • Convert values to classes : cut() :
# Reclasser en intervalles avec `cut()`
MNT <- stars::st_as_stars(MNT) # first transform in stars
MNT <- raster::cut(MNT, 
                   breaks = c(0, 6, 23.67), 
                   labels = c("0-6", "6-23.67"),
                   include.lowest = TRUE)
library(raster)
shpP16 <- rgdal::readORG(dsn = "C:/Users/BADOUARD/Desktop/Vincyane/Gaps2019/Gaps2019", layer = "Gaps2019")  
shpP16 <- raster::shapefile("C:/Users/BADOUARD/Desktop/Vincyane/Gaps2019/Gaps2019.shp")

library(sf)
shpP16 <- st_read("C:/Users/BADOUARD/Desktop/Vincyane/Gaps2019/Gaps2019.shp")

plot(shpP16)
plot(shpP16, col="#f2f2f2", bg="skyblue", lwd=0.25, border=0 )

library(ggplot2)
ggplot() + 
  geom_sf(data = shpP16, size = 3, color = "black", fill = "red") + 
  ggtitle("2019 - P16 gaps") + 
  coord_sf()

5.2.5.5 Plot raster with ggplot

# First transform raster in data.frame
TWI_df <- as.data.frame(TWI, xy = TRUE)%>%
  na.omit()

ggplot() + 
  geom_raster(data = TWI_df, aes(x = x, y = y, fill = `P13_2022_TWI`)) + # with geom_raster()
  scale_fill_viridis_c() +
  ggtitle("Paracou 2019 P13 - TWI") + 
  coord_sf()

5.2.6 Shapefiles

  • lire un shapefile : st_read() Contiennent plsrs fichiers : fichier .prj : informations du système de coordonnées

5.2.7 Systèmes de coordonnées de références (CRS)

  • Le vérifier : crs()
  • L’attribuer : obj <- st_set_crs(obj, 2972)
  • en UTM c’est en mètres

5.2.7.1 Système de coordonnées géographiques

Valeurs : Longitude/latitude

La Terre est représenté sphérique ou ellipsoïde (+juste).

Modele ellipsoïde, 2 paramètres: - rayon équatorial - rayon polaire

5.2.7.2 Système de référence de coordonnées pojetées

basé sur des coordonées catésiennes sur une surface plane. et basés sur les systèmes de coordonnées géographiques

  • une origine
  • axes x & y
  • une unité (ex: le mètre)

5.3 Functions

chercher là dedans : https://r-spatial.org/book/07-Introsf.html

Retourner des points interpolés selon une interpolation linéaire entre 2 points : approx(x= coord, y = autrecoord, xout = coord selon laquelle interpoler)

5.3.1 sf package

  • Passer d’un sf à un dataframe (supprimer la col géométrie) : st_geometry(sf) <- NULL
# from sf to df with original coordinates
XY <- st_coordinates(datasf) # stock coordinates
st_geometry(datasf) <- NULL # no more sf object
data <- cbind(datasf, XY) # bind XY coord

  • sfg to sfc : st_sfc(sfg)

  • sfc to sfg : st_point(as.numeric(unlist(sfc))) (ex ici point)

  • MULTIPOINT to POINT : st_cast(st_sfc(mutipoint), “POINT”, group_or_split = TRUE)

  • Simplification d’objet : st_cast(objet, “objt + simple” (ex : polygon to linestring)

  • Attribuer un crs à un objt sf ou sfc : st_transform(objtsanscrs, st_crs(objtaveccrs))

  • Un buffer autour d’un objet spatial: buffer <- st_buffer(pol, dist=100, endCapStyle = “SQUARE”) %>% st_union() # the square doesn’t works

  • corners of a square : st_coordinates(square)[,1:2]

  • st_union ne préserve pas l’ordre des géométries dans la géométrie finale, il faut utiliser :

objt_sf_12 <- objt_sf_1 %>% rbind(objt_sf_2)  

objt_sf_12 <- do.call(c, st_geometry(objt_sf_12))
  • Create a regular grid :
# Subplots (1ha) grid for 25 ha
grid <- data.frame(expand.grid(
  x = seq(0,500, by= 100),
  y = seq(0,500, by= 100)
)
)
plot(grid)

5.3.2 Stars package to work on rasters

voir : https://r-spatial.org/book/07-Introsf.html ### en 3D + package : plot3D for 3-D arrows, segments, polygons, boxes, rectangles 3D fixe plot3D::scatter3D()

  • rgl for visualisation in 3D 3D pivotable : plotly::plot_ly(data, x = ~x, y = ~y, z = ~z, color=var) rgl::plot3d()