Commit 3ccebd9a authored by jbferet's avatar jbferet
Browse files

- changed function split.line into split_line

- added imports (prcomp...)
- simplified 'PCA_Info.RData', now only called once in map_spectral_species
- changed dots into underscores for names of most variables
- added missing input variables
parent 81325e1c
Package: biodivMapR
Title: biodivMapR: an R package for α- and β-diversity mapping using remotely-sensed images
Version: 0.9.0
Authors@R: c( person("Jean-Baptiste", "Feret", email = "jb.feret@teledetection.fr", role = c("aut", "cre")),
person("Florian", "de Boissieu", email = "florian.deboissieu@irstea.fr", role = c("ctb"), comment = "clean code, format as package"))
Authors@R: c(person(given = "Jean-Baptiste",
family = "Feret",
email = "jb.feret@teledetection.fr",
role = c("aut", "cre"),
comment = c(ORCID = "0000-0002-0151-1334")),
person( given = "Florian",
family = "de Boissieu",
email = "florian.deboissieu@irstea.fr",
role = c("ctb"),
comment = "clean code, format as package"))
Description: this packages allows processing image data based on the method described in the following publication:
Féret, J.-B., Asner, G.P., 2014. Mapping tropical forest canopy diversity using high-fidelity imaging spectroscopy. Ecol. Appl. 24, 1289–1296. https://doi.org/10.1890/13-1824.1
It expects an image file as input, with a specific data format.
ENVI HDR image with BIL interleave required
It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
Encoding: UTF-8
License: GPL-3
License: GPL-3 + file LICENSE
LazyData: true
Imports:
dissUtils,
......@@ -19,10 +26,12 @@ Imports:
labdsv,
matlab,
matrixStats,
methods,
raster,
rgdal,
R.utils,
snow,
sp,
stringr,
tools,
vegan,
......
# Generated by roxygen2: do not edit by hand
S3method(split,line)
export(check_data)
export(diversity_from_plots)
export(get_projection)
......@@ -12,17 +11,21 @@ export(perform_PCA)
export(perform_radiometric_filtering)
export(raster2BIL)
export(select_PCA_components)
export(split_line)
import(raster)
import(tools)
importFrom(dissUtils,diss)
importFrom(ecodist,nmds)
importFrom(fields,rdist)
importFrom(future,multiprocess)
importFrom(future,plan)
importFrom(future,sequential)
importFrom(future.apply,future_lapply)
importFrom(labdsv,pco)
importFrom(matlab,ones)
importFrom(matlab,padarray)
importFrom(matrixStats,rowSds)
importFrom(methods,is)
importFrom(raster,projection)
importFrom(raster,raster)
importFrom(raster,shapefile)
......@@ -30,7 +33,13 @@ importFrom(raster,writeRaster)
importFrom(rgdal,readOGR)
importFrom(rgdal,writeOGR)
importFrom(snow,splitRows)
importFrom(sp,spTransform)
importFrom(stats,as.dist)
importFrom(stats,kmeans)
importFrom(stats,prcomp)
importFrom(stats,sd)
importFrom(stringr,str_count)
importFrom(utils,file.edit)
importFrom(utils,find)
importFrom(utils,read.table)
importFrom(vegan,fisher.alpha)
......@@ -12,158 +12,157 @@
#' converts a raster into BIL format as expected by DivMapping codes
#'
#' @param Raster.Path character. Full path for the raster to be converted
#' @param Raster_Path character. Full path for the raster to be converted
#' @param Sensor character. Name of the sensor
#' @param Convert.Integer boolean. Should data be converted into integer ?
#' @param Multiplying.Factor numeric. Multiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000).
#' @param Output.Directory character. Path to output directory.
#' @param Multiplying.Factor.Last numeric. Multiplying factor for last band.
#' @param Convert_Integer boolean. Should data be converted into integer ?
#' @param Multiplying_Factor numeric. Multiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000).
#' @param Output_Dir character. Path to output directory.
#' @param Multiplying_Factor_Last numeric. Multiplying factor for last band.
#'
#' @return Output.Path path for the image converted into ENVI BIL format
#' @return Output_Path path for the image converted into ENVI BIL format
#' @import raster
#' @import tools
#' @export
raster2BIL <- function(Raster.Path, Sensor = "unknown", Output.Directory = FALSE, Convert.Integer = TRUE, Multiplying.Factor = 1, Multiplying.Factor.Last = 1) {
raster2BIL <- function(Raster_Path, Sensor = "unknown", Output_Dir = FALSE, Convert_Integer = TRUE, Multiplying_Factor = 1, Multiplying_Factor_Last = 1) {
# get directory and file name of original image
Input.File <- basename(Raster.Path)
Input.Dir <- dirname(Raster.Path)
Input_File <- basename(Raster_Path)
Input_Dir <- dirname(Raster_Path)
# define path where data will be stored
if (Output.Directory == FALSE) {
Output.Path <- paste(Input.Dir, "/Converted_DivMAP/", file_path_sans_ext(Input.File), sep = "")
if (Output_Dir == FALSE) {
Output_Path <- paste(Input_Dir, "/Converted_DivMAP/", file_path_sans_ext(Input_File), sep = "")
} else {
dir.create(Output.Directory, showWarnings = FALSE, recursive = TRUE)
Output.Path <- paste(Output.Directory, "/", file_path_sans_ext(Input.File), sep = "")
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
Output_Path <- paste(Output_Dir, "/", file_path_sans_ext(Input_File), sep = "")
}
message("The converted file will be written in the following location:")
print(Output.Path)
Output.Dir <- dirname(Output.Path)
dir.create(Output.Dir, showWarnings = FALSE, recursive = TRUE)
print(Output_Path)
Output_Dir <- dirname(Output_Path)
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
# apply multiplying factors
message("reading initial file")
Output.Img <- Multiplying.Factor * brick(Raster.Path)
Last.Band.Name <- Output.Img@data@names[length(Output.Img@data@names)]
Output.Img[[Last.Band.Name]] <- Multiplying.Factor.Last * Output.Img[[Last.Band.Name]]
Output_Img <- Multiplying_Factor * brick(Raster_Path)
Last_Band_Name <- Output_Img@data@names[length(Output_Img@data@names)]
Output_Img[[Last_Band_Name]] <- Multiplying_Factor_Last * Output_Img[[Last_Band_Name]]
# convert into integer
if (Convert.Integer == TRUE) {
Output.Img <- round(Output.Img)
if (Convert_Integer == TRUE) {
Output_Img <- round(Output_Img)
}
# write raster
message("writing converted file")
if (Convert.Integer == TRUE) {
r <- writeRaster(Output.Img, filename = Output.Path, format = "EHdr", overwrite = TRUE, datatype = "INT2S")
if (Convert_Integer == TRUE) {
r <- writeRaster(Output_Img, filename = Output_Path, format = "EHdr", overwrite = TRUE, datatype = "INT2S")
} else {
r <- writeRaster(Output.Img, filename = Output.Path, format = "EHdr", overwrite = TRUE)
r <- writeRaster(Output_Img, filename = Output_Path, format = "EHdr", overwrite = TRUE)
}
hdr(r, format = "ENVI")
# remove unnecessary files
File2Remove <- paste(Output.Path, ".aux.xml", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".aux.xml", sep = "")
File2Remove <- paste(Output_Path, ".aux.xml", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".aux.xml", sep = "")
file.remove(File2Remove)
File2Remove <- paste(Output.Path, ".prj", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".prj", sep = "")
File2Remove <- paste(Output_Path, ".prj", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".prj", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Remove <- paste(Output.Path, ".sta", sep = "")
File2Remove <- paste(file_path_sans_ext(Output.Path), ".sta", sep = "")
File2Remove <- paste(Output_Path, ".sta", sep = "")
File2Remove <- paste(file_path_sans_ext(Output_Path), ".sta", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Remove <- paste(Output.Path, ".stx", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".stx", sep = "")
File2Remove <- paste(Output_Path, ".stx", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".stx", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Rename <- paste(file_path_sans_ext(Output.Path), ".hdr", sep = "")
File2Rename2 <- paste(Output.Path, ".hdr", sep = "")
File2Rename <- paste(file_path_sans_ext(Output_Path), ".hdr", sep = "")
File2Rename2 <- paste(Output_Path, ".hdr", sep = "")
if (file.exists(File2Rename)) {
file.rename(from = File2Rename, to = File2Rename2)
}
# change dot into underscore
Output.Path.US <- file.path(
dirname(Output.Path),
gsub(basename(Output.Path), pattern = "[.]", replacement = "_")
Output_Path_US <- file.path(
dirname(Output_Path),
gsub(basename(Output_Path), pattern = "[.]", replacement = "_")
)
if (!Output.Path.US == Output.Path) {
file.rename(from = Output.Path, to = Output.Path.US)
if (!Output_Path_US == Output_Path) {
file.rename(from = Output_Path, to = Output_Path_US)
}
Output.Path.US.HDR <- paste0(Output.Path.US, ".hdr")
if (!Output.Path.US.HDR == paste0(Output.Path, ".hdr")) {
file.rename(from = paste0(Output.Path, ".hdr"), to = Output.Path.US.HDR)
Output_Path_US_HDR <- paste0(Output_Path_US, ".hdr")
if (!Output_Path_US_HDR == paste0(Output_Path, ".hdr")) {
file.rename(from = paste0(Output_Path, ".hdr"), to = Output_Path_US_HDR)
### UTILITY?? ###
file.rename(from = Output.Path.US.HDR, to = Output.Path.US.HDR)
file.rename(from = Output_Path_US_HDR, to = Output_Path_US_HDR)
}
if (!Sensor == "unknown") {
HDR.Temp.Path <- system.file("extdata", "HDR", paste0(Sensor, ".hdr"), package = "biodivMapR")
if (file.exists(HDR.Temp.Path)) {
HDR_Temp_Path <- system.file("extdata", "HDR", paste0(Sensor, ".hdr"), package = "biodivMapR")
if (file.exists(HDR_Temp_Path)) {
message("reading header template corresponding to the sensor located here:")
print(HDR.Temp.Path)
print(HDR_Temp_Path)
# get raster template corresponding to the sensor
HDR.Template <- read_ENVI_header(HDR.Temp.Path)
HDR_Template <- read_ENVI_header(HDR_Temp_Path)
# get info to update hdr file
# read hdr
HDR.input <- read_ENVI_header(get_HDR_name(Output.Path))
if (!is.null(HDR.Template$wavelength)) {
HDR.input$wavelength <- HDR.Template$wavelength
HDR_input <- read_ENVI_header(get_HDR_name(Output_Path))
if (!is.null(HDR_Template$wavelength)) {
HDR_input$wavelength <- HDR_Template$wavelength
}
if (!is.null(HDR.Template$`sensor type`)) {
HDR.input$`sensor type` <- HDR.Template$`sensor type`
if (!is.null(HDR_Template$`sensor type`)) {
HDR_input$`sensor type` <- HDR_Template$`sensor type`
}
if (!is.null(HDR.Template$`band names`)) {
HDR.input$`band names` <- HDR.Template$`band names`
if (!is.null(HDR_Template$`band names`)) {
HDR_input$`band names` <- HDR_Template$`band names`
}
if (!is.null(HDR.Template$`wavelength units`)) {
HDR.input$`wavelength units` <- HDR.Template$`wavelength units`
if (!is.null(HDR_Template$`wavelength units`)) {
HDR_input$`wavelength units` <- HDR_Template$`wavelength units`
}
# define visual stretch in the VIS domain
HDR.input$`default stretch` <- "0 1000 linear"
HDR_input$`default stretch` <- "0 1000 linear"
# write corresponding hdr file
write_ENVI_header(HDR.input, get_HDR_name(Output.Path))
} else if (!file.exists(HDR.Temp.Path)) {
write_ENVI_header(HDR_input, get_HDR_name(Output_Path))
} else if (!file.exists(HDR_Temp_Path)) {
message("Header template corresponding to the sensor expected to be found here")
print(HDR.Temp.Path)
print(HDR_Temp_Path)
message("please provide this header template in order to write info in HDR file")
print(get_HDR_name(Output.Path))
print(get_HDR_name(Output_Path))
message("or manually add wavelength location in HDR file, if relevant")
}
} else if (Sensor == "unknown") {
message("please make sure that the follozing header file contains information required")
print(get_HDR_name(Output.Path))
print(get_HDR_name(Output_Path))
message("or manually add wavelength location in HDR file, if relevant")
}
return(Output.Path)
return(Output_Path)
}
#' Checks if the data to be processed has the format type expected
#'
#' @param Raster.Path full path for the raster to be converted
#' @param Raster_Path full path for the raster to be converted
#' @param Mask is the raster a mask?
#'
#' @export
check_data <- function(Raster.Path, Mask = FALSE) {
HDR.Path <- get_HDR_name(Raster.Path)
check_data <- function(Raster_Path, Mask = FALSE) {
HDR_Path <- get_HDR_name(Raster_Path)
# check if the hdr file exists
if (file.exists(HDR.Path)) {
HDR <- read_ENVI_header(HDR.Path)
if (file.exists(HDR_Path)) {
HDR <- read_ENVI_header(HDR_Path)
if (Mask == FALSE & (!HDR$interleave == "bil") & (!HDR$interleave == "BIL")) {
message("")
message("*********************************************************")
......@@ -217,7 +216,7 @@ check_data <- function(Raster.Path, Mask = FALSE) {
message("")
message("*********************************************************")
message("The following HDR file was expected, but could not be found:")
print(HDR.Path)
print(HDR_Path)
message("The image format may not compatible with the processing chain")
message("Image format expected:")
message("ENVI hdr file with band interleaved by line (BIL) file format")
......
......@@ -34,8 +34,8 @@ apply_continuum_removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
Spectral.Data <- splitRows(Spectral.Data, nb.CR)
# perform multithread continuum removal
plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
Schedule.Per.Thread <- ceiling(nb.CR / nbCPU)
Spectral.Data.tmp <- future_lapply(Spectral.Data, FUN = ContinuumRemoval, Spectral.Bands = Spectral$Wavelength, future.scheduling = Schedule.Per.Thread)
Schedule_Per_Thread <- ceiling(nb.CR / nbCPU)
Spectral.Data.tmp <- future_lapply(Spectral.Data, FUN = ContinuumRemoval, Spectral.Bands = Spectral$Wavelength, future.scheduling = Schedule_Per_Thread)
plan(sequential)
Spectral.Data <- do.call("rbind", Spectral.Data.tmp)
rm(Spectral.Data.tmp)
......
......@@ -11,22 +11,24 @@
#' Performs radiometric filtering based on three criteria: NDVI, NIR reflectance, Blue reflectance
#'
#' @param Image.Path Path of the image to be processed
#' @param Mask.Path Path of the mask corresponding to the image
#' @param Output.Dir output directory
#' @param Image_Path Path of the image to be processed
#' @param Mask_Path Path of the mask corresponding to the image
#' @param Output_Dir output directory
#' @param TypePCA Type of PCA: "PCA" or "SPCA"
#' @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
#' @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
#' @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
#' @param Mask.Path
#' @param NDVI_Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI_Thresh)
#' @param Blue_Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue_Thresh --> filter clouds) refl expected between 0 and 10000
#' @param NIR_Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR_Thresh) refl expected between 0 and 10000
#' @param Blue central wavelength corresponding to the blue spectral band (in nanometers)
#' @param Red central wavelength corresponding to the red spectral band (in nanometers)
#' @param NIR central wavelength corresponding to the NIR spectral band (in nanometers)
#'
#' @return ImPathShade = updated shademask file
#' @export
perform_radiometric_filtering <- function(Image.Path, Mask.Path, Output.Dir, TypePCA = "SPCA", NDVI.Thresh = 0.5, Blue.Thresh = 500, NIR.Thresh = 1500, Blue = 480, Red = 700, NIR = 835) {
perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, TypePCA = "SPCA", NDVI_Thresh = 0.5, Blue_Thresh = 500, NIR_Thresh = 1500, Blue = 480, Red = 700, NIR = 835) {
# define full output directory
Output.Dir.Full <- define_output_directory(Output.Dir, Image.Path, TypePCA)
Output_Dir_Full <- define_output_directory(Output_Dir, Image_Path, TypePCA)
# define dimensions of the image
ImPathHDR <- get_HDR_name(Image.Path)
ImPathHDR <- get_HDR_name(Image_Path)
HDR <- read_ENVI_header(ImPathHDR)
Image.Format <- ENVI_type2bytes(HDR)
ipix <- as.double(HDR$lines)
......@@ -36,14 +38,14 @@ perform_radiometric_filtering <- function(Image.Path, Mask.Path, Output.Dir, Typ
ImSizeGb <- (lenTot * Image.Format$Bytes) / (1024^3)
# Create / Update shade mask if optical data
if (Mask.Path == FALSE | Mask.Path == "") {
if (Mask_Path == FALSE | Mask_Path == "") {
print("Create mask based on NDVI, NIR and Blue threshold")
} else {
print("Update mask based on NDVI, NIR and Blue threshold")
}
Shade.Update <- paste(Output.Dir.Full, "ShadeMask_Update", sep = "")
Mask.Path <- create_mask_from_threshold(Image.Path, Mask.Path, Shade.Update, NDVI.Thresh, Blue.Thresh, NIR.Thresh, Blue, Red, NIR)
return(Mask.Path)
Shade.Update <- paste(Output_Dir_Full, "ShadeMask_Update", sep = "")
Mask_Path <- create_mask_from_threshold(Image_Path, Mask_Path, Shade.Update, NDVI_Thresh, Blue_Thresh, NIR_Thresh, Blue, Red, NIR)
return(Mask_Path)
}
# create a mask based on NDVI, Green reflectance and NIR reflectance
......@@ -55,12 +57,12 @@ perform_radiometric_filtering <- function(Image.Path, Mask.Path, Output.Dir, Typ
# @param ImPath full path of a raster file
# @param ImPathShade full path of the raster mask corresponding to the raster file
# @param ImPathShade.Update wavelength (nm) of the spectral bands to be found
# @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
# @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
# @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
# @param NDVI_Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI_Thresh)
# @param Blue_Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue_Thresh --> filter clouds) refl expected between 0 and 10000
# @param NIR_Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR_Thresh) refl expected between 0 and 10000
#
# @return ImPathShade path for the updated shademask produced
create_mask_from_threshold <- function(ImPath, ImPathShade, ImPathShade.Update, NDVI.Thresh, Blue.Thresh, NIR.Thresh, Blue = 480, Red = 700, NIR = 835) {
create_mask_from_threshold <- function(ImPath, ImPathShade, ImPathShade.Update, NDVI_Thresh, Blue_Thresh, NIR_Thresh, Blue = 480, Red = 700, NIR = 835) {
# define wavelength corresponding to the spectral domains Blue, Red and NIR
Spectral.Bands <- c(Blue, Red, NIR)
ImPathHDR <- get_HDR_name(ImPath)
......@@ -78,15 +80,15 @@ create_mask_from_threshold <- function(ImPath, ImPathShade, ImPathShade.Update,
message("Could not find the spectral bands required to compute NDVI")
}
if (Image.Bands$Distance2WL[1] > 25) {
Image.Subset[, , 1] <- Blue.Thresh + 0 * Image.Subset[, , 1]
Image.Subset[, , 1] <- Blue_Thresh + 0 * Image.Subset[, , 1]
message("Could not find a spectral band in the blue domain: will not perform filtering based on blue reflectance")
}
if (Image.Bands$Distance2WL[3] > 50) {
Image.Subset[, , 3] <- NIR.Thresh + 0 * Image.Subset[, , 3]
Image.Subset[, , 3] <- NIR_Thresh + 0 * Image.Subset[, , 3]
message("Could not find a spectral band in the NIR domain: will not perform filtering based on NIR reflectance")
}
Mask <- matrix(0, nrow = Header$lines, ncol = Header$samples)
SelPixels <- which(NDVI > NDVI.Thresh & Image.Subset[, , 1] < Blue.Thresh & Image.Subset[, , 3] > NIR.Thresh)
SelPixels <- which(NDVI > NDVI_Thresh & Image.Subset[, , 1] < Blue_Thresh & Image.Subset[, , 3] > NIR_Thresh)
Mask[SelPixels] <- 1
# update initial shade mask
ImPathShade <- update_shademask(ImPathShade, Header, Mask, ImPathShade.Update)
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -13,60 +13,63 @@
#' maps spectral species based on PCA file computed previously
#'
#' @param Input.Image.File Path and name of the image to be processed
#' @param Output.Dir output directory
#' @param Input_Image_File Path and name of the image to be processed
#' @param Output_Dir output directory
#' @param TypePCA Type of PCA: "PCA" or "SPCA"
#' @param PCA.Files path for PCA file
#' @param PCA_Files path for PCA file
#' @param ImPathShade path for shade file (= mask)
#' @param Pix_Per_Partition number of pixels for each partition
#' @param nb_partitions number of partition
#' @param nbCPU numeric. Number of CPUs to use in parallel.
#' @param MaxRAM numeric. MaxRAM maximum size of chunk in GB to limit RAM allocation when reading image file.
#' @param nbclusters number of clusters defined in k-Means
#'
#' @importFrom utils read.table
#' @export
map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePCA = "SPCA", nbclusters = 50, nbCPU = 1, MaxRAM = FALSE) {
map_spectral_species <- function(Input_Image_File, Output_Dir, PCA_Files,ImPathShade,Pix_Per_Partition,nb_partitions,TypePCA = "SPCA", nbclusters = 50, nbCPU = 1, MaxRAM = FALSE) {
# for each image
Output.Dir2 <- define_output_directory(Output.Dir, Input.Image.File, TypePCA)
Output.Dir.SS <- define_output_subdir(Output.Dir, Input.Image.File, TypePCA, "SpectralSpecies")
Output.Dir.PCA <- define_output_subdir(Output.Dir, Input.Image.File, TypePCA, "PCA")
Spectral.Species.Path <- paste(Output.Dir.SS, "SpectralSpecies", sep = "")
Output_Dir_SS <- define_output_subdir(Output_Dir, Input_Image_File, TypePCA, "SpectralSpecies")
Output_Dir_PCA <- define_output_subdir(Output_Dir, Input_Image_File, TypePCA, "PCA")
Spectral_Species_Path <- paste(Output_Dir_SS, "SpectralSpecies", sep = "")
# if no prior diversity map has been produced --> need PCA file
if (!file.exists(PCA.Files)) {
if (!file.exists(PCA_Files)) {
message("")
message("*********************************************************")
message("WARNING: PCA file expected but not found")
print(PCA.Files)
print(PCA_Files)
message("process aborted")
message("*********************************************************")
message("")
stop()
} else {
WS_Save <- paste(Output.Dir2, "PCA_Info.RData", sep = "")
WS_Save <- paste(Output_Dir_PCA, "PCA_Info.RData", sep = "")
load(file = WS_Save)
## 1- Select components used to perform clustering ##
PC.Select.Path <- paste(Output.Dir.PCA, "Selected_Components.txt", sep = "")
if (file.exists(PC.Select.Path)) {
PC.Select <- read.table(PC.Select.Path)[[1]]
dataPCA <- PCA.model$dataPCA[, PC.Select]
if (length(PC.Select) == 1) {
PC_Select_Path <- paste(Output_Dir_PCA, "Selected_Components.txt", sep = "")
if (file.exists(PC_Select_Path)) {
PC_Select <- read.table(PC_Select_Path)[[1]]
dataPCA <- PCA.model$dataPCA[, PC_Select]
if (length(PC_Select) == 1) {
dataPCA <- matrix(dataPCA, ncol = 1)
}
message("Selected components:")
print(PC.Select)
print(PC_Select)
message("Please add carriage return after last selected component if not part of the list")
message("If these do not match with your selection, please correct file following file:")
print(PC.Select.Path)
print(PC_Select_Path)
} else {
print(paste("File named ", PC.Select.Path, "needs to be created first"))
print(paste("File named ", PC_Select_Path, "needs to be created first"))
print("Image processing aborted")
exit()
}
## 2- PERFORM KMEANS FOR EACH ITERATION & DEFINE SPECTRAL SPECIES ##
print("perform k-means clustering for each subset and define centroids")
# scaling factor subPCA between 0 and 1
Kmeans.info <- init_kmeans(dataPCA, Pix.Per.Iter, nb_partitions, nbclusters, nbCPU)
Kmeans_info <- init_kmeans(dataPCA, Pix_Per_Partition, nb_partitions, nbclusters, nbCPU)
## 3- ASSIGN SPECTRAL SPECIES TO EACH PIXEL ##
print("apply Kmeans to the whole image and determine spectral species")
apply_kmeans(PCA.Files, PC.Select, ImPathShade, Kmeans.info, Spectral.Species.Path, nbCPU, MaxRAM)
apply_kmeans(PCA_Files, PC_Select, ImPathShade, Kmeans_info, Spectral_Species_Path, nbCPU, MaxRAM)
}
return()
}
......@@ -74,7 +77,7 @@ map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePC
# computes k-means from nb_partitions subsets taken from dataPCA
#
# @param dataPCA initial dataset sampled from PCA image
# @param Pix.Per.Iter number of pixels per iteration
# @param Pix_Per_Partition number of pixels per iteration
# @param nb_partitions number of k-means then averaged
# @param nbCPU
# @param nbclusters number of clusters used in kmeans
......@@ -83,18 +86,18 @@ map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePC
#' @importFrom future plan multiprocess sequential
#' @importFrom future.apply future_lapply
#' @importFrom stats kmeans
init_kmeans <- function(dataPCA, Pix.Per.Iter, nb_partitions, nbclusters, nbCPU = 1) {
init_kmeans <- function(dataPCA, Pix_Per_Partition, nb_partitions, nbclusters, nbCPU = 1) {
m0 <- apply(dataPCA, 2, function(x) min(x))
M0 <- apply(dataPCA, 2, function(x) max(x))
d0 <- M0 - m0
dataPCA <- center_reduce(dataPCA, m0, d0)
# get the dimensions of the images, and the number of subimages to process
dataPCA <- split(as.data.frame(dataPCA), rep(1:nb_partitions, each = Pix.Per.Iter))
dataPCA <- split(as.data.frame(dataPCA), rep(1:nb_partitions, each = Pix_Per_Partition))
# multiprocess kmeans clustering
plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
Schedule.Per.Thread <- ceiling(length(dataPCA) / nbCPU)
Schedule_Per_Thread <- ceiling(length(dataPCA) / nbCPU)
res <- future_lapply(dataPCA, FUN = kmeans, centers = nbclusters, iter.max = 50, nstart = 10,
algorithm = c("Hartigan-Wong"), future.scheduling = Schedule.Per.Thread)
algorithm = c("Hartigan-Wong"), future.scheduling = Schedule_Per_Thread)
plan(sequential)
Centroids <- list(nb_partitions)
for (i in (1:nb_partitions)) {
......@@ -106,61 +109,61 @@ init_kmeans <- function(dataPCA, Pix.Per.Iter, nb_partitions, nbclusters, nbCPU
# Applies Kmeans clustering to PCA image
#
# @param PCA.Path path for the PCA image
# @param PC.Select PCs selected from PCA
# @param PCA_Path path for the PCA image
# @param PC_Select PCs selected from PCA
# @param ImPathShade shade mask
# @param Kmeans.info information about kmeans computed in previous step
# @param Kmeans_info information about kmeans computed in previous step
# @param nbCPU
# @param MaxRAM
# @param Spectral.Species.Path path for spectral species file to be written
# @param Spectral_Species_Path path for spectral species file to be written
#
# @return
apply_kmeans <- function(PCA.Path, PC.Select, ImPathShade, Kmeans.info, Spectral.Species.Path, nbCPU = 1, MaxRAM = FALSE) {
ImPathHDR <- get_HDR_name(PCA.Path)
HDR.PCA <- read_ENVI_header(ImPathHDR)
PCA.Format <- ENVI_type2bytes(HDR.PCA)
HDR.Shade <- get_HDR_name(ImPathShade)
HDR.Shade <- read_ENVI_header(HDR.Shade)
# prepare for sequential processing: SeqRead.Image informs about byte location to read
apply_kmeans <- function(PCA_Path, PC_Select, ImPathShade, Kmeans_info, Spectral_Species_Path, nbCPU = 1, MaxRAM = FALSE) {
ImPathHDR <- get_HDR_name(PCA_Path)
HDR_PCA <- read_ENVI_header(ImPathHDR)
PCA_Format <- ENVI_type2bytes(HDR_PCA)
HDR_Shade <- get_HDR_name(ImPathShade)
HDR_Shade <- read_ENVI_header(HDR_Shade)
# prepare for sequential processing: SeqRead_Image informs about byte location to read
if (MaxRAM == FALSE) {