04-spatial-operations.Rmd

# Géotraitements {#spatial-operations}

## Prérequis {-}

- Ce chapitre nécessité les mêmes paquets que ceux utilisés dans le chapitre \@ref(attr): 

```{r 04-spatial-operations-1, message=FALSE, results='hide'}
library(sf)
library(terra)
library(dplyr)
library(spData)
```

- Vous devrez également charger deux jeux de données pour cette section \@ref(spatial-ras)

```{r 04-spatial-operations-1-1}
elev = rast(system.file("raster/elev.tif", package = "spData"))
grain = rast(system.file("raster/grain.tif", package = "spData"))
```

## Introduction

Les opérations spatiales, y compris les jointures spatiales entre les ensembles de données vectorielles et les opérations locales et focales sur les ensembles de données raster, constituent une partie essentielle de la géocomputation\index{geocomputation}.
Ce chapitre montre comment les objets spatiaux peuvent être modifiés d'une multitude de façons en fonction de leur emplacement et de leur forme.
De nombreuses opérations spatiales ont un équivalent non spatial (par exemple via leurs attributs), de sorte que des concepts tels que la sélection et la jonction de jeux de données démontrés dans le chapitre précédent sont applicables ici.
Cela est particulièrement vrai pour les opérations *vectorielles* : La section \@ref(vector-attribute-manipulation) sur la manipulation des tables attributaires fournit la base pour comprendre son équivalent spatial, à savoir la sélection spatial (traitée dans la section \@ref(spatial-subsetting)).
La jointure spatiale (section \@ref(spatial-joining)) et l'agrégation (section \@ref(spatial-aggr)) ont également des contreparties non spatiales, traitées dans le chapitre précédent.

Les opérations spatiales diffèrent toutefois des opérations non spatiales à plusieurs égards :
Les jointures spatiales, par exemple, peuvent être effectuées de plusieurs manières --- y compris la mise en correspondance d'entités qui se croisent ou se trouvent à une certaine distance de l'ensemble de données cible --- alors que les jointures de table attributaire abordées dans la section \@ref(vector-attribute-joining) du chapitre précédent ne peuvent être effectuées que d'une seule manière (sauf lorsqu'on utilise des jointures floues, comme décrit dans la documentation du paquet [**fuzzyjoin**](https://cran.r-project.org/package=fuzzyjoin)).
Les différents *types* de relations spatiales entre objets comme les superpositions/intersections   et les objets disjoints, sont décrits dans la section \@ref(topological-relations).
\index{spatial operations}
Un autre aspect unique des objets spatiaux est la distance : tous les objets spatiaux sont liés par l'espace et les calculs de distance peuvent être utilisés pour explorer la force de cette relation, comme décrit dans le contexte des données vectorielles à la section \@ref(relations-distance).

Les opérations spatiales sur les rasters comprennent la sélection --- traité dans la section \@ref(spatial-raster-subsetting) --- et la fusion de plusieurs " tuiles " raster en un seul objet, comme le montre la section \@ref(merging-rasters).
*L'algèbre de raster* couvre une gamme d'opérations qui modifient les valeurs des cellules, avec ou sans référence aux valeurs des cellules environnantes.
Le concept d'algèbre de raster vital pour de nombreuses applications,  est présenté dans la section \@ref(map-algebra) ; les opérations d'algèbre de raster locales, focales et zonales sont traitées respectivement dans les sections \@ref(local-operations), \@ref(focal-operations) et \@ref(zonal-operations). Les opérations d'algèbre globales, qui génèrent des statistiques synthétiques représentant l'ensemble d'un jeu de données raster, et les calculs de distance sur les données raster, sont abordés dans la section \@ref(global-operations-and-distances).
Dans la dernière section avant les exercices (\@ref(merging-rasters)), le processus de fusion de deux ensembles de données raster est abordé et démontré à l'aide d'un exemple reproductible.

```{block2 04-spatial-operations-2, type='rmdnote'}
Il est important de noter que les opérations spatiales qui utilisent deux objets spatiaux reposent sur le fait que les deux objets ont le même système  de coordonnées de référence, un sujet qui a été introduit dans la section \@ref(crs-intro) et qui sera traité plus en profondeur dans le chapitre \@ref(reproj-geo-data).
```

## Géotraitements sur des données vectorielles {#spatial-vec}

Cette section fournit une vue d'ensemble des opérations spatiales sur les données géographiques vectorielles représentées sous forme de *simple features* du package **sf**.
La section \@ref(spatial-ras) présente les opérations spatiales sur les ensembles de données raster à l'aide des classes et des fonctions du paquet **terra**.

### Sélection spatiale

La sélection spatial est le processus qui consiste à prendre un objet spatial et à renvoyer un nouvel objet contenant uniquement les caractéristiques en relation dans l'espace à un autre objet.
De manière analogue à la  *sélection d'attributs* (traité dans la section \@ref(vector-attribute-subsetting)), des sélection de jeux de données `sf` peuvent être créés avec l'opérateur de crochets (`[`) en utilisant la syntaxe `x[y, , op = st_intersects]`, où `x` est un objet `sf` à partir duquel un sous-ensemble de lignes sera retourné, `y` est l'objet de sous-ensemble et `, op = st_intersects` est un argument optionnel qui spécifie la relation topologique (également connue sous le nom de prédicat binaire) utilisée pour faire la sélection.
La relation topologique par défaut utilisée lorsqu'un argument `op` n'est pas fourni est `st_intersects()` : la commande `x[y, ]` est identique à `x[y, , op = st_intersects]` montrée ci-dessus mais pas à `x[y, , op = st_disjoint]` (la signification de ces relations topologiques et des autres est décrite dans la section suivante).
La fonction `filter()` du **tidyverse**\index{tidyverse (package)} peut également être utilisée mais cette approche est plus verbeuse, comme nous le verrons dans les exemples ci-dessous.
\index{vector!subsetting}
\index{spatial!subsetting}

To demonstrate spatial subsetting, we will use the `nz` and `nz_height` datasets in the **spData** package, which contain geographic data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure \@ref(fig:nz-subset)), in a projected coordinate system.
The following code chunk creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region:

```{r 04-spatial-operations-3}
canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
```

```{r nz-subset, echo=FALSE, warning=FALSE, fig.cap="Exemple de sélection spatiale avec des triangles rouges représentant 101 points hauts en Nouvelle-Zélande, regroupés près de la région centrale de Canterbury (à gauche). Les points dans la région de Canterbury ont été créés avec l'opérateur de sélection `[` (surligné en gris, à droite).", fig.scap="Exemple de sélection spatiale..", message=FALSE}
library(tmap)
p_hpnz1 = tm_shape(nz) + tm_polygons(col = "white") +
  tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 0.25) +
  tm_layout(main.title = "Sommets en Nouvelle Zélande", main.title.size = 1,
            bg.color = "lightblue")
p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") +
  tm_shape(canterbury) + tm_fill(col = "gray") + 
  tm_shape(canterbury_height) + tm_symbols(shape = 2, col = "red", size = 0.25) +
  tm_layout(main.title = "Sommets à Canterbury", main.title.size = 1,
            bg.color = "lightblue")
tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2)
```

Comme pour la sélection d'attributs, la commande `x[y, ]` (équivalente à `nz_height[canterbury, ]`) sélectionne les caractéristiques d'une *cible* `x` en utilisant le contenu d'un objet *source* `y`.
Cependant, au lieu que `y` soit un vecteur de classe `logical` ou `integer`, pour la sélection spatiale, `x` et `y` doivent être des objets géographiques.
Plus précisément, les objets utilisés pour la sélection spatiale de cette manière doivent avoir la classe `sf` ou `sfc` : `nz` et `nz_height` sont tous deux des jeux de données vectorielles géographiques et ont la classe `sf`, et le résultat de l'opération renvoie un autre objet `sf` représentant les caractéristiques de l'objet cible `nz_height` qui intersectent (dans ce cas, les points hauts qui sont situés dans) la région de `canterbury`.

Diverses *relations topologiques* peuvent être utilisées pour le sélection spatiale. Elles déterminent le type de relation spatiale que les caractéristiques de l'objet cible doivent avoir avec l'objet de  sélection.
Il peut s'agir de *touches* (touche), *crosses* (croisse) ou *within* (dedans), comme nous le verrons bientôt dans la section \@ref(topological-relations). 
Le paramètre par défaut `st_intersects` est une relation topologique 'attrape-tout' qui retournera les éléments de la cible qui *touchent*, *croissent* ou sont *within* (dedans) l'objet source 'sélectionnant'.
Comme indiqué ci-dessus, d'autres opérateurs spatiaux peuvent être spécifiés avec l'argument `op =`, comme le montre la commande suivante qui renvoie l'opposé de `st_intersects()`, les points qui ne sont pas en intersection avec `Canterbury` (voir la section \@ref(topological-relations)) :

```{r 04-spatial-operations-4, eval=FALSE}
nz_height[canterbury, , op = st_disjoint]
```

```{block2 04-spatial-operations-5, type='rmdnote'}
Notez que l´argument vide --- dénoté par `, ,` --- dans l´extrait de code précédent est inclus pour mettre en évidence `op`, le troisième argument dans `[` pour les objets `sf`.
On peut l´utiliser pour modifier l´opération de sélection de plusieurs façons.
`nz_height[canterbury, 2, op = st_disjoint]`, par exemple, retourne les mêmes lignes mais n´inclut que la deuxième colonne d´attributs (voir ``sf:::`[.sf`` et le `?sf`` pour plus de détails).
```

Pour de nombreuses applications, c'est tout ce que vous aurez besoin de savoir sur les sélections spatiales avec les données vectorielles !
Si vous êtes impatient d'en savoir plus sur les relations topologiques, au-delà de `st_intersects()` et `st_disjoint()`, passez à la section suivante (\@ref(topological-relations)).
Si vous êtes intéressé par les détails, y compris les autres façons de faire des sélections, c'est par ici.

Une autre façon d'effectuer une sélection spatiale est d'utiliser les objets retournés par les opérateurs topologiques.
Ces objets peuvent être utiles en soi, par exemple lors de l'exploration du réseau de relations entre des régions contiguës, mais ils peuvent également être utilisés pour sélectionner comme le montre le morceau de code ci-dessous :

```{r 04-spatial-operations-6}
sel_sgbp = st_intersects(x = nz_height, y = canterbury)
class(sel_sgbp)
sel_sgbp
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]
```

Le code ci-dessus crée un objet de classe `sgbp` (un prédicat binaire de géométrie "creuse", une liste de longueur `x` dans l'opération spatiale) et le convertit ensuite en un vecteur logique `sel_logical` (contenant seulement les valeurs `TRUE` et `FALSE`, quelque chose qui peut aussi être utilisé par la fonction filtre de **dplyr**).
\index{binary predicate|seealso {topological relations}}
La fonction `lengths()` identifie les éléments de `nz_height` qui ont une intersection avec *tout* objet de `y`.
Dans ce cas, 1 est la plus grande valeur possible, mais pour des opérations plus complexes, on peut utiliser la méthode pour sélectionner uniquement les caractéristiques qui ont une intersection avec, par exemple, 2 caractéristiques ou plus de l'objet source.

```{block2 04-spatial-operations-7, type='rmdnote'}
Note : une autre façon de retourner une sortie logique est de mettre `sparse = FALSE` (ce qui signifie retourner une matrice dense et non une matrice 'creuse') dans des opérateurs tels que `st_intersects()`. La commande `st_intersects(x = nz_height, y = canterbury, sparse = FALSE)[, 1]`, par exemple, retournerait une sortie identique à `sel_logical`.
Note : la solution impliquant les objets `sgbp` est cependant plus généralisable, car elle fonctionne pour les opérations *many-to-many* et a des besoins en mémoire plus faibles.
```

Le même résultat peut être obtenu avec la fonction de **sf** `st_filter()` qui a été [créée](https://github.com/r-spatial/sf/issues/1148) pour augmenter la compatibilité entre les objets `sf` et les manipulation de données de **dplyr** :

```{r}
canterbury_height3 = nz_height |>
  st_filter(y = canterbury, .predicate = st_intersects)
```

<!--toDo:jn-->
<!-- fix pipes -->

```{r 04-spatial-operations-7b-old, eval=FALSE, echo=FALSE}
# Additional tests of subsetting
canterbury_height4 = nz_height |>
  filter(st_intersects(x = _, y = canterbury, sparse = FALSE))
canterbury_height5 = nz_height |>
  filter(sel_logical)
identical(canterbury_height3, canterbury_height4)
identical(canterbury_height3, canterbury_height5)
identical(canterbury_height2, canterbury_height4)
identical(canterbury_height, canterbury_height4)
waldo::compare(canterbury_height2, canterbury_height4)
```

A ce stade, il y a trois versions identiques (à l'exception des noms de lignes) de `canterbury_height`, une créée en utilisant l'opérateur `[`, une créée via un objet de sélection intermédiaire, et une autre utilisant la fonction de commodité de **sf** `st_filter()`.
<!-- RL: commented out for now as old. Todo: if we ever update that vignette uncomment the next line. -->
<!-- To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). -->
La section suivante explore différents types de relations spatiales, également connues sous le nom de prédicats binaires, qui peuvent être utilisées pour identifier si deux éléments sont spatialement liés ou non.

### Relations topologiques

Les relations topologiques décrivent les relations spatiales entre les objets.
Les "relations topologiques binaires", pour leur donner leur nom complet, sont des énoncés logiques (en ce sens que la réponse ne peut être que `VRAI` ou `FAUX`) sur les relations spatiales entre deux objets définis par des ensembles ordonnés de points (formant typiquement des points, des lignes et des polygones) en deux dimensions ou plus [@egenhofer_mathematical_1990].
Cela peut sembler plutôt abstrait et, en effet, la définition et la classification des relations topologiques reposent sur des fondements mathématiques publiés pour la première fois sous forme de livre en 1966 [@spanier_algebraic_1995], le domaine de la topologie algébrique se poursuivant au 21^e^ siècle [@dieck_algebraic_2008].

Malgré leur origine mathématique, les relations topologiques peuvent être comprises intuitivement en se référant à des visualisations de fonctions couramment utilisées qui testent les types courants de relations spatiales.
La figure \@ref(fig:relations) montre une variété de paires géométriques et leurs relations associées.
Les troisième et quatrième paires de la figure \@ref(fig:relations) (de gauche à droite puis vers le bas) montrent que, pour certaines relations, l'ordre est important : alors que les relations *equals*, *intersects*, *crosses*, *touches* et *overlaps* sont symétriques, ce qui signifie que si `function(x, y)` est vraie, `function(y, x)` le sera aussi, les relations dans lesquelles l'ordre des géométries est important, comme *contains* et *within*, ne le sont pas.
Remarquez que chaque paire de géométries possède une chaîne "DE-9IM" telle que FF2F11212, décrite dans la section suivante.
\index{topological relations}

```{r relations, echo=FALSE, fig.cap="Relations topologiques entre géométries vectorielles, inspirées des figures 1 et 2 d'Egenhofer et Herring (1990). Les relations pour lesquelles la fonction(x, y) est vraie sont imprimées pour chaque paire de géométries, x étant représenté en rose et y en bleu. La nature de la relation spatiale pour chaque paire est décrite par la chaîne de caractères du Dimensionally Extended 9-Intersection Model.", fig.show='hold', message=FALSE, fig.asp=0.66, warning=FALSE}
# source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R")
source("code/de_9im.R")
library(sf)
xy2sfc = function(x, y) st_sfc(st_polygon(list(cbind(x, y))))
p1 = xy2sfc(x = c(0, 0, 1, 1,   0), y = c(0, 1, 1, 0.5, 0))
p2 = xy2sfc(x = c(0, 1, 1, 0), y = c(0, 0, 0.5, 0))
p3 = xy2sfc(x = c(0, 1, 1, 0), y = c(0, 0, 0.7, 0))
p4 = xy2sfc(x = c(0.7, 0.7, 0.9, 0.7), y = c(0.8, 0.5, 0.5, 0.8))
p5 = xy2sfc(x = c(0.6, 0.7, 1, 0.6), y = c(0.7, 0.5, 0.5, 0.7))
p6 = xy2sfc(x = c(0.1, 1, 1, 0.1), y = c(0, 0, 0.3, 0))
p7 = xy2sfc(x = c(0.05, 0.05, 0.6, 0.5, 0.05), y = c(0.4, 0.97, 0.97, 0.4, 0.4))

# todo: add 3 more with line/point relations?
tmap::tmap_arrange(de_9im(p1, p2), de_9im(p1, p3), de_9im(p1, p4),
                   de_9im(p7, p1), de_9im(p1, p5), de_9im(p1, p6), nrow = 2)
```

Dans `sf`, les fonctions testant les différents types de relations topologiques sont appelées binary predicates", comme décrit dans la vignette *Manipulating Simple Feature Geometries*, qui peut être consultée avec la commande [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), et dans la page d'aide [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html).
Pour voir comment les relations topologiques fonctionnent en pratique, créons un exemple simple et reproductible, en nous appuyant sur les relations illustrées dans la Figure \@ref(fig:relations) et en consolidant les connaissances sur la représentation des géométries vectorielles acquises dans un chapitre précédent (Section \@ref(geometry)).
Notez que pour créer des données tabulaires représentant les coordonnées (x et y) des sommets du polygone, nous utilisons la fonction R de base `cbind()` pour créer une matrice représentant les points de coordonnées, un `POLYGON`, et enfin un objet `sfc`, comme décrit au chapitre \@ref(spatial-class)) :

```{r}
polygon_matrix = cbind(
  x = c(0, 0, 1, 1,   0),
  y = c(0, 1, 1, 0.5, 0)
)
polygon_sfc = st_sfc(st_polygon(list(polygon_matrix)))
```

Nous allons créer des géométries supplémentaires pour démontrer les relations spatiales à l'aide des commandes suivantes qui, lorsqu'elles sont tracées sur le polygone créé ci-dessus, se rapportent les unes aux autres dans l'espace, comme le montre la Figure \@ref(fig:relation-objects).
Notez l'utilisation de la fonction `st_as_sf()` et de l'argument `coords` pour convertir efficacement un tableau de données contenant des colonnes représentant des coordonnées en un objet `sf` contenant des points :

```{r}
line_sfc = st_sfc(st_linestring(cbind(
  x = c(0.4, 1),
  y = c(0.2, 0.5)
)))
# créer des points
point_df = data.frame(
  x = c(0.2, 0.7, 0.4),
  y = c(0.1, 0.2, 0.8)
)
point_sf = st_as_sf(point_df, coords = c("x", "y"))
```

```{r relation-objects, echo=FALSE, fig.cap="Points (`point_df` 1 à 3), ligne et polygones arrangés pour illustrer les relations topologiques.", fig.asp=1, out.width="50%", fig.scap="Exemples des relations topologiques."}
par(pty = "s")
plot(polygon_sfc, border = "red", col = "gray", axes = TRUE)
plot(line_sfc, lwd = 5, add = TRUE)
plot(point_sf, add = TRUE, lab = 1:4, cex = 2)
text(point_df[, 1] + 0.02, point_df[, 2] + 0.04, 1:3, cex = 1.3)
```

Une première question simple pourrait être : quels sont les points de `point_sf` en intersection avec le polygone `polygon_sfc` ?
On peut répondre à cette question par inspection (les points 1 et 3 sont respectivement en contact et à l'intérieur du polygone).
On peut répondre à cette question avec le prédicat spatial `st_intersects()` comme suit:

```{r 04-spatial-operations-9, eval=FALSE}
st_intersects(point_sf, polygon_sfc)
#> Sparse geometry binary predicate... `intersects'
#>  1: 1
#>  2: (empty)
#>  3: 1
```

Le résultat devrait correspondre à votre intuition :
des résultats positifs (`1`) sont retournés pour le premier et le troisième point, et un résultat négatif (représenté par un vecteur vide "(empty)") pour le deuxième  en dehors de la frontière du polygone.
Ce qui peut être inattendu, c'est que le résultat se présente sous la forme d'une liste de vecteurs.
Cette sortie *matrice creuse* n'enregistre une relation que si elle existe, ce qui réduit les besoins en mémoire des opérations topologiques sur les objets avec de nombreuses entités.
Comme nous l'avons vu dans la section précédente, une *matrice dense* composée de valeurs `TRUE` ou `FALSE` est retournée lorsque `sparse = FALSE` :

```{r 04-spatial-operations-10}
st_intersects(point_sf, polygon_sfc, sparse = FALSE)
```

Dans la sortie ci-dessus, chaque ligne représente un élément dans l'objet cible (l'argument `x`) et chaque colonne représente un élément dans l'objet de sélection (`y`). 
Dans ce cas, il n'y a qu'un seul élément dans l'objet `y` `polygon_sfc`, donc le résultat, qui peut être utilisé pour la sélection comme nous l'avons vu dans la section \@ref(spatial-subsetting), n'a qu'une seule colonne.

`st_intersects()` renvoie `TRUE` même dans les cas où les éléments se touchent juste : *intersects* est une opération topologique "fourre-tout" qui identifie de nombreux types de relations spatiales, comme l'illustre la figure \@ref(fig:relations).
Il y a des questions plus restrictives, par exemple  : quels sont les points situés à l'intérieur du polygone, et quelles sont les caractéristiques qui sont sur ou qui contiennent une frontière partagée avec `y` ?
On peut répondre à ces questions de la manière suivante (résultats non montrés) :

```{r 04-spatial-operations-9-2, eval=FALSE}
st_within(point_sf, polygon_sfc)
st_touches(point_sf, polygon_sfc)
```

Notez que bien que le premier point *touche* la limite du polygone, il n'est pas à l'intérieur de celui-ci ; le troisième point est à l'intérieur du polygone mais ne touche aucune partie de sa frontière.
L'opposé de `st_intersects()` est `st_disjoint()`, qui retourne uniquement les objets qui n'ont aucun rapport spatial avec l'objet sélectionné (ici `[, 1]` convertit le résultat en vecteur) :

```{r 04-spatial-operations-11}
st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1]
```

La fonction `st_is_within_distance()` détecte les éléments qui  touchent *presque* l'objet de sélection. La fonction a un argument supplémentaire `dist`.
Il peut être utilisé pour définir la distance à laquelle les objets cibles doivent se trouver avant d'être sélectionnés.
Remarquez que bien que le point 2 soit à plus de 0,2 unités de distance du sommet le plus proche de `polygon_sfc`, il est quand même sélectionné lorsque la distance est fixée à 0,2. 
En effet, la distance est mesurée par rapport à l'arête la plus proche, dans ce cas la partie du polygone qui se trouve directement au-dessus du point 2 dans la figure \@ref(fig:relation-objets).
(Vous pouvez vérifier que la distance réelle entre le point 2 et le polygone est de 0,13 avec la commande `st_distance(point_sf, polygon_sfc)`).
Le prédicat spatial binaire "is within distance" (es à distance de) est démontré dans l'extrait de code ci-dessous, dont les résultats montrent que chaque point est à moins de 0,2 unité du polygone :

```{r 04-spatial-operations-14}
st_is_within_distance(point_sf, polygon_sfc, dist = 0.2, sparse = FALSE)[, 1]
```

```{r, eval=FALSE, echo=FALSE}
# verify distances to the polygon with reference to paragraph above:
st_distance(point_sf, polygon_sfc)
#           [,1]
# [1,] 0.0000000
# [2,] 0.1341641
# [3,] 0.0000000
```


```{block2 04-spatial-operations-15, type='rmdnote'}
Les fonctions de calcul des relations topologiques utilisent des indices spatiaux pour accélérer considérablement les performances des requêtes spatiales.
Elles y parviennent en utilisant l´algorithme *Sort-Tile-Recursive* (STR).
La fonction `st_join`, mentionnée dans la section suivante, utilise également l´indexation spatiale. 
Vous pouvez en savoir plus à l´adresse suivante https://www.r-spatial.org/r/2017/06/22/spatial-index.html.
```

```{r 04-spatial-operations-16, eval=FALSE, echo=FALSE}
# other tests
st_overlaps(point_sf, polygon_sfc, sparse = FALSE)
st_covers(point_sf, polygon_sfc, sparse = FALSE)
st_covered_by(point_sf, polygon_sfc, sparse = FALSE)
```

```{r 04-spatial-operations-17, eval=FALSE, echo=FALSE}
st_contains(a, p[2, ], sparse = TRUE)
```

```{r 04-spatial-operations-18, eval=FALSE, echo=FALSE}
# starting simpler so commented
a1 = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a2 = st_polygon(list(rbind(c(2, 0), c(2, 2), c(3, 2), c(3, 0), c(2, 0))))
a = st_sfc(a1, a2)

b1 = a1 * 0.5
b2 = a2 * 0.4 + c(1, 0.5)
b = st_sfc(b1, b2)

l1 = st_linestring(x = matrix(c(0, 3, -1, 1), , 2))
l2 = st_linestring(x = matrix(c(-1, -1, -0.5, 1), , 2))
l = st_sfc(l1, l2)

p = st_multipoint(x = matrix(c(0.5, 1, -1, 0, 1, 0.5), , 2))

plot(a, border = "red", axes = TRUE)
plot(b, border = "green", add = TRUE)
plot(l, add = TRUE)
plot(p, add = TRUE)
```

### Les chaines DE-9IM

Les prédicats binaires présentés dans la section précédente reposent sur le modèle *Dimensionally Extended 9-Intersection Model* (DE-9IM).
Comme le suggère son nom cryptique, ce n'est pas un sujet facile.
Y consacrer du temps peut être utile afin d'améliorer notre compréhension des relations spatiales.
En outre, les utilisations avancées de DE-9IM incluent la création de prédicats spatiaux personnalisés.
Ce modèle était à l'origine intitulé " DE + 9IM " par ses inventeurs, en référence à la " dimension des intersections des limites, des intérieurs et des extérieurs de deux entités " [@clementini_comparison_1995], mais il est désormais désigné par DE-9IM [@shen_classification_2018].
<!-- The model's workings can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). -->

```{r de-9im, echo=FALSE, eval=FALSE}
# Todo one day: revive this
b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b = st_buffer(b, dist = 1) # convert points to circles
bsf = sf::st_sf(data.frame(Object = c("a", "b")), geometry = b)
b9 = replicate(bsf, n = 9, simplify = FALSE)
b9sf = do.call(rbind, b9)
domains = c("Interior", "Boundary", "Exterior")
b9sf$domain_a = rep(rep(domains, 3), each = 2)
b9sf$domain_b = rep(rep(domains, each = 3), each = 2)
library(ggplot2)
ggplot(b9sf) +
  geom_sf() +
  facet_grid(domain_a ~ domain_b)

plot(b9sf)
tmap_arrange(
  tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Interior"),
  tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Boundary"),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  tm_shape(b) + tm_polygons(alpha = 0.5),
  nrow = 3
)

plot(b)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text
```

Pour démontrer le fonctionnement des chaînes DE-9IM, examinons les différentes façons dont la première paire de géométries peut-être reliée dans la figure \@ref(fig:relations).
La figure \@ref(fig:de9imgg) illustre le modèle à 9 intersections (9IM).  Elle montre les intersections entre chaque combinaison possible entre l'intérieur, la limite et l'extérieur de chaque objet. Chaque composant du premier objet `x` est disposé en colonnes et que chaque composant de `y` est disposé en lignes, un graphique à facettes est créé avec les intersections entre chaque élément mises en évidence.

```{r de9imgg, echo=FALSE, warning=FALSE, fig.cap="Illustration du fonctionnement du Modèle Dimensionnel Étendu à 9 Intersections (DE-9IM). Les couleurs qui ne figurent pas dans la légende représentent le chevauchement entre les différentes composantes. Les lignes épaisses mettent en évidence les intersections bidimensionnelles, par exemple entre la limite de l'objet x et l'intérieur de l'objet y, illustrées dans la facette supérieure du milieu.", message=FALSE}
p1_2 = st_as_sf(c(p1, p3))
ii = st_as_sf(st_intersection(p1, p3))
ii$Object = "Intersection"
ii$domain_a = "Intérieur"
ii$domain_b = "Intérieur"

bi = st_sf(x = st_intersection(
  st_cast(p1, "LINESTRING"),
  st_difference(p3, st_buffer(st_cast(p3, "LINESTRING"), dist = 0.01))
  ))
bi = st_buffer(bi, dist = 0.01)
bi$Object = "Intersection"
bi$domain_a = "Limite"
bi$domain_b = "Intérieur"

ei = st_sf(x = st_difference(p3, p1))
ei$Object = "Intersection"
ei$domain_a = "Extérieur"
ei$domain_b = "Intérieur"

ib = st_sf(x = st_intersection(
  st_cast(p3, "LINESTRING"),
  st_difference(p1, st_buffer(st_cast(p1, "LINESTRING"), dist = 0.005))
  ))
ib = st_buffer(ib, dist = 0.01)
ib$Object = "Intersection"
ib$domain_a = "Intérieur"
ib$domain_b = "Limite"

bb = st_cast(ii, "POINT")
bb_line = st_sf(x = st_sfc(st_linestring(matrix(c(1, 0.5, 1, 0.7), nrow = 2, byrow = TRUE))))
bb_line_buffer = st_buffer(bb_line, dist = 0.01)
bb_buffer = st_buffer(bb, dist = 0.01)
bb = st_union(bb_buffer, bb_line_buffer)
bb$Object = "Intersection"
bb$domain_a = "Limite"
bb$domain_b = "Limite"

eb = st_sf(x = st_difference(
  st_cast(p3, "LINESTRING"),
  p1
  ))
eb = st_buffer(eb, dist = 0.01)
eb$Object = "Intersection"
eb$domain_a = "Extérieur"
eb$domain_b = "Limite"

ie = st_sf(x = st_difference(p1, p3))
ie$Object = "Intersection"
ie$domain_a = "Intérieur"
ie$domain_b = "Extérieur"

be = st_sf(x = st_difference(
  st_cast(p1, "LINESTRING"),
  p3
  ))
be = st_buffer(be, dist = 0.01)
be$Object = "Intersection"
be$domain_a = "Limite"
be$domain_b = "Extérieur"

ee = st_sf(x = st_difference(
  st_buffer(st_union(p1, p3), 0.02),
  st_union(p1, p3)
  ))
ee$Object = "Intersection"
ee$domain_a = "Extérieur"
ee$domain_b = "Extérieur"

b9 = replicate(p1_2, n = 9, simplify = FALSE)
b9sf = do.call(rbind, b9)
b9sf$Object = rep(c("x", "y"), 9)
domains = c("Intérieur", "Limite", "Extérieur")
b9sf$domain_a = rep(rep(domains, 3), each = 2)
b9sf$domain_b = rep(rep(domains, each = 3), each = 2)
b9sf = rbind(b9sf, ii, bi, ei, ib, bb, eb, ie, be, ee)
b9sf$domain_a = ordered(b9sf$domain_a, levels = c("Intérieur", "Limite", "Extérieur"))
b9sf$domain_b = ordered(b9sf$domain_b, levels = c("Intérieur", "Limite", "Extérieur"))
b9sf = b9sf %>% 
  mutate(alpha = case_when(
   Object == "x" ~ 0.1, 
   Object == "y" ~ 0.1, 
   TRUE ~ 0.2 
  ))
library(ggplot2)
ggplot(b9sf) +
  geom_sf(aes(fill = Object, alpha = alpha)) +
  facet_grid(domain_b ~ domain_a) +
  scale_fill_manual(values = c("red", "lightblue", "yellow"), position = "top", name = "") +
  scale_alpha_continuous(range = c(0.3, 0.9)) +
  guides(alpha = "none") +
  theme_void() +
  theme(legend.position = "top")
```

Les chaînes DE-9IM sont dérivées de la dimension de chaque type de relation.
Dans ce cas, les intersections rouges de la figure \@ref(fig:de9imgg) ont des dimensions de 0 (points), 1 (lignes) et 2 (polygones), comme le montre le tableau \@ref(tab:de9emtable).

```{r de9emtable, echo=FALSE}
# See https://github.com/Robinlovelace/geocompr/issues/699
pattern = st_relate(p1, p3)
matrix_de_9im = function(pattern) {
    string = unlist(strsplit(pattern , ""))
    matrix_de_9im = matrix(string, nrow = 3, byrow = TRUE)
    colnames(matrix_de_9im) = c("I", "B", "E")
    row.names(matrix_de_9im) = c("I", "B", "E")
    return(matrix_de_9im)
}

m = matrix_de_9im(pattern)
colnames(m) = c("Intérieur (x)", "Limite (x)", "Extérieur (x)")
rownames(m) = c("Intérieur (y)", "Limite (y)", "Extérieur (y)")
knitr::kable(m, caption = "Tableau montrant les relations entre les intérieurs, les limites et les extérieurs des géométries x et y.")
```

En aplatissant cette matrice "ligne par ligne" (c'est-à-dire en concaténant la première ligne, puis la deuxième, puis la troisième), on obtient la chaîne `212111212`.
Un autre exemple va permettre d'expliciter ce système :
la relation représentée sur la figure \@ref(fig:relations) (la troisième paire de polygones dans la troisième colonne et la première ligne) peut être définie dans le système DE-9IM comme suit :

- Les intersections entre l'*intérieur* du grand objet `x` et l'intérieur, la limite et l'extérieur de `y` ont des dimensions respectives de 2, 1 et 2.
- Les intersections entre la *frontière* du grand objet `x` et l'intérieur, la frontière et l'extérieur de `y` ont des dimensions respectives de F, F et 1, où "F" signifie "faux", les objets sont disjoints.
- Les intersections entre l'*extérieur* de `x` et l'intérieur, la limite et l'extérieur de `y` ont des dimensions respectives de F, F et 2 : l'extérieur du plus grand objet ne touche pas l'intérieur ou la limite de `y`, mais l'extérieur du plus petit et du plus grand objet couvre la même surface.

Ces trois composants, une fois concaténés, créent la chaîne `212`, `FF1`, et `FF2`.
C'est le même résultat que celui obtenu par la fonction `st_relate()` (voir le code source de ce chapitre pour voir comment les autres géométries de la figure \@ref(fig:relations) ont été créées) :

```{r}
xy2sfc = function(x, y) st_sfc(st_polygon(list(cbind(x, y))))
x = xy2sfc(x = c(0, 0, 1, 1,   0), y = c(0, 1, 1, 0.5, 0))
y = xy2sfc(x = c(0.7, 0.7, 0.9, 0.7), y = c(0.8, 0.5, 0.5, 0.8))
st_relate(x, y)
```

La compréhension des chaînes DE-9IM permet de développer de nouveaux prédicats spatiaux binaires.
La page d'aide `?st_relate` contient des définitions de fonctions pour les relations "reine" et "tour" dans lesquelles les polygones partagent une frontière ou seulement un point, respectivement.
Les relations "reine" signifient que les relations "frontière-frontière" (la cellule de la deuxième colonne et de la deuxième ligne de la table \@ref(tab:de9emtable), ou le cinquième élément de la chaîne DE-9IM) ne doivent pas être vides, ce qui correspond au *pattern*  `F***T****`, tandis que pour les relations "tour", le même élément doit être 1 (ce qui signifie une intersection linéaire).
Ces relations sont implémentées comme suit :

```{r}
st_queen = function(x, y) st_relate(x, y, pattern = "F***T****")
st_rook = function(x, y) st_relate(x, y, pattern = "F***1****")
```

A partir de l'objet `x` créé précédemment, nous pouvons utiliser les fonctions nouvellement créées pour trouver quels éléments de la grille sont une 'reine' et une 'tour' par rapport à la case centrale de la grille comme suit :

```{r queenscode, fig.show='hide'}
grid = st_make_grid(x, n = 3)
grid_sf = st_sf(grid)
grid_sf$queens = lengths(st_queen(grid, grid[5])) > 0
plot(grid, col = grid_sf$queens)
grid_sf$rooks = lengths(st_rook(grid, grid[5])) > 0
plot(grid, col = grid_sf$rooks)
```

```{r queens, fig.cap="Démonstration de prédicats spatiaux binaires personnalisés permettant de trouver les relations 'reine' (à gauche) et 'tour' (à droite) par rapport à la case centrale dans une grille à 9 géométries.", echo=FALSE, warning=FALSE}
tm_shape(grid_sf) +
  tm_fill(col = c("queens", "rooks"), palette = c("white", "black")) +
  tm_shape(grid_sf) +
  tm_borders(col = "grey", lwd = 2) +
  tm_layout(frame = FALSE, legend.show = FALSE,
            panel.labels = c("reine", "fou"))
```


<!-- Another of a custom binary spatial predicate is 'overlapping lines' which detects lines that overlap for some or all of another line's geometry. -->
<!-- This can be implemented as follows, with the pattern signifying that the intersection between the two line interiors must be a line: -->

```{r, echo=FALSE, eval=FALSE}
st_lineoverlap = function(x, y) st_relate(x, y, pattern = "T*1******")
line1 = st_sfc(st_linestring(cbind(
  x = c(0, 0.8),
  y = c(0, 0)
)))
line2 = st_sfc(st_linestring(cbind(
  x = c(0.1, 0.5),
  y = c(0, 0)
)))
line3 = st_sfc(st_linestring(cbind(
  x = c(0, 0.5),
  y = c(0, 0.2)
)))
st_queen(line1, line2)
st_relate(line1, line2)
st_relate(line1, line3)
st_lineoverlap(line1, line2)
st_lineoverlap(line1, line3)
de_9im(line1, line2)
# test the function
rnet = pct::get_pct_rnet(region = "isle-of-wight")
osm_net = osmextract::oe_get_network(place = "isle-of-wight", mode = "driving")
sel = st_relate(rnet, osm_net, pattern = "T*1******")
summary(lengths(sel) > 0)
rnet_joined1 = st_join(rnet, osm_net, join = st_lineoverlap)
rnet_joined2 = st_join(rnet, osm_net, join = st_relate, pattern = "T*1******")
rnet_joined3 = st_join(rnet, osm_net)
summary(is.na(rnet_joined1$osm_id))
summary(is.na(rnet_joined2$osm_id))
summary(is.na(rnet_joined3$osm_id))
sel_relates = st_relate(rnet[1, ], osm_net)
dim(sel_relates)
sel_table = table(sel_relates)
sel_table
dim(sel_table)
sel_restrictive = sel_relates[1, ] == "0F1FF0102"
summary(sel_restrictive)
nrow(osm_net)
mapview::mapview(rnet[1, ]) + mapview::mapview(osm_net[sel_restrictive, ])

rnet_approx = rnet
st_precision(rnet_approx) = 100
head(st_coordinates(rnet_approx))

sel_relates = st_relate(rnet_approx[1, ], osm_net)
dim(sel_relates)
sel_table = table(sel_relates)
sel_table

```

### Jointure spatiale 

La jointure de deux jeux de données non spatiales repose sur une variable "clé" partagée, comme décrit dans la section \@ref(vector-attribute-joining).
La jointure de données spatiales applique le même concept, mais s'appuie sur les relations spatiales, décrites dans la section précédente.
Comme pour les données attributaires, la jointure ajoute de nouvelles colonnes à l'objet cible (l'argument `x` dans les fonctions de jointure), à partir d'un objet source (`y`).
\index{join!spatial}
\index{spatial!join}

Le processus est illustré par l'exemple suivant : imaginez que vous disposez de dix points répartis au hasard sur la surface de la Terre et que vous demandez, pour les points qui se trouvent sur la terre ferme, dans quels pays se trouvent-ils ?
La mise en œuvre de cette idée dans un [exemple reproductible] (https://github.com/Robinlovelace/geocompr/blob/main/code/04-spatial-join.R) renforcera vos compétences en matière de traitement des données géographiques et vous montrera comment fonctionnent les jointures spatiales.
Le point de départ consiste à créer des points dispersés de manière aléatoire sur la surface de la Terre :

```{r 04-spatial-operations-19}
set.seed(2018) # définir la seed pour la reproductibilité
(bb = st_bbox(world)) # les limites de la terre
random_df = data.frame(
  x = runif(n = 10, min = bb[1], max = bb[3]),
  y = runif(n = 10, min = bb[2], max = bb[4])
)
random_points = random_df |> 
  st_as_sf(coords = c("x", "y")) |> # définir les coordonnés
  st_set_crs("EPSG:4326") # définir le CRS
```

Le scénario illustré dans la Figure \@ref(fig:spatial-join) montre que l'objet `random_points` (en haut à gauche) n'a pas d'attributs, alors que le `world` (en haut à droite) a des attributs, y compris les noms de pays indiqués pour un échantillon de pays dans la légende.
Les jointures spatiales sont implémentées avec `st_join()`, comme illustré dans l'extrait de code ci-dessous.
La sortie est l'objet `random_joined` qui est illustré dans la Figure \@ref(fig:spatial-join) (en bas à gauche).
Avant de créer l'ensemble de données jointes, nous utilisons la sélection spatiale pour créer `world_random`, qui contient uniquement les pays qui contiennent des points aléatoires, afin de vérifier que le nombre de noms de pays retournés dans l'ensemble de données jointes doit être de quatre (cf. le panneau supérieur droit de la Figure \@ref(fig:spatial-join)).

```{r 04-spatial-operations-20, message=FALSE}
world_random = world[random_points, ]
nrow(world_random)
random_joined = st_join(random_points, world["name_long"])
```

```{r spatial-join, echo=FALSE, fig.cap="Illustration d'une jointure spatiale. Une nouvelle variable attributaire est ajoutée aux points aléatoires (en haut à gauche) de l'objet monde source (en haut à droite), ce qui donne les données représentées dans le dernier panneau.", fig.asp=0.5, warning=FALSE, message=FALSE, out.width="100%", fig.scap="Illustration d'une jointure spatiale."}
# source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-spatial-join.R")
source("code/04-spatial-join.R")
tmap_arrange(jm1, jm2, jm3, jm4, nrow = 2, ncol = 2)
```

Par défaut, `st_join()` effectue une jointure à gauche (*left join*), ce qui signifie que le résultat est un objet contenant toutes les lignes de `x`, y compris les lignes sans correspondance dans `y` (voir la section \@ref(vector-attribute-joining)), mais il peut également effectuer des jointures internes en définissant l'argument `left = FALSE`.
Comme pour les sélections spatiales, l'opérateur topologique par défaut utilisé par `st_join()` est `st_intersects()`, qui peut être modifié en définissant l'argument `join` (cf. `?st_join` pour plus de détails).
L'exemple ci-dessus montre l'ajout d'une colonne d'une couche de polygones à une couche de points, mais la même approche fonctionne indépendamment des types de géométrie.
Dans de tels cas, par exemple lorsque `x` contient des polygones, dont chacun correspond à plusieurs objets dans `y`, les jointures spatiales résulteront en des caractéristiques dupliquées, crée une nouvelle ligne pour chaque correspondance dans `y`.

<!-- Idea: demonstrate what happens when there are multiple matches with reprex (low priority, RL: 2021-12) -->

### Jointure sans chevauchement

Parfois, deux jeux de données géographiques ne se touchent pas mais ont quand même une forte relation géographique.
Les jeux de données `cycle_hire` et `cycle_hire_osm`, présent dans le paquet **spData**, en sont un bon exemple.
Leur tracé montre qu'ils sont souvent étroitement liés mais qu'ils ne se touchent pas, comme le montre la figure \@ref(fig:cycle-hire), dont une version de base est créée avec le code suivant ci-dessous :
\index{join!non-overlapping}

```{r 04-spatial-operations-21, eval=FALSE}
plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")
```

Nous pouvons vérifier si certains points se superposent  avec `st_intersects()` :

```{r 04-spatial-operations-22, message=FALSE}
any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
```

```{r 04-spatial-operations-23, echo=FALSE, eval=FALSE}
# include pour montrer d'autres façon de tester le chevauchement
sum(st_geometry(cycle_hire) %in% st_geometry(cycle_hire_osm))
sum(st_coordinates(cycle_hire)[, 1] %in% st_coordinates(cycle_hire_osm)[, 1])
```

```{r cycle-hire, fig.cap="La distribution spatiale des points de location de vélos à Londres, basée sur les données officielles (bleu) et les données OpenStreetMap (rouge).", echo=FALSE, warning=FALSE, fig.scap="La répartition spatiale des points de location de vélos à Londres."}
if (knitr::is_latex_output()){
  knitr::include_graphics("figures/cycle-hire-1.png")
} else if (knitr::is_html_output()){
  # library(tmap)
  # osm_tiles = tmaptools::read_osm(tmaptools::bb(cycle_hire, ext = 1.3), type =   "https://korona.geog.uni-heidelberg.de/tiles/roadsg/x={x}&y={y}&z={z}")
  # qtm(osm_tiles) +
    # tm_shape(cycle_hire) +
    # tm_bubbles(col = "blue", alpha = 0.5, size = 0.2) +
    # tm_shape(cycle_hire_osm) +
    # tm_bubbles(col = "red", alpha = 0.5, size = 0.2) +
    # tm_scale_bar()
  library(leaflet)
  leaflet() |>
    # addProviderTiles(providers$OpenStreetMap.BlackAndWhite) %>%
    addCircles(data = cycle_hire) %>%
    addCircles(data = cycle_hire_osm, col = "red")  
}
```

Imaginons que nous ayons besoin de joindre la variable "capacité" de `cycle_hire_osm` aux données officielles "cible" contenues dans `cycle_hire`.
Dans ce cas, une jointure sans chevauchement est nécessaire.
La méthode la plus simple est d'utiliser l'opérateur topologique `st_is_within_distance()`, comme démontré ci-dessous en utilisant une distance seuil de 20 m (notez que cela fonctionne avec des données projetées et non projetées).

```{r 04-spatial-operations-24}
sel = st_is_within_distance(cycle_hire, cycle_hire_osm, dist = 20)
summary(lengths(sel) > 0)
```

```{r 04-spatial-operations-24-without-s2-test, eval=FALSE, echo=FALSE}
sf::sf_use_s2(FALSE)
sel = st_is_within_distance(cycle_hire, cycle_hire_osm, dist = 20)
summary(lengths(sel) > 0)
# still works: must be lwgeom or some other magic!
```


```{r 04-spatial-operations-24-projected, eval=FALSE, echo=FALSE}
# This chunk contains the non-overlapping join on projected data, a step that is no longer needed:
# Note that, before performing the relation, both objects are transformed into a projected CRS.
# These projected objects are created below (note the affix `_P`, short for projected):
cycle_hire_P = st_transform(cycle_hire, 27700)
cycle_hire_osm_P = st_transform(cycle_hire_osm, 27700)
sel = st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20)
summary(lengths(sel) > 0)
```

Cela montre qu'il y a `r sum(lengths(sel) > 0)` des points dans l'objet cible `cycle_hire` dans la distance seuil (20 m) de `cycle_hire_osm`.
Comment récupérer les *valeurs* associées aux points respectifs de `cycle_hire_osm` ?
La solution est à nouveau avec `st_join()`, mais avec un argument `dist` supplémentaire (fixé à 20 m en dessous) :

```{r 04-spatial-operations-25}
z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, dist = 20)
nrow(cycle_hire)
nrow(z)
```

Remarquez que le nombre de lignes dans le résultat joint est supérieur à la cible.
Cela est dû au fait que certaines stations de location de vélos dans `cycle_hire` ont plusieurs correspondances dans `cycle_hire_osm`.
Pour agréger les valeurs des points qui se chevauchent et renvoyer la moyenne, nous pouvons utiliser les méthodes d'agrégation apprises au chapitre \@ref(attr), ce qui donne un objet avec le même nombre de lignes que la cible :

```{r 04-spatial-operations-26}
z = z %>% 
  group_by(id) %>% 
  summarize(capacity = mean(capacity))
nrow(z) == nrow(cycle_hire)
```

La capacité des stations proches peut être vérifiée en comparant les cartes de la capacité des données source `cycle_hire_osm` avec les résultats dans ce nouvel objet (cartes non montrées) :

```{r 04-spatial-operations-27, eval=FALSE}
plot(cycle_hire_osm["capacity"])
plot(z["capacity"])
```

Le résultat de cette jointure a utilisé une opération spatiale pour modifier les données attributaires associées aux entités simples ; la géométrie associée à chaque entité est restée inchangée.

### Agrégats spatiaux {#spatial-aggr}

Comme pour l'agrégation de données attributaires, l'agrégation de données spatiales *condense* les données : les sorties agrégées comportent moins de lignes que les entrées non agrégées.
Les *fonctions d'agrégation* statistiques, telles que la moyenne ou la somme, résument plusieurs valeurs \index{statistiques} d'une variable et renvoient une seule valeur par *variable de regroupement*.
La section \@ref(vector-attribute-aggregation) a montré comment `aggregate()` et `group_by() |> summarize()` condensent les données basées sur des variables d'attributs, cette section montre comment les mêmes fonctions fonctionnent avec des objets spatiaux.
\index{aggregation!spatial}

Pour revenir à l'exemple de la Nouvelle-Zélande, imaginez que vous voulez connaître la hauteur moyenne des points hauts de chaque région : c'est la géométrie de la source (`y` ou `nz` dans ce cas) qui définit comment les valeurs de l'objet cible (`x` ou `nz_height`) sont regroupées.
Ceci peut être fait en une seule ligne de code avec la méthode `aggregate()` de la base R :

```{r 04-spatial-operations-28}
nz_agg = aggregate(x = nz_height, by = nz, FUN = mean)
```

Le résultat de la commande précédente est un objet `sf` ayant la même géométrie que l'objet d'agrégation (spatiale) (`nz`), ce que vous pouvez vérifier avec la commande `identical(st_geometry(nz), st_geometry(nz_agg))`.
Le résultat de l'opération précédente est illustré dans la Figure \@ref(fig:spatial-aggregation), qui montre la valeur moyenne des caractéristiques de `nz_height` dans chacune des 16 régions de la Nouvelle-Zélande.
Le même résultat peut également être généré en passant la sortie de `st_join()` dans les fonctions 'tidy' `group_by()` et `summarize()` comme suit :

```{r spatial-aggregation, echo=FALSE, fig.cap="Hauteur moyenne des 101 points culminants par régions de la Nouvelle-Zélande.", fig.asp=1, message=FALSE, out.width="50%"}
library(tmap)
tm_shape(nz_agg) +
  tm_fill("elevation", breaks = seq(27, 30, by = 0.5) * 1e2) +
  tm_borders()
```

```{r 04-spatial-operations-29}
nz_agg2 = st_join(x = nz, y = nz_height) |>
  group_by(Name) |>
  summarize(elevation = mean(elevation, na.rm = TRUE))
```

```{r test-tidy-spatial-join, eval=FALSE, echo=FALSE}
plot(nz_agg)
plot(nz_agg2)
# aggregate looses the name of aggregating objects
```

Les entités `nz_agg` résultants ont la même géométrie que l'objet d'agrégation `nz` mais avec une nouvelle colonne résumant les valeurs de `x` dans chaque région en utilisant la fonction `mean()`.
D'autres fonctions peuvent, bien sûr, remplacer `mean()` comme la `median()`, `sd()` ou d'autres fonctions qui retournent une seule valeur par groupe.
Remarque : une différence entre les approches `aggregate()` et `group_by() |> summarize()` est que la première donne des valeurs `NA` pour les noms de régions non correspondantes, tandis que la seconde préserve les noms de régions.
L'approche "tidy" est plus flexible en termes de fonctions d'agrégation et de noms de colonnes des résultats.
Les opérations d'agrégation créant  de nouvelles géométries sont décrites dans la section \@ref(geometry-unions). 


### Jointure de couches sans superposition parfaite  {#incongruent}

La congruence spatiale\index{congruence spatiale} est un concept important lié à l'agrégation spatiale.
Un *objet d'agrégation* (que nous appellerons `y`) est *congruent* avec l'objet cible (`x`) si les deux objets ont des frontières communes.
C'est souvent le cas des données sur les limites administratives, où des unités plus grandes --- telles que les Middle Layer Super Output Areas ([MSOAs](https://www.ons.gov.uk/methodology/geography/ukgeographies/censusgeography)) au Royaume-Uni ou les districts dans de nombreux autres pays européens --- sont composées de nombreuses unités plus petites.

Les objets d'agrégation *non congruent*, en revanche, ne partagent pas de frontières communes avec la cible [@qiu_development_2012].
C'est un problème pour l'agrégation spatiale (et d'autres opérations spatiales) illustrée dans la figure \@ref(fig:areal-example) : agréger le centroïde de chaque sous-zone ne donnera pas de résultats précis.
L'interpolation surfacique résout ce problème en transférant les valeurs d'un ensemble d'unités surfaciques à un autre, à l'aide d'une gamme d'algorithmes comprenant des approches simples de pondération de la surface et des approches plus sophistiquées telles que les méthodes " pycnophylactiques " [@tobler_smooth_1979].

```{r areal-example, echo=FALSE, fig.cap="Illustration des entités surfaciques congruentes (à gauche) et non congruentes (à droite) par rapport à des zones d'agrégation plus grandes (bordures bleues translucides).", fig.asp=0.2, fig.scap="Illustration of congruent and incongruent areal units."}
source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-areal-example.R", print.eval = TRUE)
```

Le paquet **spData** contient un jeux de données nommé `incongruent` (polygones colorés avec des bordures noires dans le panneau de droite de la Figure \@ref(fig:areal-example)) et un ensemble de données nommé `aggregating_zones` (les deux polygones avec la bordure bleue translucide dans le panneau de droite de la Figure \@ref(fig:areal-example)).
Supposons que la colonne `valeur` de `incongruent` se réfère au revenu régional total en millions d'euros.
Comment pouvons-nous transférer les valeurs des neuf polygones spatiaux sous-jacents dans les deux polygones de `aggregating_zones` ?

La méthode la plus simple est l'interpolation spatiale pondérée par la surface, qui transfère les valeurs de l'objet `non congruent` vers une nouvelle colonne dans `aggregating_zones` proportionnellement à la surface de recouvrement : plus l'intersection spatiale entre les caractéristiques d'entrée et de sortie est grande, plus la valeur correspondante est grande.
Ceci est implémenté dans `st_interpolate_aw()`, comme démontré dans le morceau de code ci-dessous.

```{r 04-spatial-operations-30}
iv = incongruent["value"] # garde uniquement les valeurs à transférer
agg_aw = st_interpolate_aw(iv, aggregating_zones, extensive = TRUE)
agg_aw$value
```

Dans notre cas, il est utile d'additionner les valeurs des intersections qui se trouvent dans les zones d'agrégation, car le revenu total est une variable dite spatialement extensive (qui augmente avec la superficie), en supposant que le revenu est réparti uniformément dans les zones plus petites (d'où le message d'avertissement ci-dessus).
Il en va différemment pour les variables spatialement [intensives](https://geodacenter.github.io/workbook/3b_rates/lab3b.html#spatially-extensive-and-spatially-intensive-variables) telles que le revenu *moyen* ou les pourcentages, qui n'augmentent pas avec la superficie.
`st_interpolate_aw()` fonctionne également avec des variables spatialement intensives : mettez le paramètre `extensive` à `FALSE` et il utilisera une moyenne plutôt qu'une fonction de somme pour faire l'agrégation.

### Les relations de distance 

Alors que les relations topologiques sont binaires --- une caractéristique est soit en intersection avec une autre ou non --- les relations de distance sont continues.
La distance entre deux objets est calculée avec la fonction `st_distance()`.
Ceci est illustré dans l'extrait de code ci-dessous, qui trouve la distance entre le point le plus élevé de Nouvelle-Zélande et le centroïde géographique de la région de Canterbury, créé dans la section \@ref(spatial-subsetting) :
\index{sf!distance relations}

```{r 04-spatial-operations-31, warning=FALSE}
nz_highest = nz_height |> slice_max(n = 1, order_by = elevation)
canterbury_centroid = st_centroid(canterbury)
st_distance(nz_highest, canterbury_centroid)
```

Il y a deux choses potentiellement surprenantes dans ce résultat :

- Il a des `unités`, ce qui nous indique que la distance est de 100 000 mètres, et non de 100 000 pouces, ou toute autre mesure de distance.
- Elle est retournée sous forme de matrice, même si le résultat ne contient qu'une seule valeur.

Cette deuxième caractéristique laisse entrevoir une autre fonctionnalité utile de `st_distance()`, sa capacité à retourner des *matrices de distance* entre toutes les combinaisons de caractéristiques dans les objets `x` et `y`.
Ceci est illustré dans la commande ci-dessous, qui trouve les distances entre les trois premières caractéristiques de `nz_height` et les régions d'Otago et de Canterbury en Nouvelle-Zélande représentées par l'objet `co`.

```{r 04-spatial-operations-32}
co = filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)
```

Remarquez que le distance entre les deuxième et troisième éléments de `nz_height` et le deuxième élément de `co` est de zéro.
Cela démontre le fait que les distances entre les points et les polygones se réfèrent à la distance à *n'importe quelle partie du polygone* :
Les deuxième et troisième points dans `nz_height` sont *dans* Otago, ce qui peut être vérifié en les traçant (résultat non montré) :

```{r 04-spatial-operations-33, eval=FALSE}
plot(st_geometry(co)[2])
plot(st_geometry(nz_height)[2:3], add = TRUE)
```

## Opérations spatiales sur les données raster {#spatial-ras}

Cette section s'appuie sur la section \@ref(manipulating-raster-objects), qui met en évidence diverses méthodes de base pour manipuler des jeux de données raster. Elle illustre des opérations raster plus avancées et explicitement spatiales, et utilise les objets `elev` et `grain` créés manuellement dans la section \@ref(manipulating-raster-objects).
Pour la commodité du lecteur, ces jeux de données se trouvent également dans le paquet **spData**.

### Sélection spatiale {#spatial-raster-subsetting}

Le chapitre précédent (Section \@ref(manipulating-raster-objects)) a montré comment extraire les valeurs associées à des ID de cellules spécifiques ou à des combinaisons de lignes et de colonnes.
Les objets raster peuvent également être extraits par leur emplacement (coordonnées) et via d'autres objets spatiaux.
Pour utiliser les coordonnées pour la sélection, on peut "traduire" les coordonnées en un ID de cellule avec la fonction **terra** `cellFromXY()`.
Une alternative est d'utiliser `terra::extract()` (attention, il existe aussi une fonction appelée `extract()` dans le **tidyverse**\index{tidyverse (package)}) pour extraire des valeurs.
Les deux méthodes sont démontrées ci-dessous pour trouver la valeur de la cellule qui couvre un point situé aux coordonnées 0,1, 0,1.
\index{raster!subsetting}
\index{spatial!subsetting}

```{r 04-spatial-operations-34, eval = FALSE}
id = cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
# the same as
terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))
```

<!--jn:toDo-->
<!-- to update? -->
<!-- It is convenient that both functions also accept objects of class `Spatial* Objects`. -->
Les objets raster peuvent également être sélectionnés avec un autre objet raster, comme le montre le code ci-dessous :

```{r 04-spatial-operations-35, eval=FALSE}
clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip]
# on peut aussi utiliser extract
# terra::extract(elev, ext(clip))
```

Cela revient à récupérer les valeurs du premier objet raster (dans ce cas, `elev`) qui se trouvent dans l'étendue d'un second objet raster (ici : `clip`), comme illustré dans la figure \@ref(fig:raster-subset).

```{r raster-subset, echo = FALSE, fig.cap = "Raster originale (à gauche). Masque raster (au milieu). Résultat du clip (à droite).", fig.scap="Sélection de raster."}
knitr::include_graphics("figures/04_raster_subset.png")
```

L'exemple ci-dessus a retourné les valeurs de cellules spécifiques, mais dans de nombreux cas, ce sont des sorties spatiales qui sont nécessaires.
Cela peut être fait en utilisant l'opérateur `[`, avec `drop = FALSE`, comme indiqué dans la section \@ref(manipulating-raster-objects), qui montre également comment les objets raster peuvent être sélectionnés par divers objets.
Le code ci-dessous en est un exemple  en retournant les deux premières cellules (de la ligne supérieure) de `elev` en tant qu'objet raster  (seules les 2 premières lignes de la sortie sont montrées) :

```{r 04-spatial-operations-36, eval=FALSE}
elev[1:2, drop = FALSE]    # sélection spatiale par ID
#> class       : SpatRaster 
#> dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#> ...
```

```{r 04-spatial-operations-37, echo=FALSE, eval=FALSE}
# objectif: montrer le résultats de la sélection précédente
x = elev[1, 1:2, drop = FALSE]
plot(x)
```

Un autre cas d'utilisation courante de sélection spatiale est celle où une image raster avec des valeurs `logiques` (ou `NA`) est utilisée pour masquer une autre image raster avec la même étendue et la même résolution, comme illustré dans la Figure \@ref(fig:raster-subset).
Dans ce cas, les fonctions `[` et `mask()` peuvent être utilisées (résultats non montrés) :

```{r 04-spatial-operations-38, eval=FALSE}
# créer un masque raster
rmask = elev
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)
```

Dans le morceau de code ci-dessus, nous avons créé un objet masque appelé `rmask` avec des valeurs assignées aléatoirement à `NA` et `TRUE`.
Ensuite, nous voulons garder les valeurs de `elev` qui sont `TRUE` dans `rmask`.
En d'autres termes, nous voulons masquer `elev` avec `rmask`.

```{r 04-spatial-operations-38b, eval=FALSE}
# sélection spatiale
elev[rmask, drop = FALSE]           # avec l'opérateur [ 
mask(elev, rmask)                   # avec mask()
```

L'approche ci-dessus peut également être utilisée pour remplacer certaines valeurs (par exemple, celles qui seraient fausses) par NA.  

```{r 04-spatial-operations-38c, eval=FALSE}
elev[elev < 20] = NA
```

Ces opérations sont en fait des opérations booléennes locales puisque nous comparons deux rasters par cellule.
La sous-section suivante explore ces opérations et d'autres opérations connexes plus en détail.

### Algèbre raster

\index{map algebra}
Le terme "algèbre raster" a été inventé à la fin des années 1970 pour décrire un "ensemble de conventions, de capacités et de techniques" pour l'analyse des données géographiques raster *et* (bien que moins marquées) vectorielles [@tomlin_map_1994].
<!-- Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. -->
Dans ce contexte, nous définissons l'algèbre raster plus strictement, comme des opérations qui modifient ou résument les valeurs des cellules raster, en référence aux cellules environnantes, aux zones ou aux fonctions statistiques s'appliquant à chaque cellule.

Les opérations d'algèbre raster ont tendance à être rapides, car les jeux de données raster ne stockent qu'implicitement les coordonnées, d'où le [vieil adage](https://geozoneblog.wordpress.com/2013/04/19/raster-vs-vector/) "raster is faster but vector is corrector".
La position des cellules dans les jeux de données raster peut être calculée à l'aide de leur position matricielle, de la résolution et de l'origine du jeu de données (stockées dans l'en-tête).
Pour le traitement, cependant, la position géographique d'une cellule n'est guère pertinente tant que nous nous assurons que la position de la cellule est toujours la même après le traitement.
De plus, si deux ou plusieurs jeux de données raster partagent la même étendue, projection et résolution, on peut les traiter comme des matrices pour le traitement.

C'est de cette manière que l'algèbre raster fonctionne avec le paquet **terra**.
Premièrement, les en-têtes des jeux de données raster sont interrogés et (dans les cas où les opérations d'algèbre raster fonctionnent sur plus d'un ensemble de données) vérifiés pour s'assurer que les jeux de données sont compatibles.
Deuxièmement, l'algèbre raster conserve ce que l'on appelle la correspondance de localisation une à une, ce qui signifie que les cellules ne peuvent pas se déplacer.
Cela diffère de l'algèbre matricielle, dans laquelle les valeurs changent de position, par exemple lors de la multiplication ou de la division de matrices.

L'algèbre raster (ou modélisation cartographique avec des données raster) divise les opérations sur des rasters en quatre sous-classes [@tomlin_geographic_1990], chacune travaillant sur une ou plusieurs grilles simultanément :

1. Opérations *locales* ou par cellule
2. Les opérations *Focales* ou de voisinage.
Le plus souvent la valeur de la cellule de sortie est le résultat d'un bloc de cellules d'entrée de 3 x 3 cellules
3. Les opérations *Zonales* sont similaires aux opérations focales, mais la grille de pixels environnante sur laquelle les nouvelles valeurs sont calculées peut avoir des tailles et des formes irrégulières.
4. Les opérations *Globales* ou par-raster. 
Ici la cellule de sortie dérive potentiellement sa valeur d'un ou de plusieurs rasters entiers.

Cette typologie classe les opérations d'algèbre raster en fonction du nombre de cellules utilisées pour chaque étape de traitement des pixels et du type de sortie.
Par souci d'exhaustivité, nous devons mentionner que les opérations sur des rasters peuvent également être classées par discipline, comme le terrain, l'analyse hydrologique ou la classification des images.
Les sections suivantes expliquent comment chaque type d'opérations d'algèbre raster peut être utilisé, en se référant à des exemples documentés.

### Opérations locales

\index{map algebra!local operations}
Les opérations **locales** comprennent toutes les opérations cellule par cellule dans une ou plusieurs couches.
Elles ont un cas typique d'algèbre raster et comprennent l'ajout ou la soustraction de valeurs d'une image raster, l'élévation au carré et la multiplication d'images raster.
L'algèbre raster permet également des opérations logiques telles que la recherche de toutes les cellules qui sont supérieures à une valeur spécifique (5 dans notre exemple ci-dessous).
Le paquet **terra** prend en charge toutes ces opérations et bien plus encore, comme le montre la figure ci-dessous (\@ref(fig:04-local-operations)):

```{r 04-spatial-operations-41, eval = FALSE}
elev + elev
elev^2
log(elev)
elev > 5
```

```{r 04-local-operations, echo=FALSE, fig.cap="Exemples de différentes opérations locales de l'objet raster elev : additionner deux rasters, élever au carré, appliquer une transformation logarithmique et effectuer une opération logique."}
knitr::include_graphics("figures/04-local-operations.png")
```

Un autre bon exemple d'opérations locales est la classification d'intervalles de valeurs numériques en groupes, comme le regroupement d'un modèle numérique d'élévation en altitudes basses (classe 1), moyennes (classe 2) et hautes (classe 3).
En utilisant la commande `classify()`, nous devons d'abord construire une matrice de classification, où la première colonne correspond à l'extrémité inférieure et la deuxième colonne à l'extrémité supérieure de la classe.
La troisième colonne représente la nouvelle valeur pour les plages spécifiées dans les colonnes un et deux.

```{r 04-spatial-operations-40}
rcl = matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
rcl
```

Ici, nous affectons les valeurs du raster dans les plages 0--12, 12--24 et 24--36 et les *reclassons* pour prendre les valeurs 1, 2 et 3, respectivement.

```{r 04-spatial-operations-40b, eval = FALSE}
recl = classify(elev, rcl = rcl)
```

La fonction `classify()` peut également être utilisée lorsque nous voulons réduire le nombre de classes dans nos rasters catégorisés.
Nous effectuerons plusieurs reclassements supplémentaires dans le chapitre \@ref(location)'.

En dehors des opérateurs arithmétiques, on peut aussi utiliser les fonctions `app()`, `tapp()` et `lapp()`.
Elles sont plus efficaces, et donc préférables pour les grands jeux de données raster. 
De plus, elles vous permettent d'enregistrer directement un fichier de sortie.
La fonction `app()` applique une fonction à chaque cellule d'une couche matricielle et est utilisée pour résumer (par exemple, en calculant la somme) les valeurs de plusieurs couches en une seule couche.
`tapp()` est une extension de `app()`, nous permettant de sélectionner un sous-ensemble de couches (voir l'argument `index`) pour lesquelles nous voulons effectuer une certaine opération.
Enfin, la fonction `lapp()` permet d'appliquer une fonction à chaque cellule en utilisant les couches comme arguments -- une application de `lapp()` est présentée ci-dessous.

Le calcul du *normalized difference vegetation index* (NDVI) est une opération raster locale (pixel par pixel) bien connue.
Elle produit un raster dont les valeurs sont comprises entre -1 et 1 ; les valeurs positives indiquent la présence de plantes vivantes (le plus souvent > 0,2).
Le NDVI est calculé à partir des bandes rouge et proche infrarouge (NIR) d'images de télédétection, généralement issues de systèmes satellitaires tels que Landsat ou Sentinel.
La végétation absorbe fortement la lumière dans le spectre de la lumière visible, et surtout dans le canal rouge, tout en réfléchissant la lumière NIR, ce qui explique la formule du NDVI :

$$
\begin{split}
NDVI&= \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}\\
\end{split}
$$

Calculons le NDVI pour l'image satellite multispectrale du parc National de Zion.

```{r}
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)
```

L'objet raster comporte quatre bandes de satellite - bleu, vert, rouge et proche infrarouge (NIR).
Notre prochaine étape va être d'implémenter la formule NDVI dans une fonction R :

```{r}
ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}
```

Cette fonction accepte deux arguments numériques, `nir` et `red`, et retourne un vecteur numérique avec les valeurs NDVI.
Elle peut être utilisée comme argument `fun` de `lapp()`.
Nous devons juste nous rappeler que notre fonction n'a besoin que de deux bandes (et non de quatre comme dans le raster original), et qu'elles doivent être dans l'ordre NIR puis `red`.
C'est pourquoi nous sélection la grille d'entrée avec `multi_rast[[c(4, 3)]]` avant d'effectuer tout calcul.

```{r}
ndvi_rast = lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)
```

Le résultat, présenté sur le panneau de droite de la figure \@ref(fig:04-ndvi), peut être comparé à l'image RVB de la même zone (panneau de gauche de la même figure).
Il nous permet de constater que les plus grandes valeurs de NDVI sont liées aux zones de forêt dense dans les parties nord de la zone, tandis que les valeurs les plus faibles sont liées au lac au nord et aux crêtes montagneuses enneigées.

```{r 04-ndvi, echo=FALSE, fig.cap="Image RVB (à gauche) et valeurs NDVI (à droite) calculées pour l'exemple de l'image satellite du parc National de Zion."}
knitr::include_graphics("figures/04-ndvi.png")
```

La cartographie prédictive est une autre application intéressante des opérations raster locales.
La variable de réponse correspond à des points mesurés ou observés dans l'espace, par exemple, la richesse des espèces, la présence de glissements de terrain, les maladies des arbres ou le rendement des cultures.
Par conséquent, nous pouvons facilement récupérer des variables prédictives spatiales ou aériennes à partir de divers rasters (élévation, pH, précipitations, température, couverture végétale, classe de sol, etc.)
Ensuite, nous modélisons notre réponse en fonction de nos prédicteurs en utilisant `lm()`, `glm()`, `gam()` ou une technique d'apprentissage automatique. 
Les prédictions spatiales sur les objets raster peuvent donc être réalisées en appliquant des coefficients estimés aux valeurs de chaque cellule, et en additionnant les valeurs raster de sortie (voir chapitre \@ref(eco)).

### Operations focales

\index{map algebra!focal operations}
Alors que les fonctions locales opèrent sur une cellule, bien que pouvant provenir de plusieurs couches, les opérations **focales** prennent en compte une cellule centrale (focale) et ses voisines.
Le voisinage (également appelé noyau, filtre ou fenêtre mobile) utilisé est généralement de le taille 3 par 3 cellules (c'est-à-dire la cellule centrale et ses huit voisines), mais peut prendre toute autre forme (pas nécessairement rectangulaire) définie par l'utilisateur.
Une opération focale applique une fonction d'agrégation à toutes les cellules du voisinage spécifié, utilise la sortie correspondante comme nouvelle valeur pour la cellule centrale, puis passe à la cellule centrale suivante (figure \@ref(fig:focal-example)).
Cette opération est également appelée filtrage spatial et convolution [@burrough_principles_2015].

Dans R, nous pouvons utiliser la fonction `focal()` pour effectuer un filtrage spatial. 
Nous définissons la forme de la fenêtre mobile avec une `matrice` dont les valeurs correspondent aux poids (voir le paramètre `w` dans le code ci-dessous).
Ensuite, le paramètre `fun` nous permet de spécifier la fonction que nous souhaitons appliquer à ce voisinage.
Ici, nous choisissons le minimum, mais toute autre fonction de résumé, y compris `sum()`, `mean()`, ou `var()` peut être utilisée.

```{r 04-spatial-operations-42, eval = FALSE}
r_focal = focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)
```

Cette fonction accepte également des arguments supplémentaires, par exemple, si elle doit supprimer les NAs dans le processus (`na.rm = TRUE`) ou non (`na.rm = FALSE`).

```{r focal-example, echo = FALSE, fig.cap = "Raster d'entrée (gauche) et raster de sortie (droite) suite à une opération focale - trouver la valeur minimale dans des fenêtres mobiles 3 par 3.", fig.scap="Exemple d'une opération focale."}
knitr::include_graphics("figures/04_focal_example.png")
```

Nous pouvons rapidement vérifier si le résultat correspond à nos attentes.
En effet, dans notre exemple, la valeur minimale doit toujours se situer dans le coin supérieur gauche de la fenêtre mobile (rappelez-vous que nous avons créé le raster d'entrée en incrémentant les valeurs des cellules d'une unité par ligne, en commençant par le coin supérieur gauche).
Dans cet exemple, la matrice de pondération est composée uniquement de 1, ce qui signifie que chaque cellule a le même poids sur la sortie, mais cela peut être modifié.

Les fonctions ou filtres focaux jouent un rôle dominant dans le traitement des images.
Les filtres "passe-bas" (*low-pass*) ou de lissage utilisent la fonction moyenne pour éliminer les valeurs extrêmes.
Dans le cas de données catégorielles, on peut remplacer la moyenne par le mode, qui est la valeur la plus courante.
En revanche, les filtres "passe-haut" accentuent les entités.
Les filtres Laplace et Sobel de détection de lignes peuvent servir d'exemple ici.
Consultez la page d'aide `focal()` pour savoir comment les utiliser dans R (ils seront également utilisés dans les exercices à la fin de ce chapitre).

Le traitement des données numériques de terrain, le calcul des caractéristiques topographiques telles que la pente, l'aspect et les directions d'écoulement, repose sur des fonctions focales.
La fonction de **terra** `terrain()` peut être utilisé pour calculer ces métriques, bien que certains algorithmes de traitement des données numériques de terrain, y compris la méthode de Zevenbergen et Thorne pour calculer la pente, n'y soient pas implémentés.
De nombreux autres algorithmes --- notamment les courbures, les zones contributives et les indices d'humidité --- sont mis en œuvre dans des logiciels libres de systèmes d'information géographique (SIG) de bureau.
Le chapitre \@ref(gis) montre comment accéder à ces fonctionnalités SIG à partir de R.

### Opérations Zonales

\index{map algebra!zonal operations}
Tout comme les opérations focales, les opérations *zonales* appliquent une fonction d'agrégation à plusieurs cellules.
Cependant, un deuxième raster, généralement avec des valeurs catégorielles, définit les *filtres zonaux* (ou 'zones') dans le cas des opérations zonales, en opposition à une fenêtre de voisinage prédéfinie dans le cas de l'opération focale.
Par conséquent, les cellules  définissant le filtre zonal ne doivent pas nécessairement être voisines.
Notre raster `grain` en est un bon exemple, comme l'illustre le panneau de droite de la figure \@ref(fig:cont-raster) : différentes tailles de granulométrie sont réparties de manière irrégulière dans le raster.
Enfin, le résultat d'une opération zonale est un tableau récapitulatif groupé par zone, c'est pourquoi cette opération est également connue sous le nom de *statistiques zonales* dans le monde des SIG\index{GIS}. 
Ceci est à l'opposé des opérations focales qui retournent un objet matriciel.

Le  code suivant utilise la fonction `zonal()` pour calculer l'altitude moyenne associée à chaque classe de taille de grain.

```{r 04-spatial-operations-43}
z = zonal(elev, grain, fun = "mean")
z
```

ceci renvoie les statistiques\index{statistics} pour chaque catégorie, ici l'altitude moyenne pour chaque classe de granulométrie.
Remarque : il est aussi possible d'obtenir un raster avec les statistiques calculées pour chaque zone en mettant l'argument `as.raster` à `TRUE`.

### Opérations globales et distances

Les opérations *globales* sont un cas particulier des opérations zonales, l'ensemble des données raster représentant une seule zone.
Les opérations globales les plus courantes sont des statistiques descriptives pour l'ensemble des données raster, telles que le minimum ou le maximum, que nous avons déjà abordées dans la section \@ref(summarizing-raster-objects).

En dehors de cela, les opérations globales sont également utiles pour le calcul des rasters de distance et de poids.
Dans le premier cas, on peut calculer la distance entre chaque cellule et une cellule cible spécifique.
Par exemple, on peut vouloir calculer la distance à la côte la plus proche (voir aussi `terra::distance()`).
On peut aussi vouloir prendre en compte la topographie, c'est-à-dire qu'on n'est pas seulement intéressé par la distance pure mais qu'on voudrait aussi éviter de traverser des chaînes de montagnes en allant vers la côte.
Pour ce faire, nous pouvons pondérer la distance par l'altitude de sorte que chaque mètre d'altitude supplémentaire "prolonge" la distance euclidienne.
Les calculs de visibilité et de bassin visuel appartiennent également à la famille des opérations globales (dans les exercices du chapitre \@ref(gis), vous calculerez un raster de bassin visuel).

### Les contreparties de l'algèbre raster dans le traitement vectoriel

De nombreuses opérations d'algèbre raster ont une contrepartie dans le traitement vectoriel [@liu_essential_2009].
Le calcul d'une distance raster (opération globale) en ne considérant qu'une distance maximale (opération focale logique) est l'équivalent d'une opération de mise en mémoire tampon vectorielle (section \@ref(clipping)).
Le reclassement de données raster (fonction locale ou zonale selon l'entrée) est équivalent à la dissolution de données vectorielles (section \@ref(spatial-joining)). 
La superposition de deux données matricielles (opération locale), dont l'une contient des valeurs `NULL` ou `NA` représentant un masque, est similaire au découpage vectoriel (Section `@ref(clipping)).
L'intersection de deux couches est tout à fait similaire au détourage spatial (section \@ref(spatial-subsetting)). 
La différence est que ces deux couches (vectorielles ou matricielles) partagent simplement une zone de chevauchement (voir la figure \@ref(fig:venn-clip) pour un exemple).
Cependant, faites attention à la formulation.
Parfois, les mêmes mots ont des significations légèrement différentes pour les modèles de données rasters et vectorielles.
Dans le cas des données vectorielles, l'agrégation consiste à dissoudre les polygones, tandis que dans le cas des données rasters, elle consiste à augmenter la résolution.
En fait, on pourrait considérer que dissoudre ou agréger des polygones revient à diminuer la résolution. 
Cependant, les opérations zonales semble être le meilleur équivalent raster par rapport à la modification de la résolution des cellules. 
Les opérations zonales peuvent dissoudre les cellules d'un raster en fonction des zones (catégories) d'un autre raster en utilisant une fonction d'agrégation (voir ci-dessus).

### Fusionner des rasters

\index{raster!merge}
Supposons que nous voulions calculer le NDVI (voir la section \@ref(local-operations)), et que nous voulions en plus calculer les attributs du terrain à partir des données d'altitude pour des observations dans une zone d'étude.
Ces calculs reposent sur des informations de télédétection. 
L'imagerie correspondante est souvent divisée en tuiles couvrant une étendue spatiale spécifique, et fréquemment, une zone d'étude couvre plus d'une tuile.
Dans ce cas, nous devons fusionner les tuiles couvertes par notre zone d'étude. 
Le plus simple est alors de fusionner ces scènes, c'est-à-dire les mettre côte à côte.
Cela est possible, par exemple, avec les données numériques d'altitude (SRTM, ASTER).
Dans l'extrait de code suivant, nous téléchargeons d'abord les données d'élévation SRTM pour l'Autriche et la Suisse (pour les codes pays, voir la fonction **geodata** `country_codes()`).
Dans une seconde étape, nous fusionnons les deux rasters en un seul.

```{r 04-spatial-operations-44, eval = FALSE}
aut = geodata::elevation_30s(country = "AUT", path = tempdir())
ch = geodata::elevation_30s(country = "CHE", path = tempdir())
aut_ch = merge(aut, ch)
```

La commande `merge()` de **terra** combine deux images, et dans le cas où elles se chevauchent, elle utilise la valeur du premier raster.
<!--jn:toDo-->
<!-- gdalUtils is slower (for this files): -->
<!-- two_rast = c(terra::sources(aut)$source, terra::sources(ch)$source) -->
<!-- tf = tempfile(fileext = ".tif") -->
<!-- bench::mark({gdalUtils::mosaic_rasters(two_rast, tf)}) -->
<!-- You can do exactly the same with `gdalUtils::mosaic_rasters()` which is faster, and therefore recommended if you have to merge a multitude of large rasters stored on disk. -->

cette approche de fusion est peu utile lorsque les valeurs qui se chevauchent ne correspondent pas les unes aux autres.
C'est souvent le cas lorsque vous voulez combiner des images spectrales provenant de scènes qui ont été prises à des dates différentes.
La commande `merge()` fonctionnera toujours mais vous verrez une frontière nette dans l'image résultante.
D'autre part, la commande `mosaic()` vous permet de définir une fonction pour la zone de recouvrement. 
Par exemple, nous pourrions calculer la valeur moyenne -- cela pourrait lisser la bordure claire dans le résultat fusionné, mais ne la fera probablement pas disparaître.
<!-- The following sentences have been commented out and can be removed because the packages, and info, is now out of date -->
<!-- See https://github.com/Robinlovelace/geocompr/pull/424 for discussion -->
<!-- To do so, we need a more advanced approach.  -->
<!-- Remote sensing scientists frequently apply histogram matching or use regression techniques to align the values of the first image with those of the second image. -->
<!-- The packages **landsat** (`histmatch()`, `relnorm()`, `PIF()`), **satellite** (`calcHistMatch()`) and **RStoolbox** (`histMatch()`, `pifMatch()`) provide the corresponding functions for the **raster**'s package objects. -->
Pour une introduction plus détaillée à la télédétection avec R, il est possible de consulter @wegmann_remote_2016.
<!--jn:toDo-->
<!--update the above reference to the 2nd edition-->

## Exercises


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