- integrated stars package in order to read any file format, includin…

…g 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
jbferet · Jun 1, 2020 · 8350b41 · 8350b41
1 parent 594e8e5
commit 8350b41
Show file tree

Hide file tree

Showing 15 changed files with 556 additions and 440 deletions.
diff --git 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
@@ -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
@@ -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
@@ -4,6 +4,7 @@
 # ==============================================================================
 # PROGRAMMERS:
 # Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Florian de Boissieu <fdeboiss@gmail.com>
 # 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,47 +185,28 @@ 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()
 	  }
   } else {
     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
@@ -4,6 +4,7 @@
 # ===============================================================================
 # PROGRAMMERS:
 # Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Florian de Boissieu <fdeboiss@gmail.com>
 # 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
@@ -4,6 +4,7 @@
 # ==============================================================================
 # PROGRAMMERS:
 # Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Florian de Boissieu <fdeboiss@gmail.com>
 # Copyright 2018/07 Jean-Baptiste FERET
 # ==============================================================================
 # This library contains functions to filter raster based on radiometric criteria
@@ -24,32 +25,26 @@
 #'
 #' @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){
     check_data(Mask_Path,Mask = TRUE)
   }
   # 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()