05-geometry-operations.Rmd

# Opèrations géométriques {#geometry-operations}

## Prérequis {-}

- Ce chapitre utilise les mêmes paquets que le chapitre \@ref(spatial-operations) mais avec l'ajout de **spDataLarge**, qui a été installé dans le chapitre \@ref(spatial-class):

```{r 05-geometry-operations-1, message=FALSE}
library(sf)
library(terra)
library(dplyr)
library(spData)
library(spDataLarge)
```

## Introduction

Jusqu'à présent, ce livre a abordé la structure des jeux de données géographiques (chapitre \@ref(spatial-class)), et la manière de les manipuler en fonction de leurs attributs non géographiques (chapitre \@ref(attr)) et de leurs relations spatiales (chapitre \@ref(spatial-operations)).
Ce chapitre se concentre sur la manipulation des éléments géographiques des objets géographiques, par exemple en simplifiant et en convertissant les géométries vectorielles, en recadrant les rasters et en convertissant les objets vectoriels en rasters et les rasters en vecteurs.
Après l'avoir lu --- et avoir fait les exercices à la fin --- vous devriez comprendre et contrôler la colonne géométrique des objets `sf` ainsi que l'étendue et l'emplacement géographique des pixels représentés dans les rasters par rapport à d'autres objets géographiques.

La section \@ref(geo-vec) couvre la transformation des géométries vectorielles avec des opérations "unaires" (ou fonction avec un argument)  et "binaires" (fonction avec plus d'un argument).
Les opérations unaires portent sur une seule géométrie de manière isolée, notamment la simplification (de lignes et de polygones), la création de tampons et de centroïdes, et le déplacement/la mise à l'échelle/la rotation de géométries uniques à l'aide de " transformations affines " (sections \@ref(simplification) à \@ref(transformations affines)).
Les transformations binaires modifient une géométrie en fonction de la forme d'une autre, y compris l'écrêtage et les unions géométriques (\index{vector!union}), traités respectivement dans les sections \@ref(écrêtage) et \@ref(unions géométriques).
Les transformations de type (d'un polygone à une ligne, par exemple) sont présentées dans la section \@ref(type-trans).

La section \@ref(geo-ras) couvre les transformations géométriques sur les objets rasters.
Il s'agit de modifier la taille et le nombre des pixels, et de leur attribuer de nouvelles valeurs.
Elle enseigne comment modifier la résolution (également appelée agrégation et désagrégation), l'étendue et l'origine d'un objet matriciel.
Ces opérations sont particulièrement utiles si l'on souhaite aligner des rasters provenant de sources diverses.
Les objets rasters alignés partagent une correspondance biunivoque entre les pixels, ce qui permet de les traiter à l'aide d'opérations d'algèbre raster, décrites dans la section \@ref(map-algebra). 
L'interaction entre les objets raster et vectoriels est traitée au chapitre  \@ref(raster-vector). 
Elle montre comment les valeurs matricielles peuvent être "masquées" et "extraites" par des géométries vectorielles.
Il est important de noter qu'elle montre comment " polygoniser " les données raster et " rastériser " les veceurs, ce qui rend les deux modèles de données plus interchangeables.

## Opérations géométriques sur les données vectorielles {#geo-vec}

Cette section traite des opérations qui, d'une manière ou d'une autre, modifient la géométrie des objets vectoriels (`sf`).
Elle est plus avancée que les opérations sur les données spatiales présentées dans le chapitre précédent (dans la section \@ref(spatial-vec)), parce qu'ici nous allons plus loin dans la géométrie :
les fonctions présentées dans cette section fonctionnent sur les objets de la classe `sfc` en plus des objets de la classe `sf`.

### Simplification

\index{vector!simplification} 
La simplification est un processus de généralisation des objets vectoriels (lignes et polygones) généralement destiné à être utilisé dans des cartes à plus petite échelle.
Une autre raison de simplifier les objets est de réduire la quantité de mémoire, d'espace disque et de bande passante qu'ils consomment :
il peut être judicieux de simplifier des géométries complexes avant de les publier sous forme de cartes interactives. 
Le paquet **sf** fournit `st_simplify()`, qui utilise l'implémentation GEOS de l'algorithme de Douglas-Peucker pour réduire le nombre de sommets.
`st_simplify()` utilise la `dTolerance` pour contrôler le niveau de généralisation des unités de la carte [voir @douglas_algorithms_1973 pour plus de détails].
La figure \@ref(fig:seine-simp) illustre la simplification d'une géométrie `LINESTRING` représentant la Seine et ses affluents.
La géométrie simplifiée a été créée par la commande suivante :

```{r 05-geometry-operations-2}
seine_simp = st_simplify(seine, dTolerance = 2000)  # 2000 m
```

```{r seine-simp, echo=FALSE, fig.cap="Comparaison de la géométrie originale et simplifiée de la Seine.", warning=FALSE, fig.scap="Simplification en action.", message=FALSE}
library(tmap)
p_simp1 = tm_shape(seine) + tm_lines() +
  tm_layout(main.title = "Donnée source")
p_simp2 = tm_shape(seine_simp) + tm_lines() +
  tm_layout(main.title = "st_simplify")
tmap_arrange(p_simp1, p_simp2, ncol = 2)
```

L'objet `seine_simp` résultant est une copie de l'objet original `seine` mais avec moins de vertices.
Le résultat étant visuellement plus simple (Figure \@ref(fig:seine-simp), à droite) et consommant moins de mémoire que l'objet original, comme vérifié ci-dessous :

```{r 05-geometry-operations-3}
object.size(seine)
object.size(seine_simp)
```

La simplification est également applicable aux polygones.
Ceci est illustré par l'utilisation de `us_states`, représentant les États-Unis contigus.
Comme nous le montrons dans le chapitre \@ref(reproj-geo-data), GEOS suppose que les données sont dans un CRS projeté et cela pourrait conduire à des résultats inattendus lors de l'utilisation d'un CRS géographique.
Par conséquent, la première étape consiste à projeter les données dans un CRS projeté adéquat, tel que le US National Atlas Equal Area (epsg = 2163) (à gauche sur la figure \@ref(fig:us-simp)) :

```{r 05-geometry-operations-4}
us_states2163 = st_transform(us_states, "EPSG:2163")
us_states2163 = us_states2163
```

`st_simplify()` works equally well with projected polygons:

```{r 05-geometry-operations-5}
us_states_simp1 = st_simplify(us_states2163, dTolerance = 100000)  # 100 km
```

Une limitation de `st_simplify()` est qu'il simplifie les objets sur une base géométrique.
Cela signifie que la "topologie" est perdue, ce qui donne lieu à des polygones se superposant ou séparés par des vides, comme le montre la figure \@ref(fig:us-simp) (panneau du milieu).
`ms_simplify()` de **rmapshaper** fournit une alternative qui surmonte ce problème.
Par défaut, il utilise l'algorithme de Visvalingam, qui surmonte certaines limitations de l'algorithme de Douglas-Peucker [@visvalingam_line_1993].
<!-- https://bost.ocks.org/mike/simplify/ -->
L'extrait de code suivant utilise cette fonction pour simplifier `us_states2163`.
Le résultat n'a que 1% des sommets de l'entrée (fixée à l'aide de l'argument `keep`) mais son nombre d'objets reste intact car nous avons fixé `keep_shapes = TRUE` :^[
La simplification des objets multi-polygones peut supprimer les petits polygones internes, même si l'argument `keep_shapes` est défini à TRUE. Pour éviter cela, vous devez définir `explode = TRUE`. Cette option convertit tous les mutlipolygones en polygones séparés avant leur simplification.
]

```{r 05-geometry-operations-6, warning=FALSE, message=FALSE}
# proportion des points à garder (0-1; par defaut 0.05)
us_states_simp2 = rmapshaper::ms_simplify(us_states2163, keep = 0.01,
                                          keep_shapes = TRUE)
```


Une alternative à la simplification est le lissage des limites des géométries des polygones et des linéaires (*linestring*). Elle est implémenté dans le package **smoothr**. 
Le lissage interpole les arêtes des géométries et n'entraîne pas nécessairement une réduction du nombre de sommets, mais il peut être particulièrement utile lorsque l'on travaille avec des géométries qui résultent de la vectorisation spatiale d'un raster (un sujet traité dans le chapitre \@ref(raster-vector).
**smoothr** implémente trois techniques de lissage : une régression à noyau gaussien, l'algorithme de découpage en coins de Chaikin et l'interpolation par splines, qui sont tous décrits dans la vignette du paquetage et dans [website](https://strimas.com/smoothr/). 
Notez que, comme pour `st_simplify()`, les algorithmes de lissage ne préservent pas la 'topologie'.
La fonction phare de **smoothr** est `smooth()`, où l'argument `method` spécifie la technique de lissage à utiliser.
Vous trouverez ci-dessous un exemple d'utilisation de la régression à noyau gaussien pour lisser les frontières des états américains en utilisant `method=ksmooth`.
L'argument `smoothness` contrôle la largeur de bande de la gaussienne qui est utilisée pour lisser la géométrie et a une valeur par défaut de 1.

```{r 05-geometry-operations-6b, warning=FALSE}
us_states_simp3 = smoothr::smooth(us_states2163, method = 'ksmooth', smoothness = 6)

```

Enfin, la comparaison visuelle de l'ensemble de données originales et des deux versions simplifiées montre des différences entre les sorties des algorithmes de Douglas-Peucker (`st_simplify`), de Visvalingam (`ms_simplify`) et de régression à noyau gaussien (`smooth(method=ksmooth)`) (Figure \@ref(fig:us-simp)) :

```{r us-simp, echo=FALSE, fig.cap="Simplification des polygones, comparant la géométrie originale des États-Unis continentaux avec des versions simplifiées, générées avec les fonctions des paquets sf (haut à droite) et rmapshaper (bas à gauche) et smoothr (bas à droite).", warning=FALSE, fig.asp=0.3, fig.scap="Examples de simplification de polygones."}
library(tmap)
p_ussimp1 = tm_shape(us_states2163) + tm_polygons() + tm_layout(main.title = "Donnée source")
p_ussimp2 = tm_shape(us_states_simp1) + tm_polygons() + tm_layout(main.title = "st_simplify")
p_ussimp3 = tm_shape(us_states_simp2) + tm_polygons() + tm_layout(main.title = "ms_simplify")
p_ussimp4 = tm_shape(us_states_simp3) + tm_polygons() + tm_layout(main.title = "smooth(method=ksmooth)")
tmap_arrange(p_ussimp1, p_ussimp2, p_ussimp3, p_ussimp4, ncol = 2, nrow = 2)
```

### Centroïdes

\index{vector!centroids} 
Les opérations de centroïdes identifient le centre des objets géographiques.
Comme pour les mesures statistiques de tendance centrale (y compris les définitions de la moyenne et de la médiane), il existe de nombreuses façons de définir le centre géographique d'un objet.
Toutes créent des représentations par un point unique d'objets vectoriels plus complexes.

Le *centroïde géographique* est sans doute l'opération la plus couramment utilisée.
Ce type d'opération (souvent jute appelé "centroïde") représente le centre de masse d'un objet spatial (pensez à une assiette en équilibre sur votre doigt).
Les centroïdes géographiques ont de nombreuses utilisations, par exemple pour créer une représentation ponctuelle simple de géométries complexes, ou pour estimer les distances entre polygones.
Ils peuvent être calculés à l'aide de la fonction **sf** `st_centroid()`, comme le montre le code ci-dessous, qui génère les centroïdes géographiques de régions de Nouvelle-Zélande et d'affluents de la Seine, illustrés par des points noirs sur la figure \@ref(fig:centr).

```{r 05-geometry-operations-7, warning=FALSE}
nz_centroid = st_centroid(nz)
seine_centroid = st_centroid(seine)
```

Parfois, le centroïde géographique se trouve en dehors des limites de l'objet parent (pensez à un beignet).
Dans ce cas, les opérations dites de *point sur la surface* peuvent être utilisées pour garantir que le point se trouvera dans l'objet parent (par exemple, pour étiqueter des objets de type multipolygones irréguliers tels que des îles), comme l'illustrent les points rouges de la figure \@ref(fig:centr).
Remarquez que ces points rouges se trouvent toujours sur leurs objets parents.
Ils ont été créés avec `st_point_on_surface()` comme suit :^[
Une description du fonctionnement de `st_point_on_surface()` est fournie sur https://gis.stackexchange.com/q/76498.
]

```{r 05-geometry-operations-8, warning=FALSE}
nz_pos = st_point_on_surface(nz)
seine_pos = st_point_on_surface(seine)
```

```{r centr, warning=FALSE, echo=FALSE, fig.cap="Centroïdes (points noirs) et points sur la surface (points rouges) des ensembles de données des régions de Nouvelle-Zélande (à gauche) et de la Seine (à droite).", fig.scap="Centroïde et point sur la surface."}
p_centr1 = tm_shape(nz) + tm_borders() +
  tm_shape(nz_centroid) + tm_symbols(shape = 1, col = "black", size = 0.5) +
  tm_shape(nz_pos) + tm_symbols(shape = 1, col = "red", size = 0.5)  
p_centr2 = tm_shape(seine) + tm_lines() +
  tm_shape(seine_centroid) + tm_symbols(shape = 1, col = "black", size = 0.5) +
  tm_shape(seine_pos) + tm_symbols(shape = 1, col = "red", size = 0.5)  
tmap_arrange(p_centr1, p_centr2, ncol = 2)
```

Il existe d'autres types de centroïdes, notamment le *centre de Chebyshev* et le *centre visuel*.
Nous ne les explorerons pas ici, mais il est possible de les calculer à l'aide de R, comme nous le verrons dans le chapitre \@ref(algorithms).

### Buffers/tampons

\index{vector!buffers} 
Les buffers ou tampons sont des polygones représentant la zone située à une distance donnée d'une caractéristique géométrique :
Que le type d'origine soit un point, une ligne ou un polygone, la sortie est toujours un polygone.
Contrairement à la simplification (qui est souvent utilisée pour la visualisation et la réduction de la taille des fichiers), la mise en mémoire tampon est généralement utilisée pour l'analyse des données géographiques.
Combien de points se trouvent à une distance donnée de cette ligne ?
Quels groupes démographiques se trouvent à une distance de déplacement de ce nouveau magasin ?
Il est possible de répondre à ce genre de questions et de les visualiser en créant des tampons autour des entités géographiques d'intérêt.

La figure \@ref(fig:buffs) illustre des buffers de différentes tailles (5 et 50 km) entourant la Seine et ses affluents.
Les commandes ci-dessous, utilisées pour créer ces buffers,  montrent que la commande `st_buffer()` nécessite au moins deux arguments : une géométrie d'entrée et une distance, fournie dans les unités du SRC (dans ce cas, les mètres) :

```{r 05-geometry-operations-9}
seine_buff_5km = st_buffer(seine, dist = 5000)
seine_buff_50km = st_buffer(seine, dist = 50000)
```

```{r buffs, echo=FALSE, fig.cap="Tampons de 5 km autour du jeu de données de la Seine  (à gauche) et de 50 km (à droite). Notez les couleurs, qui reflètent le fait qu'un tampon est créé par élément géométrique.", fig.show='hold', out.width="75%", fig.scap="Tampons autour du jeu de données de la Seine."}
p_buffs1 = tm_shape(seine_buff_5km) + tm_polygons(col = "name") +
  tm_shape(seine) + tm_lines() +
  tm_layout(main.title = "Tampons de 5 km", legend.show = FALSE)
p_buffs2 = tm_shape(seine_buff_50km) + tm_polygons(col = "name") +
  tm_shape(seine) + tm_lines() +
  tm_layout(main.title = "Tampons de 50 km", legend.show = FALSE)
tmap_arrange(p_buffs1, p_buffs2, ncol = 2)
```

```{block2 05-geometry-operations-10, type = "rmdnote"}
Le troisième et dernier argument de `st_buffer()` est `nQuadSegs`, qui signifie 'nombre de segments par quadrant' et qui est fixé par défaut à 30 (ce qui signifie que les cercles créés par les buffers sont composés de $4 \times 30 = 120$ lignes).
Cet argument a rarement besoin d´être défini.
Les cas inhabituels où il peut être utile incluent lorsque la mémoire consommée par la sortie d´une opération de tampon est une préoccupation majeure (dans ce cas, il devrait être réduit) ou lorsque la très haute précision est nécessaire (dans ce cas, il devrait être augmenté).
```

```{r nQuadSegs, eval=FALSE, echo=FALSE}
# Demonstrate nQuadSegs
seine_buff_simple = st_buffer(seine, dist = 50000, nQuadSegs = 3)
plot(seine_buff_simple, key.pos = NULL, main = "Buffer de 50 km")
plot(seine, key.pos = NULL, lwd = 3, pal = rainbow, add = TRUE)
seine_points = st_cast(seine[1, ], "POINT")
buff_single = st_buffer(seine_points[1, ], 50000, 2)
buff_points = st_cast(buff_single, "POINT")
plot(st_geometry(buff_single), add = TRUE)
```

### Application affine

\index{vector!affine transformation} 
Une application affine est une transformation qui préserve les lignes et le parallélisme.
Cependant, les angles ou la longueur ne sont pas nécessairement préservés.
Les transformations affines comprennent, entre autres, le déplacement (translation), la mise à l'échelle et la rotation.
En outre, il est possible d'utiliser n'importe quelle combinaison de celles-ci.
Les applications affines sont une partie essentielle de la géocomputation.
Par exemple, le décalage est nécessaire pour le placement d'étiquettes, la mise à l'échelle est utilisée dans les cartogrammes de zones non contiguës (voir la section \@ref(other-mapping-packages)), et de nombreuses transformations affines sont appliquées lors de la reprojection ou de l'amélioration de la géométrie créée à partir d'une carte déformée ou mal projetée.
Le paquet **sf** implémente la transformation affine pour les objets des classes `sfg` et `sfc`.

```{r 05-geometry-operations-11}
nz_sfc = st_geometry(nz)
```

Le décalage déplace chaque point de la même distance en unités cartographiques.
Cela peut être fait en ajoutant un vecteur numérique à un objet vectoriel.
Par exemple, le code ci-dessous déplace toutes les coordonnées y de 100 000 mètres vers le nord, mais laisse les coordonnées x intactes (panneau gauche de la figure \@ref(fig:affine-trans)).  

```{r 05-geometry-operations-12}
nz_shift = nz_sfc + c(0, 100000)
```

La mise à l'échelle agrandit ou rétrécit les objets par un facteur.
Elle peut être appliquée de manière globale ou locale.
La mise à l'échelle globale augmente ou diminue toutes les valeurs des coordonnées par rapport aux coordonnées d'origine, tout en gardant intactes les relations topologiques de toutes les géométries.
Elle peut être effectuée par soustraction ou multiplication d'un objet `sfg` ou `sfc`.

```{r 05-geometry-operations-13, echo=FALSE,eval=FALSE}
nz_scale0 = nz_sfc * 0.5
```

Le changement à l'échelle locale traite les géométries indépendamment et nécessite des points autour desquels les géométries vont être mises à l'échelle, par exemple des centroïdes.
Dans l'exemple ci-dessous, chaque géométrie est réduite d'un facteur deux autour des centroïdes (panneau central de la figure \@ref(fig:affine-trans)).
Pour cela, chaque objet est d'abord décalé de manière à ce que son centre ait les coordonnées `0, 0` (`(nz_sfc - nz_centroid_sfc)`). 
Ensuite, les tailles des géométries sont réduites de moitié (`* 0.5`).
Enfin, le centroïde de chaque objet est ramené aux coordonnées des données d'entrée (`+ nz_centroid_sfc`). 

```{r 05-geometry-operations-14}
nz_centroid_sfc = st_centroid(nz_sfc)
nz_scale = (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc
```

La rotation de coordonnées bidimensionnelles nécessite une matrice de rotation :

$$
R =
\begin{bmatrix}
\cos \theta & -\sin \theta \\  
\sin \theta & \cos \theta \\
\end{bmatrix}
$$

Elle fait tourner les points dans le sens des aiguilles d'une montre.
La matrice de rotation peut être implémentée dans R comme suit :

```{r 05-geometry-operations-15}
rotation = function(a){
  r = a * pi / 180 #degrées en radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
} 
```

La fonction `rotation` accepte un argument `a` - un angle de rotation en degrés.
La rotation peut être effectuée autour de points sélectionnés, comme les centroïdes (panneau de droite de la figure \@ref(fig:affine-trans)).
Voir `vignette("sf3")` pour plus d'exemples.

```{r 05-geometry-operations-16}
nz_rotate = (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc
```

```{r affine-trans, echo=FALSE, fig.cap="Illustrations des transformations affines : décalage, échelle et rotation.", warning=FALSE, eval=TRUE, fig.scap="Illustrations des transformations affines."}
st_crs(nz_shift) = st_crs(nz_sfc)
st_crs(nz_scale) = st_crs(nz_sfc)
st_crs(nz_rotate) = st_crs(nz_sfc)
p_at1 = tm_shape(nz_sfc) + tm_polygons() +
  tm_shape(nz_shift) + tm_polygons(col = "red") +
  tm_layout(main.title = "Décalage")
p_at2 = tm_shape(nz_sfc) + tm_polygons() +
  tm_shape(nz_scale) + tm_polygons(col = "red") +
  tm_layout(main.title = "Chgt d'échelle")
p_at3 = tm_shape(nz_sfc) + tm_polygons() +
  tm_shape(nz_rotate) + tm_polygons(col = "red") +
  tm_layout(main.title = "Rotation")
tmap_arrange(p_at1, p_at2, p_at3, ncol = 3)
```

```{r 05-geometry-operations-17, echo=FALSE,eval=FALSE}
nz_scale_rotate = (nz_sfc - nz_centroid_sfc) * 0.25 * rotation(90) + nz_centroid_sfc
```
```{r 05-geometry-operations-18, echo=FALSE,eval=FALSE}
shearing = function(hx, hy){
  matrix(c(1, hy, hx, 1), nrow = 2, ncol = 2)
}
nz_shear = (nz_sfc - nz_centroid_sfc) * shearing(1.1, 0) + nz_centroid_sfc
```
```{r 05-geometry-operations-19, echo=FALSE,eval=FALSE}
plot(nz_sfc)
plot(nz_shear, add = TRUE, col = "red")
```

Enfin, les géométries nouvellement créées peuvent remplacer les anciennes avec la fonction `st_set_geometry()` : 

```{r 05-geometry-operations-20}
nz_scale_sf = st_set_geometry(nz, nz_scale)
```

### Découper {#clipping}

\index{vector!clipping} 
\index{spatial!subsetting} 
Le découpage spatial est une forme de sélection spatiale qui implique des changements dans les colonnes `géométriques` d'au moins certaines des entités affectées.

Le découpage ne peut s'appliquer qu'à des éléments plus complexes que des points : 
les lignes, les polygones et leurs équivalents "multi".
Pour illustrer le concept, nous allons commencer par un exemple simple :
deux cercles superposés dont le point central est distant d'une unité et dont le rayon est de un (Figure \@ref(fig:points)).

```{r points, fig.cap="cercles superposés.", fig.asp=0.4, crop = TRUE}
b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # créer 2 points
b = st_buffer(b, dist = 1) # convertir les points en cercles
plot(b, border = "grey")
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3) # ajout du texte
```

Imaginez que vous voulez sélectionner non pas un cercle ou l'autre, mais l'espace couvert par les deux `x` *et* `y`.
Cela peut être fait en utilisant la fonction `st_intersection()`\index{vecteur!intersection}, illustrée en utilisant des objets nommés `x` et `y` qui représentent les cercles de gauche et de droite (Figure \@ref(fig:circle-intersection)).

```{r circle-intersection, fig.cap="Cercles superposés avec une couleur grise pour indiquer l'intersection entre eux", fig.asp=0.4, fig.scap="Cercles superposés montrant les types d'intersection.", crop = TRUE}
x = b[1]
y = b[2]
x_and_y = st_intersection(x, y)
plot(b, border = "grey")
plot(x_and_y, col = "lightgrey", border = "grey", add = TRUE) # surface intersectée
```

Le passage de code suivant montre comment cela fonctionne pour toutes les combinaisons du diagramme de Venn représentant `x` et `y`, inspiré de la [Figure 5.1](http://r4ds.had.co.nz/transform.html#logical-operators) du livre *R for Data Science* [@grolemund_r_2016].

```{r venn-clip, echo=FALSE, fig.cap="Équivalents spatiaux des opérateurs logiques.", warning=FALSE}
source("https://github.com/Robinlovelace/geocompr/raw/main/code/05-venn-clip.R")
# source("code/05-venn-clip.R") # for testing local version, todo: remove or change
```

### Sélection et découpage

Le découpage d'objets peut modifier leur géométrie, mais il peut également sélectionner des objets, en ne renvoyant que les entités qui intersectent (ou intersectent partiellement) un objet de découpage/sélection.
Pour illustrer ce point, nous allons sélectionner les points qui incluent dans le cadre englobant (*bounding box*) des cercles `x` et `y` de la figure \@ref(fig:venn-clip).
Certains points seront à l'intérieur d'un seul cercle, d'autres à l'intérieur des deux et d'autres encore à l'intérieur d'aucun.
`st_sample()` est utilisé ci-dessous pour générer une distribution *simple et aléatoire* de points à l'intérieur de l'étendue des cercles `x` et `y`, ce qui donne le résultat illustré dans la Figure \@ref(fig:venn-subset), ce qui soulève la question suivante : comment sous-ensembler les points pour ne renvoyer que le point qui intersecte *à la fois* `x` et `y` ?

```{r venn-subset, fig.cap="Points distribués de manière aléatoire dans le cadre englobant les cercles x et y. Les points qui croisent les deux objets x et y sont mis en évidence.", fig.height=6, fig.width=9, fig.asp=0.4, fig.scap="Points distribués aléatoirement dans le cadre englobant. Notez qu'un seul point intersecte à la fois x et y, mis en évidence par un cercle rouge.", echo=FALSE}
bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
p_xy1 = p[x_and_y]
plot(box, border = "grey", lty = 2)
plot(x, add = TRUE, border = "grey")
plot(y, add = TRUE, border = "grey")
plot(p, add = TRUE)
plot(p_xy1, cex = 3, col = "red", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 2)
```


```{r venn-subset-to-show, eval=FALSE}
bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
x_and_y = st_intersection(x, y)
```

Le code ci-dessous montre trois façons d'obtenir le même résultat.
Nous pouvons utiliser directement l'intersection `index{vecteur!intersection} de `x` et `y` (représentée par `x_et_y` dans l'extrait de code précédent) comme objet de sélection, comme le montre la première ligne du morceau de code ci-dessous.
Nous pouvons également trouver l'intersection entre les points d'entrée représentés par `p` et l'objet de sélection et de découpage `x_et_y`, comme le montre la deuxième ligne du code ci-dessous.
Cette deuxième approche renvoie les entités qui ont une intersection partielle avec `x_and_y` mais avec des géométries modifiées pour les entités dont les surfaces recoupent celle de l'objet de sélection.
La troisième approche consiste à créer un objet de sélection en utilisant le prédicat spatial binaire `st_intersects()`, introduit dans le chapitre précédent.
Les résultats sont identiques (à l'exception de différences superficielles dans les noms d'attributs), mais l'implémentation diffère substantiellement :

```{r 05-geometry-operations-21}
p_xy1 = p[x_and_y]
p_xy2 = st_intersection(p, x_and_y)
sel_p_xy = st_intersects(p, x, sparse = FALSE)[, 1] &
  st_intersects(p, y, sparse = FALSE)[, 1]
p_xy3 = p[sel_p_xy]
```

```{r 05-geometry-operations-22, echo=FALSE, eval=FALSE}
# testons si les objets sont identiques :
identical(p_xy1, p_xy2)
identical(p_xy2, p_xy3)
identical(p_xy1, p_xy3)
waldo::compare(p_xy1, p_xy2) # identiques, à l'exception des noms d'attributs
waldo::compare(p_xy2, p_xy3) # identiques, à l'exception des noms d'attributs


#  Une autre façon d'échantillonner la bb
bb = st_bbox(st_union(x, y))
pmulti = st_multipoint(pmat)
box = st_convex_hull(pmulti)
```

Bien que l'exemple ci-dessus soit plutôt trivial et fourni à des fins éducatives plutôt qu'appliquées, et que nous encouragions le lecteur à reproduire les résultats pour approfondir sa compréhension de la manipulation des objets vectoriels géographiques dans R, il soulève une question importante : quelle implémentation utiliser ?
En général, les implémentations les plus concises doivent être privilégiées, ce qui signifie la première approche ci-dessus.
Nous reviendrons sur la question du choix entre différentes implémentations d'une même technique ou d'un même algorithme au chapitre \@ref(algorithmes).

### GUnions de géométries

\index{vector!union} 
\index{aggregation!spatial} 
Comme nous l'avons vu dans la section \@ref(vector-attribute-aggregation), l'agrégation spatiale peut dissoudre silencieusement les géométries des polygones se touchant dans le même groupe.
Cela est démontré dans le code ci-dessous dans lequel 49 `us_states` sont agrégés en 4 régions à l'aide des fonctions de R base et du **tidyverse**\index{tidyverse (package)} (voir les résultats dans la figure \@ref(fig:us-regions)) :

```{r 05-geometry-operations-23}
regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
regions2 = us_states |>
  group_by(REGION) |>
  summarize(pop = sum(total_pop_15, na.rm = TRUE))
```

```{r 05-geometry-operations-24, echo=FALSE}
# st_join(buff, africa[, "pop"]) |>
#   summarize(pop = sum(pop, na.rm = TRUE))
# summarize(africa[buff, "pop"], pop = sum(pop, na.rm = TRUE))
```

```{r us-regions, fig.cap="Agrégation spatiale sur des polygones contigus, illustrée par l'agrégation de la population des États américains en régions, la population étant représentée par une couleur. Notez que l'opération dissout automatiquement les frontières entre les états.", echo=FALSE, warning=FALSE, fig.asp=0.2, out.width="100%", fig.scap="Agrégation spatiale sur des polygones contigus."}
source("https://github.com/Robinlovelace/geocompr/raw/main/code/05-us-regions.R", print.eval = TRUE)
```

Que se passe-t-il au niveau des géométries ?
En coulisses, `aggregate()` et `summarize()` combinent les géométries et dissolvent les frontières entre elles en utilisant `st_union()`.
Ceci est démontré par le code ci-dessous qui crée une union des Etats-Unis de l'Ouest :  

```{r 05-geometry-operations-25}
us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)
```

La fonction peut prendre deux géométries et les unir, comme le montre ll code ci-dessous qui crée un bloc occidental uni incorporant le Texas (défi : reproduire et représenter le résultat) :

```{r 05-geometry-operations-26, message=FALSE}
texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)
```

```{r 05-geometry-operations-27, echo=FALSE, eval=FALSE}
plot(texas_union)
# objectif: expérimenter avec st_union
us_south2 = st_union(us_west[1, ], us_west[6, ])
plot(us_southhwest)
```

### Transformations de type {#type-trans}

\index{vector!geometry casting} 
La transformation d'un type de géométrie en un autre (*casting*)  est une opération puissante.
Elle est implémentée dans la fonction `st_cast()` du package **sf**.
Il est important de noter que la fonction `st_cast()` se comporte différemment selon qu'il s'agit d'un objet géométrique simple (`sfg`), d'une entité avec une colonne géométrique simple (`sfc`) ou d'un objet entité simple.

Créons un multipoint pour illustrer le fonctionnement des transformations de type géométrique sur des objets de géométrie simple (`sfg`) :

```{r 05-geometry-operations-28}
multipoint = st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))
```

Dans ce cas, `st_cast()` peut être utile pour transformer le nouvel objet en *linestring* (linéaire) ou en polygone (Figure \@ref(fig:single-cast)) :

```{r 05-geometry-operations-29}
linestring = st_cast(multipoint, "LINESTRING")
polyg = st_cast(multipoint, "POLYGON")
```

```{r single-cast, echo = FALSE, fig.cap="Exemples de lignes et de polygones créés à partir d'une géométrie multipoint", warning=FALSE, fig.asp=0.3, fig.scap="Exemples d'opérations de transformations de type de géométrie."}
p_sc1 = tm_shape(st_sfc(multipoint)) + tm_symbols(shape = 1, col = "black", size = 0.5) +
  tm_layout(main.title = "MULTIPOINT", inner.margins = c(0.05, 0.05, 0.05, 0.05))
p_sc2 = tm_shape(st_sfc(linestring)) + tm_lines() +
  tm_layout(main.title = "LIGNE")
p_sc3 = tm_shape(st_sfc(polyg)) + tm_polygons(border.col = "black") +
  tm_layout(main.title = "POLYGONE")
tmap_arrange(p_sc1, p_sc2, p_sc3, ncol = 3)
```

La conversion de multipoint en ligne est une opération courante qui crée un objet ligne à partir d'observations ponctuelles ordonnées, telles que des mesures GPS ou des sources géolocalisés.
Cela permet d'effectuer des opérations spatiales telles que la longueur du chemin parcouru.
La conversion de multipoint ou de *linestring* en polygone est souvent utilisée pour calculer une surface, par exemple à partir de l'ensemble des mesures GPS prises autour d'un lac ou des coins d'un terrain à bâtir.

Le processus de transformation peut également être inversé en utilisant `st_cast()` :

```{r 05-geometry-operations-30}
multipoint_2 = st_cast(linestring, "MULTIPOINT")
multipoint_3 = st_cast(polyg, "MULTIPOINT")
all.equal(multipoint, multipoint_2)
all.equal(multipoint, multipoint_3)
```

```{block2 05-geometry-operations-31, type='rmdnote'}
Pour les entités de géométries simples (`sfg`), `st_cast` permet également de transformer des géométries de non-multi-types vers des multi-types (par exemple, `POINT` vers `MULTIPOINT`) et de multi-types vers des non-multi-types.
Toutefois, dans le deuxième groupe de cas, seul le premier élément de l´ancien objet est conservé.
```

```{r 05-geometry-operations-32, include=FALSE}
cast_all = function(xg) {
  lapply(c("MULTIPOLYGON", "MULTILINESTRING", "MULTIPOINT", "POLYGON", "LINESTRING", "POINT"), 
         function(x) st_cast(xg, x))
}
t = cast_all(multipoint)
t2 = cast_all(polyg)
```

La transformation en différent types géométrique des d'entités de type simple colonne  (`sfc`) et des objets d'entités simples fonctionnent de la même manière que pour les entités de géométries simples (`sfg`) dans la plupart des cas. 
Une différence importante est la conversion des multi-types en non-multi-types.
À la suite de ce processus, les multi-objets, `sf` ou `sfg` sont divisés en plusieurs non-multi-objets.

Le tableau \@ref(tab:sfs-st-cast) montre les transformations de type géométrique possibles sur les objets d'entités simples.
Les géométries d'entités simples (représentées par la première colonne du tableau) peuvent être transformées en plusieurs types de géométrie, représentés par les colonnes du tableau \@ref(tab:sfs-st-cast)
Plusieurs des transformations ne sont pas possibles, par exemple, vous ne pouvez pas convertir un point unique en un multilinéaire (*multilinestring*) ou un polygone (ainsi les cellules `[1, 4:5]` dans le tableau sont NA).
Certaines transformations divisent un seul élément en plusieurs sous-éléments, en "étendant" les objets `sf` (en ajoutant de nouvelles lignes avec des valeurs d'attributs dupliquées).
Par exemple, lorsqu'une géométrie multipoint composée de cinq paires de coordonnées est transformée en une géométrie "POINT", la sortie contiendra cinq entités.

```{r sfs-st-cast, echo=FALSE}
sfs_st_cast = read.csv("extdata/sfs-st-cast.csv")
abbreviate_geomtypes = function(geomtypes) {
  geomtypes_new = gsub(pattern = "POINT", replacement = "POI", x = geomtypes)
  geomtypes_new = gsub(pattern = "POLYGON", replacement = "POL", x = geomtypes_new)
  geomtypes_new = gsub(pattern = "LINESTRING", replacement = "LIN", x = geomtypes_new)
  geomtypes_new = gsub(pattern = "MULTI", replacement = "M", x = geomtypes_new)
  geomtypes_new = gsub(pattern = "GEOMETRYCOLLECTION", replacement = "GC", x = geomtypes_new)
  geomtypes_new
}
sfs_st_cast$input_geom = abbreviate_geomtypes(sfs_st_cast$input_geom)
names(sfs_st_cast) = abbreviate_geomtypes(names(sfs_st_cast))
names(sfs_st_cast)[1] = ""
knitr::kable(sfs_st_cast, 
             caption = paste("Transformation de type de géométrie sur des entités simples", 
                             "(voir section 2.1) avec un type d'entrée par ligne et", 
                             "type de sortie par colonne"),
             caption.short = "Transformation de type de géométrique sur des entités simples.",
             booktabs = TRUE) |> 
  kableExtra::add_footnote("Note : Les valeurs comme (1) représentent le nombre d'entités ; NA signifie que l'opération n'est pas possible. Abréviations : POI, LIN, POL et GC font référence à POINT, LINESTRING, POLYGON et GEOMETRYCOLLECTION. La version MULTI de ces types de géométrie est indiquée par un M précédent, par exemple, MPOI est l'acronyme de MULTIPOINT.", notation = "none")
```

Essayons d'appliquer des transformations de type géométrique sur un nouvel objet, `multilinestring_sf`, à titre d'exemple (à gauche sur la Figure \@ref(fig:line-cast)) :

```{r 05-geometry-operations-33}
multilinestring_list = list(matrix(c(1, 4, 5, 3), ncol = 2), 
                            matrix(c(4, 4, 4, 1), ncol = 2),
                            matrix(c(2, 4, 2, 2), ncol = 2))
multilinestring = st_multilinestring(multilinestring_list)
multilinestring_sf = st_sf(geom = st_sfc(multilinestring))
multilinestring_sf
```

Vous pouvez l'imaginer comme un réseau routier ou fluvial. 
Le nouvel objet n'a qu'une seule ligne qui définit toutes les lignes.
Cela limite le nombre d'opérations qui peuvent être faites, par exemple, cela empêche d'ajouter des noms à chaque segment de ligne ou de calculer les longueurs des lignes individuelles.
La fonction `st_cast()` peut être utilisée dans cette situation, car elle sépare un mutlilinestring en trois linestrings :

```{r 05-geometry-operations-34}
linestring_sf2 = st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
```

```{r line-cast, echo=FALSE, fig.cap="Exemples de transformation de type de géométrie entre MULTILINESTRING (à gauche) et LINESTRING (à droite).", warning=FALSE, fig.scap="Exemples de transformation de type de géométrie."}
p_lc1 = tm_shape(multilinestring_sf) + tm_lines(lwd = 3) +
  tm_title("MULTILINESTRING")
linestring_sf2$name = c("Riddle Rd", "Marshall Ave", "Foulke St")
p_lc2 = tm_shape(linestring_sf2) + tm_lines(lwd = 3, col = "name", col.scale = tm_scale(values = "Set2")) +
  tm_title("LINESTRING") +
  tm_layout(legend.show = FALSE)
tmap_arrange(p_lc1, p_lc2, ncol = 2)
```

Le nouvel objet permet la création d'attributs (voir la section \@ref(vec-attr-creation)) et la mesure de la longueur :

```{r 05-geometry-operations-35}
linestring_sf2$name = c("Riddle Rd", "Marshall Ave", "Foulke St")
linestring_sf2$length = st_length(linestring_sf2)
linestring_sf2
```

## Opérations géométriques sur les données raster {#geo-ras}

\index{raster!manipulation} 
Les opérations  géométriques sur des raster comprennent le décalage, le retournement, la mise en miroir, la mise à l'échelle, la rotation ou la déformation des images.
Ces opérations sont nécessaires pour une variété d'applications, y compris le géoréférencement, utilisé pour permettre aux images d'être superposées sur une carte précise avec un CRS connu [@liu_essential_2009].
Il existe une variété de techniques de géoréférencement, notamment :

- Géorectification basée sur des [points de contrôle au sol](https://www.qgistutorials.com/en/docs/3/georeferencing_basics.html) connus 
- Orthorectification, qui tient également compte de la topographie locale.
- L'[enregistrement](https://en.wikipedia.org/wiki/Image_registration) d'images  est utilisé pour combiner des images de la même chose mais prises par différents capteurs en alignant une image sur une autre (en termes de système de coordonnées et de résolution).

R est plutôt inadapté pour les deux premiers points car ceux-ci nécessitent souvent une intervention manuelle, c'est pourquoi ils sont généralement réalisés à l'aide d'un logiciel SIG dédié (voir également le chapitre : \@ref(gis)).
En revanche, l'alignement de plusieurs images est possible dans R et cette section montre entre autres comment le faire.
Cela implique souvent de modifier l'étendue, la résolution et l'origine d'une image.
Une projection correspondante est bien sûr également nécessaire, mais elle est déjà traitée dans la section \@ref(reproj-ras).

Dans tous les cas, il existe d'autres raisons d'effectuer une opération géométrique sur une seule image raster.
Par exemple, dans le chapitre \@ref(location) nous définissons les zones métropolitaines en Allemagne comme des pixels de 20 km^2^ avec plus de 500.000 habitants. 
La trame d'habitants d'origine a cependant une résolution de 1 km^2^, c'est pourquoi nous allons diminuer (agréger) la résolution d'un facteur 20 (voir le chapitre \@ref(define-metropolitan-areas)).
Une autre raison d'agréger une image matricielle est simplement de réduire le temps d'exécution ou d'économiser de l'espace disque.
Bien entendu, cela n'est possible que si la tâche à accomplir permet une résolution plus grossière.

### Intersections géométriques

\index{raster!intersection} 
Dans la section \@ref(spatial-raster-subsetting), nous avons montré comment extraire des valeurs d'un raster superposé à d'autres objets spatiaux.
Pour récupérer une sortie spatiale, nous pouvons utiliser pratiquement la même syntaxe de sélection.
La seule différence est que nous devons préciser que nous souhaitons conserver la structure matricielle en mettant l'argument `drop` à `FALSE`.
Ceci retournera un objet raster contenant les cellules dont les points médians se chevauchent avec `clip`.

```{r 05-geometry-operations-36}
elev = rast(system.file("raster/elev.tif", package = "spData"))
clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
```

Pour la même opération, nous pouvons également utiliser les commandes `intersect()` et `crop()`.

### Étendue et origine

\index{raster!merging} 
Lors de la fusion ou de l'exécution de l'algèbre raster sur des rasters, leur résolution, leur projection, leur origine et/ou leur étendue doivent correspondre. Sinon, comment ajouter les valeurs d'un raster ayant une résolution de 0,2 degré décimal à un second raster ayant une résolution de 1 degré décimal ?
Le même problème se pose lorsque nous souhaitons fusionner des images satellite provenant de différents capteurs avec des projections et des résolutions différentes. 
Nous pouvons traiter de telles disparités en alignant les trames.

Dans le cas le plus simple, deux images ne diffèrent que par leur étendue.
Le code suivant ajoute une ligne et deux colonnes de chaque côté de l'image raster tout en fixant toutes les nouvelles valeurs à une altitude de 1000 mètres (Figure \@ref(fig:extend-example)).

```{r extend-example0}
elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_2 = extend(elev, c(1, 2))
```

```{r extend-example, fig.cap = "Trame originale (à gauche) et la même trame (à droite) agrandie d'une ligne en haut et en bas et de deux colonnes à gauche et à droite.", fig.scap="Étendre les rasters.", echo=FALSE, fig.asp=0.5}
source("https://github.com/Robinlovelace/geocompr/raw/main/code/05-extend-example.R", print.eval = TRUE)
```

Performing an algebraic operation on two objects with differing extents in R, the **terra** package returns an error.

```{r 05-geometry-operations-37, error=TRUE}
elev_3 = elev + elev_2
```

Cependant, nous pouvons aligner l'étendue de deux rasters avec `extend()`. 
Au lieu d'indiquer à la fonction le nombre de lignes ou de colonnes à ajouter (comme nous l'avons fait précédemment), nous lui permettons de le déterminer en utilisant un autre objet raster.
Ici, nous étendons l'objet `elev` à l'étendue de `elev_2`. 
Les lignes et colonnes nouvellement ajoutées prennent la valeur  `NA`.

```{r 05-geometry-operations-38}
elev_4 = extend(elev, elev_2)
```

L'origine d'un raster est le coin de la cellule le plus proche des coordonnées (0, 0).
La fonction `origin()` renvoie les coordonnées de l'origine.
Dans l'exemple ci-dessous, un coin de cellule existe avec les coordonnées (0, 0), mais ce n'est pas toujours le cas.

```{r 05-geometry-operations-39}
origin(elev_4)
```

Si deux rasters ont des origines différentes, leurs cellules ne se chevauchent pas complètement, ce qui rends l'algèbre raster impossible.
Pour changer l'origine -- utilisez `origin()`.^[
Si les origines de deux données matricielles ne sont que marginalement éloignées, il suffit parfois d'augmenter l'argument `tolerance` de `terra::terraOptions()`.
]
La figure \@ref(fig:origin-example) révèle l'effet de la modification de l'origine de cette manière.

```{r}
# changer l'origine
origin(elev_4) = c(0.25, 0.25)
```

```{r origin-example, fig.cap="Rasters avec des valeurs identiques mais des origines différentes.", echo=FALSE}
elev_poly = st_as_sf(as.polygons(elev, dissolve = FALSE))
elev4_poly = st_as_sf(as.polygons(elev_4, dissolve = FALSE))
tm_shape(elev4_poly) +
  tm_grid() +
  tm_polygons(col = "elev") +
  tm_shape(elev_poly) +
  tm_polygons(col = "elev") +
  tm_layout(frame = FALSE, legend.show = FALSE,
            inner.margins = c(0.1, 0.12, 0, 0))
# # See https://github.com/Robinlovelace/geocompr/issues/695
# knitr::include_graphics("https://user-images.githubusercontent.com/1825120/146618199-786fe3ad-9718-4dd0-a640-41180fc17e63.png")
```

Notez que le changement de résolution (section suivante) modifie souvent aussi l'origine.

### Agrégation et désagrégation

\index{raster!aggregation} 
\index{raster!disaggregation} 
Les jeux de données raster peuvent également différer en ce qui concerne leur résolution. 
Pour faire correspondre les résolutions, on peut soit diminuer (`aggregate()`) soit augmenter (`disagg()`) la résolution des rasters.^[
Nous faisons ici référence à la résolution spatiale.
En télédétection, les résolutions spectrale (bandes spectrales), temporelle (observations dans le temps de la même zone) et radiométrique (profondeur de couleur) sont également importantes.
Consultez l'exemple `tapp()` dans la documentation pour avoir une idée sur la façon de faire une agrégation de raster temporel.
]
À titre d'exemple, nous modifions ici la résolution spatiale de `dem` (trouvé dans le paquet **spDataLarge**) par un facteur 5 (Figure \@ref(fig:aggregate-example)).
De plus, la valeur de la cellule de sortie doit correspondre à la moyenne des cellules d'entrée (notez que l'on pourrait également utiliser d'autres fonctions, telles que `median()`, `sum()`, etc ):

```{r 05-geometry-operations-40}
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg = aggregate(dem, fact = 5, fun = mean)
```

```{r aggregate-example, fig.cap = "Raster original (gauche). Raster agrégé (droite).", echo=FALSE}
p_ar1 = tm_shape(dem) +
  tm_raster(style = "cont", legend.show = FALSE) +
  tm_layout(main.title = "A. Original", frame = FALSE)
p_ar2 = tm_shape(dem_agg) +
  tm_raster(style = "cont", legend.show = FALSE) +
  tm_layout(main.title = "B. Aggregated", frame = FALSE)
tmap_arrange(p_ar1, p_ar2, ncol = 2)
```

La fonction `disagg()` augmente la résolution des objets matriciels, en fournissant deux méthodes pour assigner des valeurs aux cellules nouvellement créées : la méthode par défaut (`method = "near"`) donne simplement à toutes les cellules de sortie la valeur de la cellule d'entrée, et donc duplique les valeurs, ce qui conduit à une sortie "en bloc".
La méthode `bilinear` utilise les quatre centres de pixels les plus proches de l'image d'entrée (points de couleur saumon sur la figure \@ref(fig:bilinear)) pour calculer une moyenne pondérée par la distance (flèches sur la figure \@ref(fig:bilinear).
La valeur de la cellule de sortie est représentée par un carré dans le coin supérieur gauche de la figure \@ref(fig:bilinear)).

```{r 05-geometry-operations-41}
dem_disagg = disagg(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)
```

```{r bilinear, echo = FALSE, fig.width=8, fig.height=10, fig.cap="La moyenne pondérée par la distance des quatre cellules d'entrée les plus proches détermine la sortie lors de l'utilisation de la méthode bilinéaire pour la désagrégation.", fig.scap="Bilinear disaggregation in action.", warning=FALSE}
source("https://github.com/Robinlovelace/geocompr/raw/main/code/05-bilinear.R", print.eval = TRUE)
# knitr::include_graphics("https://user-images.githubusercontent.com/1825120/146619205-3c0c2e3f-9e8b-4fda-b014-9c342a4befbb.png")
```

En comparant les valeurs de `dem` et `dem_disagg`, on constate qu'elles ne sont pas identiques (vous pouvez aussi utiliser `compareGeom()` ou `all.equal()`).
Cependant, il ne fallait pas s'y attendre, puisque la désagrégation est une simple technique d'interpolation.
Il est important de garder à l'esprit que la désagrégation permet d'obtenir une résolution plus fine ; les valeurs correspondantes, cependant, ne peuvent qu'êtres aussi précises que leur source de résolution initiale.

### Rééchantillonnage

\index{raster!resampling}
Les méthodes d'agrégation et de désagrégation ci-dessus ne conviennent que lorsque nous voulons modifier la résolution de notre raster par le facteur d'agrégation/désagrégation. 
Cependant, que faire lorsque nous avons deux ou plusieurs raster avec des résolutions et des origines différentes ?
C'est le rôle du rééchantillonnage - un processus de calcul des valeurs pour les nouveaux emplacements des pixels.
En bref, ce processus prend les valeurs de notre raster original et recalcule de nouvelles valeurs pour un raster cible avec une résolution et une origine personnalisées.

<!--toDo: jn-->
<!-- consider if adding this new figure makes sense -->

```{r, echo=FALSE, eval=FALSE}
target_rast = rast(xmin = 794600, xmax = 798200, 
                   ymin = 8931800, ymax = 8935400,
                   resolution = 150, crs = "EPSG:32717")

target_rast_p = st_as_sf(as.polygons(target_rast))
dem_resampl1 = resample(dem, target_rast, method = "near")

tm1 = tm_shape(dem) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE)

tm2 = tm_shape(dem) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE) +
  tm_shape(target_rast_p) +
  tm_borders()

tm3 = tm_shape(dem_resampl1) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE)

tmap_arrange(tm1, tm2, tm3, nrow = 1)
```


Il existe plusieurs méthodes pour estimer les valeurs d'un raster avec différentes résolutions/origines, comme le montre la figure \@ref(fig:resampl).
Ces méthodes comprennent :

- Plus proche voisin - attribue la valeur de la cellule la plus proche du raster original à la cellule du raster cible.
Cette méthode est rapide et convient généralement aux réechantillonnage de raster de catégories.
- Interpolation bilinéaire - affecte une moyenne pondérée des quatre cellules les plus proches de l'image originale à la cellule de l'image cible (Figure \@ref(fig:bilinear)). Il s'agit de la méthode la plus rapide pour les rasters continus
- Interpolation cubique - utilise les valeurs des 16 cellules les plus proches de la trame d'origine pour déterminer la valeur de la cellule de sortie, en appliquant des fonctions polynomiales du troisième ordre. Elle est aussi utilisée pour les raster continus. Elle permet d'obtenir une surface plus lissée que l'interpolation bilinéaire, mais elle est également plus exigeante en termes de calcul.
- Interpolation par spline cubique - utilise également les valeurs des 16 cellules les plus proches de la trame d'origine pour déterminer la valeur de la cellule de sortie, mais applique des splines cubiques (fonctions polynomiales du troisième ordre par morceaux) pour obtenir les résultats. Elle est utilisée pour les trames continues
- Rééchantillonnage par fenêtré de Lanczos - utilise les valeurs des 36 cellules les plus proches de la trame d'origine pour déterminer la valeur de la cellule de sortie. Il est utilisé pour les raster continues^[Une explication plus détaillée de cette méthode peut être trouvée sur https://gis.stackexchange.com/a/14361/20955.
]

Les explications ci-dessus mettent en évidence le fait que seul le rééchantillonnage par *voisin le plus proche* est adapté aux rasters contenant des catégories, alors que toutes les méthodes peuvent être utilisées (avec des résultats différents) pour les matrices continues.
En outre, chaque méthode successive nécessite plus de temps de traitement.

Pour appliquer le rééchantillonnage, le package **terra** fournit une fonction `resample()`.
Elle accepte un raster d'entrée (`x`), un raster avec des propriétés spatiales cibles (`y`), et une méthode de rééchantillonnage (`method`).

Nous avons besoin d'un raster avec des propriétés spatiales cibles pour voir comment la fonction `resample()` fonctionne.
Pour cet exemple, nous créons `target_rast`, mais vous utiliserez souvent un objet raster déjà existant.

```{r 05-geometry-operations-42}
target_rast = rast(xmin = 794600, xmax = 798200, 
                   ymin = 8931800, ymax = 8935400,
                   resolution = 150, crs = "EPSG:32717")
```

Ensuite, nous devons fournir nos deux objets rasters comme deux premiers arguments et l'une des méthodes de rééchantillonnage décrites ci-dessus.

```{r 05-geometry-operations-42b}
dem_resampl = resample(dem, y = target_rast, method = "bilinear")
```

La figure \@ref(fig:resampl) montre une comparaison de différentes méthodes de rééchantillonnage sur l'objet `dem`.

```{r resampl, echo=FALSE, fig.cap="Comparaison visuelle du raster d'entré et de cinq méthodes de rééchantillonnage différentes."}
dem_resampl1 = resample(dem, target_rast, method = "near")
dem_resampl2 = resample(dem, target_rast, method = "bilinear")
dem_resampl3 = resample(dem, target_rast, method = "cubic")
dem_resampl4 = resample(dem, target_rast, method = "cubicspline")
dem_resampl5 = resample(dem, target_rast, method = "lanczos")

library(tmap)
tm1 = tm_shape(dem) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "Raster original")
tm2 = tm_shape(dem_resampl1) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "proche voisins")
tm3 = tm_shape(dem_resampl2) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "bilinéare")
tm4 = tm_shape(dem_resampl3) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "cubique")
tm5 = tm_shape(dem_resampl4) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "cubiquespline")
tm6 = tm_shape(dem_resampl5) +
  tm_raster(breaks = seq(200, 1100, by = 150), legend.show = FALSE) +
  tm_layout(frame = FALSE, main.title = "lanczos")
tmap_arrange(tm1, tm2, tm3, tm4, tm5, tm6)
```

La fonction `resample()` dispose également de quelques méthodes de rééchantillonnage supplémentaires, dont `sum`, `min`, `q1`, `med`, `q3`, `max`, `average`, `mode`, et `rms`.
Elles calculent toutes une statistique donnée en se basant sur les valeurs de toutes les cellules de la grille (hors `NA`).
Par exemple, `sum` est utile lorsque chaque cellule de raster représente une variable étendue dans l'espace (par exemple, le nombre de personnes).
En utilisant `sum`, le raster ré-échantillonné devrait avoir le même nombre total de personnes que le raster original.

Comme vous le verrez dans la section \@ref(reproj-ras), la reprojection de raster est un cas particulier de rééchantillonnage lorsque notre raster cible a un CRS différent de la trame d'origine.

\index{GDAL}
```{block2 type='rmdnote'}
La plupart des opérations géométriques dans **terra** sont conviviales, plutôt rapides, et fonctionnent sur de grands objets rasters.
Cependant, il peut y avoir des cas où **terra** n´est pas le plus performant, que ce soit pour des objets rasters étendus ou pour de nombreux fichiers rasters, et où des alternatives doivent être envisagées.

Les alternatives les plus établies sont fournies par la bibliothèque GDAL.
Elle contient plusieurs fonctions utilitaires, dont :

- `gdalinfo` - liste diverses informations sur un fichier raster, y compris sa résolution, son CRS, sa boîte de délimitation, et plus encore.
- `gdal_translate` - convertit les données raster entre différents formats de fichiers.
- `gdal_rasterize` - Convertit les données vectorielles en fichiers raster.
- `gdalwarp` - permet le mosaïquage, le rééchantillonnage, le recadrage et la reprojection de données matricielles.

Toutes les fonctions ci-dessus sont écrites en C++, mais peuvent être appelées dans R en utilisant le paquet **gdalUtilities**.
Il est important de noter que toutes ces fonctions attendent un chemin de fichier raster en entrée et retournent souvent leur sortie sous forme de fichier raster (par exemple, `gdalUtilities::gdal_translate("mon_fichier.tif", "nouveau_fichier.tif", t_srs = "EPSG:4326")`).
Ceci est très différent de l´approche habituelle de **terra**, qui attend des objets `SpatRaster` en entrée.
```

## Exercises

```{r, echo=FALSE, results='asis'}
res = knitr::knit_child('_05-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE))
cat(res, sep = '\n')
```