diff --git a/DESCRIPTION b/DESCRIPTION index 5eb0d265..58b214d8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: biodivMapR Title: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images -Version: 1.2.1 +Version: 1.3.0 Authors@R: c(person(given = "Jean-Baptiste", family = "Feret", email = "jb.feret@teledetection.fr", @@ -19,6 +19,7 @@ License: GPL-3 + file LICENSE LazyData: true Imports: dissUtils, + data.table, ecodist, fields, future, @@ -27,11 +28,13 @@ Imports: matlab, matrixStats, methods, + mmand, raster, rgdal, R.utils, snow, sp, + stars, stringr, tools, vegan, diff --git a/NAMESPACE b/NAMESPACE index 9ae46ddf..7ad1ee1a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(Write_Big_Image) +export(Write_Image_NativeRes) export(check_data) export(diversity_from_plots) export(list_shp) @@ -13,7 +14,11 @@ export(raster2BIL) export(select_PCA_components) export(split_line) import(raster) +import(stars) import(tools) +importFrom(data.table,data.table) +importFrom(data.table,rbindlist) +importFrom(data.table,setorder) importFrom(dissUtils,diss) importFrom(ecodist,nmds) importFrom(fields,rdist) @@ -24,8 +29,10 @@ importFrom(future.apply,future_lapply) importFrom(labdsv,pco) importFrom(matlab,ones) importFrom(matlab,padarray) +importFrom(matrixStats,rowAnys) importFrom(matrixStats,rowSds) importFrom(methods,is) +importFrom(mmand,erode) importFrom(raster,projection) importFrom(raster,raster) importFrom(raster,writeRaster) diff --git a/NEWS.md b/NEWS.md index 34dc6b6f..0de2e987 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# biodivMapR v1.3.0 + +## Changes +- integrated stars package in order to read any file format, including TIFF format +- developed a generic function to write rasters +- prepared for MNF +- changed default red band for the computation of NDVI: closest band to 690 nm is now selected instead of closest band to 700 nm + # biodivMapR v1.2.1 ## Changes diff --git a/R/Lib_CheckConvertData.R b/R/Lib_CheckConvertData.R index 33acf480..1856ebff 100644 --- a/R/Lib_CheckConvertData.R +++ b/R/Lib_CheckConvertData.R @@ -4,6 +4,7 @@ # ============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2018/07 Jean-Baptiste FERET # ============================================================================== # This library checks if a raster in format expected for diversity mapping @@ -168,25 +169,13 @@ check_data <- function(Raster_Path, Mask = FALSE) { # check if the hdr file exists if (file.exists(HDR_Path)) { HDR <- read_ENVI_header(HDR_Path) - if (Mask == FALSE & (!HDR$interleave == "bil") & (!HDR$interleave == "BIL")) { - message("*********************************************************") - 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") - message("") - message("Current Image format") - print(HDR$interleave) - message("Please run the function named ") - print("raster2BIL") - message("in order to convert your raster data") - message("or use appropriate software") - message("*********************************************************") - stop() - } else if (Mask == FALSE & ((HDR$interleave == "bil") | (HDR$interleave == "BIL"))) { - if (is.null(HDR$`wavelength units`)) { + if (Mask == FALSE ) { + if (is.null(HDR$wavelength)) { message("*********************************************************") - message("Image to process is not multispectral/hyperspectral image ") - message("Format is OK, but make sure Continuum_Removal is set to FALSE") + message("No wavelength is associated to the image in the .hdr file") + message("Please add wavelengths to the .hdr file") + message("if you are processing multi or hyperspectral optical data") + message(" Otherwise, make sure Continuum_Removal is set to FALSE ") message("*********************************************************") } else { if (HDR$`wavelength units` == "Unknown") { @@ -196,33 +185,21 @@ check_data <- function(Raster_Path, Mask = FALSE) { message("if not, stop processing and convert wavelengths in nanometers in HDR file") message("*********************************************************") } - if ((!HDR$`wavelength units` == "Nanometers") & (!HDR$`wavelength units` == "nanometers")) { + if ((!HDR$`wavelength units` == "Nanometers") & (!HDR$`wavelength units` == "nanometers") & + (!HDR$`wavelength units` == "Micrometers") & (!HDR$`wavelength units` == "micrometers")) { 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("*********************************************************") - } - if (HDR$`wavelength units` == "micrometers") { - message("*********************************************************") - message("Please convert wavelengths in nanometers in HDR file") - message("*********************************************************") - stop() - } - if ((HDR$`wavelength units` == "nanometers") | (HDR$`wavelength units` == "Nanometers")) { - message("*********************************************************") - message(" Format of main raster OK for processing ") + message("Please make sure the wavelengths are in nanometers or micrometers") + message("if not, stop processing and convert wavelengths in nanometers") + message("or micrometers in HDR file") 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(" and value = 1 for pixels to keep ") message("*********************************************************") stop() } @@ -230,13 +207,6 @@ check_data <- function(Raster_Path, Mask = FALSE) { message("*********************************************************") message("The following HDR file was expected, but could not be found:") 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") - message("") - message("Please run the function named ") - print("raster2BIL") - message("in order to convert your raster data") message("*********************************************************") stop() } diff --git a/R/Lib_ContinuumRemoval.R b/R/Lib_ContinuumRemoval.R index b5422ef5..6ddc86e2 100644 --- a/R/Lib_ContinuumRemoval.R +++ b/R/Lib_ContinuumRemoval.R @@ -4,6 +4,7 @@ # =============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2018/07 Jean-Baptiste FERET # =============================================================================== # This Library is dedicated to the computation of the continuum removal @@ -53,10 +54,10 @@ apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) { # the convex hull is based on the computation of the derivative between R at a # 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 Minit numeric. initial data matrix (nb samples x nb bands) +#' @param Spectral_Bands numeric. central wavelength for the spectral bands # -# @return samples from image and updated number of pixels to sampel if necessary +#' @return samples from image and updated number of pixels to sampel if necessary ContinuumRemoval <- function(Minit, Spectral_Bands) { # Filter and prepare data prior to continuum removal @@ -154,10 +155,10 @@ ContinuumRemoval <- function(Minit, Spectral_Bands) { # - possibly remaining negative values are set to 0 # - constant spectra are eliminated # -# @param Spectral_Bands -# @param Minit initial data matrix, n rows = n samples, p cols = p spectral bands +#' @param Minit numeric. initial data matrix (nb samples x nb bands) +#' @param Spectral_Bands numeric. central wavelength for the spectral bands # -# @return updated Minit +#' @return list. updated Minit #' @importFrom matrixStats rowSds filter_prior_CR <- function(Minit, Spectral_Bands) { diff --git a/R/Lib_FilterData.R b/R/Lib_FilterData.R index 3c8907a6..3b891b07 100644 --- a/R/Lib_FilterData.R +++ b/R/Lib_FilterData.R @@ -4,6 +4,7 @@ # ============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2018/07 Jean-Baptiste FERET # ============================================================================== # This library contains functions to filter raster based on radiometric criteria @@ -24,7 +25,9 @@ #' #' @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) { +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 check_data(Image_Path) if (!Mask_Path==FALSE){ @@ -32,24 +35,16 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ } # define full output directory Output_Dir_Full <- define_output_directory(Output_Dir, Image_Path, TypePCA) - # define dimensions of the image - ImPathHDR <- get_HDR_name(Image_Path) - HDR <- read_ENVI_header(ImPathHDR) - Image_Format <- ENVI_type2bytes(HDR) - ipix <- as.double(HDR$lines) - jpix <- as.double(HDR$samples) - 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 == "") { 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) + Shade_Update <- paste(Output_Dir_Full, "ShadeMask_Update", sep = "") + Mask_Path <- create_mask_from_threshold(ImPath = Image_Path, MaskPath = Mask_Path, MaskPath_Update = Shade_Update, + NDVI_Thresh = NDVI_Thresh, Blue_Thresh = Blue_Thresh, NIR_Thresh = NIR_Thresh, + Blue = Blue, Red = Red, NIR = NIR) return(Mask_Path) } @@ -59,44 +54,57 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ # NIR (min) threshold eliminates shadows # ! only valid if Optical data!! # -# @param ImPath full path of a raster file -# @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 +#' @param ImPath character. Full path of a raster file +#' @param MaskPath character. Full path of the mask to be used with the raster file +#' @param MaskPath_Update character. Full path of the updated mask to be used with the raster file +#' @param NDVI_Thresh numeric. NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI_Thresh) +#' @param Blue_Thresh numeric. 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 numeric. NIR threshold applied to produce a mask (select pixels with NIR refl < NIR_Thresh) refl expected between 0 and 10000 +#' @param Blue numeric. spectral band corresponding to the blue channel (in nanometers) +#' @param Red numeric. spectral band corresponding to the red channel (in nanometers) +#' @param NIR numeric. spectral band corresponding to the NIR channel (in nanometers) # # @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) { +create_mask_from_threshold <- function(ImPath, MaskPath, MaskPath_Update, NDVI_Thresh, Blue_Thresh, NIR_Thresh, + Blue = 480, Red = 690, NIR = 835) { # define wavelength corresponding to the spectral domains Blue, Red and NIR Spectral_Bands <- c(Blue, Red, NIR) ImPathHDR <- get_HDR_name(ImPath) - Header <- read_ENVI_header(ImPathHDR) + HDR <- read_ENVI_header(ImPathHDR) + # distance between expected bands defining red, blue and NIR info and available band from sensor + Dist2Band <- 25 + # in case micrometers + if (max(HDR$wavelength)<100 | HDR$`wavelength units` == "micrometers"){ + Spectral_Bands <- 0.001*Spectral_Bands + Dist2Band <- 0.001*Dist2Band + } # 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, HDR$wavelength) # read band data from image - Image_Subset <- read_image_bands(ImPath, Header, Image_Bands$ImBand) + Image_Subset <- read_image_bands(ImPath = ImPath, HDR = HDR, + ImBand = Image_Bands$ImBand) # create mask # check if spectral bands required for NDVI exist - if (Image_Bands$Distance2WL[2] < 25 & Image_Bands$Distance2WL[3] < 25) { + + if (Image_Bands$Distance2WL[2] < Dist2Band & Image_Bands$Distance2WL[3] < Dist2Band) { NDVI <- ((Image_Subset[, , 3]) - (Image_Subset[, , 2])) / ((Image_Subset[, , 3]) + (Image_Subset[, , 2])) } else { - NDVI <- matrix(1, nrow = Header$lines, ncol = Header$samples) + NDVI <- matrix(1, nrow = HDR$lines, ncol = HDR$samples) message("Could not find the spectral bands required to compute NDVI") } - if (Image_Bands$Distance2WL[1] > 25) { + if (Image_Bands$Distance2WL[1] > Dist2Band) { 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) { + if (Image_Bands$Distance2WL[3] > 2*Dist2Band) { 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) + Mask <- matrix(0, nrow = HDR$lines, ncol = HDR$samples) SelPixels <- which(NDVI > NDVI_Thresh & Image_Subset[, , 1] < Blue_Thresh & Image_Subset[, , 3] > NIR_Thresh) Mask[SelPixels] <- 1 # update initial shade mask - MaskPath <- update_shademask(MaskPath, Header, Mask, MaskPath.Update) + MaskPath <- update_shademask(MaskPath, HDR, Mask, MaskPath_Update) list2Remove <- ls() rm(list=list2Remove[-which(list2Remove=='MaskPath')]) gc() diff --git a/R/Lib_ImageProcess.R b/R/Lib_ImageProcess.R index c006e5d0..d3bc74da 100644 --- a/R/Lib_ImageProcess.R +++ b/R/Lib_ImageProcess.R @@ -4,6 +4,7 @@ # ============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2018/07 Jean-Baptiste FERET # ============================================================================== # This Library contains functions to manipulate & process raster images @@ -60,9 +61,7 @@ change_resolution_HDR <- function(HDR, window_size) { return(HDR) } -# clean data matrix from inf values and identifies if some values -# (columns) are constant. variables showing constant value -# need to be eliminated before PCA +# remove constant bands # # @param DataMatrix each variable is a column # @param Spectral summary of spectral information: which spectral bands selected from initial data @@ -70,7 +69,7 @@ change_resolution_HDR <- function(HDR, window_size) { # @return updated DataMatrix and Spectral # ' @importFrom stats sd -check_invariant_bands <- function(DataMatrix, Spectral) { +rm_invariant_bands <- function(DataMatrix, Spectral) { # samples with inf value are eliminated for (i in 1:ncol(DataMatrix)) { elim <- which(DataMatrix[, i] == Inf) @@ -78,7 +77,7 @@ check_invariant_bands <- function(DataMatrix, Spectral) { DataMatrix <- DataMatrix[-elim, ] } } - # bands showing null std are eliminated + # bands showing null std are removed stdsub <- apply(DataMatrix, 2, sd) BandsNoVar <- which(stdsub == 0) # BandsNoVar = which(stdsub<=0.002) @@ -203,6 +202,12 @@ exclude_spectral_domains <- function(ImPath, Excluded_WL = FALSE) { Excluded_WL <- rbind(Excluded_WL, c(1780, 2040)) } } + message("Exclude spectral domains corresponding to water vapor absorption") + message("The following spectral domains (min and max spectral band in nanometers) will be discarded") + print(Excluded_WL) + message("Please define the input variable Excluded_WL when calling perform_PCA") + message("if you want to modify these spectral domains") + # in case a unique specrtal domain is provided as excluded domain if (length(Excluded_WL)==2){ Excluded_WL = matrix(Excluded_WL,ncol = 2) @@ -210,6 +215,9 @@ exclude_spectral_domains <- function(ImPath, Excluded_WL = FALSE) { # get image header data ImPathHDR <- get_HDR_name(ImPath) HDR <- read_ENVI_header(ImPathHDR) + if ((HDR$`wavelength units`=='Micrometers') | (HDR$`wavelength units`=='micrometers')){ + Excluded_WL <- 0.001*Excluded_WL + } nbchannels0 <- HDR$bands idOkBand <- seq(1, nbchannels0) if ("wavelength" %in% names(HDR)) { @@ -301,15 +309,15 @@ extract_samples_from_image <- function(ImPath, coordPix, MaxRAM = FALSE, Already if (dim(coordPix)[2] == 2) { if (dim(coordPix)[1] >= 2) { coordPix_tmp <- list() - coordPix_tmp$Row <- coordPix[, 1] - coordPix_tmp$Column <- coordPix[, 2] + coordPix_tmp$row <- coordPix[, 1] + coordPix_tmp$col <- coordPix[, 2] } else if (dim(coordPix)[1] == 1) { coordPix_tmp <- list() - coordPix_tmp$Row <- coordPix[1] - coordPix_tmp$Column <- coordPix[2] + coordPix_tmp$row <- coordPix[1] + coordPix_tmp$col <- coordPix[2] } } - } else if (typeof(coordPix) == "list" & length(grep("Row", names(coordPix))) > 0 & length(grep("Column", names(coordPix))) > 0) { + } else if (typeof(coordPix) == "list" & length(grep("row", names(coordPix))) > 0 & length(grep("col", names(coordPix))) > 0) { coordPix_tmp <- coordPix } # initial index value of the pixels requested in the image, following original ranking @@ -317,8 +325,8 @@ extract_samples_from_image <- function(ImPath, coordPix, MaxRAM = FALSE, Already # rank of the initial list of pixels Rank_Index_Init <- sort(Index_Init, index = TRUE) # convert - if (typeof(coordPix) == "list" & length(grep("Row", names(coordPix))) > 0 & length(grep("Column", names(coordPix))) > 0) { - coordPix <- cbind(coordPix$Row, coordPix$Column) + if (typeof(coordPix) == "list" & length(grep("row", names(coordPix))) > 0 & length(grep("col", names(coordPix))) > 0) { + coordPix <- cbind(coordPix$row, coordPix$col) } # divide data access based on the size of the image (to avoid reading 30 Gb file at once...) Image_Format <- ENVI_type2bytes(HDR) @@ -354,16 +362,16 @@ extract_samples_from_image <- function(ImPath, coordPix, MaxRAM = FALSE, Already # re-sort samples Sample_Sort <- list() - Sample_Sort$Row <- coordPix_List[, 1] - Sample_Sort$Column <- coordPix_List[, 2] + Sample_Sort$row <- coordPix_List[, 1] + Sample_Sort$col <- coordPix_List[, 2] Coord_Pixels <- sub2ind(HDR, Sample_Sort) # rank of the pixels extracted from the image Rank_Pixels <- sort(Coord_Pixels, index = TRUE) # pixel re-ordering needs to be performed in order to get back to Rank_Index_Init # first re-order pixels in order to follow increasing index value # Sample_Sort$index = Coord_Pixels[Rank_Pixels$ix] - # Sample_Sort$Row = Sample_Sort$Row[Rank_Pixels$ix] - # Sample_Sort$Column = Sample_Sort$Column[Rank_Pixels$ix] + # Sample_Sort$row = Sample_Sort$row[Rank_Pixels$ix] + # Sample_Sort$col = Sample_Sort$col[Rank_Pixels$ix] # if bug, check this line # Sample_Sel = Sample_Sel[Rank_Pixels$ix] Sample_Sel <- Sample_Sel[Rank_Pixels$ix, ] @@ -377,36 +385,106 @@ extract_samples_from_image <- function(ImPath, coordPix, MaxRAM = FALSE, Already return(Sample_Sel) } -# Perform random sampling on an image including a mask +# Extract bands of sparse pixels in big image data cube. # -# @param ImPath path for image -# @param HDR path for hdr file -# @param MaskPath path for mask -# @param nb_partitions number of k-means then averaged -# @param Pix_Per_Partition number of pixels per iteration +# Extract bands of sparse pixels in big image data cube, typically hyperspectral data cube. +# @param ImPath character. Path to the image cube +# @param rowcol matrix or data.frame with two columns: row, col. +# If columns are not named, 1st=row, 2nd=col. +# @param MaxRAM numeric. Maximum memory use at block reading. +# It constrains the maximum number rows of a block # -# @return samples from image, coordinates of these samples and updated number of pixels to sample +# @return matrix. Rows are corresponding to the samples, columns are the bands. +#' @import stars +extract.big_raster <- function(ImPath, rowcol, MaxRAM=.25){ + # ImPath = hs_file + # rowcol = arcd + # rowcol = as.data.table(rowcol) + # MaxRAM = .25 + + if(!is.data.frame(rowcol)){ + rowcol <- as.data.frame(rowcol) + } + + if(!all(c('row', 'col') %in% colnames(rowcol))){ + warning('Columns row,col not found in rowcol argument. The two first columns are considered as row, col respectively.') + colnames(rowcol)[1:2]= c('row', 'col') + } + + + r = brick(ImPath) + # nbytes = as.numeric(substring(dataType(r), 4, 4)) + # stars converts automatically values to numeric + nbytes = 8 + ImgSizeGb = prod(dim(r))*nbytes/2^30 + LineSizeGb = prod(dim(r)[2:3])*nbytes/2^30 + LinesBlock = floor(MaxRAM/LineSizeGb) + rowcol$rowInBlock = ((rowcol$row-1) %% LinesBlock)+1 # row in block number + rowcol$block=floor((rowcol$row-1)/LinesBlock)+1 # block number + rowcol$sampleIndex = 1:nrow(rowcol) # sample index to reorder result + + sampleList = lapply(unique(rowcol$block), function(iblock){ + rc = rowcol[rowcol$block==iblock,] + rr = range(rc$row) + nYSize = diff(rr)+1 + nXSize = max(rc$col) + # stars data cube dimension order is x*y*band + ipix_stars = (rc$rowInBlock-min(rc$rowInBlock))*nXSize+rc$col + values = read_stars(ImPath, RasterIO =list(nXSize=nXSize, nYOff=rr[1], nYSize=nYSize))[[1]] + values = matrix(values, nrow=nYSize*nXSize) + res = cbind(rc$sampleIndex, values[ipix_stars, ]) + rm('values') + gc() + return(res) + }) + + samples = do.call(rbind, sampleList) + samples = samples[order(samples[,1]),2:ncol(samples)] + + return(samples) +} + +# @return list, see Details +# @details +# The returned list contains: +# - DataSubset: matrix of NxP of N samples and P bands +# - nbPix2Sample: integer giving the number of pixels sampled (only central pixel of kernel) +# - coordPix: a data.table with columns 'row', 'col' of pixel in the image corresponding to each row of DataSubset, and if kernel is not NULL +# Kind (Kernel index) and 'id' the sample ID to be used with the kernel #' @importFrom matlab ones -get_random_subset_from_image <- function(ImPath, HDR, MaskPath, nb_partitions, Pix_Per_Partition) { +#' @import raster +#' @importFrom mmand erode +#' @importFrom data.table data.table rbindlist setorder +#' @importFrom matrixStats rowAnys +get_random_subset_from_image <- function(ImPath, HDR, MaskPath, nb_partitions, Pix_Per_Partition, kernel=NULL) { + r <- brick(ImPath) nbPix2Sample <- nb_partitions * Pix_Per_Partition # get total number of pixels - nbpix <- as.double(HDR$lines) * as.double(HDR$samples) + rdim <- dim(r) + nlines <- rdim[1] + nsamples <- rdim[2] + nbpix <- ncell(r) # 1- Exclude masked pixels from random subset # Read Mask if ((!MaskPath == "") & (!MaskPath == FALSE)) { - fid <- file( - description = MaskPath, open = "rb", blocking = TRUE, - encoding = getOption("encoding"), raw = FALSE - ) - ShadeMask <- readBin(fid, integer(), n = nbpix, size = 1) - close(fid) - ShadeMask <- aperm(array(ShadeMask, dim = c(HDR$samples, HDR$lines))) + mask <- matrix(raster(MaskPath),ncol= HDR$samples,nrow = HDR$lines) + if(any(dim(mask)[1:2] != dim(r)[1:2])){ + stop('Mask and Image rasters do not have the same XY diemnsions.') + } } else { - ShadeMask <- array(ones(HDR$lines * HDR$samples, 1), dim = c(HDR$lines, HDR$samples)) + mask <- array(1, dim = c(nlines, nsamples)) + } + + if(is.matrix(kernel)){ + # erode mask with kernel, to keep valid central pixels and neighbours + mask = matlab::padarray(mask, c(1,1), padval=0, direction='both') + mask = erode(mask, (kernel!=0)*1) + mask = mask[2:(nrow(mask)-1), 2:(ncol(mask)-1)] } + # get a list of samples among unmasked pixels - IndexInit <- (matrix(1:nbpix, ncol = HDR$samples)) - ValidPixels <- sample(which(ShadeMask == 1 | ShadeMask == 255)) + + ValidPixels <- which(mask > 0) NbValidPixels <- length(ValidPixels) # Check if number of valid pixels is compatible with number of pixels to be extracted # if number of pixels to sample superior to number of valid pixels, then adjust iterations @@ -416,24 +494,57 @@ get_random_subset_from_image <- function(ImPath, HDR, MaskPath, nb_partitions, P Pix_Per_Partition <- floor(NbValidPixels / nb_partitions) nbPix2Sample <- nb_partitions * Pix_Per_Partition } - # Select a subset of nbPix2Sample among pixselected - pixselected <- ValidPixels[1:nbPix2Sample] - - # define location of pixselected in binary file (avoid multiple reads) and optimize access to disk - # the file is a BIL file, which means that for each line in the image, - # we first need to define if datapoints are on the line, then read one point after the other - coordi <- ((pixselected - 1) %% HDR$lines) + 1 - coordj <- floor((pixselected - 1) / HDR$lines) + 1 - # sort based on line and col - indxPix <- order(coordi, coordj, na.last = FALSE) - coordPix <- cbind(coordi[indxPix], coordj[indxPix]) + + # Select a random subset of nbPix2Sample + seed <- sample(10000)[1] + set.seed(seed) + pixselected <- sample(ValidPixels, nbPix2Sample) + Row <- ((pixselected - 1) %% nlines) + 1 # row + Column <- floor((pixselected - 1) / nlines) + 1 # column + if(is.matrix(kernel)){ + coordPixK = list() + mesh=matlab::meshgrid(-(ncol(kernel)%/%2):(ncol(kernel)%/%2), -(nrow(kernel)%/%2):(nrow(kernel)%/%2)) + for(p in which(kernel!=0)){ + coordPixK[[p]] = data.table(row = Row+mesh$y[p], col = Column+mesh$x[p], id=1:length(Row)) + } + coordPix = rbindlist(coordPixK, idcol='Kind') + # Order along coordPix$id for further use in noise, mnf + setorder(coordPix, 'id') + }else{ + coordPix = data.table(row = Row, col = Column, id = 1:length(Row)) + } + # sort based on .bil dim order, i.e. band.x.y or band.col.row + # TODO: sorting may not be necessary anymore, neither unique coordinates + # setorder(coordPix, col, row) + + # make unique + ucoordPix <- unique(coordPix[,c('row','col')]) + ucoordPix[['sampleIndex']] = 1:nrow(ucoordPix) # 2- Extract samples from image - Sample_Sel <- extract_samples_from_image(ImPath, coordPix) - # randomize! - Sample_Sel <- Sample_Sel[sample(dim(Sample_Sel)[1]), ] - coordPix <- coordPix[sample(dim(Sample_Sel)[1]), ] - my_list <- list("DataSubset" = Sample_Sel, "nbPix2Sample" = nbPix2Sample,"coordPix"=coordPix) + # TODO: if coordPix can be a data.frame it would be easier + # Sample_Sel <- biodivMapR:::extract_samples_from_image(ImPath, ucoordPix) + Sample_Sel <- extract.big_raster(ImPath, ucoordPix[,1:2]) + samplePixIndex = merge(coordPix, ucoordPix, by=c('row', 'col'), sort = FALSE) + # randomize! It should already be random from the sample operation on mask valid pixels + Sample_Sel <- Sample_Sel[samplePixIndex$sampleIndex, ] + samplePixIndex[['sampleIndex']]=NULL + + # remove NA and Inf pixels + if(any(is.na(Sample_Sel) | is.infinite(Sample_Sel))){ + print('Removing pixels with NA values.') + rmrows <- !rowAnys(is.na(Sample_Sel) | is.infinite(Sample_Sel)) + rmpix <- unique(samplePixIndex$id[rmrows]) + Sample_Sel = Sample_Sel[!(samplePixIndex$id %in% rmpix),] + samplePixIndex = samplePixIndex[!(samplePixIndex$id %in% rmpix)] + nbPix2Sample = length(unique(samplePixIndex$id)) + } + + ### FLORIAN: REMOVED RANDOMIZATION AS IT SEEMED USELESS + # Sample_Sel <- Sample_Sel[sample(dim(Sample_Sel)[1]), ] + # coordPix <- coordPix[sample(dim(Sample_Sel)[1]), ] + # TODO: check how returned coordPix is used, as it is now a data.frame + my_list <- list("DataSubset" = Sample_Sel, "nbPix2Sample" = nbPix2Sample,"coordPix"=samplePixIndex) return(my_list) } @@ -525,18 +636,31 @@ ind2sub2 <- function(Raster, Image_Index) { # # @return rank of all spectral bands of interest in the image and corresponding wavelength #' @importFrom matlab padarray -mean_filter <- function(ImageInit, nbi, nbj, SizeFilt) { - E <- padarray(ImageInit, c(SizeFilt, SizeFilt), "symmetric", "both") - ImageSmooth <- matrix(0, nrow = nbi, ncol = nbj) - Mat2D <- MatSun <- matrix(0, nrow = ((2 * SizeFilt) + 1)^2, ncol = nbj) - spl <- split(1:nbj, 1:((2 * SizeFilt) + 1)) - mid <- ceiling((((2 * SizeFilt) + 1)^2) / 2) - for (i in (SizeFilt + 1):(nbi + SizeFilt)) { - for (j in 1:((2 * SizeFilt) + 1)) { - # create a 2D matrix - Mat2D[, spl[[j]]] <- matrix(E[(i - SizeFilt):(i + SizeFilt), (spl[[j]][1]):(tail(spl[[j]], n = 1) + 2 * SizeFilt)], nrow = ((2 * SizeFilt) + 1)^2) +mean_filter <- function(Image, SizeFilt,NA_remove = FALSE) { + nbi <- dim(Image)[1] + nbj <- dim(Image)[2] + # if matrix: convert into array + if (is.na(dim(Image)[3])){ + Image <- array(Image,dim = c(nbi,nbj,1)) + } + ImageSmooth <- array(0,c(nbi,nbj,dim(Image)[3])) + for (band in 1: dim(Image)[3]){ + E <- padarray(Image[,,band], c(SizeFilt, SizeFilt), "symmetric", "both") + ImageSmooth_tmp <- matrix(0, nrow = nbi, ncol = nbj) + Mat2D <- MatSun <- matrix(0, nrow = ((2 * SizeFilt) + 1)^2, ncol = nbj) + spl <- split(1:nbj, 1:((2 * SizeFilt) + 1)) + mid <- ceiling((((2 * SizeFilt) + 1)^2) / 2) + for (i in (SizeFilt + 1):(nbi + SizeFilt)) { + for (j in 1:((2 * SizeFilt) + 1)) { + # create a 2D matrix + Mat2D[, spl[[j]]] <- matrix(E[(i - SizeFilt):(i + SizeFilt), (spl[[j]][1]):(tail(spl[[j]], n = 1) + 2 * SizeFilt)], nrow = ((2 * SizeFilt) + 1)^2) + } + ImageSmooth_tmp[(i - SizeFilt), ] <- colMeans(Mat2D, na.rm = TRUE) + } + if (NA_remove == TRUE){ + ImageSmooth_tmp[which(is.na(Image[,,band]))] <- NA } - ImageSmooth[(i - SizeFilt), ] <- colMeans(Mat2D, na.rm = TRUE) + ImageSmooth[,,band] <- ImageSmooth_tmp } return(ImageSmooth) } @@ -645,6 +769,8 @@ read_ENVI_header <- function(HDRpath) { # @param ImBand Bands to be read # # @return Image_Subset information corresponding to ImBand +#' @import stars + read_image_bands <- function(ImPath, HDR, ImBand) { # first get image format Image_Format <- ENVI_type2bytes(HDR) @@ -652,42 +778,45 @@ read_image_bands <- function(ImPath, HDR, ImBand) { jpix <- HDR$samples nbChannels <- HDR$bands nbSubset <- length(ImBand) - # then open and read image - # depending on image size, need to read in one or multiple times - lenTot <- as.double(ipix) * as.double(jpix) * as.double(nbChannels) - lenSubset <- as.double(ipix) * as.double(jpix) * as.double(nbSubset) - ImSizeGb <- (lenTot * Image_Format$Bytes) / (1024^3) - # maximum image size read at once. If image larger, then reads in multiple pieces - LimitSizeGb <- 0.25 - if (ImSizeGb < LimitSizeGb) { - nbLinesPerCPU <- ceiling(ipix) - nbPieces <- 1 - } else { - # nb of lines corresponding to LimitSizeGb - OneLine <- as.double(jpix) * as.double(nbChannels) * Image_Format$Bytes - nbLinesPerCPU <- floor(LimitSizeGb * (1024^3) / OneLine) - # number of pieces to split the image into - nbPieces <- ceiling(ipix / nbLinesPerCPU) - } - # Define segments of image to be read - SeqRead.Image <- where_to_read(HDR, nbPieces) - # Read segments (subsets) of image - Image_Subsets <- list() - for (i in 1:nbPieces) { - # number of elements to be read - Byte_Start <- SeqRead.Image$ReadByte_Start[i] - nbLines <- SeqRead.Image$Lines_Per_Chunk[i] - lenBin <- SeqRead.Image$ReadByte_End[i] - SeqRead.Image$ReadByte_Start[i] + 1 - Image_Subsets[[i]] <- read_bin_subset(Byte_Start, nbLines, lenBin, ImPath, ImBand, jpix, nbChannels, Image_Format) + Image_Subset <- array(NA,c(ipix,jpix,length(ImBand))) + i <- 0 + for (band in ImBand){ + i <- i+1 + bndtmp = t(matrix(raster(ImPath, band = ImBand[i]),nrow = HDR$samples,ncol = HDR$lines)) + Image_Subset[,,i] <- array(bndtmp,c(ipix,jpix,1)) } - # reshape image with original size and selected bands - Image_Subset <- build_image_from_list(Image_Subsets, ipix, jpix, nbSubset) - rm(Image_Subsets) - gc() return(Image_Subset) } -# reads subset of an image +# reads subset of lines from an image +# +# @param ImPath path for the image +# @param HDR header information corresponding to the image +# @param Line_Start which line to start reading +# @param Lines_To_Read number of lines to read +# +# @return data corresponding to the subset in original 3D format +read_image_subset <- function(ImPath, HDR, Line_Start,Lines_To_Read,ImgFormat='3D'){ + # list of pixels to be extracted + ListRows <- seq(Line_Start,Line_Start+Lines_To_Read-1) + ListCols <- seq(1,HDR$samples) + # list pixels + coordPix <- expand.grid(ListRows,ListCols) + names(coordPix) <- c('row','col') + coordPix[['sampleIndex']] = 1:nrow(coordPix) + # 2- Extract samples from image + Sample_Sel <- extract.big_raster(ImPath, coordPix[,1:2]) + Sample_Sel <- array(matrix(as.numeric(unlist(Sample_Sel)),ncol = HDR$bands),c(Lines_To_Read,HDR$samples,HDR$bands)) + if (ImgFormat == "2D") { + Sample_Sel <- array(Sample_Sel, c(Lines_To_Read * HDR$samples, HDR$bands)) + } + if (ImgFormat == "Shade") { + Sample_Sel <- matrix(Sample_Sel, Lines_To_Read, HDR$samples, byrow = F) + } + return(Sample_Sel) +} + +# reads subset of an ENVI BIL image # # @param ImPath path for the image # @param HDR header information corresponding to the image @@ -698,7 +827,7 @@ read_image_bands <- function(ImPath, HDR, ImBand) { # @param ImgFormat should output be 2D or 3D (original image format)? # # @return data corresponding to the subset in original 3D format -read_image_subset <- function(ImPath, HDR, Byte_Start, lenBin, nbLines, Image_Format, ImgFormat) { +read_BIL_image_subset <- function(ImPath, HDR, Byte_Start, lenBin, nbLines, Image_Format, ImgFormat) { fidIm <- file( description = ImPath, open = "rb", blocking = TRUE, encoding = getOption("encoding"), raw = FALSE @@ -843,6 +972,7 @@ update_shademask <- function(MaskPath, HDR, Mask, MaskPath_Update) { HDR_Update$`band names` <- { "Mask" } + HDR_Update$`default stretch` <- '0 1 linear' HDR_Update$wavelength <- NULL HDR_Update$fwhm <- NULL HDR_Update$resolution <- NULL @@ -984,7 +1114,11 @@ Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){ } ImgChunk <- aperm(ImgChunk, c(2, 3, 1)) # writeBin(as.numeric(c(PCA_Chunk)), fidOUT, size = PCA_Format$Bytes,endian = .Platform$endian) - writeBin(c(ImgChunk), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) + if (!Image_Format$Bytes==1){ + writeBin(c(ImgChunk), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) + } else { + writeBin(c(as.integer(ImgChunk)), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) + } close(fidOUT) } rm(ImgWrite) @@ -993,6 +1127,102 @@ Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){ return("") } +#' write an image resulting from "window processing" at native spatial resolution +#' (assuming square windows & origin at top left corner) +#' +#' @param Image numeric. Image corresponding to a 2D matrix or 3D array +#' @param ImagePath character. Path where the image should be written +#' @param HDR list. Image header +#' @param window_size numeric. window size used to generate Image +#' +#' @return None +#' @export + +Write_Image_NativeRes <- function(Image,ImagePath,HDR,window_size){ + # update HDR: dimensions & spatial resolution + HDR_Full <- HDR + HDR_Full$samples <- HDR$samples * window_size + HDR_Full$lines <- HDR$lines * window_size + HDR_Full <- revert_resolution_HDR(HDR_Full, window_size) + # define name for native resolution image and corresponding HDR + ImagePath_FullRes <- paste(ImagePath, "_Fullres", sep = "") + headerFpath <- paste(ImagePath_FullRes, ".hdr", sep = "") + # write header + write_ENVI_header(HDR_Full, headerFpath) + Image_Format <- ENVI_type2bytes(HDR_Full) + # create image the same size as native image: convert matrix into array + if (length(dim(Image))==2){ + Image = array(Image,c(HDR$lines,HDR$samples,1)) + } + Image_FullRes <- array(NA,c(HDR_Full$lines,HDR_Full$samples,HDR_Full$bands)) + for (band in 1:HDR_Full$bands){ + for (i in 1:HDR$lines) { + for (j in 1:HDR$samples) { + Image_FullRes[((i-1)*window_size+1):(i*window_size),((j-1)*window_size+1):(j*window_size),band] <- Image[i,j,band] + } + } + } + # write image and make sure size does not matter ... + Write_Big_Image(ImgWrite = Image_FullRes,ImagePath = ImagePath_FullRes, + HDR = HDR_Full,Image_Format = Image_Format) + # zip resulting file + ZipFile(ImagePath_FullRes) + return("") +} + + +# Writes a matrix or an array into a ENVI BIL raster +# +#' @param Image numeric. matrix or array of image to be written +#' @param HDR hdr template +#' @param ImagePath path of image file to be written +#' @param window_size spatial units use dto compute diversiy (in pixels) +#' @param FullRes should full resolution image be written (original pixel size) +#' @param LowRes should low resolution image be written (one value per spatial unit) +# +#' @return None + +write_raster <- function(Image, HDR, ImagePath, window_size, FullRes = TRUE, LowRes = FALSE,SmoothImage = FALSE) { + + Image_Format <- ENVI_type2bytes(HDR) + # Write image with resolution corresponding to window_size + if (LowRes == TRUE) { + # write header + headerFpath <- paste(ImagePath, ".hdr", sep = "") + write_ENVI_header(HDR, headerFpath) + # write image + Write_Big_Image(Image,ImagePath,HDR,Image_Format) + } + + # Write image with Full native resolution + if (FullRes == TRUE) { + Write_Image_NativeRes(Image,ImagePath,HDR,window_size) + } + + # write smoothed map + if (SmoothImage == TRUE){ + SizeFilt <- 1 + if (min(c(dim(Image)[1], dim(Image)[2])) > (2 * SizeFilt + 1)) { + Image_Smooth <- mean_filter(Image, SizeFilt,NA_remove = TRUE) + ImagePath.Smooth <- paste(ImagePath, "_MeanFilter", sep = "") + if (LowRes == TRUE) { + # write header + headerFpath <- paste(ImagePath.Smooth, ".hdr", sep = "") + write_ENVI_header(HDR, headerFpath) + # write image + Write_Big_Image(Image_Smooth,ImagePath.Smooth,HDR,Image_Format) + # close(fidOUT) + } + if (FullRes == TRUE) { + # Write image with Full native resolution + Write_Image_NativeRes(Image_Smooth,ImagePath.Smooth,HDR,window_size) + } + } + } + return("") +} + + # convert image coordinates from X-Y to index # # @param HDR_Raster @@ -1000,7 +1230,7 @@ Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){ # # @return Image_Index sub2ind <- function(HDR_Raster, Pixels) { - Image_Index <- (Pixels$Column - 1) * HDR_Raster$lines + Pixels$Row + Image_Index <- (Pixels$col - 1) * HDR_Raster$lines + Pixels$row return(Image_Index) } @@ -1050,12 +1280,12 @@ revert_resolution_HDR <- function(HDR, window_size) { # @param ImagePath path for the image # @return None ZipFile <- function(ImagePath) { - PathRoot <- getwd() - ImageDir <- dirname(ImagePath) - ImageName <- basename(ImagePath) - setwd(ImageDir) - zip::zip(zipfile = paste0(ImageName, ".zip"), files = ImageName) - file.remove(ImageName) - setwd(PathRoot) + + + ImagePath <- normalizePath(ImagePath) + + zip::zipr(zipfile = paste0(ImagePath, ".zip"), files = ImagePath) + file.remove(ImagePath) + return(invisible()) } diff --git a/R/Lib_MapAlphaDiversity.R b/R/Lib_MapAlphaDiversity.R index cd94f628..c36eba26 100644 --- a/R/Lib_MapAlphaDiversity.R +++ b/R/Lib_MapAlphaDiversity.R @@ -48,36 +48,43 @@ map_alpha_div <- function(Input_Image_File, Output_Dir, window_size, if (length((grep("Fisher", Index_Alpha))) > 0) Fisher <- TRUE Output_Dir_Alpha <- define_output_subdir(Output_Dir, Input_Image_File, TypePCA, "ALPHA") + HDR <- ALPHA$HDR if (Shannon == TRUE) { Index <- "Shannon" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Shannon, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Shannon, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) if (MapSTD == TRUE) { Index <- "Shannon_SD" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Shannon.SD, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Shannon.SD, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) } } if (Fisher == TRUE) { Index <- "Fisher" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Fisher, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Fisher, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) if (MapSTD == TRUE) { Index <- "Fisher_SD" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Fisher.SD, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Fisher.SD, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) } } if (Simpson == TRUE) { Index <- "Simpson" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Simpson, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Simpson, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) if (MapSTD == TRUE) { Index <- "Simpson_SD" + HDR$`band names` <- Index Alpha_Path <- paste(Output_Dir_Alpha, Index, "_", window_size, sep = "") - write_raster_alpha(ALPHA$Simpson.SD, ALPHA$HDR, Alpha_Path, window_size, Index, FullRes = FullRes, LowRes = LowRes) + write_raster(ALPHA$Simpson.SD, HDR, Alpha_Path, window_size, FullRes = FullRes, LowRes = LowRes, SmoothImage = TRUE) } } return(invisible()) @@ -85,14 +92,14 @@ map_alpha_div <- function(Input_Image_File, Output_Dir, window_size, # Map alpha diversity metrics based on spectral species # -# @param Spectral_Species_Path path for spectral species file to be written -# @param window_size size of spatial units (in pixels) to compute diversity -# @param nbclusters number of clusters defined in k-Means -# @param pcelim -# @param nbCPU -# @param MaxRAM -# @param Index_Alpha -# @param MinSun minimum proportion of sunlit pixels required to consider plot +# @param Spectral_Species_Path character. path for spectral species file to be written +# @param window_size numeric. size of spatial units (in pixels) to compute diversity +# @param nbclusters numeric. number of clusters defined in k-Means +# @param pcelim numeric. percentage of occurence of a cluster below which cluster is eliminated +# @param nbCPU numeric. number of CPUs available +# @param MaxRAM numeric. maximum RAM available +# @param Index_Alpha list. list of alpha diversity indices to be computed from spectral species +# @param MinSun numeric. minimum proportion of sunlit pixels required to consider plot # # @return list of mean and SD of alpha diversity metrics #' @importFrom future plan multiprocess sequential @@ -203,9 +210,13 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters Shannon.SD <- do.call(rbind, Shannon_SD_Chunk) Fisher.SD <- do.call(rbind, Fisher_SD_Chunk) Simpson.SD <- do.call(rbind, Simpson_SD_Chunk) + # prepare HDR for alpha diversity + HDR <- HDR_SSD + HDR$bands <- 1 + HDR$`data type` <- 4 my_list <- list( "Shannon" = Shannon.Mean, "Fisher" = Fisher.Mean, "Simpson" = Simpson.Mean, - "Shannon.SD" = Shannon.SD, "Fisher.SD" = Fisher.SD, "Simpson.SD" = Simpson.SD, "HDR" = HDR_SSD + "Shannon.SD" = Shannon.SD, "Fisher.SD" = Fisher.SD, "Simpson.SD" = Simpson.SD, "HDR" = HDR ) return(my_list) } @@ -234,7 +245,7 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters convert_PCA_to_SSD <- function(ReadWrite, Spectral_Species_Path, HDR_SS, HDR_SSD, SS_Format, SSD_Format, ImgFormat, window_size, nbclusters, MinSun, pcelim, Index_Alpha, SSD_Path, Sunlit_Path, Sunlit_Format) { - SS_Chunk <- read_image_subset( + SS_Chunk <- read_BIL_image_subset( Spectral_Species_Path, HDR_SS, ReadWrite$RW_SS$Byte_Start, ReadWrite$RW_SS$lenBin, ReadWrite$RW_SS$nbLines, SS_Format, ImgFormat @@ -288,7 +299,7 @@ convert_PCA_to_SSD <- function(ReadWrite, Spectral_Species_Path, HDR_SS, HDR_SSD # @param window_size size of spatial units (in pixels) to compute diversity # @param nbclusters number of clusters defined in k-Means # @param MinSun minimum proportion of sunlit pixels required to consider plot -# @param Index_Alpha +# @param Index_Alpha list. list of alpha diversity indices to be computed from spectral species # @param pcelim minimum proportion for a spectral species to be included # # @return list of alpha diversity metrics for each iteration @@ -380,87 +391,3 @@ get_Simpson <- function(Distrib) { Simpson <- 1 - sum(Distrib * Distrib, na.rm = TRUE) return(Simpson) } - -# Writes image of alpha diversity indicator (1 band) and smoothed alpha diversity -# -# @param Image 2D matrix of image to be written -# @param HDR_SSD hdr template derived from SSD to modify -# @param ImagePath path of image file to be written -# @param window_size spatial units use dto compute diversiy (in pixels) -# @param Index name of the index (eg. Shannon) -# @param FullRes should full resolution image be written (original pixel size) -# @param LowRes should low resolution image be written (one value per spatial unit) -# -# @return -write_raster_alpha <- function(Image, HDR_SSD, ImagePath, window_size, Index, FullRes = TRUE, LowRes = FALSE) { - - # Write image with resolution corresponding to window_size - HDR_Alpha <- HDR_SSD - HDR_Alpha$bands <- 1 - HDR_Alpha$`data type` <- 4 - HDR_Alpha$`band names` <- Index - Image_Format <- ENVI_type2bytes(HDR_Alpha) - if (LowRes == TRUE) { - # write header - headerFpath <- paste(ImagePath, ".hdr", sep = "") - write_ENVI_header(HDR_Alpha, headerFpath) - # write image and make sure size does not matter ... - Write_Big_Image(Image,ImagePath,HDR_Alpha,Image_Format) - } - if (FullRes == TRUE) { - # Write image with Full native resolution - HDR_Full <- HDR_Alpha - HDR_Full$samples <- HDR_Alpha$samples * window_size - HDR_Full$lines <- HDR_Alpha$lines * window_size - HDR_Full <- revert_resolution_HDR(HDR_Full, window_size) - ImagePath_FullRes <- paste(ImagePath, "_Fullres", sep = "") - headerFpath <- paste(ImagePath_FullRes, ".hdr", sep = "") - write_ENVI_header(HDR_Full, headerFpath) - Image_Format <- ENVI_type2bytes(HDR_Full) - Image_FullRes <- matrix(NA, ncol = HDR_Full$samples, nrow = HDR_Full$lines) - for (i in 1:HDR_SSD$lines) { - for (j in 1:HDR_SSD$samples) { - Image_FullRes[((i - 1) * window_size + 1):(i * window_size), ((j - 1) * window_size + 1):(j * window_size)] <- Image[i, j] - } - } - # write image and make sure size does not matter ... - Write_Big_Image(Image_FullRes,ImagePath_FullRes,HDR_Full,Image_Format) - # zip resulting file - ZipFile(ImagePath_FullRes) - } - # write smoothed map - SizeFilt <- 1 - nbi <- dim(Image)[1] - nbj <- dim(Image)[2] - if (min(c(nbi, nbj)) > (2 * SizeFilt + 1)) { - Image_Smooth <- mean_filter(Image, nbi, nbj, SizeFilt) - Image_Smooth[which(is.na(Image))] <- NA - ImagePath.Smooth <- paste(ImagePath, "_MeanFilter", sep = "") - headerFpath <- paste(ImagePath.Smooth, ".hdr", sep = "") - Image_Format <- ENVI_type2bytes(HDR_Alpha) - if (LowRes == TRUE) { - write_ENVI_header(HDR_Alpha, headerFpath) - # write image and make sure size does not matter ... - Write_Big_Image(Image_Smooth,ImagePath.Smooth,HDR_Alpha,Image_Format) - # close(fidOUT) - } - if (FullRes == TRUE) { - # Write image with Full native resolution - ImagePath_FullRes <- paste(ImagePath.Smooth, "_Fullres", sep = "") - headerFpath <- paste(ImagePath_FullRes, ".hdr", sep = "") - write_ENVI_header(HDR_Full, headerFpath) - Image_Format <- ENVI_type2bytes(HDR_Full) - Image_FullRes <- matrix(NA, ncol = HDR_Full$samples, nrow = HDR_Full$lines) - for (i in 1:HDR_SSD$lines) { - for (j in 1:HDR_SSD$samples) { - Image_FullRes[((i - 1) * window_size + 1):(i * window_size), ((j - 1) * window_size + 1):(j * window_size)] <- Image_Smooth[i, j] - } - } - # write image and make sure size does not matter ... - Write_Big_Image(Image_FullRes,ImagePath_FullRes,HDR_Full,Image_Format) - # zip resulting file - ZipFile(ImagePath_FullRes) - } - } - return("") -} diff --git a/R/Lib_MapBetaDiversity.R b/R/Lib_MapBetaDiversity.R index 1ecf9bab..5661dc41 100644 --- a/R/Lib_MapBetaDiversity.R +++ b/R/Lib_MapBetaDiversity.R @@ -45,7 +45,9 @@ map_beta_div <- function(Input_Image_File, Output_Dir, window_size, print("Write beta diversity maps") Index <- paste("BetaDiversity_BCdiss_", scaling, sep = "") Beta.Path <- paste(Output_Dir_BETA, Index, "_", window_size, sep = "") - write_raster_beta(Beta$BetaDiversity, Beta$HDR, Beta.Path, window_size, FullRes = FullRes, LowRes = LowRes) + write_raster(Image = Beta$BetaDiversity, HDR = Beta$HDR, ImagePath = Beta.Path, + window_size = window_size, FullRes = FullRes, LowRes = LowRes, + SmoothImage = FALSE) return(invisible()) } @@ -226,9 +228,11 @@ compute_beta_metrics <- function(Output_Dir, MinSun, Nb_Units_Ordin, nb_partitio MatBCdist <- as.dist(MatBC, diag = FALSE, upper = FALSE) if (scaling == "NMDS") { Beta_Ordination_sel <- compute_NMDS(MatBCdist) + PCname <- 'NMDS' } else if (scaling == "PCO") { BetaPCO <- pco(MatBCdist, k = 3) Beta_Ordination_sel <- BetaPCO$points + PCname <- 'PCoA' } # Perform nearest neighbor on spatial units excluded from Ordination @@ -247,9 +251,21 @@ compute_beta_metrics <- function(Output_Dir, MinSun, Nb_Units_Ordin, nb_partitio BetaDiversityRGB[, , i] <- BetaTmp } list <- ls() - rm(list = list[-which(list == "BetaDiversityRGB" | list == "Select_Sunlit" | list == "HDR_Sunlit")]) + + # update HDR file + HDR_Beta <- HDR_Sunlit + HDR_Beta$bands <- 3 + HDR_Beta$`data type` <- 4 + PCs <- list() + for (i in 1:3) { + PCs <- c(PCs, paste(PCname,'#', i,sep = '')) + } + PCs <- paste(PCs, collapse = ", ") + HDR_Beta$`band names` <- PCs + + rm(list = list[-which(list == "BetaDiversityRGB" | list == "Select_Sunlit" | list == "HDR_Beta")]) gc() - my_list <- list("BetaDiversity" = BetaDiversityRGB, "Select_Sunlit" = Select_Sunlit, "HDR" = HDR_Sunlit) + my_list <- list("BetaDiversity" = BetaDiversityRGB, "Select_Sunlit" = Select_Sunlit, "HDR" = HDR_Beta) return(my_list) } @@ -311,74 +327,6 @@ compute_NN_from_ordination <- function(MatBC3, knn, BetaDiversity0) { return(Ordin_est) } -# Writes image of beta diversity indicator (3 bands) resulting from BC + NMDS -# -# @param Image 2D matrix of image to be written -# @param HDR_SSD hdr template derived from SSD to modify -# @param ImagePath path of image file to be written -# @param window_size spatial units use dto compute diversiy (in pixels) -# @param FullRes should full resolution image be written (original pixel size) -# @param LowRes should low resolution image be written (one value per spatial unit) -# -# @return -write_raster_beta <- function(Image, HDR_SSD, ImagePath, window_size, FullRes = TRUE, LowRes = FALSE) { - # Write image with resolution corresponding to window_size - HDR_Beta <- HDR_SSD - HDR_Beta$bands <- 3 - HDR_Beta$`data type` <- 4 - PCs <- list() - for (i in 1:3) { - PCs <- c(PCs, paste("NMDS ", i)) - } - PCs <- paste(PCs, collapse = ", ") - HDR_Beta$`band names` <- PCs - Image_Format <- ENVI_type2bytes(HDR_Beta) - if (LowRes == TRUE) { - headerFpath <- paste(ImagePath, ".hdr", sep = "") - write_ENVI_header(HDR_Beta, headerFpath) - # write image and make sure size does not matter ... - Write_Big_Image(Image,ImagePath,HDR_Beta,Image_Format) - # ImgWrite <- aperm(Image, c(2, 3, 1)) - # fidOUT <- file( - # description = ImagePath, open = "wb", blocking = TRUE, - # encoding = getOption("encoding"), raw = FALSE - # ) - # writeBin(c(ImgWrite), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) - # close(fidOUT) - } - if (FullRes == TRUE) { - # Write image with Full native resolution - HDR_Full <- HDR_Beta - HDR_Full$samples <- HDR_Beta$samples * window_size - HDR_Full$lines <- HDR_Beta$lines * window_size - HDR_Full <- revert_resolution_HDR(HDR_Full, window_size) - ImagePath_FullRes <- paste(ImagePath, "_Fullres", sep = "") - headerFpath <- paste(ImagePath_FullRes, ".hdr", sep = "") - write_ENVI_header(HDR_Full, headerFpath) - Image_Format <- ENVI_type2bytes(HDR_Full) - Image_FullRes <- array(NA, c(HDR_Full$lines, HDR_Full$samples, 3)) - for (b in 1:3) { - for (i in 1:HDR_SSD$lines) { - for (j in 1:HDR_SSD$samples) { - Image_FullRes[((i - 1) * window_size + 1):(i * window_size), ((j - 1) * window_size + 1):(j * window_size), b] <- Image[i, j, b] - } - } - } - # write image and make sure size does not matter ... - Write_Big_Image(Image_FullRes,ImagePath_FullRes,HDR_Full,Image_Format) - # ImgWrite <- aperm(Image_FullRes, c(2, 3, 1)) - # fidOUT <- file( - # description = ImagePath_FullRes, open = "wb", blocking = TRUE, - # encoding = getOption("encoding"), raw = FALSE - # ) - # writeBin(c(ImgWrite), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) - # close(fidOUT) - # zip resulting file - ZipFile(ImagePath_FullRes) - } - return("") -} - # Gets sunlit pixels from SpectralSpecies_Distribution_Sunlit # # @param ImPathSunlit path for SpectralSpecies_Distribution_Sunlit diff --git a/R/Lib_MapSpectralSpecies.R b/R/Lib_MapSpectralSpecies.R index 095d9a02..b4d899a6 100644 --- a/R/Lib_MapSpectralSpecies.R +++ b/R/Lib_MapSpectralSpecies.R @@ -4,6 +4,7 @@ # ============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2019/06 Jean-Baptiste FERET # ============================================================================== # This Library applies clustering on a selection of components stored in a PCA @@ -57,7 +58,9 @@ map_spectral_species <- function(Input_Image_File,Output_Dir,PCA_Files,PCA_model # sample data from image and perform PCA ImPathHDR <- get_HDR_name(Input_Image_File) HDR <- read_ENVI_header(ImPathHDR) - Subset <- get_random_subset_from_image(Input_Image_File, HDR, Input_Mask_File, nb_partitions, Pix_Per_Partition) + Subset <- get_random_subset_from_image(ImPath = Input_Image_File, HDR = HDR, MaskPath = Input_Mask_File, + nb_partitions = nb_partitions, Pix_Per_Partition = Pix_Per_Partition, + kernel = NULL) SubsetInit <- Subset # if needed, apply continuum removal if (Continuum_Removal == TRUE) { @@ -69,7 +72,7 @@ map_spectral_species <- function(Input_Image_File,Output_Dir,PCA_Files,PCA_model Subset$DataSubset <- Subset$DataSubset[, -SpectralFilter$BandsNoVar] } if (TypePCA == "PCA" | TypePCA == "SPCA") { - dataPCA <- t(t(PCA_model$eiV[, 1:PCA_model$Nb_PCs]) %*% t(center_reduce(Subset$DataSubset, PCA_model$mu, PCA_model$scale))) + dataPCA <- scale(Subset$DataSubset, PCA_model$center, PCA_model$scale) %*% PCA_model$rotation[, 1:PCA_model$Nb_PCs] } dataPCA <- dataPCA[, PC_Select] if (length(PC_Select) == 1) { @@ -263,14 +266,14 @@ compute_spectral_species <- function(PCA_Path, Input_Mask_File, Spectral_Species HDR_Shade <- read_ENVI_header(ShadeHDR) Shade.Format <- ENVI_type2bytes(HDR_Shade) ImgFormat <- "Shade" - Shade_Chunk <- read_image_subset(Input_Mask_File, HDR_Shade, Location_RW$Byte_Start_Shade, Location_RW$lenBin_Shade, Location_RW$nbLines, Shade.Format, ImgFormat) + Shade_Chunk <- read_BIL_image_subset(Input_Mask_File, HDR_Shade, Location_RW$Byte_Start_Shade, Location_RW$lenBin_Shade, Location_RW$nbLines, Shade.Format, ImgFormat) PCA_HDR <- get_HDR_name(PCA_Path) HDR_PCA <- read_ENVI_header(PCA_HDR) PCA_Format <- ENVI_type2bytes(HDR_PCA) # read "unfolded" (2D) PCA image ImgFormat <- "2D" - PCA_Chunk <- read_image_subset(PCA_Path, HDR_PCA, Location_RW$Byte_Start_PCA, Location_RW$lenBin_PCA, Location_RW$nbLines, PCA_Format, ImgFormat) + PCA_Chunk <- read_BIL_image_subset(PCA_Path, HDR_PCA, Location_RW$Byte_Start_PCA, Location_RW$lenBin_PCA, Location_RW$nbLines, PCA_Format, ImgFormat) PCA_Chunk <- PCA_Chunk[, PC_Select] if (length(PC_Select) == 1) { PCA_Chunk <- matrix(PCA_Chunk, ncol = 1) diff --git a/R/Lib_PerformPCA.R b/R/Lib_PerformPCA.R index 5a77e812..b1a222ef 100644 --- a/R/Lib_PerformPCA.R +++ b/R/Lib_PerformPCA.R @@ -4,6 +4,7 @@ # ============================================================================== # PROGRAMMERS: # Jean-Baptiste FERET +# Florian de Boissieu # Copyright 2018/07 Jean-Baptiste FERET # ============================================================================== # This Library is used to perform PCA on raster prior to diversity mapping @@ -25,8 +26,9 @@ #' #' @return list of paths corresponding to resulting PCA files #' @export -perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuum_Removal=TRUE, TypePCA="SPCA", NbPCs_To_Keep=FALSE, - FilterPCA = FALSE, Excluded_WL = FALSE, nb_partitions = 20, nbCPU = 1, MaxRAM = 0.25) { +perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuum_Removal=TRUE, TypePCA="SPCA", + NbPCs_To_Keep=30,FilterPCA = FALSE, Excluded_WL = FALSE, nb_partitions = 20, + nbCPU = 1, MaxRAM = 0.25) { # check if format of raster data is as expected check_data(Input_Image_File) if (!Input_Mask_File==FALSE){ @@ -47,7 +49,9 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu ImPathHDR <- get_HDR_name(Input_Image_File) HDR <- read_ENVI_header(ImPathHDR) # extract a random selection of pixels from image - Subset <- get_random_subset_from_image(Input_Image_File, HDR, Input_Mask_File, nb_partitions, Pix_Per_Partition) + Subset <- get_random_subset_from_image(ImPath = Input_Image_File, HDR = HDR, MaskPath = Input_Mask_File, + nb_partitions = nb_partitions , Pix_Per_Partition = Pix_Per_Partition, + kernel=NULL) # if needed, apply continuum removal if (Continuum_Removal == TRUE) { Subset$DataSubset <- apply_continuum_removal(Subset$DataSubset, SpectralFilter, nbCPU = nbCPU) @@ -65,7 +69,7 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu } DataSubset <- Subset$DataSubset # clean reflectance data from inf and constant values - CleanData <- check_invariant_bands(DataSubset, SpectralFilter) + CleanData <- rm_invariant_bands(DataSubset, SpectralFilter) DataSubset <- CleanData$DataMatrix SpectralFilter <- CleanData$Spectral @@ -91,16 +95,20 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu # In order to exclude these pixels, we compute mean and SD for the 3 first # components and exclude all pixels showing values ouside "mean+-3SD" range print("perform 2nd filtering: Exclude extreme PCA values") - if (dim(PCA_model$dataPCA)[2] > 5) { + if (dim(PCA_model$x)[2] > 5) { PCsel <- 1:5 } else { - PCsel <- 1:dim(PCA_model$dataPCA)[2] + PCsel <- 1:dim(PCA_model$x)[2] } Shade_Update <- paste(Output_Dir_Full, "ShadeMask_Update_PCA", sep = "") - Input_Mask_File <- filter_PCA(Input_Image_File, HDR, Input_Mask_File, Shade_Update, SpectralFilter, Continuum_Removal, PCA_model, PCsel, TypePCA, nbCPU = nbCPU, MaxRAM = MaxRAM) + Input_Mask_File <- filter_PCA(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spectral = SpectralFilter, + Continuum_Removal, PCA_model, PCsel, TypePCA, + nbCPU = nbCPU, MaxRAM = MaxRAM) ## Compute PCA 2 based on the updated shade mask ## # extract a random selection of pixels from image - Subset <- get_random_subset_from_image(Input_Image_File, HDR, Input_Mask_File, nb_partitions, Pix_Per_Partition) + Subset <- get_random_subset_from_image(ImPath = Input_Image_File, HDR = HDR, MaskPath = Input_Mask_File, + nb_partitions = nb_partitions, Pix_Per_Partition = Pix_Per_Partition, + kernel = NULL) # if needed, apply continuum removal if (Continuum_Removal == TRUE) { Subset$DataSubset <- apply_continuum_removal(Subset$DataSubset, SpectralFilter, nbCPU = nbCPU) @@ -120,7 +128,7 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu # # # assume that 1st data cleaning is enough... ## Uncommented June 5, 2019 # clean reflectance data from inf and constant values - CleanData <- check_invariant_bands(DataSubset, SpectralFilter) + CleanData <- rm_invariant_bands(DataSubset, SpectralFilter) DataSubset <- CleanData$DataMatrix SpectralFilter <- CleanData$Spectral print("perform PCA#2 on the subset image") @@ -134,16 +142,12 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu } } # Number of PCs computed and written in the PCA file: 30 if hyperspectral - if (NbPCs_To_Keep==FALSE){ - Nb_PCs <- dim(PCA_model$dataPCA)[2] - if (Nb_PCs > 30) { - Nb_PCs <- 30 - } - } else { - Nb_PCs = NbPCs_To_Keep + Nb_PCs <- dim(PCA_model$x)[2] + if (Nb_PCs > NbPCs_To_Keep){ + Nb_PCs <- NbPCs_To_Keep } PCA_model$Nb_PCs <- Nb_PCs - PCA_model$dataPCA <- NULL + PCA_model$x <- NULL # CREATE PCA FILE CONTAINING ONLY SELECTED PCs print("Apply PCA model to the whole image") Output_Dir_PCA <- define_output_subdir(Output_Dir, Input_Image_File, TypePCA, "PCA") @@ -173,12 +177,14 @@ perform_PCA <- function(Input_Image_File, Input_Mask_File, Output_Dir, Continuu # @return Shade_Update = updated shade mask # ' @importFrom stats sd # ' @importFrom matlab ones -filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spectral, CR, PCA_model, PCsel, TypePCA, nbCPU = 1, MaxRAM = 0.25) { +filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, + Spectral, Continuum_Removal, PCA_model, PCsel, TypePCA, + nbCPU = 1, MaxRAM = 0.25) { # 1- get extreme values falling outside of mean +- 3SD for PCsel first components # compute mean and sd of the 5 first components of the sampled data - MeanSub <- colMeans(PCA_model$dataPCA) - SDSub <- apply(PCA_model$dataPCA, 2, sd) + MeanSub <- colMeans(PCA_model$x) + SDSub <- apply(PCA_model$x, 2, sd) MinPC <- MeanSub - 3.0 * SDSub MaxPC <- MeanSub + 3.0 * SDSub @@ -193,6 +199,7 @@ filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spe HDR_Shade$resolution <- NULL HDR_Shade$bandwidth <- NULL HDR_Shade$purpose <- NULL + HDR_Shade$`default stretch` <- '0 1 linear' HDR_Shade$`byte order` <- get_byte_order() headerFpath <- paste(Shade_Update, ".hdr", sep = "") write_ENVI_header(HDR_Shade, headerFpath) @@ -231,28 +238,24 @@ filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spe for (i in 1:nbPieces) { print(paste("PCA Piece #", i, "/", nbPieces)) # read image and mask data - Byte_Start <- SeqRead_Image$ReadByte_Start[i] - nbLinesIMG <- SeqRead_Image$Lines_Per_Chunk[i] - lenBin <- SeqRead_Image$ReadByte_End[i] - SeqRead_Image$ReadByte_Start[i] + 1 + nbLines <- SeqRead_Image$Lines_Per_Chunk[i] ImgFormat <- "2D" - Image_Chunk <- read_image_subset(Input_Image_File, HDR, Byte_Start, lenBin, nbLinesIMG, Image_Format, ImgFormat) - - Byte_Start <- SeqRead_Shade$ReadByte_Start[i] - nbLines <- SeqRead_Shade$Lines_Per_Chunk[i] - lenBin <- SeqRead_Shade$ReadByte_End[i] - SeqRead_Shade$ReadByte_Start[i] + 1 + Image_Chunk <- read_image_subset(ImPath = Input_Image_File, HDR = HDR, + Line_Start = SeqRead_Image$Line_Start[i],Lines_To_Read = nbLines, + ImgFormat = ImgFormat) ImgFormat <- "Shade" if ((!Input_Mask_File == FALSE) & (!Input_Mask_File == "")) { - Shade_Chunk <- read_image_subset(Input_Mask_File, HDR_Shade, Byte_Start, lenBin, nbLines, Shade_Format, ImgFormat) + Shade_Chunk <- read_image_subset(ImPath = Input_Mask_File, HDR = HDR_Shade, + Line_Start = SeqRead_Image$Line_Start[i],Lines_To_Read = nbLines, + ImgFormat = ImgFormat) } else { Shade_Chunk <- ones(nbLines * HDR$samples, 1) } - - elimShade <- which(Shade_Chunk == 0) keepShade <- which(Shade_Chunk == 1) - Image_Chunk <- Image_Chunk[keepShade, ] + Image_Chunk <- Image_Chunk[keepShade,] # apply Continuum removal if needed - if (CR == TRUE) { + if (Continuum_Removal == TRUE) { Image_Chunk <- apply_continuum_removal(Image_Chunk, Spectral, nbCPU = nbCPU) } else { if (!length(Spectral$WaterVapor) == 0) { @@ -265,11 +268,11 @@ filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spe } # Apply PCA if (TypePCA == "PCA" | TypePCA == "SPCA") { - Image_Chunk <- t(t(PCA_model$eiV[, PCsel]) %*% t(center_reduce(Image_Chunk, PCA_model$mu, PCA_model$scale))) + Image_Chunk <- scale(Image_Chunk, PCA_model$center, PCA_model$scale) %*% PCA_model$rotation[, PCsel] } # get PCA of the group of line and rearrange the data to write it correctly in the output file - linetmp <- matrix(NA, ncol = ncol(Image_Chunk), nrow = (HDR$samples * nbLinesIMG)) + linetmp <- matrix(NA, ncol = ncol(Image_Chunk), nrow = (HDR$samples * nbLines)) if (length(keepShade) > 0) { linetmp[keepShade, ] <- Image_Chunk } @@ -310,16 +313,16 @@ filter_PCA <- function(Input_Image_File, HDR, Input_Mask_File, Shade_Update, Spe # @param PCA_model PCA model description # @param Spectral spectral information to be used in the image # @param Nb_PCs number of components kept in the resulting PCA raster -# @param CR Shoudl continuum removal be performed? +# @param CR boolean. If TRUE continuum removal is performed. # @param TypePCA PCA, SPCA, NLPCA # @param nbCPU number of CPUs to process data # @param MaxRAM max RAM when initial image is read (in Gb) # # @return None -write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_model, Spectral, Nb_PCs, CR, TypePCA, nbCPU = 1, MaxRAM = 0.25) { +write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_model, Spectral, Nb_PCs, Continuum_Removal, TypePCA, nbCPU = 1, MaxRAM = 0.25) { ImPathHDR <- get_HDR_name(Input_Image_File) HDR <- read_ENVI_header(ImPathHDR) - if ((!Input_Mask_File == FALSE) & (!Input_Mask_File == "")) { + if (is.character(Input_Mask_File) && (Input_Mask_File != "")) { ShadeHDR <- get_HDR_name(Input_Mask_File) HDR_Shade <- read_ENVI_header(ShadeHDR) } else { @@ -329,12 +332,7 @@ write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_mo HDR_PCA <- HDR HDR_PCA$bands <- Nb_PCs HDR_PCA$`data type` <- 4 - PCs <- list() - for (i in 1:Nb_PCs) { - PCs <- c(PCs, paste("PC", i)) - } - PCs <- paste(PCs, collapse = ", ") - HDR_PCA$`band names` <- PCs + HDR_PCA$`band names` <- paste('PC', 1:Nb_PCs, collapse = ", ") HDR_PCA$wavelength <- NULL HDR_PCA$fwhm <- NULL HDR_PCA$resolution <- NULL @@ -375,27 +373,23 @@ write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_mo for (i in 1:nbPieces) { print(paste("PCA Piece #", i, "/", nbPieces)) # read image and mask data - Byte_Start <- SeqRead_Image$ReadByte_Start[i] nbLines <- SeqRead_Image$Lines_Per_Chunk[i] - lenBin <- SeqRead_Image$ReadByte_End[i] - SeqRead_Image$ReadByte_Start[i] + 1 ImgFormat <- "2D" - Image_Chunk <- read_image_subset(Input_Image_File, HDR, Byte_Start, lenBin, nbLines, Image_Format, ImgFormat) - - Byte_Start <- SeqRead_Shade$ReadByte_Start[i] - nbLines <- SeqRead_Shade$Lines_Per_Chunk[i] - lenBin <- SeqRead_Shade$ReadByte_End[i] - SeqRead_Shade$ReadByte_Start[i] + 1 - + Image_Chunk <- read_image_subset(ImPath = Input_Image_File, HDR = HDR, + Line_Start = SeqRead_Image$Line_Start[i],Lines_To_Read = nbLines, + ImgFormat = ImgFormat) if (typeof(HDR_Shade) == 'list') { ImgFormat <- "Shade" - Shade_Chunk <- read_image_subset(Input_Mask_File, HDR_Shade, Byte_Start, lenBin, nbLines, Shade_Format, ImgFormat) - elimShade <- which(Shade_Chunk == 0) + Shade_Chunk <- read_image_subset(ImPath = Input_Mask_File, HDR = HDR_Shade, + Line_Start = SeqRead_Image$Line_Start[i],Lines_To_Read = nbLines, + ImgFormat = ImgFormat) keepShade <- which(Shade_Chunk == 1) Image_Chunk <- Image_Chunk[keepShade, ] } else { keepShade <- matrix(1,ncol = 1,nrow = nrow(Image_Chunk)) } # apply Continuum removal if needed - if (CR == TRUE) { + if (Continuum_Removal) { Image_Chunk <- apply_continuum_removal(Image_Chunk, Spectral, nbCPU = nbCPU) ## added June 5, 2019 if (length(Spectral$BandsNoVar) > 0) { @@ -412,7 +406,7 @@ write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_mo # Apply PCA if (TypePCA == "PCA" | TypePCA == "SPCA") { - Image_Chunk <- t(t(PCA_model$eiV[, 1:Nb_PCs]) %*% t(center_reduce(Image_Chunk, PCA_model$mu, PCA_model$scale))) + Image_Chunk <- scale(Image_Chunk, PCA_model$center, PCA_model$scale) %*% PCA_model$rotation[, 1:Nb_PCs] } # get PCA of the group of line and rearrange the data to write it correctly in the output file @@ -426,14 +420,11 @@ write_PCA_raster <- function(Input_Image_File, Input_Mask_File, PCA_Path, PCA_mo encoding = getOption("encoding"), raw = FALSE ) if (!SeqRead_PCA$ReadByte_Start[i] == 1) { - # nbSkip = (as.double(SeqRead_PCA$ReadByte_Start[i])*as.double(PCA_Format$Bytes))-1 - # nbSkip = (as.double(SeqRead_PCA$ReadByte_Start[i])*as.double(PCA_Format$Bytes))-as.double(1) nbSkip <- (SeqRead_PCA$ReadByte_Start[i] - 1) * PCA_Format$Bytes seek(fidOUT, where = nbSkip, origin = "start", rw = "write") } PCA_Chunk <- array(PCA_Chunk, c(nbLines, HDR_PCA$samples, HDR_PCA$bands)) PCA_Chunk <- aperm(PCA_Chunk, c(2, 3, 1)) - # writeBin(as.numeric(c(PCA_Chunk)), fidOUT, size = PCA_Format$Bytes,endian = .Platform$endian) writeBin(c(PCA_Chunk), fidOUT, size = PCA_Format$Bytes, endian = .Platform$endian, useBytes = FALSE) close(fidOUT) } @@ -454,13 +445,13 @@ pca <- function(X, type) { p <- ncol(X) if (type == "SPCA") { modPCA <- prcomp(X, scale = TRUE) - sig <- modPCA$scale + } else if (type == "PCA") { modPCA <- prcomp(X, scale = FALSE) - sig <- rep(1, p) + } - my_list <- list("dataPCA" = modPCA$x, "mu" = modPCA$center, "scale" = sig, "eiV" = modPCA$rotation, "eig" = modPCA$sdev) - return(my_list) + + return(modPCA) } # defines the number of pixels per iteration based on a trade-off between image size and sample size per iteration @@ -482,7 +473,7 @@ define_pixels_per_iter <- function(ImNames, nb_partitions = nb_partitions) { lenTot <- nbPixels * as.double(HDR$bands) ImSizeGb <- (lenTot * Image_Format$Bytes) / (1024^3) # if shade mask, update number of pixels - if ((!Input_Mask_File == FALSE) & (!Input_Mask_File == "")) { + if (is.character(Input_Mask_File) && (Input_Mask_File != "")) { # read shade mask fid <- file( description = Input_Mask_File, open = "rb", blocking = TRUE, diff --git a/man/Write_Image_NativeRes.Rd b/man/Write_Image_NativeRes.Rd new file mode 100644 index 00000000..d9ab9858 --- /dev/null +++ b/man/Write_Image_NativeRes.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Lib_ImageProcess.R +\name{Write_Image_NativeRes} +\alias{Write_Image_NativeRes} +\title{write an image resulting from "window processing" at native spatial resolution +(assuming square windows & origin at top left corner)} +\usage{ +Write_Image_NativeRes(Image, ImagePath, HDR, window_size) +} +\arguments{ +\item{Image}{numeric. Image corresponding to a 2D matrix or 3D array} + +\item{ImagePath}{character. Path where the image should be written} + +\item{HDR}{list. Image header} + +\item{window_size}{numeric. window size used to generate Image} +} +\value{ +None +} +\description{ +write an image resulting from "window processing" at native spatial resolution +(assuming square windows & origin at top left corner) +} diff --git a/man/perform_PCA.Rd b/man/perform_PCA.Rd index 01dc501a..1db2a45b 100644 --- a/man/perform_PCA.Rd +++ b/man/perform_PCA.Rd @@ -10,7 +10,7 @@ perform_PCA( Output_Dir, Continuum_Removal = TRUE, TypePCA = "SPCA", - NbPCs_To_Keep = FALSE, + NbPCs_To_Keep = 30, FilterPCA = FALSE, Excluded_WL = FALSE, nb_partitions = 20, diff --git a/vignettes/biodivMapR.Rmd b/vignettes/biodivMapR.Rmd index 84fea31a..14863526 100644 --- a/vignettes/biodivMapR.Rmd +++ b/vignettes/biodivMapR.Rmd @@ -43,29 +43,23 @@ Below is the typical flow chart of the computation of diversity maps with __biod # Define processing parameters ## Input / Output files -The input images are expected to be in ENVI HDR format, BIL interleaved. -The functions `perform_radiometric_filtering` and `perform_PCA` start with a procedure checkin if the image format is as expected. if not, the functions will retrun a message explaining the problem and will stop the process. +

**__biodivMapR__ uses the __stars__ package since version 1.3.0**

+

**This means that __biodivMapR__ now handles various input raster file formats such as tiff !!**

+

**However, you will need a ENVI style header file in order to provide spectral bands**

-If the format is not ENVI format BIL interleaved, the function `raster2BIL` allows conversion into proper format and returns updated `Input_Image_File` +The function `raster2BIL` still allows conversion into ENVI HDR format, BIL interleaved and returns updated `Input_Image_File`. Please use this function if you run into errors when processing your image. The error may be due to improper information in the header file. -[__Please read the dedicated explanation if you need to convert your image file, and keep in mind that you may need to adjust existing template `.hdr` files if you want to run `biodivMapR` on images different from 'standard' BOA Sentinel-2 data with 10 spectral bands (including B02, B03, B04, B05, B06, B07, B08, B08A, B11 and B12).__](https://jbferet.github.io/biodivMapR/articles/biodivMapR_1.html "Converting a raster file to the proper file format") +**Spectral bands should be defined if the image is multi or hyperspectral image.** -Spectral bands should be defined if the image is multi or hyperspectral image. +A mask can also be set to work on a selected part of the input image. The mask is expected to be a raster with values 0 = masked or 1 = selected. Only one band is required. If no mask is to be used set `Input_Mask_File = FALSE`. -A mask can also be set to work on a selected part of the input image. The mask is expected to be a raster in the same format as the image ([__ENVI file format__](https://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html "HDR description") with corresponding '.hdr' file), with values 0 = masked or 1 = selected. Only one band is required. If no mask is to be used set `Input_Mask_File = FALSE`. - -The output directory defined with `Output_Dir` will contain all the results. For each image processed, a subdirectory will be automatically created after its name. +The output directory defined with `Output_Dir` will contain all the results. For each image processed, a subdirectory will be automatically created after its name. The output images are writted in ENVI HDR format, BIL interleaved. ```{r Input / Output files} library(biodivMapR) Input_Image_File = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'biodivMapR') -# please visit this webpage if you need to convert your raster image file using raster2BIL -# https://jbferet.github.io/biodivMapR/articles/biodivMapR_1.html -# Input_Image_File = raster2BIL(Raster_Path = Input_Image_File, -# Sensor = 'SENTINEL_2A') - Input_Mask_File = FALSE Output_Dir = 'RESULTS' diff --git a/vignettes/biodivMapR_1.Rmd b/vignettes/biodivMapR_1.Rmd index ef746c0d..10cd9ce2 100644 --- a/vignettes/biodivMapR_1.Rmd +++ b/vignettes/biodivMapR_1.Rmd @@ -19,13 +19,14 @@ knitr::opts_chunk$set( ) ``` +

**if you are running version 1.3.0 or more recent version of __biodivMapR__ this tutorial is not needed**

-__biodivMapR__ expects a specific raster file format as input file, as the package includes dedicated readers and writers. +__biodivMapR__ was expecting a specific raster file format as input file before version 1.3.0. -* Raster files are expected to be raw flat-binary files. They should not contain any header, which are stored in a separate .hdr file. Hence format such as TIFF are not accepted. +* Raster files were expected to be raw flat-binary files. They should not contain any header, which are stored in a separate .hdr file. Hence format such as TIFF are not accepted. * This binary file should be written following a [Band-interleaved-by-line storage format](https://www.harrisgeospatial.com/docs/enviimagefiles.html "BIL definition"). * This file should have either `.bil` file extension, or no file extension. -* The header file should be named the same as the binary image file, but with `.hdr` file extension, and should be located in teh same directory. +* The header file should be named the same as the binary image file, but with `.hdr` file extension, and should be located in the same directory. * the header file should be an ASCII file and follow the [__ENVI file format__](https://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html "HDR description"). A function is dedicated to conversion of a raster file into appropriate image format, named `raster2BIL`. This function is relatively straightforward to run, but users should make sure prerequisites are met before running the function.