Skip to content

GitLab

Commit f91170f6 authored by jbferet's avatar jbferet
Browse files

updated documentation & vignettes with images

parent 0ba7dee1
Expand all Hide whitespace changes
Inline Side-by-side
doc/tutorial.R
View file @ f91170f6
...@@ -21,7 +21,7 @@ knitr::opts_chunk$set( ...@@ -21,7 +21,7 @@ knitr::opts_chunk$set(
# window_size = 10 # window_size = 10
## ----PCA filtering------------------------------------------------------- ## ----PCA filtering-------------------------------------------------------
# FilterPCA = TRUE # FilterPCA = FALSE
## ----Computing options--------------------------------------------------- ## ----Computing options---------------------------------------------------
# nbCPU = 4 # nbCPU = 4
...@@ -44,11 +44,12 @@ knitr::opts_chunk$set( ...@@ -44,11 +44,12 @@ knitr::opts_chunk$set(
# print("Select PCA components for diversity estimations") # print("Select PCA components for diversity estimations")
# select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE) # select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
## ----alpha and beta diversity maps--------------------------------------- ## ----Spectral species map------------------------------------------------
# print("MAP SPECTRAL SPECIES") # print("MAP SPECTRAL SPECIES")
# map_spectral_species(Input.Image.File, Output.Dir, PCA.Files, # map_spectral_species(Input.Image.File, Output.Dir, PCA.Files,
# nbCPU = nbCPU, MaxRAM = MaxRAM) # nbCPU = nbCPU, MaxRAM = MaxRAM)
#
## ----alpha and beta diversity maps---------------------------------------
# print("MAP ALPHA DIVERSITY") # print("MAP ALPHA DIVERSITY")
# # Index.Alpha = c('Shannon','Simpson') # # Index.Alpha = c('Shannon','Simpson')
# Index.Alpha = c('Shannon') # Index.Alpha = c('Shannon')
...@@ -76,7 +77,7 @@ knitr::opts_chunk$set( ...@@ -76,7 +77,7 @@ knitr::opts_chunk$set(
# #
# # read raster data including projection # # read raster data including projection
# RasterStack = stack(Path.Raster) # RasterStack = stack(Path.Raster)
# Projection.Raster = projection.file(Path.Raster,'raster') # Projection.Raster = get_projection(Path.Raster,'raster')
# #
# # get alpha and beta diversity indicators corresponding to shapefiles # # get alpha and beta diversity indicators corresponding to shapefiles
# Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters) # Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
...@@ -104,3 +105,52 @@ knitr::opts_chunk$set( ...@@ -104,3 +105,52 @@ knitr::opts_chunk$set(
# write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE) # write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE)
# #
## ----PCoA on Field Plots-------------------------------------------------
# # apply ordination using PCoA (same as done for map_beta_div)
# library(labdsv)
# MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
# BetaPCO = pco(MatBCdist, k = 3)
#
## ----plot PCoA & Shannon-------------------------------------------------
# # very uglily assign vegetation type to polygons in shapefiles
# nbSamples = c(6,4,7,7)
# vg = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation')
# Type_Vegetation = c()
# for (i in 1: length(nbSamples)){
# for (j in 1:nbSamples[i]){
# Type_Vegetation = c(Type_Vegetation,vg[i])
# }
# }
#
# # create data frame including alpha and beta diversity
# library(ggplot2)
# Results = data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS)
#
# # plot field data in the PCoA space, with size corresponding to shannon index
# ggplot(Results, aes(x=pco1, y=pco2, color=vgtype,size=shannon)) +
# geom_point(alpha=0.6) +
# scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
# filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA2.png')
# ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
# scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
# dpi = 600, limitsize = TRUE)
#
#
# ggplot(Results, aes(x=pco1, y=pco3, color=vgtype,size=shannon)) +
# geom_point(alpha=0.6) +
# scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
# filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA3.png')
# ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
# scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
# dpi = 600, limitsize = TRUE)
#
# ggplot(Results, aes(x=pco2, y=pco3, color=vgtype,size=shannon)) +
# geom_point(alpha=0.6) +
# scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
# filename = file.path(Path.Results,'BetaDiversity_PcoA2_vs_PcoA3.png')
# ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
# scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
# dpi = 600, limitsize = TRUE)
#
doc/tutorial.Rmd
View file @ f91170f6
...@@ -60,6 +60,8 @@ Input.Mask.File = FALSE ...@@ -60,6 +60,8 @@ Input.Mask.File = FALSE
Output.Dir = 'RESULTS' Output.Dir = 'RESULTS'
``` ```
The image provided with the package is a subset of tile T33NUD acquired by Sentinel-2A satellite over Cameroonese rainforest in January 4th, 2018.
<img align="bottom" width="100%" height="100%" src="../man/figures/01_RGB_S2A_T33NUD_20180104_Subset.png">
## Spatial resolution ## 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 `window_size`, e.g. `window_size = 10` meaning a window of 10x10 pixels. It will be the spatial resolution of the ouput rasters. 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.
...@@ -77,7 +79,7 @@ window_size = 10 ...@@ -77,7 +79,7 @@ window_size = 10
If set to `TRUE`, a second filtering based on PCA outliers is processed. If set to `TRUE`, a second filtering based on PCA outliers is processed.
```{r PCA filtering} ```{r PCA filtering}
FilterPCA = TRUE FilterPCA = FALSE
``` ```
## Computing options ## Computing options
...@@ -107,16 +109,29 @@ ImPathShade = perform_radiometric_filtering(Input.Image.File, Input.Mask.File, O ...@@ -107,16 +109,29 @@ ImPathShade = perform_radiometric_filtering(Input.Image.File, Input.Mask.File, O
``` ```
## PCA ## 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_PCA_components`. One band number by line is expected in this file. 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.
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): The resulting PCA raster are then stored in a binalry file in the results directory, which in our case corresponds to
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/PCA/OutputPCA_8_PCs`
This PCA raster file can be displayed using QGIS or any GIS / image processing software. Here, the PCs corresponding to our image look like this:
<img align="bottom" width="100%" height="100%" src="../man/figures/02_PCs_comp.png">
This PCA band selection left to user judgement, who then writes the band to be kept in a `.txt` file located in the same directory as the PCA raster file. The file is automatically created and ready to edit with function `select_PCA_components`. Each selected band should be identified per line in this file.
The main goal of PC selection is to discard PCs showing no relevant information corresponding to vegetation, or including artifacts possibly explained by sensor properties. It is somehow a subjective process, and we are currently working on automatic selection of these components.
For this example, PCA bands 1, 2, 4, 5, 6 and 8 can be kept if writing the following lines in file `selected_components.txt` opened for edition (do not forget carriage return after last value):
``` ```
1 1
2 2
4
5 5
6
8
``` ```
PC#3 and PC#7 were discarded as the main patterns observed for these components did not correspond to vegetation patterns. As a rule of thumb, between 2 and 6 selected PCs are usually sufficient to catch the main diversity patterns, but this selection strongly depends on the conditions of acquisition, the type of sensor, and the complexity of the vegetation being observed.
Here is the code to perform PCA and select PCA bands: Here is the code to perform PCA and select PCA bands:
```{r PCA} ```{r PCA}
...@@ -127,13 +142,26 @@ print("Select PCA components for diversity estimations") ...@@ -127,13 +142,26 @@ print("Select PCA components for diversity estimations")
select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE) select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
``` ```
## $\alpha$ and $\beta$ diversity maps ## $\alpha$ and $\beta$ diversity maps
```{r alpha and beta diversity maps} The first step towards \alpha and \beta diversity mapping corresponds to the computation of a `SpectralSpecies` map, which identifies the centroid a the cluster ('spectral species') assigned to each pixel in the image, after k-means clustering is performed
```{r Spectral species map}
print("MAP SPECTRAL SPECIES") print("MAP SPECTRAL SPECIES")
map_spectral_species(Input.Image.File, Output.Dir, PCA.Files, map_spectral_species(Input.Image.File, Output.Dir, PCA.Files,
nbCPU = nbCPU, MaxRAM = MaxRAM) nbCPU = nbCPU, MaxRAM = MaxRAM)
```
SpectralSpecies is then stored in a raster file located here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/SpectralSpecies`
$\alpha$ and $\beta$ diversity maps, as well as validation, are based on this `SpectralSpecies` raster.
The code to compute $\alpha$ and $\beta$ diversity maps from this file is as follows:
```{r alpha and beta diversity maps}
print("MAP ALPHA DIVERSITY") print("MAP ALPHA DIVERSITY")
# Index.Alpha = c('Shannon','Simpson') # Index.Alpha = c('Shannon','Simpson')
Index.Alpha = c('Shannon') Index.Alpha = c('Shannon')
...@@ -145,8 +173,20 @@ map_beta_div(Input.Image.File, Output.Dir, window_size, ...@@ -145,8 +173,20 @@ map_beta_div(Input.Image.File, Output.Dir, window_size,
nbCPU = nbCPU, MaxRAM = MaxRAM) nbCPU = nbCPU, MaxRAM = MaxRAM)
``` ```
$\alpha$ and $\beta$ diversity maps are then stored in raster files located here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/ALPHA`
and here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/BETA`
Different rasters can be produced and users are invited to refer to the documentation for more options.
Here, processing our example leads to the following $\alpha$ and $\beta$ diversity maps
<img align="bottom" width="100%" height="100%" src="../man/figures/03_AlphaBeta.png">
# $\alpha$ and $\beta$ diversity indices from vector layer # $\alpha$ and $\beta$ diversity indices from vector layer
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. The following code computes $\alpha$ and $\beta$ diversity from field plots and extracts the corresponding diversity indices from previouly computed `SpectralSpecies` raster in order to perform validation.
```{r alpha and beta diversity indices from vector layer} ```{r alpha and beta diversity indices from vector layer}
# location of the spectral species raster needed for validation # location of the spectral species raster needed for validation
...@@ -165,7 +205,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector)) ...@@ -165,7 +205,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector))
# read raster data including projection # read raster data including projection
RasterStack = stack(Path.Raster) RasterStack = stack(Path.Raster)
Projection.Raster = projection.file(Path.Raster,'raster') Projection.Raster = get_projection(Path.Raster,'raster')
# get alpha and beta diversity indicators corresponding to shapefiles # get alpha and beta diversity indicators corresponding to shapefiles
Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters) Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
...@@ -196,3 +236,71 @@ colnames(BC_mean) = rownames(BC_mean) = Biodiv.Indicators$Name.Plot ...@@ -196,3 +236,71 @@ 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) write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE)
``` ```
These results can then be displayed according to the need for further analysis. Here, for the purpose of illustration, we provide a code in order to visualize the differences among field plots located in the image: we first perform a PCoA on the Bray Curtis dissimilarity matrix computed from the field plots:
```{r PCoA on Field Plots}
# apply ordination using PCoA (same as done for map_beta_div)
library(labdsv)
MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
BetaPCO = pco(MatBCdist, k = 3)
```
The plots corresponding to forested areas with high, medium and low diversity, as well as low vegetation/degraded forest close tomain roads are distributed as follows:
<img align="bottom" width="100%" height="100%" src="../man/figures/04_RGB_FieldLegend.png">
Here, we produce figures in order to locate the different types of vegetation in the PCoA space:
```{r plot PCoA & Shannon}
# very uglily assign vegetation type to polygons in shapefiles
nbSamples = c(6,4,7,7)
vg = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation')
Type_Vegetation = c()
for (i in 1: length(nbSamples)){
for (j in 1:nbSamples[i]){
Type_Vegetation = c(Type_Vegetation,vg[i])
}
}
# create data frame including alpha and beta diversity
library(ggplot2)
Results = data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS)
# plot field data in the PCoA space, with size corresponding to shannon index
ggplot(Results, aes(x=pco1, y=pco2, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA2.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
ggplot(Results, aes(x=pco1, y=pco3, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA3.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
ggplot(Results, aes(x=pco2, y=pco3, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA2_vs_PcoA3.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
```
The resulting figures are displayed here:
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA1_vs_PcoA2.png">
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA1_vs_PcoA3.png">
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA2_vs_PcoA3.png">
doc/tutorial.html
View file @ f91170f6
This diff is collapsed.
examples/example_script.R
View file @ f91170f6
...@@ -55,7 +55,7 @@ Output.Dir = 'RESULTS' ...@@ -55,7 +55,7 @@ Output.Dir = 'RESULTS'
window_size = 10 window_size = 10
# PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower # PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower
FilterPCA = TRUE FilterPCA = FALSE
# type of PCA: # type of PCA:
# PCA: no rescaling of the data # PCA: no rescaling of the data
...@@ -70,7 +70,7 @@ check_data(Input.Image.File) ...@@ -70,7 +70,7 @@ check_data(Input.Image.File)
################################################################################ ################################################################################
## DEFINE PARAMETERS FOR METHOD ## ## DEFINE PARAMETERS FOR METHOD ##
################################################################################ ################################################################################
nbCPU = 4 nbCPU = 2
MaxRAM = 0.5 MaxRAM = 0.5
nbclusters = 50 nbclusters = 50
...@@ -128,7 +128,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector)) ...@@ -128,7 +128,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector))
# read raster data including projection # read raster data including projection
RasterStack = stack(Path.Raster) RasterStack = stack(Path.Raster)
Projection.Raster = projection.file(Path.Raster,'raster') Projection.Raster = get_projection(Path.Raster,'raster')
# get alpha and beta diversity indicators corresponding to shapefiles # get alpha and beta diversity indicators corresponding to shapefiles
Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters) Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
...@@ -157,7 +157,11 @@ write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t ...@@ -157,7 +157,11 @@ write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t
#################################################### ####################################################
# illustrate results # illustrate results
#################################################### ####################################################
# apply ordination using PCoA (same as done for map_beta_div)
library(labdsv) library(labdsv)
MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
BetaPCO = pco(MatBCdist, k = 3)
# very uglily assign vegetation type to polygons in shapefiles # very uglily assign vegetation type to polygons in shapefiles
nbSamples = c(6,4,7,7) nbSamples = c(6,4,7,7)
vg = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation') vg = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation')
...@@ -168,10 +172,6 @@ for (i in 1: length(nbSamples)){ ...@@ -168,10 +172,6 @@ for (i in 1: length(nbSamples)){
} }
} }
# apply ordination unsing PCoA (same as done for map_beta_div)
MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
BetaPCO = pco(MatBCdist, k = 3)
# create data frame including alpha and beta diversity # create data frame including alpha and beta diversity
library(ggplot2) library(ggplot2)
Results = data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS) Results = data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS)
......
vignettes/tutorial.Rmd
View file @ f91170f6
...@@ -60,6 +60,8 @@ Input.Mask.File = FALSE ...@@ -60,6 +60,8 @@ Input.Mask.File = FALSE
Output.Dir = 'RESULTS' Output.Dir = 'RESULTS'
``` ```
The image provided with the package is a subset of tile T33NUD acquired by Sentinel-2A satellite over Cameroonese rainforest in January 4th, 2018.
<img align="bottom" width="100%" height="100%" src="../man/figures/01_RGB_S2A_T33NUD_20180104_Subset.png">
## Spatial resolution ## 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 `window_size`, e.g. `window_size = 10` meaning a window of 10x10 pixels. It will be the spatial resolution of the ouput rasters. 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.
...@@ -77,7 +79,7 @@ window_size = 10 ...@@ -77,7 +79,7 @@ window_size = 10
If set to `TRUE`, a second filtering based on PCA outliers is processed. If set to `TRUE`, a second filtering based on PCA outliers is processed.
```{r PCA filtering} ```{r PCA filtering}
FilterPCA = TRUE FilterPCA = FALSE
``` ```
## Computing options ## Computing options
...@@ -107,16 +109,29 @@ ImPathShade = perform_radiometric_filtering(Input.Image.File, Input.Mask.File, O ...@@ -107,16 +109,29 @@ ImPathShade = perform_radiometric_filtering(Input.Image.File, Input.Mask.File, O
``` ```
## PCA ## 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_PCA_components`. One band number by line is expected in this file. 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.
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): The resulting PCA raster are then stored in a binalry file in the results directory, which in our case corresponds to
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/PCA/OutputPCA_8_PCs`
This PCA raster file can be displayed using QGIS or any GIS / image processing software. Here, the PCs corresponding to our image look like this:
<img align="bottom" width="100%" height="100%" src="../man/figures/02_PCs_comp.png">
This PCA band selection left to user judgement, who then writes the band to be kept in a `.txt` file located in the same directory as the PCA raster file. The file is automatically created and ready to edit with function `select_PCA_components`. Each selected band should be identified per line in this file.
The main goal of PC selection is to discard PCs showing no relevant information corresponding to vegetation, or including artifacts possibly explained by sensor properties. It is somehow a subjective process, and we are currently working on automatic selection of these components.
For this example, PCA bands 1, 2, 4, 5, 6 and 8 can be kept if writing the following lines in file `selected_components.txt` opened for edition (do not forget carriage return after last value):
``` ```
1 1
2 2
4
5 5
6
8
``` ```
PC#3 and PC#7 were discarded as the main patterns observed for these components did not correspond to vegetation patterns. As a rule of thumb, between 2 and 6 selected PCs are usually sufficient to catch the main diversity patterns, but this selection strongly depends on the conditions of acquisition, the type of sensor, and the complexity of the vegetation being observed.
Here is the code to perform PCA and select PCA bands: Here is the code to perform PCA and select PCA bands:
```{r PCA} ```{r PCA}
...@@ -127,13 +142,26 @@ print("Select PCA components for diversity estimations") ...@@ -127,13 +142,26 @@ print("Select PCA components for diversity estimations")
select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE) select_PCA_components(Input.Image.File, Output.Dir, PCA.Files, File.Open = TRUE)
``` ```
## $\alpha$ and $\beta$ diversity maps ## $\alpha$ and $\beta$ diversity maps
```{r alpha and beta diversity maps} The first step towards \alpha and \beta diversity mapping corresponds to the computation of a `SpectralSpecies` map, which identifies the centroid a the cluster ('spectral species') assigned to each pixel in the image, after k-means clustering is performed
```{r Spectral species map}
print("MAP SPECTRAL SPECIES") print("MAP SPECTRAL SPECIES")
map_spectral_species(Input.Image.File, Output.Dir, PCA.Files, map_spectral_species(Input.Image.File, Output.Dir, PCA.Files,
nbCPU = nbCPU, MaxRAM = MaxRAM) nbCPU = nbCPU, MaxRAM = MaxRAM)
```
SpectralSpecies is then stored in a raster file located here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/SpectralSpecies`
$\alpha$ and $\beta$ diversity maps, as well as validation, are based on this `SpectralSpecies` raster.
The code to compute $\alpha$ and $\beta$ diversity maps from this file is as follows:
```{r alpha and beta diversity maps}
print("MAP ALPHA DIVERSITY") print("MAP ALPHA DIVERSITY")
# Index.Alpha = c('Shannon','Simpson') # Index.Alpha = c('Shannon','Simpson')
Index.Alpha = c('Shannon') Index.Alpha = c('Shannon')
...@@ -145,8 +173,20 @@ map_beta_div(Input.Image.File, Output.Dir, window_size, ...@@ -145,8 +173,20 @@ map_beta_div(Input.Image.File, Output.Dir, window_size,
nbCPU = nbCPU, MaxRAM = MaxRAM) nbCPU = nbCPU, MaxRAM = MaxRAM)
``` ```
$\alpha$ and $\beta$ diversity maps are then stored in raster files located here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/ALPHA`
and here:
`RESULTS/S2A_T33NUD_20180104_Subset/SPCA/BETA`
Different rasters can be produced and users are invited to refer to the documentation for more options.
Here, processing our example leads to the following $\alpha$ and $\beta$ diversity maps
<img align="bottom" width="100%" height="100%" src="../man/figures/03_AlphaBeta.png">
# $\alpha$ and $\beta$ diversity indices from vector layer # $\alpha$ and $\beta$ diversity indices from vector layer
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. The following code computes $\alpha$ and $\beta$ diversity from field plots and extracts the corresponding diversity indices from previouly computed `SpectralSpecies` raster in order to perform validation.
```{r alpha and beta diversity indices from vector layer} ```{r alpha and beta diversity indices from vector layer}
# location of the spectral species raster needed for validation # location of the spectral species raster needed for validation
...@@ -165,7 +205,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector)) ...@@ -165,7 +205,7 @@ Name.Vector = tools::file_path_sans_ext(basename(Path.Vector))
# read raster data including projection # read raster data including projection
RasterStack = stack(Path.Raster) RasterStack = stack(Path.Raster)
Projection.Raster = projection.file(Path.Raster,'raster') Projection.Raster = get_projection(Path.Raster,'raster')
# get alpha and beta diversity indicators corresponding to shapefiles # get alpha and beta diversity indicators corresponding to shapefiles
Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters) Biodiv.Indicators = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
...@@ -196,3 +236,71 @@ colnames(BC_mean) = rownames(BC_mean) = Biodiv.Indicators$Name.Plot ...@@ -196,3 +236,71 @@ 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) write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE)
``` ```
These results can then be displayed according to the need for further analysis. Here, for the purpose of illustration, we provide a code in order to visualize the differences among field plots located in the image: we first perform a PCoA on the Bray Curtis dissimilarity matrix computed from the field plots:
```{r PCoA on Field Plots}
# apply ordination using PCoA (same as done for map_beta_div)
library(labdsv)
MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
BetaPCO = pco(MatBCdist, k = 3)
```
The plots corresponding to forested areas with high, medium and low diversity, as well as low vegetation/degraded forest close tomain roads are distributed as follows:
<img align="bottom" width="100%" height="100%" src="../man/figures/04_RGB_FieldLegend.png">
Here, we produce figures in order to locate the different types of vegetation in the PCoA space:
```{r plot PCoA & Shannon}
# very uglily assign vegetation type to polygons in shapefiles
nbSamples = c(6,4,7,7)
vg = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation')
Type_Vegetation = c()
for (i in 1: length(nbSamples)){
for (j in 1:nbSamples[i]){
Type_Vegetation = c(Type_Vegetation,vg[i])
}
}
# create data frame including alpha and beta diversity
library(ggplot2)
Results = data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS)
# plot field data in the PCoA space, with size corresponding to shannon index
ggplot(Results, aes(x=pco1, y=pco2, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA2.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
ggplot(Results, aes(x=pco1, y=pco3, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA3.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
ggplot(Results, aes(x=pco2, y=pco3, color=vgtype,size=shannon)) +
geom_point(alpha=0.6) +
scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
filename = file.path(Path.Results,'BetaDiversity_PcoA2_vs_PcoA3.png')
ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
dpi = 600, limitsize = TRUE)
```
The resulting figures are displayed here:
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA1_vs_PcoA2.png">
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA1_vs_PcoA3.png">
<img align="bottom" width="100%" height="100%" src="../man/figures/BetaDiversity_PcoA2_vs_PcoA3.png">
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