Mercurial > repos > ecology > srs_pca

comparison prosail-master/R/Lib_PROSAIL_HybridInversion.R @ 0:0931c98270c9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author ecology
date Mon, 09 Jan 2023 13:38:09 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0931c98270c9
1 # ============================================================================= =
2 # prosail
3 # Lib_PROSAIL_HybridInversion.R
4 # ============================================================================= =
5 # PROGRAMMERS:
6 # Jean-Baptiste FERET <jb.feret@teledetection.fr>
7 # Florian de BOISSIEU <fdeboiss@gmail.com>
8 # Copyright 2019 / 11 Jean-Baptiste FERET
9 # ============================================================================= =
10 # This Library includes functions dedicated to PROSAIL inversion using hybrid
11 # approach based on SVM regression
12 # ============================================================================= =
13
14
15 #" This function applies SVR model on raster data in order to estimate
16 #" vegetation biophysical properties
17 #"
18 #" @param raster_path character. path for a raster file
19 #" @param hybridmodel list. hybrid models produced from train_prosail_inversion
20 #" each element of the list corresponds to a set of hybrid models for a vegetation parameter
21 #" @param pathout character. path for directory where results are written
22 #" @param selectedbands list. list of spectral bands to be selected from raster (identified by name of vegetation parameter)
23 #" @param bandname character. spectral bands corresponding to the raster
24 #" @param maskraster character. path for binary mask defining ON (1) and OFF (0) pixels in the raster
25 #" @param multiplyingfactor numeric. multiplying factor used to write reflectance in the raster
26 #" --> PROSAIL simulates reflectance between 0 and 1, and raster data expected in the same range
27 #"
28 #" @return None
29 #" @importFrom progress progress_bar
30 #" @importFrom stars read_stars
31 #" @importFrom raster raster brick blockSize readStart readStop getValues writeStart writeStop writeValues
32 #" @import rgdal
33 #" @export
34 apply_prosail_inversion <- function(raster_path, hybridmodel, pathout,
35 selectedbands, bandname,
36 maskraster = FALSE, multiplyingfactor = 10000) {
37
38 # explain which biophysical variables will be computed
39 bpvar <- names(hybridmodel)
40 print("The following biophysical variables will be computed")
41 print(bpvar)
42
43 # get image dimensions
44 if (attr(rgdal::GDALinfo(raster_path, returnStats = FALSE), "driver") == "ENVI") {
45 hdr <- read_envi_header(get_hdr_name(raster_path))
46 dimsraster <- list("rows" = hdr$lines, "cols" = hdr$samples, "bands" = hdr$bands)
47 } else {
48 dimsraster <- dim(read_stars(raster_path))
49 dimsraster <- list("rows" = as.numeric(dimsraster[2]), "cols" = as.numeric(dimsraster[1]), "bands" = as.numeric(dimsraster[3]))
50 }
51
52 # Produce a map for each biophysical property
53 for (parm in bpvar) {
54 print(paste("Computing", parm, sep = " "))
55 # read by chunk to avoid memory problem
56 blk <- blockSize(brick(raster_path))
57 # reflectance file
58 r_in <- readStart(brick(raster_path))
59 # mask file
60 r_inmask <- FALSE
61 if (maskraster == FALSE) {
62 selectpixels <- "ALL"
63 } else if (!maskraster == FALSE) {
64 if (file.exists(maskraster)) {
65 r_inmask <- readStart(raster(maskraster))
66 } else if (!file.exists(maskraster)) {
67 message("WARNING: Mask file does not exist:")
68 print(maskraster)
69 message("Processing all image")
70 selectpixels <- "ALL"
71 }
72 }
73 # initiate progress bar
74 pgbarlength <- length(hybridmodel[[parm]]) * blk$n
75 pb <- progress_bar$new(
76 format = "Hybrid inversion on raster [:bar] :percent in :elapsedfull, estimated time remaining :eta",
77 total = pgbarlength, clear = FALSE, width = 100)
78
79 # output files
80 bpvarpath <- file.path(pathout, paste(basename(raster_path), parm, sep = "_"))
81 bpvarsdpath <- file.path(pathout, paste(basename(raster_path), parm, "STD", sep = "_"))
82 r_outmean <- writeStart(raster(raster_path), filename = bpvarpath, format = "ENVI", overwrite = TRUE)
83 r_outsd <- writeStart(raster(raster_path), filename = bpvarsdpath, format = "ENVI", overwrite = TRUE)
84 selbands <- match(selectedbands[[parm]], bandname)
85
86 # loop over blocks
87 for (i in seq_along(blk$row)) {
88 # read values for block
89 # format is a matrix with rows the cells values and columns the layers
90 blockval <- getValues(r_in, row = blk$row[i], nrows = blk$nrows[i])
91 fulllength <- dim(blockval)[1]
92
93 if (typeof(r_inmask) == "logical") {
94 blockval <- blockval[, selbands]
95 # automatically filter pixels corresponding to negative values
96 selectpixels <- which(blockval[, 1] > 0)
97 blockval <- blockval[selectpixels, ]
98 } else if (typeof(r_inmask) == "S4") {
99 maskval <- getValues(r_inmask, row = blk$row[i], nrows = blk$nrows[i])
100 selectpixels <- which(maskval == 1)
101 blockval <- blockval[selectpixels, selbands]
102 }
103 mean_estimatefull <- NA * vector(length = fulllength)
104 std_estimatefull <- NA * vector(length = fulllength)
105 if (length(selectpixels) > 0) {
106 blockval <- blockval / multiplyingfactor
107 modelsvr_estimate <- list()
108 for (modind in 1:seq_along(hybridmodel[[parm]])) {
109 pb$tick()
110 modelsvr_estimate[[modind]] <- predict(hybridmodel[[parm]][[modind]], blockval)
111 }
112 modelsvr_estimate <- do.call(cbind, modelsvr_estimate)
113 # final estimated value = mean parm value for all models
114 mean_estimate <- rowMeans(modelsvr_estimate)
115 # "uncertainty" = STD value for all models
116 std_estimate <- rowSds(modelsvr_estimate)
117 mean_estimatefull[selectpixels] <- mean_estimate
118 std_estimatefull[selectpixels] <- std_estimate
119 } else {
120 for (modind in 1:seq_along(hybridmodel[[parm]])) {
121 pb$tick()
122 }
123 }
124 r_outmean <- writeValues(r_outmean, mean_estimatefull, blk$row[i], format = "ENVI", overwrite = TRUE)
125 r_outsd <- writeValues(r_outsd, std_estimatefull, blk$row[i], format = "ENVI", overwrite = TRUE)
126 }
127 # close files
128 r_in <- readStop(r_in)
129 if (typeof(r_inmask) == "S4") {
130 r_inmask <- readStop(r_inmask)
131 }
132 r_outmean <- writeStop(r_outmean)
133 r_outsd <- writeStop(r_outsd)
134 # write biophysical variable name in headers
135 hdr <- read_envi_header(get_hdr_name(bpvarpath))
136 hdr$`band names` <- paste("{", parm, "}", sep = "")
137 write_envi_header(hdr, get_hdr_name(bpvarpath))
138 }
139 print("processing completed")
140 return(invisible())
141 }
142
143 #" get hdr name from image file name, assuming it is BIL format
144 #"
145 #" @param impath path of the image
146 #"
147 #" @return corresponding hdr
148 #" @importFrom tools file_ext file_path_sans_ext
149 #" @export
150 get_hdr_name <- function(impath) {
151 if (tools::file_ext(impath) == "") {
152 impathhdr <- paste(impath, ".hdr", sep = "")
153 } else if (tools::file_ext(impath) == "bil") {
154 impathhdr <- gsub(".bil", ".hdr", impath)
155 } else if (tools::file_ext(impath) == "zip") {
156 impathhdr <- gsub(".zip", ".hdr", impath)
157 } else {
158 impathhdr <- paste(tools::file_path_sans_ext(impath), ".hdr", sep = "")
159 }
160
161 if (!file.exists(impathhdr)) {
162 message("WARNING : COULD NOT FIND hdr FILE")
163 print(impathhdr)
164 message("Process may stop")
165 }
166 return(impathhdr)
167 }
168
169 #" This function applies the regression models trained with prosail_hybrid_train
170 #"
171 #" @param regressionmodels list. List of regression models produced by prosail_hybrid_train
172 #" @param refl numeric. LUT of bidirectional reflectances factors used for training
173 #"
174 #" @return hybridres list. Estimated values corresponding to refl. Includes
175 #" - meanestimate = mean value for the ensemble regression model
176 #" - stdestimate = std value for the ensemble regression model
177 #" @importFrom stats predict
178 #" @importFrom matrixStats rowSds
179 #" @importFrom progress progress_bar
180 #" @export
181
182 prosail_hybrid_apply <- function(regressionmodels, refl) {
183
184 # make sure refl is right dimensions
185 refl <- t(refl)
186 nbfeatures <- regressionmodels[[1]]$dim
187 if (!ncol(refl) == nbfeatures && nrow(refl) == nbfeatures) {
188 refl <- t(refl)
189 }
190 nbensemble <- length(regressionmodels)
191 estimatedval <- list()
192 pb <- progress_bar$new(
193 format = "Applying SVR models [:bar] :percent in :elapsed",
194 total = nbensemble, clear = FALSE, width = 100)
195 for (i in 1:nbensemble) {
196 pb$tick()
197 estimatedval[[i]] <- predict(regressionmodels[[i]], refl)
198 }
199 estimatedval <- do.call(cbind, estimatedval)
200 meanestimate <- rowMeans(estimatedval)
201 stdestimate <- rowSds(estimatedval)
202 hybridres <- list("meanestimate" = meanestimate, "stdestimate" = stdestimate)
203 return(hybridres)
204 }
205
206 #" This function trains a suppot vector regression for a set of variables based on spectral data
207 #"
208 #" @param brf_lut numeric. LUT of bidirectional reflectances factors used for training
209 #" @param inputvar numeric. biophysical parameter corresponding to the reflectance
210 #" @param figplot Boolean. Set to TRUE if you want a scatterplot
211 #" @param nbensemble numeric. Number of individual subsets should be generated from brf_lut
212 #" @param withreplacement Boolean. should subsets be generated with or without replacement?
213 #"
214 #" @return modelssvr list. regression models trained for the retrieval of inputvar based on brf_lut
215 #" @importFrom liquidSVM svmRegression
216 #" @importFrom stats predict
217 #" @importFrom progress progress_bar
218 #" @importFrom graphics par
219 #" @importFrom expandFunctions reset.warnings
220 #" @importFrom stringr str_split
221 #" @importFrom simsalapar tryCatch.W.E
222 #" @import dplyr
223 #" @import ggplot2
224 # @" @import caret
225 #" @export
226
227 prosail_hybrid_train <- function(brf_lut, inputvar, figplot = FALSE, nbensemble = 20, withreplacement = FALSE) {
228
229 x <- y <- ymean <- ystdmin <- ystdmax <- NULL
230 # split the LUT into nbensemble subsets
231 nbsamples <- length(inputvar)
232 if (dim(brf_lut)[2] == nbsamples) {
233 brf_lut <- t(brf_lut)
234 }
235
236 # if subsets are generated from brf_lut with replacement
237 if (withreplacement == TRUE) {
238 subsets <- list()
239 samples_per_run <- round(nbsamples / nbensemble)
240 for (run in (1:nbensemble)) {
241 subsets[[run]] <- sample(seq(1, nbsamples), samples_per_run, replace = TRUE)
242 }
243 # if subsets are generated from brf_lut without replacement
244 } else if (withreplacement == FALSE) {
245 subsets <- split(sample(seq(1, nbsamples, by = 1)), seq(1, nbensemble, by = 1))
246 }
247
248 # run training for each subset
249 modelssvr <- list()
250 predictedyall <- list()
251 tunedmodelyall <- list()
252 pb <- progress_bar$new(
253 format = "Training SVR on subsets [:bar] :percent in :elapsed",
254 total = nbensemble, clear = FALSE, width = 100)
255 for (i in 1:nbensemble) {
256 pb$tick()
257 Sys.sleep(1 / 100)
258 trainingset <- list()
259 trainingset$X <- brf_lut[subsets[i][[1]], ]
260 trainingset$Y <- inputvar[subsets[i][[1]]]
261 # liquidSVM
262 r1 <- tryCatch.W.E(tunedmodel <- liquidSVM::svmRegression(trainingset$X, trainingset$Y))
263 if (!is.null(r1$warning)) {
264 msg <- r1$warning$message
265 valgamma <- str_split(string = msg, pattern = "gamma=")[[1]][2]
266 vallambda <- str_split(string = msg, pattern = "lambda=")[[1]][2]
267 if (!is.na(as.numeric(valgamma))) {
268 message("Adjusting Gamma accordingly")
269 valgamma <- as.numeric(valgamma)
270 tunedmodel <- liquidSVM::svmRegression(trainingset$X, trainingset$Y, min_gamma = valgamma)
271 }
272 if (!is.na(as.numeric(vallambda))) {
273 message("Adjusting Lambda accordingly")
274 vallambda <- as.numeric(vallambda)
275 tunedmodel <- liquidSVM::svmRegression(trainingset$X, trainingset$Y, min_lambda = vallambda)
276 }
277 }
278 modelssvr[[i]] <- tunedmodel
279 }
280
281 # if scatterplots needed
282 if (figplot == TRUE) {
283 # predict for full brf_lut
284 for (i in 1:nbensemble) {
285 tunedmodely <- stats::predict(modelssvr[[i]], brf_lut)
286 tunedmodelyall <- cbind(tunedmodelyall, matrix(tunedmodely, ncol = 1))
287 }
288 # plot prediction
289 df <- data.frame(x = rep(1:nbsamples, nbensemble), y = as.numeric(matrix(tunedmodelyall, ncol = 1)))
290 df_summary <- df %>%
291 dplyr::group_by(x) %>%
292 summarize(ymin = min(y), ystdmin = mean(y) - sd(y),
293 ymax = max(y), ystdmax = mean(y) + sd(y),
294 ymean = mean(y))
295 par(mar = rep(.1, 4))
296 p <- ggplot(df_summary, aes(x = inputvar, y = ymean)) +
297 geom_point(size = 2) +
298 geom_errorbar(aes(ymin = ystdmin, ymax = ystdmax))
299 meanpredict <- rowMeans(matrix(as.numeric(tunedmodelyall), ncol = nbensemble))
300 print(p)
301 }
302 return(modelssvr)
303 }
304
305 #" Reads ENVI hdr file
306 #"
307 #" @param hdrpath Path of the hdr file
308 #"
309 #" @return list of the content of the hdr file
310 #" @export
311 read_envi_header <- function(hdrpath) {
312 if (!grepl(".hdr$", hdrpath)) {
313 stop("File extension should be .hdr")
314 }
315 hdr <- readLines(hdrpath)
316 ## check ENVI at beginning of file
317 if (!grepl("ENVI", hdr[1])) {
318 stop("Not an ENVI header (ENVI keyword missing)")
319 } else {
320 hdr <- hdr [-1]
321 }
322 ## remove curly braces and put multi-line key-value-pairs into one line
323 hdr <- gsub("\\{([^}]*)\\}", "\\1", hdr)
324 l <- grep("\\{", hdr)
325 r <- grep("\\}", hdr)
326
327 if (length(l) != length(r)) {
328 stop("Error matching curly braces in header (differing numbers).")
329 }
330
331 if (any(r <= l)) {
332 stop("Mismatch of curly braces in header.")
333 }
334
335 hdr[l] <- sub("\\{", "", hdr[l])
336 hdr[r] <- sub("\\}", "", hdr[r])
337
338 for (i in rev(seq_along(l))) {
339 hdr <- c(
340 hdr [seq_len(l [i] - 1)],
341 paste(hdr [l [i]:r [i]], collapse = "\n"),
342 hdr [-seq_len(r [i])]
343 )
344 }
345
346 ## split key = value constructs into list with keys as names
347 hdr <- sapply(hdr, split_line, " = ", USE.NAMES = FALSE)
348 names(hdr) <- tolower(names(hdr))
349
350 ## process numeric values
351 tmp <- names(hdr) %in% c(
352 "samples", "lines", "bands", "header offset", "data type",
353 "byte order", "default bands", "data ignore value",
354 "wavelength", "fwhm", "data gain values"
355 )
356 hdr [tmp] <- lapply(hdr [tmp], function(x) {
357 as.numeric(unlist(strsplit(x, ",")))
358 })
359
360 return(hdr)
361 }
362
363 #" ENVI functions
364 #"
365 #" based on https: / / github.com / cran / hyperSpec / blob / master / R / read.ENVI.R
366 #" added wavelength, fwhm, ... to header reading
367 #" Title
368 #"
369 #" @param x character.
370 #" @param separator character
371 #" @param trim_blank boolean.
372 #"
373 #" @return list.
374 #" @export
375 split_line <- function(x, separator, trim_blank = TRUE) {
376 tmp <- regexpr(separator, x)
377 key <- substr(x, 1, tmp - 1)
378 value <- substr(x, tmp + 1, nchar(x))
379 if (trim_blank) {
380 blank_pattern <- "^[[:blank:]]*([^[:blank:]] + .*[^[:blank:]] + )[[:blank:]]*$"
381 key <- sub(blank_pattern, "\\1", key)
382 value <- sub(blank_pattern, "\\1", value)
383 }
384 value <- as.list(value)
385 names(value) <- key
386 return(value)
387 }
388
389 #" This function performs full training for hybrid invrsion using SVR with
390 #" values for default parameters
391 #"
392 #" @param minval list. minimum value for input parameters sampled to produce a training LUT
393 #" @param maxval list. maximum value for input parameters sampled to produce a training LUT
394 #" @param typedistrib list. Type of distribution. Either "Uniform" or "Gaussian"
395 #" @param gaussiandistrib list. Mean value and STD corresponding to the parameters sampled with gaussian distribution
396 #" @param parmset list. list of input parameters set to a specific value
397 #" @param nbsamples numeric. number of samples in training LUT
398 #" @param nbsamplesperrun numeric. number of training sample per individual regression model
399 #" @param nbmodels numeric. number of individual models to be run for ensemble
400 #" @param replacement bolean. is there replacement in subsampling?
401 #" @param sailversion character. Either 4SAIL or 4SAIL2
402 #" @param parms2estimate list. list of input parameters to be estimated
403 #" @param bands2select list. list of bands used for regression for each input parameter
404 #" @param noiselevel list. list of noise value added to reflectance (defined per input parm)
405 #" @param specprospect list. Includes optical constants required for PROSPECT
406 #" @param specsoil list. Includes either dry soil and wet soil, or a unique soil sample if the psoil parameter is not inverted
407 #" @param specatm list. Includes direct and diffuse radiation for clear conditions
408 #" @param path_results character. path for results
409 #" @param figplot boolean. Set TRUE to get scatterplot of estimated biophysical variable during training step
410 #" @param force4lowlai boolean. Set TRUE to artificially reduce leaf chemical constituent content for low LAI
411 #"
412 #"
413 #" @return modelssvr list. regression models trained for the retrieval of inputvar based on brf_lut
414 #" @export
415
416 train_prosail_inversion <- function(minval = NULL, maxval = NULL,
417 typedistrib = NULL, gaussiandistrib = NULL, parmset = NULL,
418 nbsamples = 2000, nbsamplesperrun = 100, nbmodels = 20, replacement = TRUE,
419 sailversion = "4SAIL",
420 parms2estimate = "lai", bands2select = NULL, noiselevel = NULL,
421 specprospect = NULL, specsoil = NULL, specatm = NULL,
422 path_results = "./", figplot = FALSE, force4lowlai = TRUE) {
423
424 ###===================================================================###
425 ### 1- PRODUCE A LUT TO TRAIN THE HYBRID INVERSION ###
426 ###===================================================================###
427 # Define sensor characteristics
428 if (is.null(specprospect)) {
429 specprospect <- prosail::specprospect
430 }
431 if (is.null(specsoil)) {
432 specsoil <- prosail::specsoil
433 }
434 if (is.null(specprospect)) {
435 specatm <- prosail::specatm
436 }
437 # define distribution for parameters to be sampled
438 if (is.null(typedistrib)) {
439 typedistrib <- data.frame("CHL" = "Uniform", "CAR" = "Uniform", "EWT" = "Uniform", "ANT" = "Uniform", "LMA" = "Uniform", "N" = "Uniform", "BROWN" = "Uniform",
440 "psoil" = "Uniform", "LIDFa" = "Uniform", "lai" = "Uniform", "q" = "Uniform", "tto" = "Uniform", "tts" = "Uniform", "psi" = "Uniform")
441 }
442 if (is.null(gaussiandistrib)) {
443 gaussiandistrib <- list("Mean" = NULL, "Std" = NULL)
444 }
445 if (is.null(minval)) {
446 minval <- data.frame("CHL" = 10, "CAR" = 0, "EWT" = 0.01, "ANT" = 0, "LMA" = 0.005, "N" = 1.0, "psoil" = 0.0, "BROWN" = 0.0,
447 "LIDFa" = 20, "lai" = 0.5, "q" = 0.1, "tto" = 0, "tts" = 20, "psi" = 80)
448 }
449 if (is.null(maxval)) {
450 maxval <- data.frame("CHL" = 75, "CAR" = 15, "EWT" = 0.03, "ANT" = 2, "LMA" = 0.03, "N" = 2.0, "psoil" = 1.0, "BROWN" = 0.5,
451 "LIDFa" = 70, "lai" = 7, "q" = 0.2, "tto" = 5, "tts" = 30, "psi" = 110)
452 }
453 # define min and max values
454 # fixed parameters
455 if (is.null(parmset)) {
456 parmset <- data.frame("TypeLidf" = 2, "alpha" = 40)
457 }
458 # produce input parameters distribution
459 if (sailversion == "4SAIL") {
460 inputprosail <- get_distribution_input_prosail(minval, maxval, parmset, nbsamples,
461 typedistrib = typedistrib,
462 Mean = gaussiandistrib$Mean, Std = gaussiandistrib$Std,
463 force4lowlai = force4lowlai)
464 } else if (sailversion == "4SAIL2") {
465 inputprosail <- get_distribution_input_prosail2(minval, maxval, parmset, nbsamples,
466 typedistrib = typedistrib,
467 Mean = gaussiandistrib$Mean, Std = gaussiandistrib$Std,
468 force4lowlai = force4lowlai)
469 }
470 if (sailversion == "4SAIL2") {
471 # Definition of Cv && update LAI
472 maxlai <- min(c(maxval$lai), 4)
473 inputprosail$Cv <- NA * inputprosail$lai
474 inputprosail$Cv[which(inputprosail$lai > maxlai)] <- 1
475 inputprosail$Cv[which(inputprosail$lai <= maxlai)] <- (1 / maxlai) + inputprosail$lai[which(inputprosail$lai <= maxlai)] / (maxlai + 1)
476 inputprosail$Cv <- inputprosail$Cv * matrix(rnorm(length(inputprosail$Cv), mean = 1, sd = 0.1))
477 inputprosail$Cv[which(inputprosail$Cv < 0)] <- 0
478 inputprosail$Cv[which(inputprosail$Cv > 1)] <- 1
479 inputprosail$Cv[which(inputprosail$lai > maxlai)] <- 1
480 inputprosail$fraction_brown <- 0 + 0 * inputprosail$lai
481 inputprosail$diss <- 0 + 0 * inputprosail$lai
482 inputprosail$Zeta <- 0.2 + 0 * inputprosail$lai
483 inputprosail$lai <- inputprosail$lai * inputprosail$Cv
484 }
485
486 # generate LUT of BRF corresponding to inputprosail, for a sensor
487 brf_lut <- Generate_LUT_BRF(sailversion = sailversion, inputprosail = inputprosail,
488 specprospect = specprospect, specsoil = specsoil, specatm = specatm)
489
490 # write parameters LUT
491 output <- matrix(unlist(inputprosail), ncol = length(inputprosail), byrow = FALSE)
492 filename <- file.path(path_results, "PROSAIL_LUT_InputParms.txt")
493 write.table(x = format(output, digits = 3), file = filename, append = FALSE, quote = FALSE,
494 col.names = names(inputprosail), row.names = FALSE, sep = "\t")
495 # Write BRF LUT corresponding to parameters LUT
496 filename <- file.path(path_results, "PROSAIL_LUT_reflectance.txt")
497 write.table(x = format(t(brf_lut), digits = 5), file = filename, append = FALSE, quote = FALSE,
498 col.names = specprospect$lambda, row.names = FALSE, sep = "\t")
499
500 # Which bands will be used for inversion?
501 if (is.null(bands2select)) {
502 bands2select <- list()
503 for (parm in parms2estimate) {
504 bands2select[[parm]] <- seq(1, length(specprospect$lambda))
505 }
506 }
507 # Add gaussian noise to reflectance LUT: one specific LUT per parameter
508 if (is.null(noiselevel)) {
509 noiselevel <- list()
510 for (parm in parms2estimate) {
511 noiselevel[[parm]] <- 0.01
512 }
513 }
514
515 # produce LIT with noise
516 brf_lut_noise <- list()
517 for (parm in parms2estimate) {
518 brf_lut_noise[[parm]] <- brf_lut[bands2select[[parm]], ] + brf_lut[bands2select[[parm]], ] * matrix(rnorm(nrow(brf_lut[bands2select[[parm]], ]) * ncol(brf_lut[bands2select[[parm]], ]),
519 0, noiselevel[[parm]]), nrow = nrow(brf_lut[bands2select[[parm]], ]))
520 }
521
522 ###===================================================================###
523 ### PERFORM HYBRID INVERSION ###
524 ###===================================================================###
525 # train SVR for each variable and each run
526 modelsvr <- list()
527 for (parm in parms2estimate) {
528 colparm <- which(parm == names(inputprosail))
529 inputvar <- inputprosail[[colparm]]
530 modelsvr[[parm]] <- prosail_hybrid_train(brf_lut = brf_lut_noise[[parm]], inputvar = inputvar,
531 figplot = figplot, nbensemble = nbmodels, withreplacement = replacement)
532 }
533 return(modelsvr)
534 }
535
536 #" writes ENVI hdr file
537 #"
538 #" @param hdr content to be written
539 #" @param hdrpath Path of the hdr file
540 #"
541 #" @return None
542 #" @importFrom stringr str_count
543 #" @export
544 write_envi_header <- function(hdr, hdrpath) {
545 h <- lapply(hdr, function(x) {
546 if (length(x) > 1 || (is.character(x) && stringr::str_count(x, "\\w + ") > 1)) {
547 x <- paste0("{", paste(x, collapse = ","), "}")
548 }
549 # convert last numerics
550 x <- as.character(x)
551 })
552 writeLines(c("ENVI", paste(names(hdr), h, sep = "=")), con = hdrpath)
553 return(invisible())
554 }