Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@Marie59
Marie59 committed Dec 16, 2022
1 parent bbfda7d commit 1bde5e1
Show file tree
Hide file tree
Showing 12 changed files with 1,057 additions and 531 deletions.
141 changes: 141 additions & 0 deletions alpha_beta.r
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,141 @@
#Rscript

###########################################
## Mapping alpha and beta diversity ##
###########################################

#####Packages : stars
# utils
# biodivmapr
# raster
# sf
# mapview
# leafpop
# RColorBrewer
# labdsv
# rgdal
# ggplot2
# gridExtra

#####Load arguments

args <- commandArgs(trailingOnly = TRUE)

#####Import the S2 data

if (length(args) < 1) {
stop("This tool needs at least 1 argument")
}else {
data_raster <- args[1]
rasterheader <- args[2]
data <- args[3]
# type of PCA:
# PCA: no rescaling of the data
# SPCA: rescaling of the data
typepca <- as.character(args[4])
alpha <- as.logical(args[5])
beta <- as.logical(args[6])
funct <- as.logical(args[7])
all <- as.logical(args[8])
source(args[9])
}

################################################################################
## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ##
################################################################################
if (data_raster == "") {
#Create a directory where to unzip your folder of data
dir.create("data_dir")
unzip(data, exdir = "data_dir")
# Path to raster
data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])

} else {
input_image_file <- file.path(getwd(), data_raster, fsep = "/")
input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
}

################################################################################
## PROCESS IMAGE ##
################################################################################
# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels

print("PERFORM PCA ON RASTER")
pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)

pca_files <- pca_output$PCA_Files
pix_per_partition <- pca_output$Pix_Per_Partition
nb_partitions <- pca_output$nb_partitions
# path for the updated mask
input_mask_file <- pca_output$MaskPath


selected_pcs <- seq(1, dim(raster::stack(input_image_file))[3])

selected_pcs <- all(selected_pcs)
################################################################################
## MAP ALPHA AND BETA DIVERSITY ##
################################################################################
print("MAP SPECTRAL SPECIES")

kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = input_image_file, Output_Dir = output_dir, PCA_Files = pca_files, Input_Mask_File = input_mask_file, SelectedPCs = selected_pcs, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, TypePCA = typepca, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters)

image_name <- tools::file_path_sans_ext(basename(input_image_file))
if (alpha == TRUE || beta == TRUE || all == TRUE) {
## alpha
print("MAP ALPHA DIVERSITY")
index_alpha <- c("Shannon")
alpha_div <- biodivMapR::map_alpha_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, Index_Alpha = index_alpha, nbclusters = nbclusters, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE)

alpha_zip <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres.zip")
alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA")
unzip(alpha_zip, exdir = alpha_path)
alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres")
alpha_raster <- raster::raster(alpha_path)
get_alpha <- convert_raster(alpha_raster)

if (alpha == TRUE || all == TRUE) {
colnames(get_alpha) <- c("Alpha", "longitude", "latitude")
plot_indices(get_alpha, titre = "Alpha")

write.table(get_alpha, file = "alpha.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
}
if (beta == TRUE || all == TRUE) {
## beta
print("MAP BETA DIVERSITY")
beta_div <- biodivMapR::map_beta_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters)

beta_path <- file.path(output_dir, image_name, typepca, "BETA", "BetaDiversity_BCdiss_PCO_10")
beta_raster <- raster::raster(beta_path)
get_beta <- convert_raster(beta_raster)

colnames(get_beta) <- c("Beta", "longitude", "latitude")
plot_indices(get_beta, titre = "Beta")

write.table(get_beta, file = "beta.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
}
}


################################################################################
## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ##
################################################################################

if (funct == TRUE || all == TRUE) {
mapper <- biodivMapR::map_functional_div(Original_Image_File = input_image_file, Functional_File = pca_files, Selected_Features = selected_pcs, Output_Dir = output_dir, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, TypePCA = typepca, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE)

funct_zip <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres.zip")
funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL")
unzip(funct_zip, exdir = funct_path)
funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres")
funct_raster <- raster::raster(funct_path)
get_funct <- convert_raster(funct_raster)

colnames(get_funct) <- c("Functionnal", "longitude", "latitude")
plot_indices(get_funct, titre = "Functionnal")

write.table(get_funct, file = "Functionnal.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
}
161 changes: 161 additions & 0 deletions alpha_beta.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
<tool id="srs_diversity_maps" name="Map diversity" version="@VERSION@" profile="20.01">
<description>from remote sensing data</description>
<macros>
<import>macrobis.xml</import>
</macros>
<expand macro="SRS_requirements"/>
<command detect_errors="exit_code"><![CDATA[
#import re
#if $method.origin == 'envi_bil':
#set input_raster = $method.input_raster
#set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier))
#set input_header = $method.input_header
#set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier))
ln '${input_raster}' '${input_raster_identifier}' &&
ln '${input_header}' '${input_header_identifier}' &&
#end if
Rscript
'$__tool_directory__/alpha_beta.r'
#if $method.origin == 'envi_bil':
'$input_raster_identifier'
'$input_header_identifier'
''
#else:
''
''
'$method.input'
#end if
'$typepca'
#if $type == 'alpha':
'TRUE'
'FALSE'
'FALSE'
'FALSE'
#else if $type == 'beta':
'FALSE'
'TRUE'
'FALSE'
'FALSE'
#else if $type == 'funct':
'FALSE'
'FALSE'
'TRUE'
'FALSE'
#else:
'FALSE'
'FALSE'
'FALSE'
'TRUE'
#end if
'$__tool_directory__/functions.r'
'$output_alpha'
'$output_beta'
'$output_funct'
'$plots'
]]>
</command>
<inputs>
<conditional name="method">
<param name="origin" type="select" label="In which format are your data ?">
<option value="zipper">The data you are using are in a zip folder Reflectance</option>
<option value="envi_bil">Your already have the files in ENVI BIL format</option>
</param>
<when value="zipper">
<param name="input" type="data" format="zip" multiple="true" label="Input data"/>
</when>
<when value="envi_bil">
<param name="input_raster" type="data" format="bil" label="Input raster" help="It can be the raw data in bil or the PCA raster layer in bil"/>
<param name="input_header" type="data" format="hdr" label="Input header"/>
</when>
</conditional>
<param name="typepca" type="select" label="Do you want to do a PCA or a SPCA ?" display="radio" help="If you choose PCA there is no rescaling of the data as oppposed as if you choose SPCA">
<option value="SPCA">SPCA</option>
<option value="PCA">PCA</option>
</param>
<param name="type" type="select" label="Alpha, beta, functional diversity and comparison plot and map" display="radio">
<option value="alpha">Alpha diversity map</option>
<option value="beta">Beta diversity map</option>
<option value="funct">Functional diversity map</option>
<option value="all">All of the above</option>
</param>
</inputs>
<outputs>
<data name="output_alpha" from_work_dir="alpha.tabular" format="tabular" label="Alpha diversity">
<filter> type == 'alpha' or type == 'all'</filter>
</data>
<data name="output_beta" from_work_dir="beta.tabular" format="tabular" label="Beta diversity">
<filter> type == 'beta' or type == 'all'</filter>
</data>
<data name="output_funct" from_work_dir="Functionnal.tabular" format="tabular" label="Functionnal diversity">
<filter> type == 'funct' or type == 'all'</filter>
</data>
<collection type="list" name="plots" label="${type} plot">
<discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
</collection>
</outputs>
<tests>
<test>
<param name="origin" value="envi_bil"/>
<param name="input_raster" value="S2A_Subset"/>
<param name="input_header" value="S2A_Subset.hdr"/>
<param name="type" value="all"/>
<output name="output_alpha">
<assert_contents>
<has_n_columns n="3"/>
</assert_contents>
</output>
<output name="output_beta">
<assert_contents>
<has_n_columns n="3"/>
</assert_contents>
</output>
<output name="output_funct">
<assert_contents>
<has_n_columns n="3"/>
</assert_contents>
</output>
<output_collection name="plots" type="list" count="3"/>
</test>
</tests>
<help><![CDATA[
========================================================================
Process satellite remote sensing data to produce biodiversity indicators
========================================================================
**What it does**
Féret and Asner (2014) developed a method for **tropical forest** diversity mapping based on very high spatial resolution airborne imaging spectroscopy.
The goal of this tool using the package biodivMapR is to produce (spectral) diversity maps based on (optical) images.
**Input description**
It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
The image is an ENVI raster including :
- A binary file (which has no extension here).
- A header file (with .hdr extension).
The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
In order to get such input we advise to use the tool preprocessing sentinel 2 data.
+--------------+----------+
| BIL | ENVI HDR |
+==============+==========+
| raster stack | Metadata |
+--------------+----------+
| ... | ... |
+--------------+----------+
**Output**
- Three tabulars : alpha, beta, functionnal each of them with 3 colomns latitude, longitude and the indice.
- Three png graph for each indice
]]></help>
<expand macro="SRS_BDMRref"/>
</tool>
Loading

0 comments on commit 1bde5e1

Please sign in to comment.