0
|
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 }
|