Analyse d’images raster

class: center, middle, inverse, title-slide

# Analyse d’images raster
## (et télédétection)
### Malika Madelin <a href="mailto:malika.madelin@u-paris.fr" class="email">malika.madelin@u-paris.fr</a>
### 1/7/2021

---

class:  center middle
background-image: url("fig/Couv.png")
background-size: cover

# Les images raster

---

![](fig/Landsat.png)
[Landsat NASA](https://landsat.gsfc.nasa.gov/sites/landsat/)

---

.pull-left[
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-1-1.png" width="100%" />
]

.pull-right[
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-2-1.png" width="100%" />
]

--
class:center

**qq mots sur le poids des images...**

--
<center><img src="fig/datatype.png" height="200px" /></center>
[dataType {raster}](https://search.r-project.org/CRAN/refmans/raster/html/dataType.html)

---
## Intérêts/avantages
.pull-left[
- structure des données ~ simple, topologie implicite
]

---
## Intérêts/avantages
.pull-left[
- structure des données ~ simple, topologie implicite
- information spatiale **continue**, 
souvent récoltée **à distance** et **globalement**
]
<center><img src="fig/Deforestation.png"></center>

---
## Intérêts/avantages
![](fig/Diapo_CollHokkaido.png)

---
## Intérêts/avantages

![](fig/Diapo_CollHokkaido2.png)

---
## Intérêts/avantages

.pull-left[
- structure des données ~ simple, topologie implicite
- information spatiale **continue**, 
souvent récoltée **à distance** et **globalement**
- information figée dans le **temps** + répétitivité
- caractère **multispectral** (bandes/couches),
- parfois **multidimensionnel** (hauteur, profondeur / temps)
]
.pull-right[
<center>Composition colorée, Sentinel 2A, 28/5/21</center>
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-3-1.png" width="100%" />
]

---
## Intérêts/avantages

.pull-left[
- structure des données ~ simple, topologie implicite
- information spatiale **continue**, 
souvent récoltée **à distance** et **globalement**
- information figée dans le **temps** + répétitivité
- caractère **multispectral** (bandes/couches),
- parfois **multidimensionnel** (hauteur, profondeur / temps)
- plusieurs thématiques possibles pour une même image
- relative facilité des traitements
]
.pull-right[
![](fig/Diapo_CollICUC.png)
Extrait d'une présentation de Madelin et Dupuis, 2018 
]

---
## Intérêts/avantages

.pull-left[
- structure des données ~ simple, topologie implicite
- information spatiale **continue**, 
souvent récoltée **à distance** et **globalement**
- information figée dans le **temps** + répétitivité
- caractère **multispectral** (bandes/couches),
- parfois **multidimensionnel** (hauteur, profondeur / temps)
- plusieurs thématiques possibles pour une même image
- relative facilité des traitements

]
.pull-right[
![](fig/Diapo_CollICUC.png)
Extrait d'une présentation de Madelin et Dupuis, 2018 
]

... mais **volume** des données non négligeable

---

## Par rapport à la reproductibilité

- souvent homogénéité dans la production des données (produits globaux), harmonisation "minime"

--
<center><img src="fig/Diapo_CollCIST.png" height="370px" /></center>
<center>Extrait d'une présentation de Pécout <it>et al</it>, 2016</center>
<center>Projet Grandes métropoles</center>
---

## Par rapport à la reproductibilité

- souvent homogénéité dans la production des données (produits globaux), harmonisation "minime"

- données souvent accessibles 
_via_ des **plateformes** ([Earth Explorer](), [Sentinel Hub](), [GEE](), ...)

---

## Par rapport à la reproductibilité

- souvent homogénéité dans la production des données (produits globaux), harmonisation "minime"

- données souvent accessibles 
_via_ des **plateformes** ([Earth Explorer](), [Sentinel Hub](), [GEE](), ...) 
_via_ des **outils d'extraction** pour certaints formats (ncdf) 
_via_ certains packages R (`{raster}`, `{modisTSP}`, `{sen2r}`...)

--
- de nombreux formats ouverts 
[geotiff](https://www.ogc.org/standards/geotiff) / [jp2](https://fr.wikipedia.org/wiki/JP2) / asc / [gpkg-raster](https://gdal.org/drivers/raster/gpkg.html) / [coverageJSON](https://covjson.org/...) 
([ncdf](https://www.unidata.ucar.edu/software/netcdf/) / hdf /grib...)

---

![](fig/Isere_ascii.png)
(BD Alti 25m)

---

## Par rapport à la reproductibilité

- souvent homogénéité dans la production des données (produits globaux), harmonisation "minime"

- de nombreux formats ouverts 
[geotiff](https://www.ogc.org/standards/geotiff) / [jp2](https://fr.wikipedia.org/wiki/JP2) / asc / [gpkg-raster](https://gdal.org/drivers/raster/gpkg.html) / [coverageJSON](https://covjson.org/...) 
([ncdf](https://www.unidata.ucar.edu/software/netcdf/) / hdf /grib...)

--
<center><img src="fig/ncdf.png" height="150px"></center>
[UCAR](http://www.narccap.ucar.edu)

---
## Par rapport à la reproductibilité

.pull-left[
<center>Modèle numérique d'élévation SRTM - Argentine</center>
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-4-1.png" width="100%" />

]

.pull-right[

```r
dim(srtm)
```

```
## [1] 6000 6000    1
```
Export `writeRaster()` selon 3 formats 
et calcul de la taille des fichiers résultants `file.size()`

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> Type du fichier </th>
 <th style="text-align:right;"> Taille en Mo </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> ascii </td>
 <td style="text-align:right;"> 155.7 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> tif </td>
 <td style="text-align:right;"> 60.6 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> gpkg </td>
 <td style="text-align:right;"> 48.7 </td>
 </tr>
</tbody>
</table>
]

---
## Par rapport à la reproductibilité

- souvent homogénéité dans la production des données (produits globaux), harmonisation "minime"

- des libraires C, communes à plusieurs logiciels/langages (GDAL, ...)

---

## Sous R, 3 packages "phares"

.left-column[
[`{raster}`](https://cran.r-project.org/web/packages/raster/)

![](fig/Hijmans.png)
]

.right-column[
 .pull-left[
 <center><img src="fig/raster.png" /></center>
 ]
 .pull-right[
 Documentation de 249p. (18/6/2018), 
 avec en préambule, une organisation par blocs : 
 ![](fig/raster2.png) 
 ]

]

---

## Sous R, 3 packages "phares"

.left-column[
[`{stars}`](https://cran.r-project.org/web/packages/stars/index.html)

![](fig/EdzerPebesma.png)
]
.right-column[

![](fig/stars2.png)
  ]

---

## Sous R, 3 packages "phares"

.left-column[
[`{terra}`](https://cran.r-project.org/web/packages/terra/)

![](fig/Hijmans.png)
]
.right-column[
 .pull-left[
 <center><img src="fig/terra.png" /></center>
 ]
 .pull-right[
 Documentation de 218p. (20/6/2018), 
 avec en préambule, une organisation par blocs : 
 ![](fig/terra2.png)
 ]

]

---
## Sous R, 3 packages "phares"

La "popularité" des 3 packages, à partir du nombre de téléchargements

---
## Sous R, 3 packages "phares"

Comparatif du [blog de seascape](https://www.seascapemodels.org/rstats/2021/06/01/STARS.html) : "So terra was about 8.5 times faster than raster and stars was about 3 times faster than terra (so about 26 times faster than raster)."

Pour les opérations réalisées...

![](fig/Benchmark.png)
[Krzysztof Dyba - kadyb](https://github.com/kadyb/raster-benchmark)

---
class: inverse center middle

# Les traitements sur les images raster

---

# Lecture des fichiers raster sous R

.left-column[
![](fig/Img_zoom_sen2b.png)
]
.right-column[

```r
file.size("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif") * 10^-3 #  (en Ko)
```

```
## [1] 183.038
```

```r
GDALinfo("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif")
```

```
## rows        218 
## columns     279 
## bands       3 
## lower left origin.x        632850.6 
## lower left origin.y        5093260 
## res.x       9.984642 
## res.y       9.984154 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## file        data/Sentinel/Sen2_zoom_TCI_280521_10m.tif 
## apparent band summary:
##   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1   Byte          FALSE           0          9        279
## 2   Byte          FALSE           0          9        279
## 3   Byte          FALSE           0          9        279
## apparent band statistics:
##   Bmin Bmax Bmean Bsd
## 1    0  255    NA  NA
## 2    0  255    NA  NA
## 3    0  255    NA  NA
## Metadata:
## AREA_OR_POINT=Area
```
]
---

# Lecture des fichiers raster sous R

Avec [`{raster}`](https://cran.r-project.org/web/packages/raster/)

```r
r <- raster("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif")
r
```

```
## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 218, 279, 60822  (nrow, ncol, ncell)
## resolution : 9.984642, 9.984154  (x, y)
## extent     : 632850.6, 635636.3, 5093260, 5095437  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : Sen2_zoom_TCI_280521_10m.tif 
## names      : Sen2_zoom_TCI_280521_10m 
## values     : 0, 255  (min, max)
```

```r
object.size(r) ; inMemory(r)
```

```
## 14608 bytes
```

```
## [1] FALSE
```

---

# Lecture des fichiers raster sous R

Avec [`{raster}`](https://cran.r-project.org/web/packages/raster/)

```r
r_brick <- brick("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif")
r_brick
```

```
## class      : RasterBrick 
## dimensions : 218, 279, 60822, 3  (nrow, ncol, ncell, nlayers)
## resolution : 9.984642, 9.984154  (x, y)
## extent     : 632850.6, 635636.3, 5093260, 5095437  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : Sen2_zoom_TCI_280521_10m.tif 
## names      : Sen2_zoom_TCI_280521_10m.1, Sen2_zoom_TCI_280521_10m.2, Sen2_zoom_TCI_280521_10m.3 
## min values :                          0,                          0,                          0 
## max values :                        255,                        255,                        255
```

```r
object.size(r_brick)
```

```
## 15192 bytes
```

---

# Lecture des fichiers raster sous R

Avec [`{terra}`](https://cran.r-project.org/web/packages/terra/)

```r
t <- rast("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif")
t
```

```
## class       : SpatRaster 
## dimensions  : 218, 279, 3  (nrow, ncol, nlyr)
## resolution  : 9.984642, 9.984154  (x, y)
## extent      : 632850.6, 635636.3, 5093260, 5095437  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source      : Sen2_zoom_TCI_280521_10m.tif 
## red-grn-blue: 1, 2, 3 
## names       : Sen_1, Sen_2, Sen_3
```

```r
object.size(t)
```

```
## 1304 bytes
```

---

# Lecture des fichiers raster sous R

Avec [`{stars}`](https://cran.r-project.org/web/packages/stars/index.html)

```r
s <- read_stars("data/Sentinel/Sen2_zoom_TCI_280521_10m.tif")
s
```

```
## stars object with 3 dimensions and 1 attribute
## attribute(s):
##                               Min. 1st Qu. Median     Mean 3rd Qu. Max.
## Sen2_zoom_TCI_280521_10m.tif    10      30     51 57.15503      64  255
## dimension(s):
##      from  to  offset    delta                refsys point values x/y
## x       1 279  632851  9.98464 WGS 84 / UTM zone 30N FALSE   NULL [x]
## y       1 218 5095437 -9.98415 WGS 84 / UTM zone 30N FALSE   NULL [y]
## band    1   3      NA       NA                    NA    NA   NULL
```

```r
object.size(s)
```

```
## 1471360 bytes
```

---
# Lecture des fichiers raster sous R

** Opérations classiques après lecture / création d'un raster : **

.left-column[
- dimensions
- étendue
- résolution(s)
- projection
]

---
# Lecture des fichiers raster sous R

** Opérations classiques après lecture / création d'un raster :**

.left-column[
- dimensions
- étendue
- résolution(s)
- projection
- visualisation (rapide)
- histogramme des valeurs
- parfois améliorations des contrastes ( _stretching_ )
]

.right-column[
.pull-left[
![](fig/Img_zoom_sen2b.png)
]

.pull-right[

<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-14-1.png" width="100%" />
À partir de Sentinel B08
]

]

---
# Modifications de la zone d'étude

--
### • (re)projection

![](fig/FigCanada.png)
Source : Centre Canadien de Télédétection

---
# Modifications de la zone d'étude

### • (re)projection
`raster:projectRaster()` / `terra:project()` / `stars:st_transform()`

--
<center><img src="fig/stars_projection.png" height="300px"></center>
[avec `{stars}`](https://cran.r-project.org/web/packages/stars/vignettes/stars5.html)
---
# Modifications de la zone d'étude

### • découpage et masque ( _crop & mask_ )

```r
# on reprend la compo colorée de Sentinel 2 (sen)
# ajout du buffer 30 km autour d'ici (cf. présentation de Ronan)
buf <- st_read(dsn = paste0("data/KeynoteRonan/oleron.gpkg"), 
 layer = "cnrsBuf", quiet = TRUE)
buf <- st_transform(buf, crs(sen)) # pour avoir le même CRS
buf <- vect(buf) # transformation de l'objet sf en object SpatVector (terra)

# découpage (selon l'étendue)
sen_Oleron_crop <- crop(sen, buf)

# masque
sen_Oleron_mask <- mask(sen, buf)
```

---
# Modifications de la zone d'étude

### • découpage et masque ( _crop & mask_ )

.left-column[
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-16-1.png" width="100%" />
]

---
# Modifications de la zone d'étude

### • découpage et masque ( _crop & mask_ )

.left-column[
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-17-1.png" width="100%" />
]

.right-column[

.pull-left[
<center> _CROP_ <center> 
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-18-1.png" width="100%" />
]

.pull-right[
<center> _MASK_ </center> 
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-19-1.png" width="100%" />
]

]

---
# Modifications de la zone d'étude

### • agrégation / désagrégation

![](fig/agg_disag.png)

---
# Modifications de la zone d'étude

### • Fusion _merge_ et mosaïque

![](fig/2SRTM.png)

---
# Opérations locales

![](fig/OpLocales.png)

Source : J. Mennis

---
# Opérations locales

## • Reclassification ( _reclassify_ )

Des 44 postes les plus fins de la nomenclature CLC2018 aux 5 principaux 
![](fig/Reclass.png)

---
# Opérations locales

## • Création de "masque raster" (1/0)

De la lecture des valeurs spectrales à une typologie eau / pas eau
<img src="SIGR2021_raster_MM_files/figure-html/unnamed-chunk-20-1.png" width="100%" />

Autre possibilité : `r[r>500] <- NA`

---
# Opérations locales