Lib_MapBetaDiversity.R 20.1 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_MapBetaDiversity.R
# ==============================================================================
# PROGRAMMERS:
# Jean-Baptiste FERET <jb.feret@irstea.fr>
# Copyright 2018/07 Jean-Baptiste FERET
# ==============================================================================
# This Library computes bray curtis dissimilarity among spatial units based on
# spectral species distribution file and generates a RGB representation with NMDS
# ==============================================================================

#' maps beta diversity indicator based on spectral species distribution
#'
15
16
17
18
19
20
#' @param Input.Image.File character. Path and name of the image to be processed.
#' @param Output.Dir character. Output directory.
#' @param Spatial.Unit numeric. Dimensions of the spatial unit.
#' @param TypePCA character. Type of PCA (PCA, SPCA, NLPCA...).
#' @param nbclusters numeric. Number of clusters defined in k-Means.
#' @param Nb.Units.Ordin numeric. Maximum number of spatial units to be processed in NMDS.
jbferet's avatar
jbferet committed
21
22
#' --> 1000 will be fast but may not capture important patterns if large area
#' --> 4000 will be slow but may show better ability to capture landscape patterns
23
24
25
26
27
28
29
#' @param MinSun numeric. Minimum proportion of sunlit pixels required to consider plot.
#' @param scaling character. "PCO" or "NMDS".
#' @param FullRes boolean.
#' @param LowRes boolean.
#' @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 pcelim numeric. Minimum contribution (in \%) required for a spectral species
jbferet's avatar
jbferet committed
30
31
#'
#' @export
32
33
34
35
36
37
38
39
40
Map.Beta.Diversity <- function(Input.Image.File, Output.Dir, Spatial.Unit,
                               TypePCA = "SPCA", nbclusters = 50,
                               Nb.Units.Ordin = 2000, MinSun = 0.25,
                               pcelim = 0.02, scaling = "PCO", FullRes = TRUE,
                               LowRes = FALSE, nbCPU = 1, MaxRAM = 0.25) {
  Output.Dir2 <- Define.Output.Dir(Output.Dir, Input.Image.File, TypePCA)
  WS_Save <- paste(Output.Dir2, "PCA_Info.RData", sep = "")
  Output.Dir.SS <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "SpectralSpecies")
  Output.Dir.BETA <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "BETA")
jbferet's avatar
jbferet committed
41
  load(file = WS_Save)
42
43
44
  Beta <- Compute.BetaDiversity(Output.Dir.SS, MinSun, Nb.Units.Ordin, NbIter,
                                nbclusters, pcelim, scaling = scaling,
                                nbCPU = nbCPU, MaxRAM = MaxRAM)
jbferet's avatar
jbferet committed
45
  # Create images corresponding to Beta-diversity
46
47
48
49
  print("Write beta diversity maps")
  Index <- paste("BetaDiversity_BCdiss_", scaling, sep = "")
  Beta.Path <- paste(Output.Dir.BETA, Index, "_", Spatial.Unit, sep = "")
  Write.Image.Beta(Beta$BetaDiversity, Beta$HDR, Beta.Path, Spatial.Unit, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
50
51
52
  return()
}

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
# maps beta diversity indicator based on spectral species distribution
# and saves in directories specifying the number of clusters
#
# @param Input.Image.File Path and name of the image 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 Nb.Units.Ordin maximum number of spatial units to be processed in NMDS
# @param MinSun minimum proportion of sunlit pixels required to consider plot
# @param pcelim minimum contribution (in %) required for a spectral species
# --> 1000 will be fast but may not capture important patterns if large area
# @param scaling
# @param FullRes
# @param LowRes
# @param nbCPU
# @param MaxRAM
# --> 4000 will be slow but may show better ability to capture landscape patterns
# @return
Map.Beta.Diversity.TestnbCluster <- function(Input.Image.File, Output.Dir, Spatial.Unit, TypePCA = "SPCA", nbclusters = 50, Nb.Units.Ordin = 2000, MinSun = 0.25, pcelim = 0.02, scaling = "PCO", FullRes = TRUE, LowRes = FALSE, nbCPU = FALSE, MaxRAM = FALSE) {
  Output.Dir2 <- Define.Output.Dir(Output.Dir, Input.Image.File, TypePCA)
  WS_Save <- paste(Output.Dir, "PCA_Info.RData", sep = "")
  Output.Dir.SS <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, paste("SpectralSpecies_", nbclusters, sep = ""))
  Output.Dir.BETA <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, paste("BETA_", nbclusters, sep = ""))
jbferet's avatar
jbferet committed
77
  load(file = WS_Save)
78
  Beta <- Compute.BetaDiversity(Output.Dir.SS, MinSun, Nb.Units.Ordin, NbIter, nbclusters, pcelim, scaling = scaling, nbCPU = nbCPU, MaxRAM = MaxRAM)
jbferet's avatar
jbferet committed
79
  # Create images corresponding to Beta-diversity
80
81
82
83
  print("Write beta diversity maps")
  Index <- paste("BetaDiversity_BCdiss_", scaling, sep = "")
  Beta.Path <- paste(Output.Dir.BETA, Index, "_", Spatial.Unit, sep = "")
  Write.Image.Beta(Beta$BetaDiversity, Beta$HDR, Beta.Path, Spatial.Unit, FullRes = FullRes, LowRes = LowRes)
jbferet's avatar
jbferet committed
84
85
86
  return()
}

87
88
89
90
91
92
93
# Gets sunlit pixels from SpectralSpecies_Distribution_Sunlit
#
# @param ImPathSunlit path for SpectralSpecies_Distribution_Sunlit
# @param MinSun minimum proportion of sunlit pixels required
#
# @return
Get.Sunlit.Pixels <- function(ImPathSunlit, MinSun) {
jbferet's avatar
jbferet committed
94
95

  # Filter out spatial units showing poor illumination
96
97
98
99
100
101
102
103
  Sunlit.HDR <- Get.HDR.Name(ImPathSunlit)
  HDR.Sunlit <- read.ENVI.header(Sunlit.HDR)
  nbpix <- as.double(HDR.Sunlit$lines) * as.double(HDR.Sunlit$samples)
  fid <- file(
    description = ImPathSunlit, open = "rb", blocking = TRUE,
    encoding = getOption("encoding"), raw = FALSE
  )
  Sunlit <- readBin(fid, numeric(), n = nbpix, size = 4)
jbferet's avatar
jbferet committed
104
  close(fid)
105
  Sunlit <- aperm(array(Sunlit, dim = c(HDR.Sunlit$samples, HDR.Sunlit$lines)))
jbferet's avatar
jbferet committed
106
  # define sunlit spatial units
107
  Select.Sunlit <- which(Sunlit > MinSun)
jbferet's avatar
jbferet committed
108
  # define where to extract each datapoint in the file
109
110
  coordi <- ((Select.Sunlit - 1) %% HDR.Sunlit$lines) + 1
  coordj <- floor((Select.Sunlit - 1) / HDR.Sunlit$lines) + 1
jbferet's avatar
jbferet committed
111
  # sort based on line and col (important for optimal scan of large files)
112
  coordTot <- cbind(coordi, coordj)
jbferet's avatar
jbferet committed
113
114
  # sort samples from top to bottom in order to optimize read/write of the image
  # indxTot saves the order of the data for reconstruction later
115
116
117
  indxTot <- order(coordTot[, 1], coordTot[, 2], na.last = FALSE)
  coordTotSort <- coordTot[indxTot, ]
  Select.Sunlit <- Select.Sunlit[indxTot]
jbferet's avatar
jbferet committed
118

119
  my_list <- list("Select.Sunlit" = Select.Sunlit, "coordTotSort" = coordTotSort)
jbferet's avatar
jbferet committed
120
121
122
  return(my_list)
}

123
124
125
126
127
# computes NMDS
#
# @param MatBCdist BC dissimilarity matrix
#
# @return BetaNMDS.sel
De Boissieu Florian's avatar
De Boissieu Florian committed
128
#' @importFrom future plan multiprocess sequential
129
130
131
132
133
134
135
#' @importFrom future.apply future_lapply
Compute.NMDS <- function(MatBCdist) {
  nbiterNMDS <- 4
  if (Sys.info()["sysname"] == "Windows") {
    nbCoresNMDS <- 2
  } else if (Sys.info()["sysname"] == "Linux") {
    nbCoresNMDS <- 4
jbferet's avatar
jbferet committed
136
137
138
  }
  # multiprocess of spectral species distribution and alpha diversity metrics
  plan(multiprocess, workers = nbCoresNMDS) ## Parallelize using four cores
139
  BetaNMDS <- future_lapply(MatBCdist, FUN = NMDS, mindim = 3, maxdim = 3, nits = 1, future.packages = c("ecodist"))
jbferet's avatar
jbferet committed
140
141
  plan(sequential)
  # find iteration with minimum stress
142
143
144
  Stress <- vector(length = nbiterNMDS)
  for (i in 1:nbiterNMDS) {
    Stress[i] <- BetaNMDS[[i]]$stress
jbferet's avatar
jbferet committed
145
146
147
148
149
150
151
152
  }
  print("Stress obtained for NMDS iterations:")
  print(Stress)
  print("Rule of thumb")
  print("stress < 0.05 provides an excellent represention in reduced dimensions")
  print("stress < 0.1 is great")
  print("stress < 0.2 is good")
  print("stress > 0.3 provides a poor representation")
153
154
155
  MinStress <- find(Stress == min(Stress))
  BetaNMDS.sel <- BetaNMDS[[MinStress]]$conf
  BetaNMDS.sel <- data.frame(BetaNMDS.sel[[1]])
jbferet's avatar
jbferet committed
156
157
158
  return(BetaNMDS.sel)
}

159
160
161
162
163
164
165
166
167
168
169
170
171
# identifies ordination coordinates based on nearest neighbors
#
# @param Beta.Ordination.sel ordination
# @param SSD.Path ath for spectral species distribution file
# @param Sample.Sel Samples selected during ordination
# @param coordTotSort coordinates of sunlit spatial units
# @param NbIter number of iterations
# @param nbclusters number of clusters
# @param pcelim number of CPUs available
# @param nbCPU number of CPUs available
#
# @return Ordination.est coordinates of each spatial unit in ordination space
#' @importFrom snow splitRows
De Boissieu Florian's avatar
De Boissieu Florian committed
172
#' @importFrom future plan multiprocess sequential
173
174
175
#' @importFrom future.apply future_lapply
Ordination.to.NN <- function(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, NbIter, nbclusters, pcelim, nbCPU = FALSE) {
  nb.Sunlit <- dim(coordTotSort)[1]
jbferet's avatar
jbferet committed
176
  # define number of samples to be sampled each time during paralle processing
177
  nb.samples.per.sub <- round(1e7 / dim(Sample.Sel)[1])
jbferet's avatar
jbferet committed
178
  # number of paralle processes to run
179
180
181
  nb.sub <- round(nb.Sunlit / nb.samples.per.sub)
  if (nb.sub == 0) nb.sub <- 1
  id.sub <- splitRows(as.matrix(seq(1, nb.Sunlit, by = 1), ncol = 1), ncl = nb.sub)
jbferet's avatar
jbferet committed
182
  # compute ordination coordinates from each subpart
183
  Nb.Units.Ordin <- dim(Sample.Sel)[1]
jbferet's avatar
jbferet committed
184
  plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
185
186
187
188
189
190
191
  Schedule.Per.Thread <- ceiling(nb.sub / nbCPU)
  OutPut <- future_lapply(id.sub,
    FUN = Ordination.Parallel, coordTotSort = coordTotSort, SSD.Path = SSD.Path,
    Sample.Sel = Sample.Sel, Beta.Ordination.sel = Beta.Ordination.sel, Nb.Units.Ordin = Nb.Units.Ordin,
    NbIter = NbIter, nbclusters = nbclusters, pcelim = pcelim, future.scheduling = Schedule.Per.Thread,
    future.packages = c("vegan", "dissUtils", "R.utils", "tools", "snow", "matlab")
  )
jbferet's avatar
jbferet committed
192
  plan(sequential)
193
  Ordination.est <- do.call("rbind", OutPut)
jbferet's avatar
jbferet committed
194
195
196
197
  gc()
  return(Ordination.est)
}

198
199
200
201
202
203
204
205
206
207
208
209
210
211
# applies results of ordination to full image based on nearest neighbors
#
# @param id.sub
# @param coordTotSort
# @param SSD.Path
# @param Sample.Sel
# @param Beta.Ordination.sel
# @param Nb.Units.Ordin
# @param NbIter
# @param nbclusters
# @param pcelim
#
# @return list of mean and SD of alpha diversity metrics
Ordination.Parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta.Ordination.sel, Nb.Units.Ordin, NbIter, nbclusters, pcelim) {
jbferet's avatar
jbferet committed
212
213

  # get Spectral species distribution
214
215
  coordPix <- coordTotSort[id.sub, ]
  SSD.NN <- Extract.Samples.From.Image(SSD.Path, coordPix)
jbferet's avatar
jbferet committed
216
  # compute the mean BC dissimilarity sequentially for each iteration
217
218
219
220
221
222
223
224
225
  MatBCtmp <- matrix(0, nrow = nrow(id.sub), ncol = Nb.Units.Ordin)
  SSDList <- list()
  for (i in 1:NbIter) {
    lb <- 1 + (i - 1) * nbclusters
    ub <- i * nbclusters
    SSDList[[1]] <- SSD.NN[, lb:ub]
    SSDList[[2]] <- Sample.Sel[, lb:ub]
    MatBCtmp0 <- Compute.BCdiss(SSDList, pcelim)
    MatBCtmp <- MatBCtmp + MatBCtmp0
jbferet's avatar
jbferet committed
226
  }
227
  MatBCtmp <- MatBCtmp / NbIter
jbferet's avatar
jbferet committed
228
  # get the knn closest neighbors from each kernel
229
230
  knn <- 3
  OutPut <- Compute.NN.From.Ordination(MatBCtmp, knn, Beta.Ordination.sel)
jbferet's avatar
jbferet committed
231
232
233
  return(OutPut)
}

234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
# computes beta diversity
#
# @param Output.Dir directory where spectral species are stored
# @param MinSun minimum proportion of sunlit pixels required to consider plot
# @param Nb.Units.Ordin maximum number of spatial units to be processed in Ordination
# @param NbIter number of iterations
# @param nbclusters number of clusters defined in k-Means
# @param scaling
# @param nbCPU
# @param MaxRAM
# @param pcelim min proprtion for a spectral species in a spatial unit to be considerd
#
# @return
#' @importFrom labdsv pco
Compute.BetaDiversity <- function(Output.Dir, MinSun, Nb.Units.Ordin, NbIter, nbclusters, pcelim, scaling = "PCO", nbCPU = FALSE, MaxRAM = FALSE) {
jbferet's avatar
jbferet committed
249
  # Define path for images to be used
250
251
  SSD.Path <- paste(Output.Dir, "SpectralSpecies_Distribution", sep = "")
  ImPathSunlit <- paste(Output.Dir, "SpectralSpecies_Distribution_Sunlit", sep = "")
jbferet's avatar
jbferet committed
252
  # Get illuminated pixels based on  SpectralSpecies_Distribution_Sunlit
253
254
255
  Sunlit.Pixels <- Get.Sunlit.Pixels(ImPathSunlit, MinSun)
  Select.Sunlit <- Sunlit.Pixels$Select.Sunlit
  nb.Sunlit <- length(Select.Sunlit)
jbferet's avatar
jbferet committed
256
257
  # Define spatial units processed through ordination and those processed through
  # Nearest neighbor based on the first ordination batch
258
259
260
261
  print("Read Spectral Species distribution")
  RandPermKernels <- sample(seq(1, nb.Sunlit, by = 1))
  if (Nb.Units.Ordin <= nb.Sunlit) {
    Kernels.NN <- RandPermKernels[(Nb.Units.Ordin + 1):nb.Sunlit]
jbferet's avatar
jbferet committed
262
  } else {
263
264
    Nb.Units.Ordin <- nb.Sunlit
    Kernels.NN <- c()
jbferet's avatar
jbferet committed
265
266
  }
  # read spectral species distribution file
267
  SSD.All <- Extract.Samples.From.Image(SSD.Path, Sunlit.Pixels$coordTotSort)
jbferet's avatar
jbferet committed
268
  # define kernels used for Ordination
269
270
  Kernels.Ordination <- RandPermKernels[1:Nb.Units.Ordin]
  Sample.Sel <- SSD.All[Kernels.Ordination, ]
jbferet's avatar
jbferet committed
271
272
273
274
  rm(SSD.All)
  gc()

  # create a Bray curtis dissimilarity matrix for each iteration
275
  print("compute BC dissimilarity for selected kernels")
jbferet's avatar
jbferet committed
276
  # create a list in with each element is an iteration
277
278
279
280
281
282
283
284
285
286
  MatBC <- matrix(0, ncol = Nb.Units.Ordin, nrow = Nb.Units.Ordin)
  SSDList <- list()
  BC.from.SSD <- list()
  for (i in 1:NbIter) {
    lb <- 1 + (i - 1) * nbclusters
    ub <- i * nbclusters
    SSDList[[1]] <- Sample.Sel[, lb:ub]
    SSDList[[2]] <- Sample.Sel[, lb:ub]
    BC.from.SSD <- Compute.BCdiss(SSDList, pcelim)
    MatBC <- MatBC + BC.from.SSD
jbferet's avatar
jbferet committed
287
  }
288
  MatBC <- MatBC / NbIter
jbferet's avatar
jbferet committed
289
290
291
292
293

  # Perform Ordination based on BC dissimilarity matrix
  print("perform Ordination on the BC dissimilarity matrix averaged from all iterations")
  # parallel computing of Ordination can be run on 2 cores on Windows.
  # core management seems better on linux --> 4 cores possible
294
295
296
297
298
299
  MatBCdist <- as.dist(MatBC, diag = FALSE, upper = FALSE)
  if (scaling == "NMDS") {
    Beta.Ordination.sel <- Compute.NMDS(MatBCdist)
  } else if (scaling == "PCO") {
    BetaPCO <- pco(MatBCdist, k = 3)
    Beta.Ordination.sel <- BetaPCO$points
jbferet's avatar
jbferet committed
300
301
302
  }

  # Perform nearest neighbor on spatial units excluded from Ordination
303
304
305
  print("BC dissimilarity between samples selected for Ordination and remaining")
  coordTotSort <- Sunlit.Pixels$coordTotSort
  Ordination.est <- Ordination.to.NN(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, NbIter, nbclusters, pcelim, nbCPU = nbCPU)
jbferet's avatar
jbferet committed
306
307

  # Reconstuct spatialized beta diversity map from previous computation
308
309
310
311
312
313
314
315
  Sunlit.HDR <- Get.HDR.Name(ImPathSunlit)
  HDR.Sunlit <- read.ENVI.header(Sunlit.HDR)
  BetaDiversity <- as.matrix(Ordination.est, ncol = 3)
  BetaDiversityRGB <- array(NA, c(as.double(HDR.Sunlit$lines), as.double(HDR.Sunlit$samples), 3))
  BetaTmp <- matrix(NA, nrow = as.double(HDR.Sunlit$lines), ncol = as.double(HDR.Sunlit$samples))
  for (i in 1:3) {
    BetaTmp[Select.Sunlit] <- BetaDiversity[, i]
    BetaDiversityRGB[, , i] <- BetaTmp
jbferet's avatar
jbferet committed
316
  }
317
318
  list <- ls()
  rm(list = list[-which(list == "BetaDiversityRGB" | list == "Select.Sunlit" | list == "HDR.Sunlit")])
jbferet's avatar
jbferet committed
319
  gc()
320
  my_list <- list("BetaDiversity" = BetaDiversityRGB, "Select.Sunlit" = Select.Sunlit, "HDR" = HDR.Sunlit)
jbferet's avatar
jbferet committed
321
322
323
  return(my_list)
}

324
325
326
327
328
329
330
331
332
333
334
# compute the bray curtis dissimilarity matrix corresponding to a list of kernels
# (rows) defined by their spectral species (columns)
# SSDList is a list containing spectral species distribution for two sets of kernels ([[1]] and [[2]])
# pcelim is the threshold for minimum contributin of a spctral species to be kept
#
# @param SSDList list of 2 groups to compute BC dissimilarity from
# @param pcelim minimum proportion required for a species to be included (currently deactivated)
#
# @return MatBC matrix of bray curtis dissimilarity matrix
#' @importFrom dissUtils diss
Compute.BCdiss <- function(SSDList, pcelim) {
jbferet's avatar
jbferet committed
335
336
337
338
339
  # compute the proportion of each spectral species
  # Here, the proportion is computed with regards to the total number of sunlit pixels
  # One may want to determine if the results are similar when the proportion is computed
  # with regards to the total number of pixels (se*se)
  # however it would increase dissimilarity betwen kernels with different number of sunlit pixels
340
341
  SSD <- list()
  for (i in 1:length(SSDList)) {
jbferet's avatar
jbferet committed
342
    # get the total number of sunlit pixels in spatial unit
343
    SumSpecies <- rowSums(SSDList[[i]])
jbferet's avatar
jbferet committed
344
345
346
347
348
349
350
351
352
    # EDIT 02-Feb-2019: changed this
    # check if some kernels contain no data and eliminate if it is the case
    # elim          = which(SumSpecies==0)
    # if (length(elim)>0){
    #   SumSpecies  = SumSpecies[-elim]
    #   SSDList[[i]]= SSDList[[i]][-elim,]
    # }
    # SSD[[i]]      = apply(SSDList[[i]],2,function(x,c1) x/c1,"c1"=SumSpecies)
    # SSD[[i]][which(SSD[[i]]<pcelim)]  = 0
353
354
355
356
    elim <- which(SumSpecies == 0)
    if (length(elim) > 0) {
      SumSpecies[elim] <- 1
      SSDList[[i]][elim, ] <- 0
jbferet's avatar
jbferet committed
357
    }
358
359
    SSD[[i]] <- apply(SSDList[[i]], 2, function(x, c1) x / c1, "c1" = SumSpecies)
    SSD[[i]][which(SSD[[i]] < pcelim)] <- 0
jbferet's avatar
jbferet committed
360
361
362
  }
  # matrix of bray curtis dissimilarity (size = nb kernels x nb kernels)
  # Here use the package "dissUtils" to compute dissimilarity matrices sequentially
363
  MatBC <- diss(SSD[[1]], SSD[[2]], method = "braycurtis")
jbferet's avatar
jbferet committed
364
  # EDIT 06-Feb-2019: added this to fix problem when empty kernels occur, leading to NA BC value
365
366
  if (length(which(is.na(MatBC) == TRUE)) > 0) {
    MatBC[which(is.na(MatBC) == TRUE)] <- 1
jbferet's avatar
jbferet committed
367
368
369
370
  }
  return(MatBC)
}

371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
# compute the nearest neighbors among kernels used in NMDS
#
# @param MatBC3 matrix of BC dissimilarity between the kernels excluded from Ordination (rows)
# @param knn number of neighbors
# @param BetaDiversity0
#
# @return estimated NMDS position based on nearest neighbors from NMDS
Compute.NN.From.Ordination <- function(MatBC3, knn, BetaDiversity0) {
  Ordin.est <- matrix(0, ncol = 3, nrow = nrow(MatBC3))
  for (i in 1:nrow(MatBC3)) {
    NNtmp <- sort(MatBC3[i, ], decreasing = FALSE, index.return = TRUE)
    Dist.Tot <- sum(NNtmp$x[1:knn])
    aa <- as.numeric(((Dist.Tot - NNtmp$x[1]) / (2 * Dist.Tot)) * BetaDiversity0[NNtmp$ix[1], ]
      + ((Dist.Tot - NNtmp$x[2]) / (2 * Dist.Tot)) * (BetaDiversity0[NNtmp$ix[2], ])
      + ((Dist.Tot - NNtmp$x[3]) / (2 * Dist.Tot)) * BetaDiversity0[NNtmp$ix[3], ])
    Ordin.est[i, 1:3] <- aa
jbferet's avatar
jbferet committed
387
388
389
390
  }
  return(Ordin.est)
}

391
392
393
394
395
396
397
398
399
400
401
# 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 Spatial.Unit 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.Image.Beta <- function(Image, HDR.SSD, ImagePath, Spatial.Unit, FullRes = TRUE, LowRes = FALSE) {
jbferet's avatar
jbferet committed
402
  # Write image with resolution corresponding to Spatial.Unit
403
404
405
406
407
408
  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))
jbferet's avatar
jbferet committed
409
  }
410
411
412
413
414
  PCs <- paste(PCs, collapse = ", ")
  HDR.Beta$`band names` <- PCs
  Image.Format <- ENVI.Type2Bytes(HDR.Beta)
  if (LowRes == TRUE) {
    headerFpath <- paste(ImagePath, ".hdr", sep = "")
jbferet's avatar
jbferet committed
415
    write.ENVI.header(HDR.Beta, headerFpath)
416
417
418
419
420
421
    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)
jbferet's avatar
jbferet committed
422
423
    close(fidOUT)
  }
424
  if (FullRes == TRUE) {
jbferet's avatar
jbferet committed
425
    # Write image with Full native resolution
426
427
428
429
430
431
    HDR.Full <- HDR.Beta
    HDR.Full$samples <- HDR.Beta$samples * Spatial.Unit
    HDR.Full$lines <- HDR.Beta$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
432
    write.ENVI.header(HDR.Full, headerFpath)
433
434
435
436
437
438
    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) * Spatial.Unit + 1):(i * Spatial.Unit), ((j - 1) * Spatial.Unit + 1):(j * Spatial.Unit), b] <- Image[i, j, b]
jbferet's avatar
jbferet committed
439
440
441
        }
      }
    }
442
443
444
445
446
447
    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)
jbferet's avatar
jbferet committed
448
449
450
451
452
453
    close(fidOUT)
    # zip resulting file
    ZipFile(ImagePath.FullRes)
  }
  return("")
}