Lib_MapAlphaDiversity.R 23.4 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
# Lib_MapAlphaDiversity.R
# ==============================================================================
# PROGRAMMERS:
# Jean-Baptiste FERET <jb.feret@irstea.fr>
# Copyright 2018/07 Jean-Baptiste FERET
# ==============================================================================
# This Library produces maps of alpha diversity indicators (Shannon, Simpson,
# Fischer...) based on spectral species file
# ==============================================================================

#' maps alpha diversity indicators  based on prior selection of PCs
#'
15
16
17
18
19
20
21
22
23
24
25
26
#' @param Input.Image.File character. Path and name of the image to be processed.
#' @param Output.Dir character. Output directory.
#' @param Spatial.Unit numeric. Size of spatial units (in pixels) to compute diversity.
#' @param TypePCA character. Type of PCA (PCA, SPCA, NLPCA...).
#' @param nbclusters numeric. Number of clusters defined in k-Means.
#' @param MinSun numeric. Minimum proportion of sunlit pixels required to consider plot.
#' @param pcelim numeric. Minimum contribution (in \%) required for a spectral species.
#' @param FullRes boolean. Full resolution.
#' @param LowRes boolean. Low resolution.
#' @param nbCPU numeric. Number of CPUs to use in parallel.
#' @param MaxRAM numeric. MaxRAM maximum size of chunk in GB to limit RAM allocation when reading image file.
#' @param Index.Alpha character. Either 'Shannon', 'Simpson' or 'Fisher'.
jbferet's avatar
jbferet committed
27
28
#'
#' @export
De Boissieu Florian's avatar
De Boissieu Florian committed
29
alpha_div <- function(Input.Image.File, Output.Dir, Spatial.Unit,
30
31
32
33
34
35
36
37
                                TypePCA = "SPCA", nbclusters = 50,
                                MinSun = 0.25, pcelim = 0.02,
                                Index.Alpha = "Shannon", FullRes = TRUE,
                                LowRes = FALSE, MapSTD = FALSE,
                                nbCPU = FALSE, MaxRAM = FALSE) {
  Output.Dir.SS <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "SpectralSpecies")
  Output.Dir.PCA <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "PCA")
  Spectral.Species.Path <- paste(Output.Dir.SS, "SpectralSpecies", sep = "")
jbferet's avatar
jbferet committed
38
  # 1- COMPUTE ALPHA DIVERSITY
39
  ALPHA <- Compute.Alpha.Diversity(Spectral.Species.Path, Spatial.Unit, nbclusters, MinSun, pcelim, nbCPU = nbCPU, MaxRAM = MaxRAM, Index.Alpha = Index.Alpha)
jbferet's avatar
jbferet committed
40
  # 2- SAVE ALPHA DIVERSITY MAPS
41
  print("Write alpha diversity maps")
jbferet's avatar
jbferet committed
42
  # which spectral indices will be computed
43
44
45
46
  Shannon <- Simpson <- Fisher <- FALSE
  if (length((grep("Shannon", Index.Alpha))) > 0) Shannon <- TRUE
  if (length((grep("Simpson", Index.Alpha))) > 0) Simpson <- TRUE
  if (length((grep("Fisher", Index.Alpha))) > 0) Fisher <- TRUE
jbferet's avatar
jbferet committed
47

48
49
50
51
52
53
54
55
56
  Output.Dir.Alpha <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "ALPHA")
  if (Shannon == TRUE) {
    Index <- "Shannon"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Shannon, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    if (MapSTD == TRUE) {
      Index <- "Shannon_SD"
      Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
      Write.Image.Alpha(ALPHA$Shannon.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
57
58
59
    }
  }

60
61
62
63
64
65
66
67
  if (Fisher == TRUE) {
    Index <- "Fisher"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Fisher, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    if (MapSTD == TRUE) {
      Index <- "Fisher_SD"
      Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
      Write.Image.Alpha(ALPHA$Fisher.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
68
69
70
    }
  }

71
72
73
74
75
76
77
78
  if (Simpson == TRUE) {
    Index <- "Simpson"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Simpson, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    if (MapSTD == TRUE) {
      Index <- "Simpson_SD"
      Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
      Write.Image.Alpha(ALPHA$Simpson.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
79
80
81
82
83
    }
  }
  return()
}

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# maps alpha diversity indicators  based on prior selection of PCs
#
# @param ImNames Path and name of the images to be processed
# @param Output.Dir output directory
# @param Spatial.Unit dimensions of the spatial unit
# @param TypePCA Type of PCA (PCA, SPCA, NLPCA...)
# @param nbclusters number of clusters defined in k-Means
# @param MinSun minimum proportion of sunlit pixels required to consider plot
# @param pcelim minimum contribution (in %) required for a spectral species
# @param FullRes
# @param LowRes
# @param nbCPU
# @param MaxRAM
# @param Index.Alpha 'Shannon', 'Simpson, 'Fisher'
#
# @return
Map.Alpha.Diversity.TestnbCluster <- function(ImNames, Output.Dir, Spatial.Unit, TypePCA = "SPCA", nbclusters = 50, MinSun = 0.25, pcelim = 0.02, Index.Alpha = "Shannon", FullRes = TRUE, LowRes = FALSE, nbCPU = FALSE, MaxRAM = FALSE) {
  Output.Dir.SS <- Define.Output.SubDir(Output.Dir, ImNames$Input.Image, TypePCA, paste("SpectralSpecies_", nbclusters, sep = ""))
  Output.Dir.PCA <- Define.Output.SubDir(Output.Dir, ImNames$Input.Image, TypePCA, "PCA")
  Spectral.Species.Path <- paste(Output.Dir.SS, "SpectralSpecies", sep = "")
jbferet's avatar
jbferet committed
104
  # 1- COMPUTE ALPHA DIVERSITY
105
  ALPHA <- Compute.Alpha.Diversity(Spectral.Species.Path, Spatial.Unit, nbclusters, MinSun, pcelim, nbCPU = nbCPU, MaxRAM = MaxRAM, Index.Alpha = Index.Alpha)
jbferet's avatar
jbferet committed
106
  # 2- SAVE ALPHA DIVERSITY MAPS
107
  print("Write alpha diversity maps")
jbferet's avatar
jbferet committed
108
  # which spectral indices will be computed
109
110
111
112
  Shannon <- Simpson <- Fisher <- FALSE
  if (length((grep("Shannon", Index.Alpha))) > 0) Shannon <- TRUE
  if (length((grep("Simpson", Index.Alpha))) > 0) Simpson <- TRUE
  if (length((grep("Fisher", Index.Alpha))) > 0) Fisher <- TRUE
jbferet's avatar
jbferet committed
113

114
115
116
117
118
119
120
121
  Output.Dir.Alpha <- Define.Output.SubDir(Output.Dir, ImNames$Input.Image, TypePCA, paste("ALPHA_", nbclusters, sep = ""))
  if (Shannon == TRUE) {
    Index <- "Shannon"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Shannon, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    Index <- "Shannon_SD"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Shannon.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
122
123
  }

124
125
126
127
128
129
130
  if (Fisher == TRUE) {
    Index <- "Fisher"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Fisher, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    Index <- "Fisher_SD"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Fisher.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
131
132
  }

133
134
135
136
137
138
139
  if (Simpson == TRUE) {
    Index <- "Simpson"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Simpson, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
    Index <- "Simpson_SD"
    Alpha.Path <- paste(Output.Dir.Alpha, Index, "_", Spatial.Unit, sep = "")
    Write.Image.Alpha(ALPHA$Simpson.SD, ALPHA$HDR, Alpha.Path, Spatial.Unit, Index, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
140
141
142
143
  }
  return()
}

144
145
146
147
148
149
150
151
152
153
154
155
# Map alpha diversity metrics based on spectral species
#
# @param Spectral.Species.Path path for spectral species file to be written
# @param Spatial.Unit 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
#
# @return list of mean and SD of alpha diversity metrics
De Boissieu Florian's avatar
De Boissieu Florian committed
156
#' @importFrom future plan multiprocess sequential
157
158
#' @importFrom future.apply future_lapply
Compute.Alpha.Diversity <- function(Spectral.Species.Path, Spatial.Unit, nbclusters, MinSun, pcelim, nbCPU = FALSE, MaxRAM = FALSE, Index.Alpha = "Shannon") {
jbferet's avatar
jbferet committed
159
  ##      read SpectralSpecies file and write distribution per spatial unit   ##
160
161
162
163
  SS.HDR <- Get.HDR.Name(Spectral.Species.Path)
  HDR.SS <- read.ENVI.header(SS.HDR)
  if (MaxRAM == FALSE) {
    MaxRAM <- 0.25
jbferet's avatar
jbferet committed
164
  }
165
166
167
  nbPieces.Min <- Split.Image(HDR.SS, MaxRAM)
  if (nbCPU == FALSE) {
    nbCPU <- detectCores()
jbferet's avatar
jbferet committed
168
  }
169
170
  if (nbPieces.Min < nbCPU) {
    nbPieces.Min <- nbCPU
jbferet's avatar
jbferet committed
171
  }
172
  SeqRead.SS <- Where.To.Read.Kernel(HDR.SS, nbPieces.Min, Spatial.Unit)
jbferet's avatar
jbferet committed
173
174
175

  ##          prepare SS distribution map and corresponding sunlit map        ##
  # prepare SS distribution map
176
177
  SSD.Path <- paste(Spectral.Species.Path, "_Distribution", sep = "")
  HDR.SSD <- HDR.SS
jbferet's avatar
jbferet committed
178
  # define number of bands
179
  HDR.SSD$bands <- HDR.SS$bands * nbclusters
jbferet's avatar
jbferet committed
180
  # define image size
181
182
  HDR.SSD$samples <- floor(HDR.SS$samples / Spatial.Unit)
  HDR.SSD$lines <- floor(HDR.SS$lines / Spatial.Unit)
jbferet's avatar
jbferet committed
183
  # change resolution
184
185
  HDR.SSD <- Change.Resolution.HDR(HDR.SSD, Spatial.Unit)
  HDR.SSD$`band names` <- NULL
jbferet's avatar
jbferet committed
186
  # create SSD file
187
188
189
190
  fidSSD <- file(
    description = SSD.Path, open = "wb", blocking = TRUE,
    encoding = getOption("encoding"), raw = FALSE
  )
jbferet's avatar
jbferet committed
191
  close(fidSSD)
192
  headerFpath <- paste(SSD.Path, ".hdr", sep = "")
jbferet's avatar
jbferet committed
193
  write.ENVI.header(HDR.SSD, headerFpath)
194
  SeqWrite.SSD <- Where.To.Write.Kernel(HDR.SS, HDR.SSD, nbPieces.Min, Spatial.Unit)
jbferet's avatar
jbferet committed
195
196

  # prepare proportion of sunlit pixels from each spatial unit
197
198
  Sunlit.Path <- paste(SSD.Path, "_Sunlit", sep = "")
  HDR.Sunlit <- HDR.SSD
jbferet's avatar
jbferet committed
199
  # define number of bands
200
  HDR.Sunlit$bands <- 1
jbferet's avatar
jbferet committed
201
  # define number of bands
202
  HDR.Sunlit$`data type` <- 4
jbferet's avatar
jbferet committed
203
  # create SSD Sunlit mask
204
205
206
207
  fidSunlit <- file(
    description = Sunlit.Path, open = "wb", blocking = TRUE,
    encoding = getOption("encoding"), raw = FALSE
  )
jbferet's avatar
jbferet committed
208
  close(fidSunlit)
209
  headerFpath <- paste(Sunlit.Path, ".hdr", sep = "")
jbferet's avatar
jbferet committed
210
  write.ENVI.header(HDR.Sunlit, headerFpath)
211
  SeqWrite.Sunlit <- Where.To.Write.Kernel(HDR.SS, HDR.Sunlit, nbPieces.Min, Spatial.Unit)
jbferet's avatar
jbferet committed
212
213

  # for each piece of image
214
215
216
217
218
219
220
  ReadWrite <- list()
  for (i in 1:nbPieces.Min) {
    ReadWrite[[i]] <- list()
    ReadWrite[[i]]$RW.SS <- ReadWrite[[i]]$RW.SSD <- ReadWrite[[i]]$RW.Sunlit <- list()
    ReadWrite[[i]]$RW.SS$Byte.Start <- SeqRead.SS$ReadByte.Start[i]
    ReadWrite[[i]]$RW.SS$nbLines <- SeqRead.SS$Lines.Per.Chunk[i]
    ReadWrite[[i]]$RW.SS$lenBin <- SeqRead.SS$ReadByte.End[i] - SeqRead.SS$ReadByte.Start[i] + 1
jbferet's avatar
jbferet committed
221

222
223
224
    ReadWrite[[i]]$RW.SSD$Byte.Start <- SeqWrite.SSD$ReadByte.Start[i]
    ReadWrite[[i]]$RW.SSD$nbLines <- SeqWrite.SSD$Lines.Per.Chunk[i]
    ReadWrite[[i]]$RW.SSD$lenBin <- SeqWrite.SSD$ReadByte.End[i] - SeqWrite.SSD$ReadByte.Start[i] + 1
jbferet's avatar
jbferet committed
225

226
227
228
    ReadWrite[[i]]$RW.Sunlit$Byte.Start <- SeqWrite.Sunlit$ReadByte.Start[i]
    ReadWrite[[i]]$RW.Sunlit$nbLines <- SeqWrite.Sunlit$Lines.Per.Chunk[i]
    ReadWrite[[i]]$RW.Sunlit$lenBin <- SeqWrite.Sunlit$ReadByte.End[i] - SeqWrite.Sunlit$ReadByte.Start[i] + 1
jbferet's avatar
jbferet committed
229
  }
230
231
232
233
  ImgFormat <- "3D"
  SSD.Format <- ENVI.Type2Bytes(HDR.SSD)
  SS.Format <- ENVI.Type2Bytes(HDR.SS)
  Sunlit.Format <- ENVI.Type2Bytes(HDR.Sunlit)
jbferet's avatar
jbferet committed
234
235
236

  # multiprocess of spectral species distribution and alpha diversity metrics
  plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
237
238
239
240
241
242
243
244
  Schedule.Per.Thread <- ceiling(nbPieces.Min / nbCPU)
  ALPHA <- future_lapply(ReadWrite,
    FUN = Convert.PCA.to.SSD, Spectral.Species.Path = Spectral.Species.Path,
    HDR.SS = HDR.SS, HDR.SSD = HDR.SSD, SS.Format = SS.Format, SSD.Format = SSD.Format,
    ImgFormat = ImgFormat, Spatial.Unit = Spatial.Unit, nbclusters = nbclusters, MinSun = MinSun,
    pcelim = pcelim, Index.Alpha = Index.Alpha, SSD.Path = SSD.Path, Sunlit.Path = Sunlit.Path,
    Sunlit.Format = Sunlit.Format, future.scheduling = Schedule.Per.Thread
  )
jbferet's avatar
jbferet committed
245
246
  plan(sequential)
  # create ful map from chunks
247
248
249
250
251
252
253
254
255
  Shannon.Mean.Chunk <- Fisher.Mean.Chunk <- Simpson.Mean.Chunk <- list()
  Shannon.SD.Chunk <- Fisher.SD.Chunk <- Simpson.SD.Chunk <- list()
  for (i in 1:length(ALPHA)) {
    Shannon.Mean.Chunk[[i]] <- ALPHA[[i]]$Shannon
    Fisher.Mean.Chunk[[i]] <- ALPHA[[i]]$Fisher
    Simpson.Mean.Chunk[[i]] <- ALPHA[[i]]$Simpson
    Shannon.SD.Chunk[[i]] <- ALPHA[[i]]$Shannon.SD
    Fisher.SD.Chunk[[i]] <- ALPHA[[i]]$Fisher.SD
    Simpson.SD.Chunk[[i]] <- ALPHA[[i]]$Simpson.SD
jbferet's avatar
jbferet committed
256
  }
257
258
259
260
261
262
263
264
265
266
  Shannon.Mean <- do.call(rbind, Shannon.Mean.Chunk)
  Fisher.Mean <- do.call(rbind, Fisher.Mean.Chunk)
  Simpson.Mean <- do.call(rbind, Simpson.Mean.Chunk)
  Shannon.SD <- do.call(rbind, Shannon.SD.Chunk)
  Fisher.SD <- do.call(rbind, Fisher.SD.Chunk)
  Simpson.SD <- do.call(rbind, Simpson.SD.Chunk)
  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
  )
jbferet's avatar
jbferet committed
267
268
269
  return(my_list)
}

270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
# Convert PCA into SSD based on previous clustering
#
# @param ReadWrite
# @param Spectral.Species.Path
# @param HDR.SS
# @param HDR.SSD
# @param SS.Format
# @param SSD.Format
# @param ImgFormat
# @param Spatial.Unit
# @param nbclusters
# @param MinSun
# @param pcelim
# @param Index.Alpha
# @param SSD.Path
# @param Sunlit.Path
# @param Sunlit.Format
#
# @param
# @param
# @return
Convert.PCA.to.SSD <- function(ReadWrite, Spectral.Species.Path, HDR.SS, HDR.SSD,
                               SS.Format, SSD.Format, ImgFormat, Spatial.Unit, nbclusters,
                               MinSun, pcelim, Index.Alpha, SSD.Path, Sunlit.Path, Sunlit.Format) {
  SS.Chunk <- Read.Image.Subset(
    Spectral.Species.Path, HDR.SS,
    ReadWrite$RW.SS$Byte.Start, ReadWrite$RW.SS$lenBin,
    ReadWrite$RW.SS$nbLines, SS.Format, ImgFormat
  )
  SSD.Alpha <- Compute.SSD(SS.Chunk, Spatial.Unit, nbclusters, MinSun, pcelim, Index.Alpha = Index.Alpha)
jbferet's avatar
jbferet committed
300
  # write spectral Species ditribution file
301
302
303
304
305
306
  fidSSD <- file(
    description = SSD.Path, open = "r+b", blocking = TRUE,
    encoding = getOption("encoding"), raw = FALSE
  )
  if (!ReadWrite$RW.SSD$Byte.Start == 1) {
    seek(fidSSD, where = ReadWrite$RW.SSD$Byte.Start - 1, origin = "start", rw = "write")
jbferet's avatar
jbferet committed
307
  }
308
309
  SSD.Chunk <- aperm(array(SSD.Alpha$SSD, c(ReadWrite$RW.SSD$nbLines, HDR.SSD$samples, HDR.SSD$bands)), c(2, 3, 1))
  writeBin(c(SSD.Chunk), fidSSD, size = SSD.Format$Bytes, endian = .Platform$endian, useBytes = FALSE)
jbferet's avatar
jbferet committed
310
311
312
  close(fidSSD)

  # write PCsunlit pixels corresponding to SSD file
313
314
315
316
317
318
  fidSunlit <- file(
    description = Sunlit.Path, open = "r+b", blocking = TRUE,
    encoding = getOption("encoding"), raw = FALSE
  )
  if (!ReadWrite$RW.Sunlit$Byte.Start == 1) {
    seek(fidSunlit, where = ReadWrite$RW.Sunlit$Byte.Start - 1, origin = "start", rw = "write")
jbferet's avatar
jbferet committed
319
  }
320
321
  Sunlit.Chunk <- t(SSD.Alpha$PCsun)
  writeBin(c(Sunlit.Chunk), fidSunlit, size = Sunlit.Format$Bytes, endian = .Platform$endian, useBytes = FALSE)
jbferet's avatar
jbferet committed
322
323
324
325
326
  close(fidSunlit)

  rm(SSD.Chunk)
  rm(Sunlit.Chunk)
  gc()
327
328
329
330
331
332
  Shannon.Mean.Chunk <- apply(SSD.Alpha$Shannon, 1:2, mean)
  Fisher.Mean.Chunk <- apply(SSD.Alpha$Fisher, 1:2, mean)
  Simpson.Mean.Chunk <- apply(SSD.Alpha$Simpson, 1:2, mean)
  Shannon.SD.Chunk <- apply(SSD.Alpha$Shannon, 1:2, sd)
  Fisher.SD.Chunk <- apply(SSD.Alpha$Fisher, 1:2, sd)
  Simpson.SD.Chunk <- apply(SSD.Alpha$Simpson, 1:2, sd)
jbferet's avatar
jbferet committed
333
334
  rm(SSD.Alpha)
  gc()
335
336
337
338
  my_list <- list(
    "Shannon" = Shannon.Mean.Chunk, "Fisher" = Fisher.Mean.Chunk, "Simpson" = Simpson.Mean.Chunk,
    "Shannon.SD" = Shannon.SD.Chunk, "Fisher.SD" = Fisher.SD.Chunk, "Simpson.SD" = Simpson.SD.Chunk
  )
jbferet's avatar
jbferet committed
339
340
341
  return(my_list)
}

342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
# compute spectral species distribution from original spectral species map
#
# @param Image.Chunk 3D image chunk of spectral species
# @param Spatial.Unit 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 pcelim minimum proportion for a spectral species to be included
#
# @return list of alpha diversity metrics for each iteration
Compute.SSD <- function(Image.Chunk, Spatial.Unit, nbclusters, MinSun, pcelim, Index.Alpha = "Shannon") {
  nbi <- floor(dim(Image.Chunk)[1] / Spatial.Unit)
  nbj <- floor(dim(Image.Chunk)[2] / Spatial.Unit)
  nbIter <- dim(Image.Chunk)[3]
  SSDMap <- array(NA, c(nbi, nbj, nbIter * nbclusters))
  shannonIter <- FisherAlpha <- SimpsonAlpha <- array(NA, dim = c(nbi, nbj, nbIter))
  PCsun <- matrix(NA, nrow = nbi, ncol = nbj)
jbferet's avatar
jbferet committed
359
360

  # which spectral indices will be computed
361
362
363
364
  Shannon <- Simpson <- Fisher <- FALSE
  if (length((grep("Shannon", Index.Alpha))) > 0) Shannon <- TRUE
  if (length((grep("Simpson", Index.Alpha))) > 0) Simpson <- TRUE
  if (length((grep("Fisher", Index.Alpha))) > 0) Fisher <- TRUE
jbferet's avatar
jbferet committed
365
366

  # for each kernel in the line
367
  for (ii in 1:nbi) {
jbferet's avatar
jbferet committed
368
    # for each kernel in the column
369
370
371
372
373
    for (jj in 1:nbj) {
      li <- ((ii - 1) * Spatial.Unit) + 1
      ui <- ii * Spatial.Unit
      lj <- ((jj - 1) * Spatial.Unit) + 1
      uj <- jj * Spatial.Unit
jbferet's avatar
jbferet committed
374
      # put all iterations in a 2D matrix shape
375
      ijit <- t(matrix(Image.Chunk[li:ui, lj:uj, ], ncol = nbIter))
jbferet's avatar
jbferet committed
376
      # keep non zero values
377
378
379
380
      ijit <- matrix(ijit[, which(!ijit[1, ] == 0)], nrow = nbIter)
      nb.Pix.Sunlit <- dim(ijit)[2]
      PCsun[ii, jj] <- nb.Pix.Sunlit / Spatial.Unit**2
      if (PCsun[ii, jj] > MinSun) {
jbferet's avatar
jbferet committed
381
        # for each iteration
382
383
384
385
386
387
388
389
        for (it in 1:nbIter) {
          lbk <- (it - 1) * nbclusters
          SSD <- as.vector(table(ijit[it, ]))
          ClusterID <- sort(unique(ijit[it, ]))
          if (pcelim > 0) {
            KeepSS <- which(SSD >= pcelim * nb.Pix.Sunlit)
            ClusterID <- ClusterID[KeepSS]
            SSD <- SSD[KeepSS]
jbferet's avatar
jbferet committed
390
          }
391
392
393
          SSDMap[ii, jj, (lbk + ClusterID)] <- SSD
          if (Shannon == TRUE) {
            shannonIter[ii, jj, it] <- Get.Shannon(SSD)
jbferet's avatar
jbferet committed
394
          }
395
396
          if (Simpson == TRUE) {
            SimpsonAlpha[ii, jj, it] <- Get.Simpson(SSD)
jbferet's avatar
jbferet committed
397
          }
398
399
400
          if (Fisher == TRUE) {
            if (length(SSD) > 2) {
              FisherAlpha[ii, jj, it] <- fisher.alpha(SSD)
jbferet's avatar
jbferet committed
401
            } else {
402
              FisherAlpha[ii, jj, it] <- 0
jbferet's avatar
jbferet committed
403
404
405
406
            }
          }
        }
      } else {
407
408
409
        shannonIter[ii, jj, ] <- NA
        FisherAlpha[ii, jj, ] <- NA
        SimpsonAlpha[ii, jj, ] <- NA
jbferet's avatar
jbferet committed
410
411
412
      }
    }
  }
413
  my_list <- list("Shannon" = shannonIter, "Fisher" = FisherAlpha, "Simpson" = SimpsonAlpha, "SSD" = SSDMap, "PCsun" = PCsun)
jbferet's avatar
jbferet committed
414
415
416
  return(my_list)
}

417
418
419
420
421
422
423
424
425
# computes shannon index from a distribution
#
# @param Distrib Distribution
#
# @return Shannon index correspnding to the distribution
Get.Shannon <- function(Distrib) {
  Distrib <- Distrib / sum(Distrib, na.rm = TRUE)
  Distrib <- Distrib[which(!Distrib == 0)]
  shannon <- -1 * sum(Distrib * log(Distrib), na.rm = TRUE)
jbferet's avatar
jbferet committed
426
427
428
  return(shannon)
}

429
430
431
432
433
434
435
436
# computes Simpson index from a distribution
#
# @param Distrib Distribution
#
# @return Simpson index correspnding to the distribution
Get.Simpson <- function(Distrib) {
  Distrib <- Distrib / sum(Distrib, na.rm = TRUE)
  Simpson <- 1 - sum(Distrib * Distrib, na.rm = TRUE)
jbferet's avatar
jbferet committed
437
438
439
  return(Simpson)
}

440
441
442
443
444
445
446
447
448
449
450
451
# 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 Spatial.Unit 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.Image.Alpha <- function(Image, HDR.SSD, ImagePath, Spatial.Unit, Index, FullRes = TRUE, LowRes = FALSE) {
jbferet's avatar
jbferet committed
452
453

  # Write image with resolution corresponding to Spatial.Unit
454
455
456
457
458
459
460
  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) {
    headerFpath <- paste(ImagePath, ".hdr", sep = "")
jbferet's avatar
jbferet committed
461
    write.ENVI.header(HDR.Alpha, headerFpath)
462
463
464
465
466
467
468
    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)
jbferet's avatar
jbferet committed
469
470
    close(fidOUT)
  }
471
  if (FullRes == TRUE) {
jbferet's avatar
jbferet committed
472
    # Write image with Full native resolution
473
474
475
476
477
478
    HDR.Full <- HDR.Alpha
    HDR.Full$samples <- HDR.Alpha$samples * Spatial.Unit
    HDR.Full$lines <- HDR.Alpha$lines * Spatial.Unit
    HDR.Full <- Revert.Resolution.HDR(HDR.Full, Spatial.Unit)
    ImagePath.FullRes <- paste(ImagePath, "_Fullres", sep = "")
    headerFpath <- paste(ImagePath.FullRes, ".hdr", sep = "")
jbferet's avatar
jbferet committed
479
    write.ENVI.header(HDR.Full, headerFpath)
480
481
482
483
484
    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) * Spatial.Unit + 1):(i * Spatial.Unit), ((j - 1) * Spatial.Unit + 1):(j * Spatial.Unit)] <- Image[i, j]
jbferet's avatar
jbferet committed
485
486
      }
    }
487
488
489
490
491
492
493
    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)
jbferet's avatar
jbferet committed
494
495
496
497
498
    close(fidOUT)
    # zip resulting file
    ZipFile(ImagePath.FullRes)
  }
  # write smoothed map
499
500
501
502
503
504
505
506
507
508
  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) {
jbferet's avatar
jbferet committed
509
      write.ENVI.header(HDR.Alpha, headerFpath)
510
511
512
513
514
515
516
      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)
jbferet's avatar
jbferet committed
517
518
      close(fidOUT)
    }
519
    if (FullRes == TRUE) {
jbferet's avatar
jbferet committed
520
      # Write image with Full native resolution
521
522
      ImagePath.FullRes <- paste(ImagePath.Smooth, "_Fullres", sep = "")
      headerFpath <- paste(ImagePath.FullRes, ".hdr", sep = "")
jbferet's avatar
jbferet committed
523
      write.ENVI.header(HDR.Full, headerFpath)
524
525
526
527
528
      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) * Spatial.Unit + 1):(i * Spatial.Unit), ((j - 1) * Spatial.Unit + 1):(j * Spatial.Unit)] <- Image.Smooth[i, j]
jbferet's avatar
jbferet committed
529
530
        }
      }
531
532
533
534
535
536
537
      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)
jbferet's avatar
jbferet committed
538
539
540
541
542
543
544
      close(fidOUT)
      # zip resulting file
      ZipFile(ImagePath.FullRes)
    }
  }
  return("")
}