Commit 93023eed authored by Florian de Boissieu's avatar Florian de Boissieu
Browse files

update vignette and README

parent 3da63885
......@@ -9,5 +9,14 @@ devtools::install_git('https://gitlab.irstea.fr/jean-baptiste.feret/diversitymap
```
# 2 Tutorial
A tutorial script can be found in `Main_DiversityMapping.R`, showing the main steps of the processing. A tutorial vignette should soon be available.
A tutorial vignette showing the main steps of the processing can be visualised with the following command line:
```
rstudioapi::viewer(system.file('doc', 'tutorial.html', package='DiversityMappR'))
```
or for non Rstudio session:
```
vignette('tutorial', package='DiversityMappR')
```
The corresponding script is available in file `examples/tutorial.R`.
## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval=FALSE
)
## ----Input / Output files------------------------------------------------
# Input.Image.File = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'DiversityMappR')
# Check.Data.Format(Input.Image.File)
#
# Input.Image.File = Convert.Raster2BIL(Raster.Path = Input.Image.File,
# Sensor = 'SENTINEL_2A',
# Convert.Integer = TRUE,
# Output.Directory = '~/test')
# Input.Mask.File = FALSE
#
# Output.Dir = 'RESULTS'
## ----Spatial resolution--------------------------------------------------
# Spatial.Res = 10
## ----PCA filtering-------------------------------------------------------
# FilterPCA = TRUE
## ----Computing options---------------------------------------------------
# nbCPU = 4
# MaxRAM = 0.5
# nbclusters = 50
## ----Mask non vegetated / shaded / cloudy pixels-------------------------
# NDVI.Thresh = 0.5
# Blue.Thresh = 500
# NIR.Thresh = 1500
# print("PERFORM RADIOMETRIC FILTERING")
# ImPathShade = Perform.Radiometric.Filtering(Input.Image.File, Input.Mask.File, Output.Dir,
# NDVI.Thresh = NDVI.Thresh, Blue.Thresh = Blue.Thresh,
# NIR.Thresh = NIR.Thresh)
## ----PCA-----------------------------------------------------------------
# print("PERFORM PCA ON RASTER")
# PCA.Files = Perform.PCA.Image(Input.Image.File, ImPathShade, Output.Dir,
# FilterPCA = TRUE, nbCPU = nbCPU, MaxRAM = MaxRAM)
# print("Select PCA components for diversity estimations")
# Select.Components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
## ----alpha and beta diversity maps---------------------------------------
# print("MAP SPECTRAL SPECIES")
# Map.Spectral.Species(Input.Image.File, Output.Dir, PCA.Files,
# nbCPU = nbCPU, MaxRAM = MaxRAM)
#
# print("MAP ALPHA DIVERSITY")
# # Index.Alpha = c('Shannon','Simpson')
# Index.Alpha = c('Shannon')
# Map.Alpha.Diversity(Input.Image.File, Output.Dir, Spatial.Res,
# nbCPU = nbCPU, MaxRAM = MaxRAM, Index.Alpha = Index.Alpha)
#
# print("MAP BETA DIVERSITY")
# Map.Beta.Diversity(Input.Image.File, Output.Dir, Spatial.Res,
# nbCPU = nbCPU, MaxRAM = MaxRAM)
## ----Validation----------------------------------------------------------
# # 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
# vect = system.file('extdata', 'VECTOR', package = 'DiversityMappR')
# Shannon.All = list() # ??
#
# # list vector data
# Path.Vector = Get.List.Shp(vect)
# Name.Vector = tools::file_path_sans_ext(basename(Path.Vector))
#
# # read raster data including projection
# RasterStack = stack(Path.Raster)
# Projection.Raster = Get.Projection(Path.Raster,'raster')
#
# # get alpha and beta diversity indicators corresponding to shapefiles
# Biodiv.Indicators = Get.Diversity.From.Plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
# # if no name
# Biodiv.Indicators$Name.Plot = seq(1,length(Biodiv.Indicators$Shannon[[1]]),by = 1)
# Shannon.RS = c(Biodiv.Indicators$Shannon)[[1]]
## ----Write validation----------------------------------------------------
# # 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)
#
......@@ -11,6 +11,15 @@ vignette: >
\usepackage[utf8]{inputenc}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval=FALSE
)
```
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
......@@ -35,8 +44,10 @@ If not they should be converted with function `Convert.Raster2BIL`.
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`.
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.
```
```{r Input / Output files}
Input.Image.File = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'DiversityMappR')
Check.Data.Format(Input.Image.File)
......@@ -45,14 +56,11 @@ Input.Image.File = Convert.Raster2BIL(Raster.Path = Input.Image.File,
Convert.Integer = TRUE,
Output.Directory = '~/test')
Input.Mask.File = FALSE
```
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.
```
Output.Dir = 'RESULTS'
```
## Spatial resolution
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 `Spatial.Res`, e.g. `Spatial.Res = 10` meaning a window of 10x10 pixels. It will be the spatial resolution of the ouput rasters.
......@@ -61,14 +69,14 @@ A Spatial.Res too small results in low number of pixels per spatial unit, hence
In this example, the spatial resolution of the input raster. Setting `Spatial.Res = 10` will result in diversity maps of spatial resolution 100x100m.
```
```{r Spatial resolution}
Spatial.Res = 10
```
## PCA filtering
If set to `TRUE`, a second filtering based on PCA outliers is processed.
```
```{r PCA filtering}
FilterPCA = TRUE
```
......@@ -79,7 +87,7 @@ The use of computing ressources can be controled with the following parameters:
* `MaxRAM` controls the size in GB of the input image chunks processed by each thread,
* `nbclusters` controls the number of clusters (or centroids) used in kmeans of each thread. The larger the value the longer the computation time.
```
```{r Computing options}
nbCPU = 4
MaxRAM = 0.5
nbclusters = 50
......@@ -88,7 +96,7 @@ nbclusters = 50
# Main processing worflow
## Mask non vegetated / shaded / cloudy pixels
```
```{r Mask non vegetated / shaded / cloudy pixels}
NDVI.Thresh = 0.5
Blue.Thresh = 500
NIR.Thresh = 1500
......@@ -110,7 +118,7 @@ For this example PCA bands 1, 2 and 5 should be kept writting the following line
Here is the code to perform PCA and select PCA bands:
```
```{r PCA}
print("PERFORM PCA ON RASTER")
PCA.Files = Perform.PCA.Image(Input.Image.File, ImPathShade, Output.Dir,
FilterPCA = TRUE, nbCPU = nbCPU, MaxRAM = MaxRAM)
......@@ -120,7 +128,7 @@ Select.Components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
## $\alpha$ and $\beta$ diversity maps
```
```{r alpha and beta diversity maps}
print("MAP SPECTRAL SPECIES")
Map.Spectral.Species(Input.Image.File, Output.Dir, PCA.Files,
nbCPU = nbCPU, MaxRAM = MaxRAM)
......@@ -139,7 +147,7 @@ Map.Beta.Diversity(Input.Image.File, Output.Dir, Spatial.Res,
# Validation
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.
```
```{r Validation}
# location of the spectral species raster needed for validation
TypePCA = 'SPCA'
Dir.Raster = file.path(Output.Dir,basename(Input.Image.File),TypePCA,'SpectralSpecies')
......@@ -167,7 +175,7 @@ Shannon.RS = c(Biodiv.Indicators$Shannon)[[1]]
The tables are then written to tab-seperated files.
```
```{r Write validation}
# write RS indicators
####################################################
# write indicators for alpha diversity
......
......@@ -15,15 +15,15 @@ Perform.PCA.Image(ImPath, ImPathShade, Output.Dir,
\item{TypePCA}{Type of PCA (PCA, SPCA, NLPCA...)}
\item{FilterPCA}{2nd filtering based on PCA}
\item{FilterPCA}{boolean. If TRUE 2nd filtering based on PCA}
\item{Excluded.WL}{Water Vapor Absorption domains (in nanometers). Can also be used to exclude spectific domains}
\item{Excluded.WL}{boolean. Water Vapor Absorption domains (in nanometers). Can also be used to exclude spectific domains}
\item{NbIter}{number of iteration to estimate diversity from the raster}
\item{NbIter}{numeric. Number of iteration to estimate diversity from the raster.}
\item{nbCPU}{number fo CPUs in use.}
\item{nbCPU}{numeric. Number fo CPUs in use.}
\item{MaxRAM}{maximum size of chunk in GB to limit RAM allocation when reading image file.}
\item{MaxRAM}{numeric. Maximum size of chunk in GB to limit RAM allocation when reading image file.}
\item{ImNames}{Path and name of the images to be processed}
}
......
---
title: "Tutorial for DiversityMappR"
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}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval=FALSE
)
```
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
The input images are expected to be in ENVI HDR format, BIL interleaved. To check if the flormat is good use fucntion `Check.Data.Format`.
If not they should be converted with function `Convert.Raster2BIL`.
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`.
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.
```{r Input / Output files}
Input.Image.File = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'DiversityMappR')
Check.Data.Format(Input.Image.File)
Input.Image.File = Convert.Raster2BIL(Raster.Path = Input.Image.File,
Sensor = 'SENTINEL_2A',
Convert.Integer = TRUE,
Output.Directory = '~/test')
Input.Mask.File = FALSE
Output.Dir = 'RESULTS'
```
## Spatial resolution
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 `Spatial.Res`, e.g. `Spatial.Res = 10` meaning a window of 10x10 pixels. It will be the spatial resolution of the ouput rasters.
As a rule of thumb, spatial units between 0.25 and 4 ha usually match with ground data.
A Spatial.Res too small results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image.
In this example, the spatial resolution of the input raster. Setting `Spatial.Res = 10` will result in diversity maps of spatial resolution 100x100m.
```{r Spatial resolution}
Spatial.Res = 10
```
## PCA filtering
If set to `TRUE`, a second filtering based on PCA outliers is processed.
```{r PCA filtering}
FilterPCA = TRUE
```
## Computing options
The use of computing ressources can be controled with the following parameters:
* `nbCPU` controls the parallelisation of the processing,
* `MaxRAM` controls the size in GB of the input image chunks processed by each thread,
* `nbclusters` controls the number of clusters (or centroids) used in kmeans of each thread. The larger the value the longer the computation time.
```{r Computing options}
nbCPU = 4
MaxRAM = 0.5
nbclusters = 50
```
# Main processing worflow
## Mask non vegetated / shaded / cloudy pixels
```{r Mask non vegetated / shaded / cloudy pixels}
NDVI.Thresh = 0.5
Blue.Thresh = 500
NIR.Thresh = 1500
print("PERFORM RADIOMETRIC FILTERING")
ImPathShade = Perform.Radiometric.Filtering(Input.Image.File, Input.Mask.File, Output.Dir,
NDVI.Thresh = NDVI.Thresh, Blue.Thresh = Blue.Thresh,
NIR.Thresh = NIR.Thresh)
```
## PCA
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.Components`. One band number by line is expected in this file.
For this example PCA bands 1, 2 and 5 should be kept writting the following lines in file `selected_components.txt` opened for edition:
```
1
2
5
```
Here is the code to perform PCA and select PCA bands:
```{r PCA}
print("PERFORM PCA ON RASTER")
PCA.Files = Perform.PCA.Image(Input.Image.File, ImPathShade, Output.Dir,
FilterPCA = TRUE, nbCPU = nbCPU, MaxRAM = MaxRAM)
print("Select PCA components for diversity estimations")
Select.Components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
```
## $\alpha$ and $\beta$ diversity maps
```{r alpha and beta diversity maps}
print("MAP SPECTRAL SPECIES")
Map.Spectral.Species(Input.Image.File, Output.Dir, PCA.Files,
nbCPU = nbCPU, MaxRAM = MaxRAM)
print("MAP ALPHA DIVERSITY")
# Index.Alpha = c('Shannon','Simpson')
Index.Alpha = c('Shannon')
Map.Alpha.Diversity(Input.Image.File, Output.Dir, Spatial.Res,
nbCPU = nbCPU, MaxRAM = MaxRAM, Index.Alpha = Index.Alpha)
print("MAP BETA DIVERSITY")
Map.Beta.Diversity(Input.Image.File, Output.Dir, Spatial.Res,
nbCPU = nbCPU, MaxRAM = MaxRAM)
```
# Validation
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.
```{r Validation}
# 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
vect = system.file('extdata', 'VECTOR', package = 'DiversityMappR')
Shannon.All = list() # ??
# list vector data
Path.Vector = Get.List.Shp(vect)
Name.Vector = tools::file_path_sans_ext(basename(Path.Vector))
# read raster data including projection
RasterStack = stack(Path.Raster)
Projection.Raster = Get.Projection(Path.Raster,'raster')
# get alpha and beta diversity indicators corresponding to shapefiles
Biodiv.Indicators = Get.Diversity.From.Plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
# 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.
```{r Write validation}
# 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)
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment