Merge branch 'dev'

jbferet · Apr 27, 2020 · 954e0c8 · 954e0c8
2 parents 6f405e6 + 386fb00
commit 954e0c8
Show file tree

Hide file tree

Showing 10 changed files with 187 additions and 51 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: biodivMapR
 Title: biodivMapR: an R package for α- and β-diversity mapping using remotely-sensed images
-Version: 1.1.0
+Version: 1.2.0
 Authors@R: c(person(given = "Jean-Baptiste", 
 					family = "Feret", 
 					email = "jb.feret@teledetection.fr", 

diff --git a/NAMESPACE b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(Write_Big_Image)
 export(check_data)
 export(diversity_from_plots)
 export(list_shp)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,13 @@
+# biodivMapR v1.2.0
+
+## Fixes
+- fixed bug raised when processing data over large number of pixels (image products >2^31 bits)
+- fixed bug by adding elimination of bands disturbed by water vapor even when ContinuumRemoval set to FALSE
+
+## Changes
+- Changed name of default directory when saving image after calling raster2BIL
+- added documentation for raster data conversion using raster2BIL
+
 # biodivMapR v1.1.0
 
 ## Fixes

diff --git a/R/Lib_CheckConvertData.R b/R/Lib_CheckConvertData.R
@@ -145,7 +145,7 @@ raster2BIL <- function(Raster_Path, Sensor = "unknown", Output_Dir = FALSE, Conv
       message("or manually add wavelength location in HDR file, if relevant")
     }
   } else if (Sensor == "unknown") {
-    message("please make sure that the follozing header file contains information required")
+    message("please make sure that the following header file contains information required")
     print(get_HDR_name(Output_Path))
     message("or manually add wavelength location in HDR file, if relevant")
   }

diff --git a/R/Lib_ImageProcess.R b/R/Lib_ImageProcess.R
@@ -870,7 +870,7 @@ where_to_read <- function(HDR, nbPieces) {
   ReadByte_Start <- lb[1:nbPieces]
   ReadByte_End <- ub[2:(nbPieces + 1)]
   Lines_Per_Chunk <- diff(Start_Per_Chunk)
-  my_list <- list("ReadByte_Start" = ReadByte_Start, "ReadByte_End" = ReadByte_End, "Lines_Per_Chunk" = Lines_Per_Chunk)
+  my_list <- list("ReadByte_Start" = ReadByte_Start, "ReadByte_End" = ReadByte_End, "Lines_Per_Chunk" = Lines_Per_Chunk,'Line_Start'=Start_Per_Chunk)
   return(my_list)
 }
 
@@ -942,6 +942,57 @@ write_ENVI_header <- function(HDR, HDRpath) {
   writeLines(c("ENVI", paste(names(HDR), h, sep = " = ")), con = HDRpath)
 }
 
+#' write an image which size is > 2**31-1
+#'
+#' @param ImgWrite numeric. Image as array
+#' @param ImagePath character. Path where the image should be written
+#' @param HDR list. Image header
+#' @param Image_Format list. description of data format corresponding to ENVI type
+#'
+#' @return None
+#' @export
+
+Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){
+  nbPieces <- split_image(HDR, LimitSizeGb = 1.8)
+  SeqRead_Image <- where_to_read(HDR, nbPieces)
+  fidOUT <- file(
+    description = ImagePath, open = "wb", blocking = TRUE,
+    encoding = getOption("encoding"), raw = FALSE
+  )
+  close(fidOUT)
+  # for each piece of image
+  for (i in 1:nbPieces) {
+    print(paste("Writing Image, piece #", i, "/", nbPieces))
+    # read image and mask data
+    Byte_Start <- SeqRead_Image$ReadByte_Start[i]
+    Line_Start <- SeqRead_Image$Line_Start[i]
+    nbLines <- SeqRead_Image$Lines_Per_Chunk[i]
+    lenBin <- SeqRead_Image$ReadByte_End[i] - SeqRead_Image$ReadByte_Start[i] + 1
+    # files to write in
+    fidOUT <- file(
+      description = ImagePath, open = "r+b", blocking = TRUE,
+      encoding = getOption("encoding"), raw = FALSE
+    )
+    if (!SeqRead_Image$ReadByte_Start[i] == 1) {
+      nbSkip <- (SeqRead_Image$ReadByte_Start[i] - 1) * Image_Format$Bytes
+      seek(fidOUT, where = nbSkip, origin = "start", rw = "write")
+    }
+    if (is.na(dim(ImgWrite)[3])){
+      ImgChunk <- array(ImgWrite[Line_Start:(nbLines+Line_Start-1),], c(nbLines, HDR$samples, HDR$bands))
+    } else {
+      ImgChunk <- array(ImgWrite[Line_Start:(nbLines+Line_Start-1),,], c(nbLines, HDR$samples, HDR$bands))
+    }
+    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)
+    close(fidOUT)
+  }
+  rm(ImgWrite)
+  rm(ImgChunk)
+  gc()
+  return("")
+}
+
 # convert image coordinates from X-Y to index
 #
 # @param HDR_Raster

diff --git a/R/Lib_MapAlphaDiversity.R b/R/Lib_MapAlphaDiversity.R
@@ -401,16 +401,11 @@ write_raster_alpha <- function(Image, HDR_SSD, ImagePath, window_size, Index, Fu
   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)
-    ImgWrite <- array(Image, c(HDR_Alpha$lines, HDR_Alpha$samples, 1))
-    ImgWrite <- aperm(ImgWrite, 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)
+    # 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
@@ -428,14 +423,8 @@ write_raster_alpha <- function(Image, HDR_SSD, ImagePath, window_size, Index, Fu
         Image_FullRes[((i - 1) * window_size + 1):(i * window_size), ((j - 1) * window_size + 1):(j * window_size)] <- Image[i, j]
       }
     }
-    ImgWrite <- array(Image_FullRes, c(HDR_Full$lines, HDR_Full$samples, 1))
-    ImgWrite <- aperm(ImgWrite, 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)
+    # 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)
   }
@@ -451,14 +440,9 @@ write_raster_alpha <- function(Image, HDR_SSD, ImagePath, window_size, Index, Fu
     Image_Format <- ENVI_type2bytes(HDR_Alpha)
     if (LowRes == TRUE) {
       write_ENVI_header(HDR_Alpha, headerFpath)
-      ImgWrite <- array(Image_Smooth, c(HDR_Alpha$lines, HDR_Alpha$samples, 1))
-      ImgWrite <- aperm(ImgWrite, c(2, 3, 1))
-      fidOUT <- file(
-        description = ImagePath.Smooth, open = "wb", blocking = TRUE,
-        encoding = getOption("encoding"), raw = FALSE
-      )
-      writeBin(c(ImgWrite), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE)
-      close(fidOUT)
+      # 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
@@ -472,14 +456,8 @@ write_raster_alpha <- function(Image, HDR_SSD, ImagePath, window_size, Index, Fu
           Image_FullRes[((i - 1) * window_size + 1):(i * window_size), ((j - 1) * window_size + 1):(j * window_size)] <- Image_Smooth[i, j]
         }
       }
-      ImgWrite <- array(Image_FullRes, c(HDR_Full$lines, HDR_Full$samples, 1))
-      ImgWrite <- aperm(ImgWrite, 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)
+      # 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)
     }

diff --git a/R/Lib_MapBetaDiversity.R b/R/Lib_MapBetaDiversity.R
@@ -49,7 +49,6 @@ map_beta_div <- function(Input_Image_File, Output_Dir, window_size,
   return(invisible())
 }
 
-
 # computes NMDS
 #
 # @param MatBCdist BC dissimilarity matrix
@@ -337,13 +336,15 @@ write_raster_beta <- function(Image, HDR_SSD, ImagePath, window_size, FullRes =
   if (LowRes == TRUE) {
     headerFpath <- paste(ImagePath, ".hdr", sep = "")
     write_ENVI_header(HDR_Beta, headerFpath)
-    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)
+    # 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
@@ -363,13 +364,15 @@ write_raster_beta <- function(Image, HDR_SSD, ImagePath, window_size, FullRes =
         }
       }
     }
-    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)
+    # 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)
   }

diff --git a/man/Write_Big_Image.Rd b/man/Write_Big_Image.Rd
diff --git a/vignettes/biodivMapR.Rmd b/vignettes/biodivMapR.Rmd
@@ -48,9 +48,11 @@ The functions `perform_radiometric_filtering` and `perform_PCA` start with a pro
 
 If the format is not ENVI format BIL interleaved, the function `raster2BIL` allows conversion into proper format and returns updated `Input_Image_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. 
 
-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 HDR), 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.
 

diff --git a/vignettes/biodivMapR_1.Rmd b/vignettes/biodivMapR_1.Rmd
@@ -0,0 +1,68 @@
+---
+title: "Converting a raster file to the proper file format"
+author: "Jean-Baptiste Féret, Florian de Boissieu"
+date: "`r Sys.Date()`"
+output: 
+  html_vignette:
+    number_sections: true
+vignette: >
+  %\VignetteIndexEntry{Tutorial}
+  %\VignetteEngine{knitr::rmarkdown}
+  \usepackage[utf8]{inputenc}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  eval=FALSE
+)
+```
+
+
+__biodivMapR__ expects a specific raster file format as input file, as the package includes dedicated readers and writers. 
+
+* 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. 
+* 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 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.
+
+# How to produce raster data with proper format using `raster2BIL`
+
+`raster2BIL` can be run as follows: 
+
+```{r raster2BIL}
+library(biodivMapR)
+
+Input_Image_File  = raster2BIL(Raster_Path = Path_Raster_File,
+                                       Sensor = Path_Template_HDR,
+                                       Convert_Integer = TRUE,
+                                       Output_Dir = Path_Output_Directory)
+
+```
+
+## Input variables
+
+The following information is expected as input information:  
+
+* `Raster_Path`: a string containing the full path for the raster to be converted.
+* `Sensor`: a string containing the name of the sensor, which should actually be linked to a template for the corresponding header file. 
+* `Convert_Integer` [optional]: a boolean stating if the raster data should be converted into integer
+* `Multiplying_Factor` [optional]:  a multiplying factor for the raster data. This is useful when converting real reflectance values between 0 and 1 into integer between 0 and 10000. Set `Multiplying_Factor =10000` in this case.
+* `Output_Dir` [optional]: a string defining the path for the output directory where converted file will be stored.
+* `Multiplying_Factor_Last` [optional]: a multiplying factor for last band. This option is useful if working with dry vegetation acquired with Sentinel-2 images, in order to artifically increase the reflectance of the B12 band (SWIR 2, centered at 2200 nm), so that SWIR 1 band remains below the convex hull defined during Continuum removal. 
+
+Please note that 
+i) the coordinate system and minimal metadata (spatial resolution, geographic information...) are expected to be provided in the original image or image + header defined in `Raster_Path` if user want them to be included in the final image products.
+
+ii) `Sensor` should refer to the name of a `.hdr` file stored in the installation directory of `biodivMapR`
+Please type `system.file(package = "biodivMapR")` in order to get the location of the install directory and create a proper `.hdr` file to be stored in `extdata/HDR/`, containing proper spectral bands __defined in nanometers__ if using optical data and using `perform_radiometric_filtering` and/or `Continuum_Removal = TRUE` in your process. Make sure that the __spectral bands in the binary file are stored with wavelengths following ascending order__. 
+
+iii) `Convert_Integer=TRUE` saves 50% space if original image is stored in real values (32bits per value). However, user need to make sure hat their original data are compatible with conversion: if reflectance values are stored as  real value between 0 and 1, the final image file will be unusable, unless proper multiplying factor is applied to the full image data (for example `Multiplying_Factor =10000`). Sentinel-2 reflectance data downloaded from ESA hub are stored as integer (16 bits) values between 0 and 10 000 instead of real (32 bits) values between 0 and 1, so integer conversion is unnecessary.
+
+## output variable
+
+`raster2BIL` returns a string containing the full path for the raster once converted.