Lib_FilterData.R 4.97 KB
Newer Older
jbferet's avatar
jbferet committed
1
# ==============================================================================
Florian de Boissieu's avatar
Florian de Boissieu committed
2
# biodivMapR
jbferet's avatar
jbferet committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Lib_FilterData.R
# ==============================================================================
# PROGRAMMERS:
# Jean-Baptiste FERET <jb.feret@irstea.fr>
# Copyright 2018/07 Jean-Baptiste FERET
# ==============================================================================
# This library contains functions to filter raster based on radiometric criteria
# ==============================================================================

#' Performs radiometric filtering based on three criteria: NDVI, NIR reflectance, Blue reflectance
#'
#' @param Image.Path Path of the image to be processed
#' @param Mask.Path Path of the mask corresponding to the image
#' @param Output.Dir output directory
De Boissieu Florian's avatar
De Boissieu Florian committed
17
#' @param TypePCA Type of PCA: "PCA" or "SPCA"
jbferet's avatar
jbferet committed
18
19
20
21
22
23
24
#' @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
#' @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
#' @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
#' @param Mask.Path
#'
#' @return ImPathShade = updated shademask file
#' @export
25
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) {
jbferet's avatar
jbferet committed
26
  # define full output directory
27
  Output.Dir.Full <- Define.Output.Dir(Output.Dir, Image.Path, TypePCA)
jbferet's avatar
jbferet committed
28
  # define dimensions of the image
29
30
31
32
33
34
35
36
  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)
  Nb.Pixels <- ipix * jpix
  lenTot <- Nb.Pixels * as.double(HDR$bands)
  ImSizeGb <- (lenTot * Image.Format$Bytes) / (1024^3)
jbferet's avatar
jbferet committed
37
38

  # Create / Update shade mask if optical data
39
40
  if (Mask.Path == FALSE | Mask.Path == "") {
    print("Create mask based on NDVI, NIR and Blue threshold")
jbferet's avatar
jbferet committed
41
  } else {
42
    print("Update mask based on NDVI, NIR and Blue threshold")
jbferet's avatar
jbferet committed
43
  }
44
45
  Shade.Update <- paste(Output.Dir.Full, "ShadeMask_Update", sep = "")
  Mask.Path <- Create.Mask.Optical(Image.Path, Mask.Path, Shade.Update, NDVI.Thresh, Blue.Thresh, NIR.Thresh, Blue, Red, NIR)
jbferet's avatar
jbferet committed
46
47
48
  return(Mask.Path)
}

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
# create a mask based on NDVI, Green reflectance and NIR reflectance
# NDVI (min) threshold eliminates non vegetated pixels
# Blue (max) threshold eliminates Clouds
# NIR (min) threshold eliminates shadows
# ! only valid if Optical data!!
#
# @param ImPath full path of a raster file
# @param ImPathShade full path of the raster mask corresponding to the raster file
# @param ImPathShade.Update wavelength (nm) of the spectral bands to be found
# @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
# @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
# @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
#
# @return ImPathShade path for the updated shademask produced
Create.Mask.Optical <- function(ImPath, ImPathShade, ImPathShade.Update, NDVI.Thresh, Blue.Thresh, NIR.Thresh, Blue = 480, Red = 700, NIR = 835) {
jbferet's avatar
jbferet committed
64
  # define wavelength corresponding to the spectral domains Blue, Red and NIR
65
66
67
  Spectral.Bands <- c(Blue, Red, NIR)
  ImPathHDR <- Get.HDR.Name(ImPath)
  Header <- read.ENVI.header(ImPathHDR)
jbferet's avatar
jbferet committed
68
  # get image bands correponding to spectral bands of interest
69
  Image.Bands <- Get.Image.Bands(Spectral.Bands, Header$wavelength)
jbferet's avatar
jbferet committed
70
  # read band data from image
71
  Image.Subset <- Read.Image.Bands(ImPath, Header, Image.Bands$ImBand)
jbferet's avatar
jbferet committed
72
73
  # create mask
  # check if spectral bands required for NDVI exist
74
75
  if (Image.Bands$Distance2WL[2] < 25 & Image.Bands$Distance2WL[3] < 25) {
    NDVI <- ((Image.Subset[, , 3]) - (Image.Subset[, , 2])) / ((Image.Subset[, , 3]) + (Image.Subset[, , 2]))
jbferet's avatar
jbferet committed
76
  } else {
77
78
    NDVI <- matrix(1, nrow = Header$lines, ncol = Header$samples)
    message("Could not find the spectral bands required to compute NDVI")
jbferet's avatar
jbferet committed
79
  }
80
81
82
  if (Image.Bands$Distance2WL[1] > 25) {
    Image.Subset[, , 1] <- Blue.Thresh + 0 * Image.Subset[, , 1]
    message("Could not find a spectral band in the blue domain: will not perform filtering based on blue reflectance")
jbferet's avatar
jbferet committed
83
  }
84
85
86
  if (Image.Bands$Distance2WL[3] > 50) {
    Image.Subset[, , 3] <- NIR.Thresh + 0 * Image.Subset[, , 3]
    message("Could not find a spectral band in the NIR domain: will not perform filtering based on NIR reflectance")
jbferet's avatar
jbferet committed
87
  }
88
89
90
  Mask <- matrix(0, nrow = Header$lines, ncol = Header$samples)
  SelPixels <- which(NDVI > NDVI.Thresh & Image.Subset[, , 1] < Blue.Thresh & Image.Subset[, , 3] > NIR.Thresh)
  Mask[SelPixels] <- 1
jbferet's avatar
jbferet committed
91
  # update initial shade mask
92
  ImPathShade <- Update.Shademask(ImPathShade, Header, Mask, ImPathShade.Update)
jbferet's avatar
jbferet committed
93
94
  return(ImPathShade)
}