tutorial.Rmd 8.31 KB
Newer Older
Florian de Boissieu's avatar
Florian de Boissieu committed
1
---
Florian de Boissieu's avatar
Florian de Boissieu committed
2
title: "Tutorial for biodivMapR"
Florian de Boissieu's avatar
Florian de Boissieu committed
3
4
5
6
7
8
9
10
11
12
13
author: "Jean-Baptiste Féret, Florian de Boissieu"
date: "`r Sys.Date()`"
output: 
  html_vignette:
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

14
15
16
17
18
19
20
21
22
```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval=FALSE
)
```


Florian de Boissieu's avatar
Florian de Boissieu committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
This tutorial aims at describing the processing workflow and giving the associated code to compute $\alpha$ and $\beta$ diversity maps on an extraction of Sentinel-2 image taken over Cameroun forest. The workflow is composed of the following steps:

* define the processing parameters

    * input / output files paths
    * output spatial resolution
    * computing options
    
* compute the $\alpha$ and $\beta$ diversity maps:
    * compute PCA and select best components
    * compute the diversity maps
    
* validate results comparing to field plots measurements



# Processing parameters

## Input / Output files
jbferet's avatar
jbferet committed
42
43
The input images are expected to be in ENVI HDR format, BIL interleaved. To check if the flormat is good use fucntion `check_data`.
If not they should be converted with function `raster2BIL`.
Florian de Boissieu's avatar
Florian de Boissieu committed
44
45
46

A mask can also be set to work on a selected part of the input image. The mask is expected to be a raster in the same format as the image (ENVI HDR), with values 0 = masked or 1 = selected. If no mask is to be used set `Input.Mask.File = FALSE`.

47
The output directory defined with `Output.Dir` will contain all the results. For each image processed, a subdirectory will be automatically created after its name.
Florian de Boissieu's avatar
Florian de Boissieu committed
48

49
50

```{r Input / Output files}
Florian de Boissieu's avatar
Florian de Boissieu committed
51
Input.Image.File  = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'biodivMapR')
jbferet's avatar
jbferet committed
52
check_data(Input.Image.File)
Florian de Boissieu's avatar
Florian de Boissieu committed
53

jbferet's avatar
jbferet committed
54
Input.Image.File  = raster2BIL(Raster.Path = Input.Image.File,
Florian de Boissieu's avatar
Florian de Boissieu committed
55
56
57
58
59
60
61
62
                                       Sensor = 'SENTINEL_2A',
                                       Convert.Integer = TRUE,
                                       Output.Directory = '~/test')
Input.Mask.File   = FALSE

Output.Dir        = 'RESULTS'
```

63

Florian de Boissieu's avatar
Florian de Boissieu committed
64
## Spatial resolution
jbferet's avatar
jbferet committed
65
The algorithm estimates \alpha and \beta diversity within a window, that is also the output spatial resolution. It is defined in number of pixel s of the input image with parameter `window_size`, e.g. `window_size = 10` meaning a window of 10x10 pixels. It will be the spatial resolution of the ouput rasters.
Florian de Boissieu's avatar
Florian de Boissieu committed
66
67

As a rule of thumb, spatial units between 0.25 and 4 ha usually match with ground data.
jbferet's avatar
jbferet committed
68
A window_size too small results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image.
Florian de Boissieu's avatar
Florian de Boissieu committed
69

jbferet's avatar
jbferet committed
70
In this example, the spatial resolution of the input raster. Setting `window_size = 10` will result in diversity maps of spatial resolution 100x100m.
Florian de Boissieu's avatar
Florian de Boissieu committed
71

72
```{r Spatial resolution}
jbferet's avatar
jbferet committed
73
window_size = 10
Florian de Boissieu's avatar
Florian de Boissieu committed
74
75
76
77
78
```

## PCA filtering
If set to `TRUE`, a second filtering based on PCA outliers is processed.

79
```{r PCA filtering}
Florian de Boissieu's avatar
Florian de Boissieu committed
80
81
82
83
84
85
FilterPCA = TRUE
```

## Computing options
The use of computing ressources can be controled with the following parameters: 

jbferet's avatar
jbferet committed
86
87
88
* `nbCPU` controls the parallelisation of the processing: how many CPUs will be asssigned for multithreading, 
* `MaxRAM` controls the size in GB of the input image chunks processed by each thread (this does not correspond to the max amount of RAM allocated),
* `nbclusters` controls the number of clusters defined by k-means clustering for each repetition. Images showing moderate diversity (temperate forests) may need lower number of clusters (20), whereas highly diverse tropical forest require more (50). The larger the value the longer the computation time. 
Florian de Boissieu's avatar
Florian de Boissieu committed
89

90
```{r Computing options}
Florian de Boissieu's avatar
Florian de Boissieu committed
91
92
93
94
95
96
97
98
nbCPU         = 4
MaxRAM        = 0.5
nbclusters    = 50
```

# Main processing worflow
## Mask non vegetated / shaded / cloudy pixels

99
```{r Mask non vegetated / shaded / cloudy pixels}
Florian de Boissieu's avatar
Florian de Boissieu committed
100
101
102
103
NDVI.Thresh = 0.5
Blue.Thresh = 500
NIR.Thresh  = 1500
print("PERFORM RADIOMETRIC FILTERING")
jbferet's avatar
jbferet committed
104
ImPathShade = perform_radiometric_filtering(Input.Image.File, Input.Mask.File, Output.Dir,
Florian de Boissieu's avatar
Florian de Boissieu committed
105
106
107
108
109
                                            NDVI.Thresh = NDVI.Thresh, Blue.Thresh = Blue.Thresh,
                                            NIR.Thresh = NIR.Thresh)
```

## PCA
jbferet's avatar
jbferet committed
110
A pixel-based PCA is run on the input image across the spectral bands to select the most interesting spectral information relative to spectral diversity and remove shaded pixels, spatial noise and sensor artefacts. This PCA band selection left to user judgement, wrinting to a file the bands to keep. The file is automatically created and ready to edit with function `select_PCA_components`. One band number by line is expected in this file. 
Florian de Boissieu's avatar
Florian de Boissieu committed
111

jbferet's avatar
jbferet committed
112
For this example PCA bands 1, 2 and 5 should be kept if writing the following lines in file `selected_components.txt` opened for edition (do not forget carriage return after last value):
Florian de Boissieu's avatar
Florian de Boissieu committed
113
114
115
116
```
1
2
5
jbferet's avatar
jbferet committed
117

Florian de Boissieu's avatar
Florian de Boissieu committed
118
119
120
121
```

Here is the code to perform PCA and select PCA bands:

122
```{r PCA}
Florian de Boissieu's avatar
Florian de Boissieu committed
123
print("PERFORM PCA ON RASTER")
jbferet's avatar
jbferet committed
124
PCA.Files  = perform_PCA(Input.Image.File, ImPathShade, Output.Dir,
Florian de Boissieu's avatar
Florian de Boissieu committed
125
126
                               FilterPCA = TRUE, nbCPU = nbCPU, MaxRAM = MaxRAM)
print("Select PCA components for diversity estimations")
jbferet's avatar
jbferet committed
127
select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
Florian de Boissieu's avatar
Florian de Boissieu committed
128
129
130
131
```

## $\alpha$ and $\beta$ diversity maps

132
```{r alpha and beta diversity maps}
Florian de Boissieu's avatar
Florian de Boissieu committed
133
print("MAP SPECTRAL SPECIES")
jbferet's avatar
jbferet committed
134
map_spectral_species(Input.Image.File, Output.Dir, PCA.Files,
Florian de Boissieu's avatar
Florian de Boissieu committed
135
136
137
138
139
                     nbCPU = nbCPU, MaxRAM = MaxRAM)

print("MAP ALPHA DIVERSITY")
# Index.Alpha   = c('Shannon','Simpson')
Index.Alpha   = c('Shannon')
jbferet's avatar
jbferet committed
140
map_alpha_div(Input.Image.File, Output.Dir, window_size,
Florian de Boissieu's avatar
Florian de Boissieu committed
141
142
143
                    nbCPU = nbCPU, MaxRAM = MaxRAM, Index.Alpha = Index.Alpha)

print("MAP BETA DIVERSITY")
jbferet's avatar
jbferet committed
144
map_beta_div(Input.Image.File, Output.Dir, window_size,
Florian de Boissieu's avatar
Florian de Boissieu committed
145
146
147
                   nbCPU = nbCPU, MaxRAM = MaxRAM)
```

Florian de Boissieu's avatar
Florian de Boissieu committed
148
# $\alpha$ and $\beta$ diversity indices from vector layer
Florian de Boissieu's avatar
Florian de Boissieu committed
149
150
The folowing code computes $\alpha$ and $\beta$ diversity from field plots and extracts the corresponding diversity index from previouly computed rasters in order to have a validation analysis.

Florian de Boissieu's avatar
Florian de Boissieu committed
151
```{r alpha and beta diversity indices from vector layer}
Florian de Boissieu's avatar
Florian de Boissieu committed
152
153
154
155
156
157
158
# location of the spectral species raster needed for validation
TypePCA     = 'SPCA'
Dir.Raster  = file.path(Output.Dir,basename(Input.Image.File),TypePCA,'SpectralSpecies')
Name.Raster = 'SpectralSpecies'
Path.Raster = file.path(Dir.Raster,Name.Raster)

# location of the directory where shapefiles used for validation are saved
Florian de Boissieu's avatar
Florian de Boissieu committed
159
vect        = system.file('extdata', 'VECTOR', package = 'biodivMapR')
Florian de Boissieu's avatar
Florian de Boissieu committed
160
161
162
Shannon.All = list() # ??

# list vector data
jbferet's avatar
jbferet committed
163
Path.Vector         = list.shp(vect)
Florian de Boissieu's avatar
Florian de Boissieu committed
164
165
166
167
Name.Vector         = tools::file_path_sans_ext(basename(Path.Vector))

# read raster data including projection
RasterStack         = stack(Path.Raster)
jbferet's avatar
jbferet committed
168
Projection.Raster   = projection.file(Path.Raster,'raster')
Florian de Boissieu's avatar
Florian de Boissieu committed
169
170

# get alpha and beta diversity indicators corresponding to shapefiles
jbferet's avatar
jbferet committed
171
Biodiv.Indicators           = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
Florian de Boissieu's avatar
Florian de Boissieu committed
172
173
174
175
176
177
178
# if no name
Biodiv.Indicators$Name.Plot = seq(1,length(Biodiv.Indicators$Shannon[[1]]),by = 1)
Shannon.RS                  = c(Biodiv.Indicators$Shannon)[[1]]
```

The tables are then written to tab-seperated files.

179
```{r Write validation}
Florian de Boissieu's avatar
Florian de Boissieu committed
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
# write RS indicators
####################################################
# write indicators for alpha diversity
Path.Results = file.path(Output.Dir, basename(Input.Image.File), TypePCA, 'VALIDATION')
dir.create(Path.Results, showWarnings = FALSE, recursive = TRUE)
ShannonIndexFile <- file.path(Path.Results, "ShannonIndex.tab")
write.table(Shannon.RS, file = ShannonIndexFile, sep = "\t", dec = ".", na = " ", 
            row.names = Biodiv.Indicators$Name.Plot, col.names= F, quote=FALSE)

Results =  data.frame(Name.Vector, Biodiv.Indicators$Richness, Biodiv.Indicators$Fisher,                                Biodiv.Indicators$Shannon, Biodiv.Indicators$Simpson)
names(Results)  = c("ID_Plot", "Species_Richness", "Fisher", "Shannon", "Simpson")
write.table(Results, file = paste(Path.Results,"AlphaDiversity.tab",sep=''), sep="\t", dec=".",               na=" ", row.names = F, col.names= T,quote=FALSE)

# write indicators for beta diversity
BC_mean = Biodiv.Indicators$BCdiss
colnames(BC_mean) = rownames(BC_mean) = Biodiv.Indicators$Name.Plot
write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE)

```