Le package raster

class: center, middle, inverse, title-slide

# Le package raster - TP4
### CEFE-CNRS
### 20 mars 2019

---

# TP4. Rastérisation de données vecteur.

Chargement des packages pour le TP 4

```r
library(sp)
library(rgdal)
library(raster)
library(gdalUtils)
#library(sf)
#library(fasterize) 
```

---

## 4.2. Rastériser un shapefile

Rastériser = transformer les données vectorielles en raster.

Il existe 3 moyens de rastériser des données issues d'un shapefile avec R :

* avec le package `raster`, fonction rasterize. Très lent en cas de données volumineuses.
* avec le package `gdalUtils`, fonction gdal_rasterize. Plus efficace.
* avec le package `fasterize` : très perfomant, mais uniquement pour les données vectorielles de type polygone avec `sf`.

---

## 4.3. Exemple avec Corine Land Cover
Dans les exemples suivants, nous transformerons en raster le fichier CLC12_D080_RGF.shp.

Les valeurs à rastériser (Code CLC niveau 3) se trouvent dans le champ CODE_12 (type String). Dans le raster en sortie, les codes CLC devront être enregistrés en tant qu'entier (INT2U).

Le raster obtenu, pour pouvoir être croisé avec les données satellitaires Sentinel, devra partager avec le raster Sentinel2 les caratéristiques suivantes :

* Système de coordonnées = UTM31N WGS84 (EPSG:32631)
* ULX, ULY = 395000,5575000
* Résolution = 10,10
* Format = GeoTIFF

Dans la foulée nous générerons également 2 autres raster avec les niveaux 1 et 2 de nomenclature.

---

## 4.4. Conversion en UTM31N
Le fichier de départ CLC12_D080_RGF.shp est en Lambert 93 (EPSG:2154). La première tâche est de le convertir en UTM31N WGS84 (EPSG:32631).

```r
shp_clc80_l93 <- readOGR(dsn="data2/CLC", layer="CLC12_D080_RGF")
```

```
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GitRepo\RSpatial\formation_2018\data2\CLC", layer: "CLC12_D080_RGF"
## with 3267 features
## It has 3 fields
```

```r
shp_clc80_utm31n <- spTransform(shp_clc80_l93, CRS("+init=EPSG:32631"))
```

---

## 4.5. Préparer le champ dont les valeurs seront enregistrées dans le raster
Il faut également ajouter 3 champs NIV3_12, NIV2_12 et NIV1_12 qui contiendront les codes en entier pour les 3 niveaux de détails de la nomenclature.

```r
code_12 <- shp_clc80_utm31n$CODE_12 # attention factor !
shp_clc80_utm31n$NIV3_12 <- as.numeric(levels(code_12))[code_12]
shp_clc80_utm31n$NIV2_12 <- as.numeric(substr(levels(code_12),1,2))[code_12]
shp_clc80_utm31n$NIV1_12 <- as.numeric(substr(levels(code_12),1,1))[code_12]
head(shp_clc80_utm31n@data)
```

```
##          ID CODE_12      AREA_HA NIV3_12 NIV2_12 NIV1_12
## 0 FR-271955     423   130.829138     423      42       4
## 1  FR-62263     211 24285.214749     211      21       2
## 2  FR-24807     112    76.021063     112      11       1
## 3  FR-24804     112   107.040152     112      11       1
## 4 FR-104366     231  1854.517282     231      23       2
## 5 FR-172279     243     4.530164     243      24       2
```

```r
writeOGR(shp_clc80_utm31n, dsn="data2/CLC", layer="CLC12_D080_UTM31N", 
         driver="ESRI Shapefile", overwrite_layer=T)
```

---

## 4.6. Rasterisation avec le package raster
Rastérisons les codes de niveau 3.

```r
f_img_sat <- "data2/S2_Marquenterre_20180504.tif "
rast_template <- raster(f_img_sat, band=1)

start_time <- Sys.time()
rast_clc80_niv3 <- rasterize(shp_clc80_utm31n, rast_template, "NIV3_12", 
                             filename="data2/CLC/CLC12_D080_NIV3.tif", datatype="INT2S", 
                             format="GTiff", overwrite=T, NAflag=-99)
Sys.time() - start_time
```

```
## Time difference of 20.45717 secs
```

---

## 4.7. Rasterisation avec le package gdalUtils
Rastérisons les codes de niveau 2.

```r
f_in <- "data2/CLC/CLC12_D080_UTM31N.shp"
f_out <- "data2/CLC/CLC12_D080_NIV2.tif"
ext_in <- extent(rast_template)
te_in <- ext_in[c(1,3,2,4)]
start_time <- Sys.time()
rast_niv2 <- gdal_rasterize(src_datasource=f_in,
                            dst_filename=f_out,
                            a="NIV2_12",
                            of="GTiff",
                            a_srs="EPSG:32631",
                            a_nodata=-99,
                            te=te_in,
                            tr=c(10,10),
                            ot="Int16"
                            )
Sys.time() - start_time
```

```
## Time difference of 15.7909 secs
```

---

## 4.8. Rasterisation avec sf et fasterize
Rastérisons les codes de niveau 3.

Utilisation du package fasterize :
* uniquement pour les shapefile de type "polygone"" !
* en conjonction avec sf

Voir https://github.com/ecohealthalliance/fasterize

---

## 4.8. Rasterisation avec sf et fasterize

```r
library(sf)
```

```
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
```

```
## 
## Attaching package: 'sf'
```

```
## The following object is masked from 'package:gdalUtils':
## 
##     gdal_rasterize
```

```r
library(fasterize) 
```

```
## 
## Attaching package: 'fasterize'
```

```
## The following object is masked from 'package:graphics':
## 
##     plot
```

```r
f_in <- "data2/CLC/CLC12_D080_UTM31N.shp"
f_out <- "data2/CLC/CLC12_D080_NIV1.tif"
shp_in <- st_read(f_in)
```

```
## Reading layer `CLC12_D080_UTM31N' from data source `D:\GitRepo\RSpatial\formation_2018\data2\CLC\CLC12_D080_UTM31N.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3267 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 379543.1 ymin: 5486582 xmax: 519049.6 ymax: 5584755
## epsg (SRID):    32631
## proj4string:    +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs
```

```r
start_time <- Sys.time()
rast_clc80_niv1 <- fasterize(shp_in, rast_template, field = "NIV1_12")
writeRaster(rast_clc80_niv1, filename=f_out, 
            datatype="INT2S", format="GTiff", overwrite=T, NAflag=-99)
Sys.time() - start_time
```

```
## Time difference of 0.4870279 secs
```