Commit 0ae3f0a7 authored by jbferet's avatar jbferet
Browse files

- updated name of variables (. into _)

- changed check_data function to handle mask files
parent 7d29405f
......@@ -3,5 +3,4 @@ doc
.Rproj.user
.Rbuildignore
.Rhistory
RESULTS
02_REDACTION
TODO.md
......@@ -165,7 +165,6 @@ check_data <- function(Raster_Path, Mask = FALSE) {
if (file.exists(HDR_Path)) {
HDR <- read_ENVI_header(HDR_Path)
if (Mask == FALSE & (!HDR$interleave == "bil") & (!HDR$interleave == "BIL")) {
message("")
message("*********************************************************")
message("The image format may not compatible with the processing chain")
message("Image format expected:")
......@@ -178,43 +177,45 @@ check_data <- function(Raster_Path, Mask = FALSE) {
message("in order to convert your raster data")
message("or use appropriate software")
message("*********************************************************")
message("")
stop()
} else if (Mask == FALSE & ((HDR$interleave == "bil") | (HDR$interleave == "BIL"))) {
if (HDR$`wavelength units` == "Unknown") {
message("")
message("*********************************************************")
message("IF MULTI / HYPERSPECTRAL DATA: ")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("*********************************************************")
message("")
}
if ((!HDR$`wavelength units` == "Nanometers") & (!HDR$`wavelength units` == "nanometers")) {
message("")
message("*********************************************************")
message("IF MULTI / HYPERSPECTRAL DATA: ")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("*********************************************************")
message("")
}
if (HDR$`wavelength units` == "micrometers") {
message("")
message("*********************************************************")
message("Please convert wavelengths in nanometers in HDR file")
message("*********************************************************")
message("")
stop()
}
if ((HDR$`wavelength units` == "nanometers") | (HDR$`wavelength units` == "Nanometers")) {
message("")
message("*********************************************************")
message(" Raster format OK for processing ")
message(" Format of main raster OK for processing ")
message("*********************************************************")
message("")
}
} else if (Mask == TRUE & HDR$bands == 1 & ((HDR$interleave == "bil") | (HDR$interleave == "BIL") | (HDR$interleave == "bsq") | (HDR$interleave == "BSQ"))) {
message("*********************************************************")
message(" Format of mask raster OK for processing ")
message("*********************************************************")
} else if (Mask == TRUE & HDR$bands > 1) {
message("*********************************************************")
message(" Mask raster should contain only one layer ")
message(" Please produce a binary mask with a unique layer ")
message("*********************************************************")
stop()
}
}
} else {
message("")
message("*********************************************************")
message("The following HDR file was expected, but could not be found:")
print(HDR_Path)
......@@ -226,7 +227,6 @@ check_data <- function(Raster_Path, Mask = FALSE) {
print("raster2BIL")
message("in order to convert your raster data")
message("*********************************************************")
message("")
stop()
}
return(invisible())
......
......@@ -11,7 +11,7 @@
# prepares data to run multithreaded continuum removal
#
# @param Spectral.Data initial data matrix (nb samples x nb bands)
# @param Spectral_Data initial data matrix (nb samples x nb bands)
# @param nbCPU
# @param Spectral information about spectral bands
#
......@@ -19,33 +19,33 @@
#' @importFrom snow splitRows
#' @importFrom future plan multiprocess sequential
#' @importFrom future.apply future_lapply
apply_continuum_removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) {
if (!length(Spectral$WaterVapor) == 0) {
Spectral.Data <- Spectral.Data[, -Spectral$WaterVapor]
Spectral_Data <- Spectral_Data[, -Spectral$WaterVapor]
}
# split data to perform continuum removal on into reasonable amount of data
nb.Values <- dim(Spectral.Data)[1] * dim(Spectral.Data)[2]
nb.Values <- dim(Spectral_Data)[1] * dim(Spectral_Data)[2]
if (nb.Values > 0) {
# corresponds to ~ 40 Mb data, but CR tends to requires ~ 10 times memory
# avoids memory crash
Max.nb.Values <- 2e6
nb.CR <- ceiling(nb.Values / Max.nb.Values)
Spectral.Data <- splitRows(Spectral.Data, nb.CR)
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)
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)
Spectral_Data <- do.call("rbind", Spectral_Data_tmp)
rm(Spectral_Data_tmp)
} else {
# edit 31-jan-2018
# resize to delete first and last band as in continuum removal
Spectral.Data <- Spectral.Data[, -c(1, 2)]
Spectral_Data <- Spectral_Data[, -c(1, 2)]
}
gc()
return(Spectral.Data)
return(Spectral_Data)
}
# Computes continuum removal for matrix shaped data: more efficient than
......@@ -54,47 +54,47 @@ apply_continuum_removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
# given spectral band and R at the following bands
#
# @param Minit initial data matrix (nb samples x nb bands)
# @param Spectral.Bands information about spectral bands
# @param Spectral_Bands information about spectral bands
#
# @return samples from image and updated number of pixels to sampel if necessary
ContinuumRemoval <- function(Minit, Spectral.Bands) {
ContinuumRemoval <- function(Minit, Spectral_Bands) {
# Filter and prepare data prior to continuum removal
CR.data <- filter_prior_CR(Minit, Spectral.Bands)
Minit <- CR.data$Minit
nb.Bands <- dim(Minit)[2]
CR.data$Minit <- c()
Spectral.Bands <- CR.data$Spectral.Bands
nb.Samples <- CR.data$nb.Samples
CR_data <- filter_prior_CR(Minit, Spectral_Bands)
Minit <- CR_data$Minit
nbBands <- dim(Minit)[2]
CR_data$Minit <- c()
Spectral_Bands <- CR_data$Spectral_Bands
nbSamples <- CR_data$nbSamples
# if samples to be considered
if (nb.Samples > 0) {
if (nbSamples > 0) {
# initialization:
# spectral band corresponding to each element of the data matrix
Lambda <- repmat(matrix(Spectral.Bands, nrow = 1), nb.Samples, 1)
Lambda <- repmat(matrix(Spectral_Bands, nrow = 1), nbSamples, 1)
# prepare matrices used to check evolution of the CR process
# - elements still not processed through continuum removal: initialization to 1
Still.Need.CR <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
Still.Need.CR <- matrix(1, nrow = nbSamples, ncol = nbBands)
# - value of the convex hull: initially set to 0
Convex.Hull <- matrix(0, nrow = nb.Samples, ncol = nb.Bands)
Convex_Hull <- matrix(0, nrow = nbSamples, ncol = nbBands)
# - reflectance value for latest interception with convex hull:
# initialization to value of the first reflectance measurement
Intercept.Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nb.Bands)
Intercept_Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nbBands)
# - spectral band of latest interception
Latest.Intercept <- repmat(matrix(Spectral.Bands[1], ncol = 1), nb.Samples, nb.Bands)
Latest.Intercept <- repmat(matrix(Spectral_Bands[1], ncol = 1), nbSamples, nbBands)
# number of spectral bands found as intercept
nb.Intercept <- 0
# continues until arbitrary stopping criterion:
# stops when reach last spectral band (all values before last = 0)
while (max(Still.Need.CR[, 1:(nb.Bands - 2)]) == 1 & (nb.Intercept <= (nb.Bands / 2))) {
while (max(Still.Need.CR[, 1:(nbBands - 2)]) == 1 & (nb.Intercept <= (nbBands / 2))) {
nb.Intercept <- nb.Intercept + 1
# Mstep give the position of the values to be updated
Update.Data <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
Update.Data[, nb.Bands] <- 0
Update_Data <- matrix(1, nrow = nbSamples, ncol = nbBands)
Update_Data[, nbBands] <- 0
# initial step: first column set to 0; following steps: all bands below
# max of the convex hull are set to 0
Update.Data[which((Lambda - Latest.Intercept) < 0)] <- 0
Update_Data[which((Lambda - Latest.Intercept) < 0)] <- 0
# compute slope for each coordinate
Slope <- (Minit - Intercept.Hull) / (Lambda - Latest.Intercept) * Still.Need.CR
Slope <- (Minit - Intercept_Hull) / (Lambda - Latest.Intercept) * Still.Need.CR
# set current spectral band and previous bands to -9999
if (!length(which(Still.Need.CR == 0)) == 0) {
Slope[which(Still.Need.CR == 0)] <- -9999
......@@ -103,25 +103,25 @@ ContinuumRemoval <- function(Minit, Spectral.Bands) {
Slope[which(is.na(Slope))] <- -9999
}
# get max index for each row and convert into linear index
Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nb.Samples, nb.Bands)
Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nbSamples, nbBands)
# !!!! OPTIM: replace repmat with column operation
# update coordinates of latest intercept
Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nb.Bands)
Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nbBands)
# update latest intercept
Intercept.Hull <- repmat(matrix(Minit[Index.Max.Slope], ncol = 1), 1, nb.Bands)
Intercept_Hull <- repmat(matrix(Minit[Index.Max.Slope], ncol = 1), 1, nbBands)
# values corresponding to the domain between the two continuum maxima
Update.Data[which((Lambda - Latest.Intercept) >= 0 | Latest.Intercept == Spectral.Bands[nb.Bands])] <- 0
Update_Data[which((Lambda - Latest.Intercept) >= 0 | Latest.Intercept == Spectral_Bands[nbBands])] <- 0
# values to eliminate for the next analysis: all spectral bands before latest intercept
Still.Need.CR[which((Lambda - Latest.Intercept) < 0)] <- 0
# the max slope is known, as well as the coordinates of the beginning and ending
# a matrix now has to be built
Convex.Hull <- Convex.Hull + Update.Data * (Intercept.Hull + sweep((Lambda - Latest.Intercept), 1, Slope[Index.Max.Slope], "*"))
Convex_Hull <- Convex_Hull + Update_Data * (Intercept_Hull + sweep((Lambda - Latest.Intercept), 1, Slope[Index.Max.Slope], "*"))
}
CR_Results0 <- Minit[, 2:(nb.Bands - 2)] / Convex.Hull[, 2:(nb.Bands - 2)]
CR_Results <- matrix(0, ncol = (nb.Bands - 3), nrow = nb.Samples)
CR_Results[CR.data$Samples.To.Keep, ] <- CR_Results0
CR_Results0 <- Minit[, 2:(nbBands - 2)] / Convex_Hull[, 2:(nbBands - 2)]
CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
CR_Results[CR_data$SamplesToKeep, ] <- CR_Results0
} else {
CR_Results <- matrix(0, ncol = (nb.Bands - 3), nrow = nb.Samples)
CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
}
list <- ls()
rm(list = list[-which(list == "CR_Results")])
......@@ -136,15 +136,15 @@ ContinuumRemoval <- function(Minit, Spectral.Bands) {
# - possibly remaining negative values are set to 0
# - constant spectra are eliminated
#
# @param Spectral.Bands
# @param Spectral_Bands
# @param Minit initial data matrix, n rows = n samples, p cols = p spectral bands
#
# @return updated Minit
#' @importFrom matrixStats rowSds
filter_prior_CR <- function(Minit, Spectral.Bands) {
filter_prior_CR <- function(Minit, Spectral_Bands) {
# number of samples to be processed
nb.Samples <- nrow(Minit)
nbSamples <- nrow(Minit)
# make sure there is no negative values
Minit <- Minit + 100.0
Minit[which(Minit < 0)] <- 0
......@@ -164,10 +164,10 @@ filter_prior_CR <- function(Minit, Spectral.Bands) {
Minit <- matrix(Minit, nrow = 1)
}
# add negative values to the last column and update spectral bands
Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nb.Samples))
nb.Bands <- ncol(Minit)
Spectral.Bands[nb.Bands] <- Spectral.Bands[nb.Bands - 1] + 10
my_list <- list("Minit" = Minit, "Spectral.Bands" = Spectral.Bands, "nb.Samples" = nb.Samples, "Samples.To.Keep" = keep)
Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nbSamples))
nbBands <- ncol(Minit)
Spectral_Bands[nbBands] <- Spectral_Bands[nbBands - 1] + 10
my_list <- list("Minit" = Minit, "Spectral_Bands" = Spectral_Bands, "nbSamples" = nbSamples, "SamplesToKeep" = keep)
return(my_list)
}
......
......@@ -22,7 +22,7 @@
#' @param Red numeric. central wavelength corresponding to the red spectral band (in nanometers)
#' @param NIR numeric. central wavelength corresponding to the NIR spectral band (in nanometers)
#'
#' @return ImPathShade = updated shademask file
#' @return MaskPath = updated mask 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) {
# check if format of raster data is as expected
......@@ -35,12 +35,12 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ
# define dimensions of the image
ImPathHDR <- get_HDR_name(Image_Path)
HDR <- read_ENVI_header(ImPathHDR)
Image.Format <- ENVI_type2bytes(HDR)
Image_Format <- ENVI_type2bytes(HDR)
ipix <- as.double(HDR$lines)
jpix <- as.double(HDR$samples)
Nb.Pixels <- ipix * jpix
lenTot <- Nb.Pixels * as.double(HDR$bands)
ImSizeGb <- (lenTot * Image.Format$Bytes) / (1024^3)
nbPixels <- ipix * jpix
lenTot <- nbPixels * as.double(HDR$bands)
ImSizeGb <- (lenTot * Image_Format$Bytes) / (1024^3)
# Create / Update shade mask if optical data
if (Mask_Path == FALSE | Mask_Path == "") {
......@@ -60,42 +60,42 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ
# ! only valid if Optical data!!
#
# @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 MaskPath full path of the raster mask corresponding to the raster file
# @param MaskPath.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
#
# @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) {
# @return MaskPath path for the updated shademask produced
create_mask_from_threshold <- function(ImPath, MaskPath, MaskPath.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)
Spectral_Bands <- c(Blue, Red, NIR)
ImPathHDR <- get_HDR_name(ImPath)
Header <- read_ENVI_header(ImPathHDR)
# get image bands correponding to spectral bands of interest
Image.Bands <- get_image_bands(Spectral.Bands, Header$wavelength)
Image_Bands <- get_image_bands(Spectral_Bands, Header$wavelength)
# read band data from image
Image.Subset <- read_image_bands(ImPath, Header, Image.Bands$ImBand)
Image_Subset <- read_image_bands(ImPath, Header, Image_Bands$ImBand)
# create mask
# check if spectral bands required for NDVI exist
if (Image.Bands$Distance2WL[2] < 25 & Image.Bands$Distance2WL[3] < 25) {
NDVI <- ((Image.Subset[, , 3]) - (Image.Subset[, , 2])) / ((Image.Subset[, , 3]) + (Image.Subset[, , 2]))
if (Image_Bands$Distance2WL[2] < 25 & Image_Bands$Distance2WL[3] < 25) {
NDVI <- ((Image_Subset[, , 3]) - (Image_Subset[, , 2])) / ((Image_Subset[, , 3]) + (Image_Subset[, , 2]))
} else {
NDVI <- matrix(1, nrow = Header$lines, ncol = Header$samples)
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]
if (Image_Bands$Distance2WL[1] > 25) {
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]
if (Image_Bands$Distance2WL[3] > 50) {
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)
return(ImPathShade)
MaskPath <- update_shademask(MaskPath, Header, Mask, MaskPath.Update)
return(MaskPath)
}
......@@ -94,42 +94,42 @@ check_invariant_bands <- function(DataMatrix, Spectral) {
# define output directory and create it if necessary
#
# @param Output.Dir output directory
# @param Output_Dir output directory
# @param ImPath image path
# @param TypePCA Type of PCA (PCA, SPCA, NLPCA...)
#
# @return path of the output directory
define_output_directory <- function(Output.Dir, ImPath, TypePCA) {
Image.Name <- strsplit(basename(ImPath), "\\.")[[1]][1]
Output.Dir <- paste(Output.Dir, "/", Image.Name, "/", TypePCA, "/", sep = "")
dir.create(Output.Dir, showWarnings = FALSE, recursive = TRUE)
return(Output.Dir)
define_output_directory <- function(Output_Dir, ImPath, TypePCA) {
Image_Name <- strsplit(basename(ImPath), "\\.")[[1]][1]
Output_Dir <- paste(Output_Dir, "/", Image_Name, "/", TypePCA, "/", sep = "")
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
return(Output_Dir)
}
# define output directory and subdirectory and create it if necessary
#
# @param Output.Dir output directory
# @param Output_Dir output directory
# @param ImPath image path
# @param TypePCA Type of PCA (PCA, SPCA, NLPCA...)
# @param Sub subdirectory
#
# @return path of the output directory
define_output_subdir <- function(Output.Dir, ImPath, TypePCA, Sub) {
Image.Name <- strsplit(basename(ImPath), "\\.")[[1]][1]
Output.Dir <- paste(Output.Dir, "/", Image.Name, "/", TypePCA, "/", Sub, "/", sep = "")
dir.create(Output.Dir, showWarnings = FALSE, recursive = TRUE)
return(Output.Dir)
define_output_subdir <- function(Output_Dir, ImPath, TypePCA, Sub) {
Image_Name <- strsplit(basename(ImPath), "\\.")[[1]][1]
Output_Dir <- paste(Output_Dir, "/", Image_Name, "/", TypePCA, "/", Sub, "/", sep = "")
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
return(Output_Dir)
}
# get information corresponding to data type defined in ENVI
#
# @param Header Path of the hdr file
# @param HDR header file
#
# @return description of data format corresponding to ENVI type
ENVI_type2bytes <- function(Header) {
ENVI_type2bytes <- function(HDR) {
# http://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html
DataTypeImage <- Header$`data type`
DataTypeImage <- HDR$`data type`
if (DataTypeImage == 1) {
# 1 = Byte: 8-bit unsigned integer
nbBytes <- 1
......@@ -176,9 +176,9 @@ ENVI_type2bytes <- function(Header) {
Type <- "INT"
is_Signed <- FALSE
}
if (Header$`byte order` == 0) {
if (HDR$`byte order` == 0) {
ByteOrder <- "little"
} else if (Header$`byte order` == 1) {
} else if (HDR$`byte order` == 1) {
ByteOrder <- "big"
}
my_list <- list("Bytes" = nbBytes, "Type" = Type, "Signed" = is_Signed, "ByteOrder" = ByteOrder)
......@@ -369,25 +369,25 @@ extract_samples_from_image <- function(ImPath, coordPix, MaxRAM = FALSE, Already
return(Sample_Sel)
}
# Perform random sampling on an image including a shade mask
# Perform random sampling on an image including a mask
#
# @param ImPath path for image
# @param HDR path for hdr file
# @param ImPathShade path for shade mask
# @param MaskPath path for mask
# @param nb_partitions number of k-means then averaged
# @param Pix_Per_Partition number of pixels per iteration
#
# @return samples from image and updated number of pixels to sampel if necessary
#' @importFrom matlab ones
get_random_subset_from_image <- function(ImPath, HDR, ImPathShade, nb_partitions, Pix_Per_Partition) {
get_random_subset_from_image <- function(ImPath, HDR, MaskPath, nb_partitions, Pix_Per_Partition) {
nbPix2Sample <- nb_partitions * Pix_Per_Partition
# get total number of pixels
nbpix <- as.double(HDR$lines) * as.double(HDR$samples)
# 1- Exclude masked pixels from random subset
# Read Mask
if ((!ImPathShade == "") & (!ImPathShade == FALSE)) {
if ((!MaskPath == "") & (!MaskPath == FALSE)) {
fid <- file(
description = ImPathShade, open = "rb", blocking = TRUE,
description = MaskPath, open = "rb", blocking = TRUE,
encoding = getOption("encoding"), raw = FALSE
)
ShadeMask <- readBin(fid, integer(), n = nbpix, size = 1)
......@@ -468,17 +468,17 @@ get_HDR_name <- function(ImPath) {
# gets rank of spectral bands in an image
#
# @param Spectral.Bands wavelength (nm) of the spectral bands to be found
# @param Spectral_Bands wavelength (nm) of the spectral bands to be found
# @param wavelength wavelength (nm) of all wavelengths in the image
#
# @return rank of all spectral bands of interest in the image and corresponding wavelength
get_image_bands <- function(Spectral.Bands, wavelength) {
get_image_bands <- function(Spectral_Bands, wavelength) {
ImBand <- c()
Distance2WL <- c()
for (band in Spectral.Bands) {
Closest.Band <- order(abs(wavelength - band))[1]
ImBand <- c(ImBand, Closest.Band)
Distance2WL <- c(Distance2WL, abs(wavelength[Closest.Band] - band))
for (band in Spectral_Bands) {
Closest_Band <- order(abs(wavelength - band))[1]
ImBand <- c(ImBand, Closest_Band)
Distance2WL <- c(Distance2WL, abs(wavelength[Closest_Band] - band))
}
my_list <- list("ImBand" = ImBand, "Distance2WL" = Distance2WL)
return(my_list)
......@@ -487,10 +487,10 @@ get_image_bands <- function(Spectral.Bands, wavelength) {
# convert image coordinates from index to X-Y
#
# @param Raster image raster object
# @param Image.Index coordinates corresponding to the raster
ind2sub <- function(Raster, Image.Index) {
c <- ((Image.Index - 1) %% Raster@ncols) + 1
r <- floor((Image.Index - 1) / Raster@ncols) + 1
# @param Image_Index coordinates corresponding to the raster
ind2sub <- function(Raster, Image_Index) {
c <- ((Image_Index - 1) %% Raster@ncols) + 1
r <- floor((Image_Index - 1) / Raster@ncols) + 1
my_list <- list("Column" = c, "Row" = r)
return(my_list)
}
......@@ -499,10 +499,10 @@ ind2sub <- function(Raster, Image.Index) {
# image coordinates are given as index = (ID.col-1) * total.lines + ID.row
#
# @param Raster image raster object
# @param Image.Index coordinates corresponding to the raster
ind2sub2 <- function(Raster, Image.Index) {
r <- ((Image.Index - 1) %% Raster@nrows) + 1
c <- floor((Image.Index - 1) / Raster@nrows) + 1
# @param Image_Index coordinates corresponding to the raster
ind2sub2 <- function(Raster, Image_Index) {
r <- ((Image_Index - 1) %% Raster@nrows) + 1
c <- floor((Image_Index - 1) / Raster@nrows) + 1
my_list <- list("Column" = c, "Row" = r)
return(my_list)
}
......@@ -573,25 +573,25 @@ read_bin_subset <- function(Byte_Start, nbLines, lenBin, ImPath, ImBand, jpix, n
# Reads ENVI hdr file
#
# @param headerFpath Path of the hdr file
# @param HDRpath Path of the hdr file
#
# @return list of the content of the hdr file
read_ENVI_header <- function(headerFpath) {
read_ENVI_header <- function(HDRpath) {
# header <- paste(header, collapse = "\n")
if (!grepl(".hdr$", headerFpath)) {
if (!grepl(".hdr$", HDRpath)) {
stop("File extension should be .hdr")
}
header <- readLines(headerFpath)
HDR <- readLines(HDRpath)
## check ENVI at beginning of file
if (!grepl("ENVI", header[1])) {
if (!grepl("ENVI", HDR[1])) {
stop("Not an ENVI header (ENVI keyword missing)")
} else {
header <- header [-1]
HDR <- HDR [-1]
}
## remove curly braces and put multi-line key-value-pairs into one line
header <- gsub("\\{([^}]*)\\}", "\\1", header)
l <- grep("\\{", header)
r <- grep("\\}", header)
HDR <- gsub("\\{([^}]*)\\}", "\\1", HDR)
l <- grep("\\{", HDR)
r <- grep("\\}", HDR)
if (length(l) != length(r)) {
stop("Error matching curly braces in header (differing numbers).")
......@@ -601,47 +601,47 @@ read_ENVI_header <- function(headerFpath) {
stop("Mismatch of curly braces in header.")
}
header[l] <- sub("\\{", "", header[l])
header[r] <- sub("\\}", "", header[r])
HDR[l] <- sub("\\{", "", HDR[l])
HDR[r] <- sub("\\}", "", HDR[r])
for (i in rev(seq_along(l))) {
header <- c(
header [seq_len(l [i] - 1)],
paste(header [l [i]:r [i]], collapse = "\n"),
header [-seq_len(r [i])]
HDR <- c(
HDR [seq_len(l [i] - 1)],
paste(HDR [l [i]:r [i]], collapse = "\n"),
HDR [-seq_len(r [i])]
)
}
## split key = value constructs into list with keys as names
header <- sapply(header, split_line, "=", USE.NAMES = FALSE)
names(header) <- tolower(names(header))
HDR <- sapply(HDR, split_line, "=", USE.NAMES = FALSE)
names(HDR) <- tolower(names(HDR))