0
|
1 # ============================================================================= =
|
|
2 # prosail
|
|
3 # Lib_PROSAIL.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 simulation
|
|
11 # SAIL versions available are 4SAIL and 4SAIL2
|
|
12 # ============================================================================= =
|
|
13
|
|
14 #" computes bidirectional reflectance factor based on outputs from PROSAIL and sun position
|
|
15 #"
|
|
16 #" The direct and diffuse light are taken into account as proposed by:
|
|
17 #" Francois et al. (2002) conversion of 400-1100 nm vegetation albedo
|
|
18 #" measurements into total shortwave broadband albedo using a canopy
|
|
19 #" radiative transfer model, Agronomie
|
|
20 #" Es = direct
|
|
21 #" Ed = diffuse
|
|
22 #"
|
|
23 #" @param rdot numeric. Hemispherical-directional reflectance factor in viewing direction
|
|
24 #" @param rsot numeric. Bi-directional reflectance factor
|
|
25 #" @param tts numeric. Solar zenith angle
|
|
26 #" @param specatm_sensor list. direct and diffuse radiation for clear conditions
|
|
27 #" @return brf numeric. Bidirectional reflectance factor
|
|
28 #" @export
|
|
29 compute_brf <- function(rdot, rsot, tts, specatm_sensor) {
|
|
30
|
|
31 ############################## #
|
|
32 ## direct / diffuse light ##
|
|
33 ############################## #
|
|
34 es <- specatm_sensor$Direct_Light
|
|
35 ed <- specatm_sensor$Diffuse_Light
|
|
36 rd <- pi / 180
|
|
37 skyl <- 0.847 - 1.61 * sin((90 - tts) * rd) + 1.04 * sin((90 - tts) * rd) * sin((90 - tts) * rd) # diffuse radiation (Francois et al., 2002)
|
|
38 pardiro <- (1 - skyl) * es
|
|
39 pardifo <- skyl * ed
|
|
40 brf <- (rdot * pardifo + rsot * pardiro) / (pardiro + pardifo)
|
|
41 return(brf)
|
|
42 }
|
|
43
|
|
44 #" Performs PROSAIL simulation based on a set of combinations of input parameters
|
|
45 #" @param spec_sensor list. Includes optical constants required for PROSPECT
|
|
46 #" refractive index, specific absorption coefficients and corresponding spectral bands
|
|
47 #" @param input_prospect list. PROSPECT input variables
|
|
48 #" @param n numeric. Leaf structure parameter
|
|
49 #" @param chl numeric. chlorophyll content (microg.cm-2)
|
|
50 #" @param car numeric. carotenoid content (microg.cm-2)
|
|
51 #" @param ant numeric. anthocyain content (microg.cm-2)
|
|
52 #" @param brown numeric. brown pigment content (Arbitrary units)
|
|
53 #" @param ewt numeric. Equivalent Water Thickness (g.cm-2)
|
|
54 #" @param lma numeric. Leaf Mass per Area (g.cm-2)
|
|
55 #" @param prot numeric. protein content (g.cm-2)
|
|
56 #" @param cbc numeric. nonprotcarbon-based constituent content (g.cm-2)
|
|
57 #" @param alpha numeric. Solid angle for incident light at surface of leaf (simulation of roughness)
|
|
58 #" @param typelidf numeric. Type of leaf inclination distribution function
|
|
59 #" @param lidfa numeric.
|
|
60 #" if typelidf == 1, controls the average leaf slope
|
|
61 #" if typelidf == 2, corresponds to average leaf angle
|
|
62 #" @param lidfb numeric.
|
|
63 #" if typelidf == 1, unused
|
|
64 #" if typelidf == 2, controls the distribution"s bimodality
|
|
65 #" @param lai numeric. Leaf Area Index
|
|
66 #" @param q numeric. Hot Spot parameter
|
|
67 #" @param tts numeric. Sun zeith angle
|
|
68 #" @param tto numeric. Observer zeith angle
|
|
69 #" @param psi numeric. Azimuth Sun / Observer
|
|
70 #" @param rsoil numeric. Soil reflectance
|
|
71 #" @param fraction_brown numeric. Fraction of brown leaf area
|
|
72 #" @param diss numeric. Layer dissociation factor
|
|
73 #" @param cv numeric. vertical crown cover percentage
|
|
74 #" = % ground area covered with crowns as seen from nadir direction
|
|
75 #" @param zeta numeric. Tree shape factor
|
|
76 #" = ratio of crown diameter to crown height
|
|
77 #" @param sailversion character. choose between 4SAIL and 4SAIL2
|
|
78 #" @param brownvegetation list. Defines optical properties for brown vegetation, if not nULL
|
|
79 #" - WVL, reflectance, Transmittance
|
|
80 #" - Set to nULL if use PROSPECT to generate it
|
|
81 #"
|
|
82 #" @return list. rdot, rsot, rddt, rsdt
|
|
83 #" rdot: hemispherical-directional reflectance factor in viewing direction
|
|
84 #" rsot: bi-directional reflectance factor
|
|
85 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
|
|
86 #" rddt: bi-hemispherical reflectance factor
|
|
87 #" @import prospect
|
|
88 #" @export
|
|
89 pro4sail <- function(spec_sensor, input_prospect = nULL, n = 1.5, chl = 40.0,
|
|
90 car = 8.0, ant = 0.0, brown = 0.0, ewt = 0.01,
|
|
91 lma = 0.008, prot = 0.0, cbc = 0.0, alpha = 40.0,
|
|
92 typelidf = 2, lidfa = nULL, lidfb = nULL, lai = nULL,
|
|
93 q = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
|
|
94 fraction_brown = 0.0, diss = 0.0, cv = 1, zeta = 1,
|
|
95 sailversion = "4SAIL", brownvegetation = nULL) {
|
|
96
|
|
97 ############################ #
|
|
98 # LEAF OPTICAL PROPERTIES ##
|
|
99 ############################ #
|
|
100 if (is.null(input_prospect)) {
|
|
101 input_prospect <- data.frame("chl" = chl, "car" = car, "ant" = ant, "brown" = brown, "ewt" = ewt,
|
|
102 "lma" = lma, "prot" = prot, "cbc" = cbc, "n" = n, "alpha" = alpha)
|
|
103 }
|
|
104 greenvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
|
|
105 n = input_prospect$n[1],
|
|
106 chl = input_prospect$chl[1],
|
|
107 car = input_prospect$car[1],
|
|
108 ant = input_prospect$ant[1],
|
|
109 brown = input_prospect$brown[1],
|
|
110 ewt = input_prospect$ewt[1],
|
|
111 lma = input_prospect$lma[1],
|
|
112 prot = input_prospect$prot[1],
|
|
113 cbc = input_prospect$cbc[1],
|
|
114 alpha = input_prospect$alpha[1])
|
|
115
|
|
116 if (sailversion == "4SAIL2") {
|
|
117 # 4SAIL2 requires one of the following combination of input parameters
|
|
118 # Case #1: valid optical properties for brown vegetation
|
|
119 if (!is.null(brownvegetation)) {
|
|
120 # need to define reflectance and Transmittance for brownvegetation
|
|
121 if (length(grep("reflectance", names(brownvegetation))) == 0 || length(grep("Transmittance", names(brownvegetation))) == 0) {
|
|
122 message("Please define brownvegetation as a list including reflectance and Transmittance")
|
|
123 stop()
|
|
124 }
|
|
125 # check if spectral domain for optical properties of brown vegetation match
|
|
126 # with spectral domain for optical properties of green vegetation
|
|
127 if (length(setdiff(spec_sensor$lambda, brownvegetation$wvl)) > 0) {
|
|
128 message("Please define same spectral domain for brownvegetation and SpecPROSPECT")
|
|
129 stop()
|
|
130 }
|
|
131 if (length(unique(lengths(input_prospect))) == 1) {
|
|
132 if (!unique(lengths(input_prospect)) == 1) {
|
|
133 message("brownvegetation defined along with multiple leaf chemical properties")
|
|
134 message("Only first set of leaf chemical properties will be used to simulate green vegetation")
|
|
135 }
|
|
136 }
|
|
137 # if no leaf optical properties brown vegetation defined
|
|
138 } else if (is.null(brownvegetation)) {
|
|
139 # if all PROSPECT input parameters have the same length
|
|
140 if (length(unique(lengths(input_prospect))) == 1) {
|
|
141 # if all PROSPECT input parameters are unique (no possibility to simulate 2 types of leaf optics)
|
|
142 if (unique(lengths(input_prospect)) == 1) {
|
|
143 # if fraction_brown set to 0, then assign green vegetation optics to brown vegetation optics
|
|
144 if (fraction_brown == 0) {
|
|
145 brownvegetation <- greenvegetation
|
|
146 # else run 4SAIL
|
|
147 } else {
|
|
148 message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
|
|
149 message("Currently one set is defined. will run 4SAIL instead of 4SAIL2")
|
|
150 sailversion <- "4SAIL"
|
|
151 }
|
|
152 # if all PROSPECT parameters have at least 2 elements
|
|
153 } else if (unique(lengths(input_prospect)) >= 2) {
|
|
154 # compute leaf optical properties
|
|
155 brownvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
|
|
156 n = input_prospect$n[2],
|
|
157 chl = input_prospect$chl[2],
|
|
158 car = input_prospect$car[2],
|
|
159 ant = input_prospect$ant[2],
|
|
160 brown = input_prospect$brown[2],
|
|
161 ewt = input_prospect$ewt[2],
|
|
162 lma = input_prospect$lma[2],
|
|
163 prot = input_prospect$prot[2],
|
|
164 cbc = input_prospect$cbc[2],
|
|
165 alpha = input_prospect$alpha[2])
|
|
166 if (unique(lengths(input_prospect)) > 2) {
|
|
167 message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
|
|
168 message("Currently more than 2 sets are defined. will only use the first 2")
|
|
169 }
|
|
170 }
|
|
171 }
|
|
172 }
|
|
173 }
|
|
174 if (sailversion == "4SAIL") {
|
|
175 if (length(unique(lengths(input_prospect))) == 1) {
|
|
176 if (unique(lengths(input_prospect)) > 1) {
|
|
177 message("4SAIL needs only one set of optical properties")
|
|
178 message("Currently more than one set of leaf chemical constituents is defined.")
|
|
179 message("Will run 4SAIL with the first set of leaf chemical constituents")
|
|
180 }
|
|
181 }
|
|
182 }
|
|
183
|
|
184 if (sailversion == "4SAIL") {
|
|
185 # run 4SAIL
|
|
186 ref <- foursail(leafoptics = greenvegetation,
|
|
187 typelidf, lidfa, lidfb, lai, q, tts, tto, psi, rsoil)
|
|
188 } else if (sailversion == "4SAIL2") {
|
|
189 # run 4SAIL2
|
|
190 ref <- foursail2(leafgreen = greenvegetation, leafbrown = brownvegetation,
|
|
191 typelidf, lidfa, lidfb, lai, q, tts, tto, psi, rsoil,
|
|
192 fraction_brown, diss, cv, zeta)
|
|
193 }
|
|
194 return(ref)
|
|
195 }
|
|
196
|
|
197 #" Performs PROSAIL simulation based on a set of combinations of input parameters
|
|
198 #" @param leafoptics list. Includes leaf optical properties (reflectance and transmittance)
|
|
199 #" and corresponding spectral bands
|
|
200 #" @param typelidf numeric. Type of leaf inclination distribution function
|
|
201 #" @param lidfa numeric.
|
|
202 #" if typelidf == 1, controls the average leaf slope
|
|
203 #" if typelidf == 2, corresponds to average leaf angle
|
|
204 #" @param lidfb numeric.
|
|
205 #" if typelidf == 1, unused
|
|
206 #" if typelidf == 2, controls the distribution"s bimodality
|
|
207 #" @param lai numeric. Leaf Area Index
|
|
208 #" @param q numeric. Hot Spot parameter
|
|
209 #" @param tts numeric. Sun zeith angle
|
|
210 #" @param tto numeric. Observer zeith angle
|
|
211 #" @param psi numeric. Azimuth Sun / Observer
|
|
212 #" @param rsoil numeric. Soil reflectance
|
|
213 #"
|
|
214 #" @return list. rdot, rsot, rddt, rsdt
|
|
215 #" rdot: hemispherical-directional reflectance factor in viewing direction
|
|
216 #" rsot: bi-directional reflectance factor
|
|
217 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
|
|
218 #" rddt: bi-hemispherical reflectance factor
|
|
219 #" @export
|
|
220
|
|
221 foursail <- function(leafoptics, typelidf = 2, lidfa = nULL, lidfb = nULL, lai = nULL,
|
|
222 q = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL) {
|
|
223
|
|
224 ############################## #
|
|
225 # LEAF OPTICAL PROPERTIES ##
|
|
226 ############################## #
|
|
227 rho <- leafoptics$Reflectance
|
|
228 tau <- leafoptics$Transmittance
|
|
229
|
|
230 # Geometric quantities
|
|
231 rd <- pi / 180
|
|
232 cts <- cos(rd * tts)
|
|
233 cto <- cos(rd * tto)
|
|
234 ctscto <- cts * cto
|
|
235 tants <- tan(rd * tts)
|
|
236 tanto <- tan(rd * tto)
|
|
237 cospsi <- cos(rd * psi)
|
|
238 dso <- sqrt(tants * tants + tanto * tanto - 2. * tants * tanto * cospsi)
|
|
239
|
|
240 # Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
|
|
241 if (typelidf == 1) {
|
|
242 foliar_distrib <- dladgen(lidfa, lidfb)
|
|
243 lidf <- foliar_distrib$lidf
|
|
244 litab <- foliar_distrib$litab
|
|
245
|
|
246 } else if (typelidf == 2) {
|
|
247 foliar_distrib <- campbell(lidfa)
|
|
248 lidf <- foliar_distrib$lidf
|
|
249 litab <- foliar_distrib$litab
|
|
250 }
|
|
251
|
|
252 # angular distance, compensation of shadow length
|
|
253 # Calculate geometric factors associated with extinction and scattering
|
|
254 # Initialise sums
|
|
255 ks <- 0
|
|
256 ko <- 0
|
|
257 bf <- 0
|
|
258 sob <- 0
|
|
259 sof <- 0
|
|
260
|
|
261 # Weighted sums over LIDF
|
|
262 na <- length(litab)
|
|
263 for (i in 1:na) {
|
|
264 ttl <- litab[i] # leaf inclination discrete values
|
|
265 ctl <- cos(rd * ttl)
|
|
266 # SAIL volume scattering phase function gives interception and portions to be
|
|
267 # multiplied by rho and tau
|
|
268 resvolscatt <- volscatt(tts, tto, psi, ttl)
|
|
269 chi_s <- resvolscatt$chi_s
|
|
270 chi_o <- resvolscatt$chi_o
|
|
271 frho <- resvolscatt$frho
|
|
272 ftau <- resvolscatt$ftau
|
|
273
|
|
274 #********************************************************************************
|
|
275 #* SUITS SYSTEM coEFFICIEnTS
|
|
276 #*
|
|
277 #* ks : Extinction coefficient for direct solar flux
|
|
278 #* ko : Extinction coefficient for direct observed flux
|
|
279 #* att : Attenuation coefficient for diffuse flux
|
|
280 #* sigb : Backscattering coefficient of the diffuse downward flux
|
|
281 #* sigf : Forwardscattering coefficient of the diffuse upward flux
|
|
282 #* sf : Scattering coefficient of the direct solar flux for downward diffuse flux
|
|
283 #* sb : Scattering coefficient of the direct solar flux for upward diffuse flux
|
|
284 #* vf : Scattering coefficient of upward diffuse flux in the observed direction
|
|
285 #* vb : Scattering coefficient of downward diffuse flux in the observed direction
|
|
286 #* w : Bidirectional scattering coefficient
|
|
287 #********************************************************************************
|
|
288
|
|
289 # Extinction coefficients
|
|
290 ksli <- chi_s / cts
|
|
291 koli <- chi_o / cto
|
|
292
|
|
293 # Area scattering coefficient fractions
|
|
294 sobli <- frho * pi / ctscto
|
|
295 sofli <- ftau * pi / ctscto
|
|
296 bfli <- ctl * ctl
|
|
297 ks <- ks + ksli * lidf[i]
|
|
298 ko <- ko + koli * lidf[i]
|
|
299 bf <- bf + bfli * lidf[i]
|
|
300 sob <- sob + sobli * lidf[i]
|
|
301 sof <- sof + sofli * lidf[i]
|
|
302 }
|
|
303
|
|
304 # Geometric factors to be used later with rho and tau
|
|
305 sdb <- 0.5 * (ks + bf)
|
|
306 sdf <- 0.5 * (ks - bf)
|
|
307 dob <- 0.5 * (ko + bf)
|
|
308 dof <- 0.5 * (ko - bf)
|
|
309 ddb <- 0.5 * (1. + bf)
|
|
310 ddf <- 0.5 * (1. - bf)
|
|
311
|
|
312 # Here rho and tau come in
|
|
313 sigb <- ddb * rho + ddf * tau
|
|
314 sigf <- ddf * rho + ddb * tau
|
|
315 att <- 1 - sigf
|
|
316 m2 <- (att + sigb) * (att - sigb)
|
|
317 m2[which(m2 <= 0)] <- 0
|
|
318 m <- sqrt(m2)
|
|
319
|
|
320 sb <- sdb * rho + sdf * tau
|
|
321 sf <- sdf * rho + sdb * tau
|
|
322 vb <- dob * rho + dof * tau
|
|
323 vf <- dof * rho + dob * tau
|
|
324 w <- sob * rho + sof * tau
|
|
325
|
|
326 # Here the LAI comes in
|
|
327 # Outputs for the case LAI = 0
|
|
328 if (lai < 0) {
|
|
329 tss <- 1
|
|
330 too <- 1
|
|
331 tsstoo <- 1
|
|
332 rdd <- 0
|
|
333 tdd <- 1
|
|
334 rsd <- 0
|
|
335 tsd <- 0
|
|
336 rdo <- 0
|
|
337 tdo <- 0
|
|
338 rso <- 0
|
|
339 rsos <- 0
|
|
340 rsod <- 0
|
|
341
|
|
342 rddt <- rsoil
|
|
343 rsdt <- rsoil
|
|
344 rdot <- rsoil
|
|
345 rsodt <- 0 * rsoil
|
|
346 rsost <- rsoil
|
|
347 rsot <- rsoil
|
|
348 } else {
|
|
349 # Other cases (LAI > 0)
|
|
350 e1 <- exp(-m * lai)
|
|
351 e2 <- e1 * e1
|
|
352 rinf <- (att - m) / sigb
|
|
353 rinf2 <- rinf * rinf
|
|
354 re <- rinf * e1
|
|
355 denom <- 1. - rinf2 * e2
|
|
356
|
|
357 j1ks <- jfunc1(ks, m, lai)
|
|
358 j2ks <- jfunc2(ks, m, lai)
|
|
359 j1ko <- jfunc1(ko, m, lai)
|
|
360 j2ko <- jfunc2(ko, m, lai)
|
|
361
|
|
362 ps <- (sf + sb * rinf) * j1ks
|
|
363 qs <- (sf * rinf + sb) * j2ks
|
|
364 pv <- (vf + vb * rinf) * j1ko
|
|
365 qv <- (vf * rinf + vb) * j2ko
|
|
366
|
|
367 rdd <- rinf * (1. - e2) / denom
|
|
368 tdd <- (1. - rinf2) * e1 / denom
|
|
369 tsd <- (ps - re * qs) / denom
|
|
370 rsd <- (qs - re * ps) / denom
|
|
371 tdo <- (pv - re * qv) / denom
|
|
372 rdo <- (qv - re * pv) / denom
|
|
373
|
|
374 tss <- exp(-ks * lai)
|
|
375 too <- exp(-ko * lai)
|
|
376 z <- jfunc3(ks, ko, lai)
|
|
377 g1 <- (z - j1ks * too) / (ko + m)
|
|
378 g2 <- (z - j1ko * tss) / (ks + m)
|
|
379
|
|
380 tv1 <- (vf * rinf + vb) * g1
|
|
381 tv2 <- (vf + vb * rinf) * g2
|
|
382 t1 <- tv1 * (sf + sb * rinf)
|
|
383 t2 <- tv2 * (sf * rinf + sb)
|
|
384 t3 <- (rdo * qs + tdo * ps) * rinf
|
|
385
|
|
386 # Multiple scattering contribution to bidirectional canopy reflectance
|
|
387 rsod <- (t1 + t2 - t3) / (1. - rinf2)
|
|
388
|
|
389 # Treatment of the hotspot-effect
|
|
390 alf <- 1e6
|
|
391 # Apply correction 2 / (K + k) suggested by F.-M. Breon
|
|
392 if (q > 0) {
|
|
393 alf <- (dso / q) * 2. / (ks + ko)
|
|
394 }
|
|
395 if (alf > 200) {
|
|
396 # inserted H. Bach 1 / 3 / 04
|
|
397 alf <- 200
|
|
398 }
|
|
399 if (alf == 0) {
|
|
400 # The pure hotspot - no shadow
|
|
401 tsstoo <- tss
|
|
402 sumint <- (1 - tss) / (ks * lai)
|
|
403 } else {
|
|
404 # Outside the hotspot
|
|
405 fhot <- lai * sqrt(ko * ks)
|
|
406 # Integrate by exponential Simpson method in 20 steps
|
|
407 # the steps are arranged according to equal partitioning
|
|
408 # of the slope of the joint probability function
|
|
409 x1 <- 0
|
|
410 y1 <- 0
|
|
411 f1 <- 1
|
|
412 fint <- (1. - exp(-alf)) * 0.05
|
|
413 sumint <- 0
|
|
414 for (i in 1:20) {
|
|
415 if (i < 20) {
|
|
416 x2 <- -log(1. - i * fint) / alf
|
|
417 } else {
|
|
418 x2 <- 1
|
|
419 }
|
|
420 y2 <- -(ko + ks) * lai * x2 + fhot * (1. - exp(-alf * x2)) / alf
|
|
421 f2 <- exp(y2)
|
|
422 sumint <- sumint + (f2 - f1) * (x2 - x1) / (y2 - y1)
|
|
423 x1 <- x2
|
|
424 y1 <- y2
|
|
425 f1 <- f2
|
|
426 }
|
|
427 tsstoo <- f1
|
|
428 }
|
|
429 # Bidirectional reflectance
|
|
430 # Single scattering contribution
|
|
431 rsos <- w * lai * sumint
|
|
432 # Total canopy contribution
|
|
433 rso <- rsos + rsod
|
|
434 # Interaction with the soil
|
|
435 dn <- 1. - rsoil * rdd
|
|
436 # rddt: bi-hemispherical reflectance factor
|
|
437 rddt <- rdd + tdd * rsoil * tdd / dn
|
|
438 # rsdt: directional-hemispherical reflectance factor for solar incident flux
|
|
439 rsdt <- rsd + (tsd + tss) * rsoil * tdd / dn
|
|
440 # rdot: hemispherical-directional reflectance factor in viewing direction
|
|
441 rdot <- rdo + tdd * rsoil * (tdo + too) / dn
|
|
442 # rsot: bi-directional reflectance factor
|
|
443 rsodt <- rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
|
|
444 rsost <- rsos + tsstoo * rsoil
|
|
445 rsot <- rsost + rsodt
|
|
446 }
|
|
447 my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt)
|
|
448 return(my_list)
|
|
449 }
|
|
450
|
|
451 #" Performs pro4sail2 simulation based on a set of combinations of input parameters
|
|
452 #" @param leafgreen list. includes relfectance and transmittance for vegetation #1 (e.g. green vegetation)
|
|
453 #" @param leafbrown list. includes relfectance and transmittance for vegetation #2 (e.g. brown vegetation)
|
|
454 #" @param typelidf numeric. Type of leaf inclination distribution function
|
|
455 #" @param lidfa numeric.
|
|
456 #" if typelidf == 1, controls the average leaf slope
|
|
457 #" if typelidf == 2, corresponds to average leaf angle
|
|
458 #" @param lidfb numeric.
|
|
459 #" if typelidf == 1, unused
|
|
460 #" if typelidf == 2, controls the distribution"s bimodality
|
|
461 #" @param lai numeric. Leaf Area Index
|
|
462 #" @param hot numeric. Hot Spot parameter = ratio of the correlation length of leaf projections in the horizontal plane and the canopy height (doi:10.1016 / j.rse.2006.12.013)
|
|
463 #" @param tts numeric. Sun zeith angle
|
|
464 #" @param tto numeric. Observer zeith angle
|
|
465 #" @param psi numeric. Azimuth Sun / Observer
|
|
466 #" @param rsoil numeric. Soil reflectance
|
|
467 #" @param fraction_brown numeric. Fraction of brown leaf area
|
|
468 #" @param diss numeric. Layer dissociation factor
|
|
469 #" @param cv numeric. vertical crown cover percentage
|
|
470 #" = % ground area covered with crowns as seen from nadir direction
|
|
471 #" @param zeta numeric. Tree shape factor
|
|
472 #" = ratio of crown diameter to crown height
|
|
473 #"
|
|
474 #" @return list. rdot, rsot, rddt, rsdt
|
|
475 #" rdot: hemispherical-directional reflectance factor in viewing direction
|
|
476 #" rsot: bi-directional reflectance factor
|
|
477 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
|
|
478 #" rddt: bi-hemispherical reflectance factor
|
|
479 #" alfast: canopy absorptance for direct solar incident flux
|
|
480 #" alfadt: canopy absorptance for hemispherical diffuse incident flux
|
|
481 #" @export
|
|
482
|
|
483 foursail2 <- function(leafgreen, leafbrown,
|
|
484 typelidf = 2, lidfa = nULL, lidfb = nULL,
|
|
485 lai = nULL, hot = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
|
|
486 fraction_brown = 0.5, diss = 0.5, cv = 1, zeta = 1) {
|
|
487
|
|
488 # This version does not include non-Lambertian soil properties.
|
|
489 # original codes do, and only need to add the following variables as input
|
|
490 rddsoil <- rdosoil <- rsdsoil <- rsosoil <- rsoil
|
|
491
|
|
492 # Geometric quantities
|
|
493 rd <- pi / 180
|
|
494
|
|
495 # Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
|
|
496 if (typelidf == 1) {
|
|
497 foliar_distrib <- dladgen(lidfa, lidfb)
|
|
498 lidf <- foliar_distrib$lidf
|
|
499 litab <- foliar_distrib$litab
|
|
500
|
|
501 } else if (typelidf == 2) {
|
|
502 foliar_distrib <- campbell(lidfa)
|
|
503 lidf <- foliar_distrib$lidf
|
|
504 litab <- foliar_distrib$litab
|
|
505 }
|
|
506
|
|
507 if (lai < 0) {
|
|
508 message("Please define positive LAI value")
|
|
509 rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
|
|
510 alfast <- alfadt <- 0 * rsoil
|
|
511 } else if (lai == 0) {
|
|
512 tss <- too <- tsstoo <- tdd <- 1.0
|
|
513 rdd <- rsd <- tsd <- rdo <- tdo <- 0.0
|
|
514 rso <- rsos <- rsod <- rsodt <- 0.0
|
|
515 rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
|
|
516 alfast <- alfadt <- 0 * rsoil
|
|
517 } else if (lai > 0) {
|
|
518 cts <- cos(rd * tts)
|
|
519 cto <- cos(rd * tto)
|
|
520 ctscto <- cts * cto
|
|
521 tants <- tan(rd * tts)
|
|
522 tanto <- tan(rd * tto)
|
|
523 cospsi <- cos(rd * psi)
|
|
524 dso <- sqrt(tants * tants + tanto * tanto - 2.0 * tants * tanto * cospsi)
|
|
525
|
|
526 # Clumping effects
|
|
527 cs <- co <- 1.0
|
|
528 if (cv <= 1.0) {
|
|
529 cs <- 1.0 - (1.0 - cv)^(1.0 / cts)
|
|
530 co <- 1.0 - (1.0 - cv)^(1.0 / cto)
|
|
531 }
|
|
532 overlap <- 0.0
|
|
533 if (zeta > 0.0) {
|
|
534 overlap <- min(cs * (1.0 - co), co * (1.0 - cs)) * exp(-dso / zeta)
|
|
535 }
|
|
536 fcd <- cs * co + overlap
|
|
537 fcs <- (1.0 - cs) * co - overlap
|
|
538 fod <- cs * (1.0 - co) - overlap
|
|
539 fos <- (1.0 - cs) * (1.0 - co) + overlap
|
|
540 fcdc <- 1.0 - (1.0 - fcd)^(0.5 / cts + 0.5 / cto)
|
|
541
|
|
542 # Part depending on diss, fraction_brown, and leaf optical properties
|
|
543 # First save the input fraction_brown as the old fraction_brown, as the following change is only artificial
|
|
544 # Better define an fraction_brown that is actually used: fb, so that the input is not modified!
|
|
545
|
|
546 fb <- fraction_brown
|
|
547 # if only green leaves
|
|
548 if (fraction_brown == 0.0) {
|
|
549 fb <- 0.5
|
|
550 leafbrown$Reflectance <- leafgreen$Reflectance
|
|
551 leafbrown$Transmittance <- leafgreen$Transmittance
|
|
552 }
|
|
553 if (fraction_brown == 1.0) {
|
|
554 fb <- 0.5
|
|
555 leafgreen$Reflectance <- leafbrown$Reflectance
|
|
556 leafgreen$Transmittance <- leafbrown$Transmittance
|
|
557 }
|
|
558 s <- (1.0 - diss) * fb * (1.0 - fb)
|
|
559 # rho1 && tau1 : green foliage
|
|
560 # rho2 && tau2 : brown foliage (bottom layer)
|
|
561 rho1 <- ((1 - fb - s) * leafgreen$Reflectance + s * leafbrown$Reflectance) / (1 - fb)
|
|
562 tau1 <- ((1 - fb - s) * leafgreen$Transmittance + s * leafbrown$Transmittance) / (1 - fb)
|
|
563 rho2 <- (s * leafgreen$Reflectance + (fb - s) * leafbrown$Reflectance) / fb
|
|
564 tau2 <- (s * leafgreen$Transmittance + (fb - s) * leafbrown$Transmittance) / fb
|
|
565
|
|
566 # angular distance, compensation of shadow length
|
|
567 # Calculate geometric factors associated with extinction and scattering
|
|
568 # Initialise sums
|
|
569 ks <- ko <- bf <- sob <- sof <- 0
|
|
570
|
|
571 # Weighted sums over LIDF
|
|
572
|
|
573 for (i in 1:seq_along(litab)) {
|
|
574 ttl <- litab[i]
|
|
575 ctl <- cos(rd * ttl)
|
|
576 # SAIL volscatt function gives interception coefficients
|
|
577 # and two portions of the volume scattering phase function to be
|
|
578 # multiplied by rho and tau, respectively
|
|
579 resvolscatt <- volscatt(tts, tto, psi, ttl)
|
|
580 chi_s <- resvolscatt$chi_s
|
|
581 chi_o <- resvolscatt$chi_o
|
|
582 frho <- resvolscatt$frho
|
|
583 ftau <- resvolscatt$ftau
|
|
584 # Extinction coefficients
|
|
585 ksli <- chi_s / cts
|
|
586 koli <- chi_o / cto
|
|
587 # Area scattering coefficient fractions
|
|
588 sobli <- frho * pi / ctscto
|
|
589 sofli <- ftau * pi / ctscto
|
|
590 bfli <- ctl * ctl
|
|
591 ks <- ks + ksli * lidf[i]
|
|
592 ko <- ko + koli * lidf[i]
|
|
593 bf <- bf + bfli * lidf[i]
|
|
594 sob <- sob + sobli * lidf[i]
|
|
595 sof <- sof + sofli * lidf[i]
|
|
596 }
|
|
597 # Geometric factors to be used later in combination with rho and tau
|
|
598 sdb <- 0.5 * (ks + bf)
|
|
599 sdf <- 0.5 * (ks - bf)
|
|
600 dob <- 0.5 * (ko + bf)
|
|
601 dof <- 0.5 * (ko - bf)
|
|
602 ddb <- 0.5 * (1. + bf)
|
|
603 ddf <- 0.5 * (1. - bf)
|
|
604
|
|
605 # LAIs in two layers
|
|
606 lai1 <- (1 - fb) * lai
|
|
607 lai2 <- fb * lai
|
|
608
|
|
609 tss <- exp(-ks * lai)
|
|
610 ck <- exp(-ks * lai1)
|
|
611 alf <- 1e6
|
|
612 if (hot > 0.0) {
|
|
613 alf <- (dso / hot) * 2.0 / (ks + ko)
|
|
614 }
|
|
615 if (alf > 200.0) {
|
|
616 alf <- 200.0 # inserted H. Bach 1 / 3 / 04
|
|
617 }
|
|
618 if (alf == 0.0) {
|
|
619 # The pure hotspot
|
|
620 tsstoo <- tss
|
|
621 s1 <- (1 - ck) / (ks * lai)
|
|
622 s2 <- (ck - tss) / (ks * lai)
|
|
623 } else {
|
|
624 # Outside the hotspot
|
|
625 fhot <- lai * sqrt(ko * ks)
|
|
626 # Integrate 2 layers by exponential simpson method in 20 steps
|
|
627 # the steps are arranged according to equal partitioning
|
|
628 # of the derivative of the joint probability function
|
|
629 x1 <- y1 <- 0.0
|
|
630 f1 <- 1.0
|
|
631 ca <- exp(alf * (fb - 1.0))
|
|
632 fint <- (1.0 - ca) * .05
|
|
633 s1 <- 0.0
|
|
634 for (istep in 1:20) {
|
|
635 if (istep < 20) {
|
|
636 x2 <- -log(1. - istep * fint) / alf
|
|
637 } else {
|
|
638 x2 <- 1. - fb
|
|
639 }
|
|
640 y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
|
|
641 f2 <- exp(y2)
|
|
642 s1 <- s1 + (f2 - f1) * (x2 - x1) / (y2 - y1)
|
|
643 x1 <- x2
|
|
644 y1 <- y2
|
|
645 f1 <- f2
|
|
646 }
|
|
647 fint <- (ca - exp(-alf)) * .05
|
|
648 s2 <- 0.0
|
|
649 for (istep in 1:20) {
|
|
650 if (istep < 20) {
|
|
651 x2 <- -log(ca - istep * fint) / alf
|
|
652 } else {
|
|
653 x2 <- 1.0
|
|
654 }
|
|
655 y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
|
|
656 f2 <- exp(y2)
|
|
657 s2 <- s2 + (f2 - f1) * (x2 - x1) / (y2 - y1)
|
|
658 x1 <- x2
|
|
659 y1 <- y2
|
|
660 f1 <- f2
|
|
661 }
|
|
662 tsstoo <- f1
|
|
663 }
|
|
664
|
|
665 # Calculate reflectances and transmittances
|
|
666 # Bottom layer
|
|
667 tss <- exp(-ks * lai2)
|
|
668 too <- exp(-ko * lai2)
|
|
669 sb <- sdb * rho2 + sdf * tau2
|
|
670 sf <- sdf * rho2 + sdb * tau2
|
|
671
|
|
672 vb <- dob * rho2 + dof * tau2
|
|
673 vf <- dof * rho2 + dob * tau2
|
|
674
|
|
675 w2 <- sob * rho2 + sof * tau2
|
|
676
|
|
677 sigb <- ddb * rho2 + ddf * tau2
|
|
678 sigf <- ddf * rho2 + ddb * tau2
|
|
679 att <- 1.0 - sigf
|
|
680 m2 <- (att + sigb) * (att - sigb)
|
|
681 m2[m2 < 0] <- 0
|
|
682 m <- sqrt(m2)
|
|
683 which_ncs <- which(m > 0.01)
|
|
684 which_cs <- which(m <= 0.01)
|
|
685
|
|
686 tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
|
|
687 rsod <- 0 * m
|
|
688 if (length(which_ncs) > 0) {
|
|
689 resncs <- nonconservativescattering(m[which_ncs], lai2, att[which_ncs], sigb[which_ncs],
|
|
690 ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
|
|
691 tdd[which_ncs] <- resncs$tdd
|
|
692 rdd[which_ncs] <- resncs$rdd
|
|
693 tsd[which_ncs] <- resncs$tsd
|
|
694 rsd[which_ncs] <- resncs$rsd
|
|
695 tdo[which_ncs] <- resncs$tdo
|
|
696 rdo[which_ncs] <- resncs$rdo
|
|
697 rsod[which_ncs] <- resncs$rsod
|
|
698 }
|
|
699 if (length(which_cs) > 0) {
|
|
700 rescs <- conservativescattering(m[which_cs], lai2, att[which_cs], sigb[which_cs],
|
|
701 ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
|
|
702 tdd[which_cs] <- rescs$tdd
|
|
703 rdd[which_cs] <- rescs$rdd
|
|
704 tsd[which_cs] <- rescs$tsd
|
|
705 rsd[which_cs] <- rescs$rsd
|
|
706 tdo[which_cs] <- rescs$tdo
|
|
707 rdo[which_cs] <- rescs$rdo
|
|
708 rsod[which_cs] <- rescs$rsod
|
|
709 }
|
|
710
|
|
711 # Set background properties equal to those of the bottom layer on a black soil
|
|
712 rddb <- rdd
|
|
713 rsdb <- rsd
|
|
714 rdob <- rdo
|
|
715 rsodb <- rsod
|
|
716 tddb <- tdd
|
|
717 tsdb <- tsd
|
|
718 tdob <- tdo
|
|
719 toob <- too
|
|
720 tssb <- tss
|
|
721 # Top layer
|
|
722 tss <- exp(-ks * lai1)
|
|
723 too <- exp(-ko * lai1)
|
|
724
|
|
725 sb <- sdb * rho1 + sdf * tau1
|
|
726 sf <- sdf * rho1 + sdb * tau1
|
|
727
|
|
728 vb <- dob * rho1 + dof * tau1
|
|
729 vf <- dof * rho1 + dob * tau1
|
|
730
|
|
731 w1 <- sob * rho1 + sof * tau1
|
|
732
|
|
733 sigb <- ddb * rho1 + ddf * tau1
|
|
734 sigf <- ddf * rho1 + ddb * tau1
|
|
735 att <- 1.0 - sigf
|
|
736
|
|
737 m2 <- (att + sigb) * (att - sigb)
|
|
738 m2[m2 < 0] <- 0
|
|
739 m <- sqrt(m2)
|
|
740 which_ncs <- which(m > 0.01)
|
|
741 which_cs <- which(m <= 0.01)
|
|
742
|
|
743 tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
|
|
744 rsod <- 0 * m
|
|
745 if (length(which_ncs) > 0) {
|
|
746 resncs <- nonconservativescattering(m[which_ncs], lai1, att[which_ncs], sigb[which_ncs],
|
|
747 ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
|
|
748 tdd[which_ncs] <- resncs$tdd
|
|
749 rdd[which_ncs] <- resncs$rdd
|
|
750 tsd[which_ncs] <- resncs$tsd
|
|
751 rsd[which_ncs] <- resncs$rsd
|
|
752 tdo[which_ncs] <- resncs$tdo
|
|
753 rdo[which_ncs] <- resncs$rdo
|
|
754 rsod[which_ncs] <- resncs$rsod
|
|
755 }
|
|
756 if (length(which_cs) > 0) {
|
|
757 rescs <- conservativescattering(m[which_cs], lai1, att[which_cs], sigb[which_cs],
|
|
758 ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
|
|
759 tdd[which_cs] <- rescs$tdd
|
|
760 rdd[which_cs] <- rescs$rdd
|
|
761 tsd[which_cs] <- rescs$tsd
|
|
762 rsd[which_cs] <- rescs$rsd
|
|
763 tdo[which_cs] <- rescs$tdo
|
|
764 rdo[which_cs] <- rescs$rdo
|
|
765 rsod[which_cs] <- rescs$rsod
|
|
766 }
|
|
767
|
|
768 # combine with bottom layer reflectances and transmittances (adding method)
|
|
769 rn <- 1.0 - rdd * rddb
|
|
770 tup <- (tss * rsdb + tsd * rddb) / rn
|
|
771 tdn <- (tsd + tss * rsdb * rdd) / rn
|
|
772 rsdt <- rsd + tup * tdd
|
|
773 rdot <- rdo + tdd * (rddb * tdo + rdob * too) / rn
|
|
774 rsodt <- rsod + (tss * rsodb + tdn * rdob) * too + tup * tdo
|
|
775
|
|
776 rsost <- (w1 * s1 + w2 * s2) * lai
|
|
777
|
|
778 rsot <- rsost + rsodt
|
|
779
|
|
780 # Diffuse reflectances at the top and the bottom are now different
|
|
781 rddt_t <- rdd + tdd * rddb * tdd / rn
|
|
782 rddt_b <- rddb + tddb * rdd * tddb / rn
|
|
783
|
|
784 # Transmittances of the combined canopy layers
|
|
785 tsst <- tss * tssb
|
|
786 toot <- too * toob
|
|
787 tsdt <- tss * tsdb + tdn * tddb
|
|
788 tdot <- tdob * too + tddb * (tdo + rdd * rdob * too) / rn
|
|
789 tddt <- tdd * tddb / rn
|
|
790
|
|
791 # Apply clumping effects to vegetation layer
|
|
792 rddcb <- cv * rddt_b
|
|
793 rddct <- cv * rddt_t
|
|
794 tddc <- 1 - cv + cv * tddt
|
|
795 rsdc <- cs * rsdt
|
|
796 tsdc <- cs * tsdt
|
|
797 rdoc <- co * rdot
|
|
798 tdoc <- co * tdot
|
|
799 tssc <- 1 - cs + cs * tsst
|
|
800 tooc <- 1 - co + co * toot
|
|
801
|
|
802 # new weight function fcdc for crown contribution (W. Verhoef, 22-05-08)
|
|
803 rsoc <- fcdc * rsot
|
|
804 tssooc <- fcd * tsstoo + fcs * toot + fod * tsst + fos
|
|
805 # Canopy absorptance for black background (W. Verhoef, 02-03-04)
|
|
806 alfas <- 1. - tssc - tsdc - rsdc
|
|
807 alfad <- 1. - tddc - rddct
|
|
808 # Add the soil background
|
|
809 rn <- 1 - rddcb * rddsoil
|
|
810 tup <- (tssc * rsdsoil + tsdc * rddsoil) / rn
|
|
811 tdn <- (tsdc + tssc * rsdsoil * rddcb) / rn
|
|
812
|
|
813 rddt <- rddct + tddc * rddsoil * tddc / rn
|
|
814 rsdt <- rsdc + tup * tddc
|
|
815 rdot <- rdoc + tddc * (rddsoil * tdoc + rdosoil * tooc) / rn
|
|
816 rsot <- rsoc + tssooc * rsosoil + tdn * rdosoil * tooc + tup * tdoc
|
|
817
|
|
818 # Effect of soil background on canopy absorptances (W. Verhoef, 02-03-04)
|
|
819 alfast <- alfas + tup * alfad
|
|
820 alfadt <- alfad * (1. + tddc * rddsoil / rn)
|
|
821 }
|
|
822 my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt,
|
|
823 "alfast" = alfast, "alfadt" = alfadt)
|
|
824 return(my_list)
|
|
825 }
|
|
826
|
|
827
|
|
828
|
|
829 #" computes non conservative scattering conditions
|
|
830 #" @param m numeric.
|
|
831 #" @param lai numeric. Leaf Area Index
|
|
832 #" @param att numeric.
|
|
833 #" @param sigb numeric.
|
|
834 #" @param ks numeric.
|
|
835 #" @param ko numeric.
|
|
836 #" @param sf numeric.
|
|
837 #" @param sb numeric.
|
|
838 #" @param vf numeric.
|
|
839 #" @param vb numeric.
|
|
840 #" @param tss numeric.
|
|
841 #" @param too numeric.
|
|
842 #"
|
|
843 #" @return list. tdd, rdd, tsd, rsd, tdo, rdo, rsod
|
|
844 #"
|
|
845 #" @export
|
|
846 nonconservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
|
|
847
|
|
848 e1 <- exp(-m * lai)
|
|
849 e2 <- e1 * e1
|
|
850 rinf <- (att - m) / sigb
|
|
851 rinf2 <- rinf * rinf
|
|
852 re <- rinf * e1
|
|
853 denom <- 1. - rinf2 * e2
|
|
854
|
|
855 j1ks <- jfunc1(ks, m, lai)
|
|
856 j2ks <- jfunc2(ks, m, lai)
|
|
857 j1ko <- jfunc1(ko, m, lai)
|
|
858 j2ko <- jfunc2(ko, m, lai)
|
|
859
|
|
860 ps <- (sf + sb * rinf) * j1ks
|
|
861 qs <- (sf * rinf + sb) * j2ks
|
|
862 pv <- (vf + vb * rinf) * j1ko
|
|
863 qv <- (vf * rinf + vb) * j2ko
|
|
864
|
|
865 tdd <- (1. - rinf2) * e1 / denom
|
|
866 rdd <- rinf * (1. - e2) / denom
|
|
867 tsd <- (ps - re * qs) / denom
|
|
868 rsd <- (qs - re * ps) / denom
|
|
869 tdo <- (pv - re * qv) / denom
|
|
870 rdo <- (qv - re * pv) / denom
|
|
871
|
|
872 z <- jfunc2(ks, ko, lai)
|
|
873 g1 <- (z - j1ks * too) / (ko + m)
|
|
874 g2 <- (z - j1ko * tss) / (ks + m)
|
|
875
|
|
876 tv1 <- (vf * rinf + vb) * g1
|
|
877 tv2 <- (vf + vb * rinf) * g2
|
|
878
|
|
879 t1 <- tv1 * (sf + sb * rinf)
|
|
880 t2 <- tv2 * (sf * rinf + sb)
|
|
881 t3 <- (rdo * qs + tdo * ps) * rinf
|
|
882
|
|
883 # Multiple scattering contribution to bidirectional canopy reflectance
|
|
884 rsod <- (t1 + t2 - t3) / (1. - rinf2)
|
|
885 my_list <- list("tdd" = tdd, "rdd" = rdd, "tsd" = tsd,
|
|
886 "rsd" = rsd, "tdo" = tdo, "rdo" = rdo,
|
|
887 "rsod" = rsod)
|
|
888 return(my_list)
|
|
889 }
|
|
890
|
|
891 #" computes conservative scattering conditions
|
|
892 #" @param m numeric.
|
|
893 #" @param lai numeric. Leaf Area Index
|
|
894 #" @param att numeric.
|
|
895 #" @param sigb numeric.
|
|
896 #" @param ks numeric.
|
|
897 #" @param ko numeric.
|
|
898 #" @param sf numeric.
|
|
899 #" @param sb numeric.
|
|
900 #" @param vf numeric.
|
|
901 #" @param vb numeric.
|
|
902 #" @param tss numeric.
|
|
903 #" @param too numeric.
|
|
904 #"
|
|
905 #" @return list. tdd, rdd, tsd, rsd, tdo, rdo, rsod
|
|
906 #"
|
|
907 #" @export
|
|
908 conservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
|
|
909
|
|
910 # near or complete conservative scattering
|
|
911 j4 <- jfunc4(m, lai)
|
|
912 amsig <- att - sigb
|
|
913 apsig <- att + sigb
|
|
914 rtp <- (1 - amsig * j4) / (1 + amsig * j4)
|
|
915 rtm <- (-1 + apsig * j4) / (1 + apsig * j4)
|
|
916 rdd <- 0.5 * (rtp + rtm)
|
|
917 tdd <- 0.5 * (rtp - rtm)
|
|
918
|
|
919 dns <- ks * ks - m * m
|
|
920 dno <- ko * ko - m * m
|
|
921 cks <- (sb * (ks - att) - sf * sigb) / dns
|
|
922 cko <- (vb * (ko - att) - vf * sigb) / dno
|
|
923 dks <- (-sf * (ks + att) - sb * sigb) / dns
|
|
924 dko <- (-vf * (ko + att) - vb * sigb) / dno
|
|
925 ho <- (sf * cko + sb * dko) / (ko + ks)
|
|
926
|
|
927 rsd <- cks * (1 - tss * tdd) - dks * rdd
|
|
928 rdo <- cko * (1 - too * tdd) - dko * rdd
|
|
929 tsd <- dks * (tss - tdd) - cks * tss * rdd
|
|
930 tdo <- dko * (too - tdd) - cko * too * rdd
|
|
931 # Multiple scattering contribution to bidirectional canopy reflectance
|
|
932 rsod <- ho * (1 - tss * too) - cko * tsd * too - dko * rsd
|
|
933
|
|
934 my_list <- list("tdd" = tdd, "rdd" = rdd, "tsd" = tsd,
|
|
935 "rsd" = rsd, "tdo" = tdo, "rdo" = rdo,
|
|
936 "rsod" = rsod)
|
|
937 return(my_list)
|
|
938 }
|
|
939
|
|
940
|
|
941
|
|
942
|
|
943
|
|
944
|
|
945 #" computes the leaf angle distribution function value (freq)
|
|
946 #"
|
|
947 #" Ellipsoidal distribution function characterised by the average leaf
|
|
948 #" inclination angle in degree (ala)
|
|
949 #" Campbell 1986
|
|
950 #" @param ala average leaf angle
|
|
951 #" @return foliar_distrib list. lidf and litab
|
|
952 #" @export
|
|
953 campbell <- function(ala) {
|
|
954
|
|
955 tx1 <- c(10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88., 90.)
|
|
956 tx2 <- c(0., 10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88.)
|
|
957 litab <- (tx2 + tx1) / 2
|
|
958 n <- length(litab)
|
|
959 tl1 <- tx1 * (pi / 180)
|
|
960 tl2 <- tx2 * (pi / 180)
|
|
961 excent <- exp(-1.6184e-5 * ala**3 + 2.1145e-3 * ala**2 - 1.2390e-1 * ala + 3.2491)
|
|
962 sum0 <- 0
|
|
963
|
|
964 freq <- c()
|
|
965 for (i in 1:n) {
|
|
966 x1 <- excent / (sqrt(1. + excent**2. * tan(tl1[i])**2))
|
|
967 x2 <- excent / (sqrt(1. + excent**2. * tan(tl2[i])**2))
|
|
968 if (excent == 1) {
|
|
969 freq[i] <- abs(cos(tl1[i]) - cos(tl2[i]))
|
|
970 } else {
|
|
971 alpha <- excent / sqrt(abs(1 - excent**2))
|
|
972 alpha2 <- alpha**2
|
|
973 x12 <- x1**2
|
|
974 x22 <- x2**2
|
|
975 alpx1 <- 0 * alpha2
|
|
976 alpx2 <- 0 * alpha2
|
|
977 almx1 <- 0 * alpha2
|
|
978 almx2 <- 0 * alpha2
|
|
979 if (excent > 1) {
|
|
980 alpx1 <- sqrt(alpha2[excent > 1] + x12[excent > 1])
|
|
981 alpx2[excent > 1] <- sqrt(alpha2[excent > 1] + x22[excent > 1])
|
|
982 dum <- x1 * alpx1 + alpha2 * log(x1 + alpx1)
|
|
983 freq[i] <- abs(dum - (x2 * alpx2 + alpha2 * log(x2 + alpx2)))
|
|
984 } else {
|
|
985 almx1 <- sqrt(alpha2 - x12)
|
|
986 almx2 <- sqrt(alpha2 - x22)
|
|
987 dum <- x1 * almx1 + alpha2 * asin(x1 / alpha)
|
|
988 freq[i] <- abs(dum - (x2 * almx2 + alpha2 * asin(x2 / alpha)))
|
|
989 }
|
|
990 }
|
|
991 }
|
|
992 sum0 <- sum(freq)
|
|
993 freq0 <- freq / sum0
|
|
994 foliar_distrib <- list("lidf" = freq0, "litab" = litab)
|
|
995 return(foliar_distrib)
|
|
996 }
|
|
997
|
|
998 #" computes the leaf angle distribution function value (freq)
|
|
999 #"
|
|
1000 #" Using the original bimodal distribution function initially proposed in SAIL
|
|
1001 #" references
|
|
1002 #" ----------
|
|
1003 #" (Verhoef1998) Verhoef, Wout. Theory of radiative transfer models applied
|
|
1004 #" in optical remote sensing of vegetation canopies.
|
|
1005 #" nationaal Lucht en Ruimtevaartlaboratorium, 1998.
|
|
1006 #" http: / / library.wur.nl / WebQuery / clc / 945481.
|
|
1007 #" @param a controls the average leaf slope
|
|
1008 #" @param b controls the distribution"s bimodality
|
|
1009 #" LIDF type a b
|
|
1010 #" Planophile 1 0
|
|
1011 #" Erectophile -1 0
|
|
1012 #" Plagiophile 0 -1
|
|
1013 #" Extremophile 0 1
|
|
1014 #" Spherical -0.35 -0.15
|
|
1015 #" Uniform 0 0
|
|
1016 #" requirement: ||lidfa|| + ||lidfb|| < 1
|
|
1017 #"
|
|
1018 #" @return foliar_distrib list. lidf and litab
|
|
1019 #" @export
|
|
1020 dladgen <- function(a, b) {
|
|
1021 litab <- c(5., 15., 25., 35., 45., 55., 65., 75., 81., 83., 85., 87., 89.)
|
|
1022 freq <- c()
|
|
1023 for (i1 in 1:8) {
|
|
1024 t <- i1 * 10
|
|
1025 freq[i1] <- dcum(a, b, t)
|
|
1026 }
|
|
1027 for (i2 in 9:12) {
|
|
1028 t <- 80. + (i2 - 8) * 2.
|
|
1029 freq[i2] <- dcum(a, b, t)
|
|
1030 }
|
|
1031 freq[13] <- 1
|
|
1032 for (i in 13:2) {
|
|
1033 freq[i] <- freq[i] - freq[i - 1]
|
|
1034 }
|
|
1035 foliar_distrib <- list("lidf" = freq, "litab" = litab)
|
|
1036 return(foliar_distrib)
|
|
1037 }
|
|
1038
|
|
1039 #" dcum function
|
|
1040 #" @param a numeric. controls the average leaf slope
|
|
1041 #" @param b numeric. controls the distribution"s bimodality
|
|
1042 #" @param t numeric. angle
|
|
1043 #" @return f
|
|
1044 #" @export
|
|
1045 dcum <- function(a, b, t) {
|
|
1046 rd <- pi / 180
|
|
1047 if (a >= 1) {
|
|
1048 f <- 1 - cos(rd * t)
|
|
1049 } else {
|
|
1050 eps <- 1e-8
|
|
1051 delx <- 1
|
|
1052 x <- 2 * rd * t
|
|
1053 p <- x
|
|
1054 while (delx >= eps) {
|
|
1055 y <- a * sin(x) + .5 * b * sin(2. * x)
|
|
1056 dx <- .5 * (y - x + p)
|
|
1057 x <- x + dx
|
|
1058 delx <- abs(dx)
|
|
1059 }
|
|
1060 f <- (2. * y + p) / pi
|
|
1061 }
|
|
1062 return(f)
|
|
1063 }
|
|
1064
|
|
1065 #" J1 function with avoidance of singularity problem
|
|
1066 #"
|
|
1067 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
|
|
1068 #" @param l numeric.
|
|
1069 #" @param t numeric. Leaf Area Index
|
|
1070 #" @return jout numeric.
|
|
1071 #" @export
|
|
1072 jfunc1 <- function(k, l, t) {
|
|
1073 # J1 function with avoidance of singularity problem
|
|
1074 del <- (k - l) * t
|
|
1075 jout <- 0 * l
|
|
1076 jout[which(abs(del) > 1e-3)] <- (exp(-l[which(abs(del) > 1e-3)] * t) - exp(-k * t)) / (k - l[which(abs(del) > 1e-3)])
|
|
1077 jout[which(abs(del) <= 1e-3)] <- 0.5 * t * (exp(-k * t) + exp(-l[which(abs(del) <= 1e-3)] * t)) * (1 - del[which(abs(del) <= 1e-3)] * del[which(abs(del) <= 1e-3)] / 12)
|
|
1078 return(jout)
|
|
1079 }
|
|
1080
|
|
1081 #" J2 function with avoidance of singularity problem
|
|
1082 #"
|
|
1083 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
|
|
1084 #" @param l numeric.
|
|
1085 #" @param t numeric. Leaf Area Index
|
|
1086 #" @return jout numeric.
|
|
1087 #" @export
|
|
1088 jfunc2 <- function(k, l, t) {
|
|
1089 # J2 function
|
|
1090 jout <- (1. - exp(-(k + l) * t)) / (k + l)
|
|
1091 return(jout)
|
|
1092 }
|
|
1093
|
|
1094 #" J3 function with avoidance of singularity problem
|
|
1095 #"
|
|
1096 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
|
|
1097 #" @param l numeric.
|
|
1098 #" @param t numeric. Leaf Area Index
|
|
1099 #" @return jout numeric.
|
|
1100 #" @export
|
|
1101 jfunc3 <- function(k, l, t) {
|
|
1102 out <- (1. - exp(-(k + l) * t)) / (k + l)
|
|
1103 return(out)
|
|
1104 }
|
|
1105
|
|
1106
|
|
1107 #" j4 function for treating (near) conservative scattering
|
|
1108 #"
|
|
1109 #" @param m numeric. Extinction coefficient for direct (solar or observer) flux
|
|
1110 #" @param t numeric. Leaf Area Index
|
|
1111 #" @return jout numeric.
|
|
1112 #" @export
|
|
1113 jfunc4 <- function(m, t) {
|
|
1114
|
|
1115 del <- m * t
|
|
1116 out <- 0 * del
|
|
1117 out[del > 1e-3] <- (1 - exp(-del)) / (m * (1 + exp(-del)))
|
|
1118 out[del <= 1e-3] <- 0.5 * t * (1. - del * del / 12.)
|
|
1119 return(out)
|
|
1120 }
|
|
1121
|
|
1122
|
|
1123 #" compute volume scattering functions and interception coefficients
|
|
1124 #" for given solar zenith, viewing zenith, azimuth and leaf inclination angle.
|
|
1125 #"
|
|
1126 #" @param tts numeric. solar zenith
|
|
1127 #" @param tto numeric. viewing zenith
|
|
1128 #" @param psi numeric. azimuth
|
|
1129 #" @param ttl numeric. leaf inclination angle
|
|
1130 #" @return res list. includes chi_s, chi_o, frho, ftau
|
|
1131 #" @export
|
|
1132 volscatt <- function(tts, tto, psi, ttl) {
|
|
1133 #********************************************************************************
|
|
1134 #* chi_s = interception functions
|
|
1135 #* chi_o = interception functions
|
|
1136 #* frho = function to be multiplied by leaf reflectance rho
|
|
1137 #* ftau = functions to be multiplied by leaf transmittance tau
|
|
1138 #********************************************************************************
|
|
1139 # Wout Verhoef, april 2001, for CROMA
|
|
1140
|
|
1141 rd <- pi / 180
|
|
1142 costs <- cos(rd * tts)
|
|
1143 costo <- cos(rd * tto)
|
|
1144 sints <- sin(rd * tts)
|
|
1145 sinto <- sin(rd * tto)
|
|
1146 cospsi <- cos(rd * psi)
|
|
1147 psir <- rd * psi
|
|
1148 costl <- cos(rd * ttl)
|
|
1149 sintl <- sin(rd * ttl)
|
|
1150 cs <- costl * costs
|
|
1151 co <- costl * costo
|
|
1152 ss <- sintl * sints
|
|
1153 so <- sintl * sinto
|
|
1154
|
|
1155 #c ..............................................................................
|
|
1156 #c betas -bts- and betao -bto- computation
|
|
1157 #c Transition angles (beta) for solar (betas) and view (betao) directions
|
|
1158 #c if thetav + thetal > pi / 2, bottom side of the leaves is observed for leaf azimut
|
|
1159 #c interval betao + phi<leaf azimut<2pi-betao + phi.
|
|
1160 #c if thetav + thetal<pi / 2, top side of the leaves is always observed, betao=pi
|
|
1161 #c same consideration for solar direction to compute betas
|
|
1162 #c ..............................................................................
|
|
1163
|
|
1164 cosbts <- 5
|
|
1165 if (abs(ss) > 1e-6) {
|
|
1166 cosbts <- -cs / ss
|
|
1167 }
|
|
1168 cosbto <- 5
|
|
1169 if (abs(so) > 1e-6) {
|
|
1170 cosbto <- -co / so
|
|
1171 }
|
|
1172
|
|
1173 if (abs(cosbts) < 1) {
|
|
1174 bts <- acos(cosbts)
|
|
1175 ds <- ss
|
|
1176 } else {
|
|
1177 bts <- pi
|
|
1178 ds <- cs
|
|
1179 }
|
|
1180 chi_s <- 2. / pi * ((bts - pi * .5) * cs + sin(bts) * ss)
|
|
1181 if (abs(cosbto) < 1) {
|
|
1182 bto <- acos(cosbto)
|
|
1183 doo <- so
|
|
1184 } else if (tto < 90) {
|
|
1185 bto <- pi
|
|
1186 doo <- co
|
|
1187 } else {
|
|
1188 bto <- 0
|
|
1189 doo <- -co
|
|
1190 }
|
|
1191 chi_o <- 2. / pi * ((bto - pi * .5) * co + sin(bto) * so)
|
|
1192
|
|
1193 #c ..............................................................................
|
|
1194 #c computation of auxiliary azimut angles bt1, bt2, bt3 used
|
|
1195 #c for the computation of the bidirectional scattering coefficient w
|
|
1196 #c .............................................................................
|
|
1197
|
|
1198 btran1 <- abs(bts - bto)
|
|
1199 btran2 <- pi - abs(bts + bto - pi)
|
|
1200
|
|
1201 if (psir <= btran1) {
|
|
1202 bt1 <- psir
|
|
1203 bt2 <- btran1
|
|
1204 bt3 <- btran2
|
|
1205 } else {
|
|
1206 bt1 <- btran1
|
|
1207 if (psir <= btran2) {
|
|
1208 bt2 <- psir
|
|
1209 bt3 <- btran2
|
|
1210 } else {
|
|
1211 bt2 <- btran2
|
|
1212 bt3 <- psir
|
|
1213 }
|
|
1214 }
|
|
1215 t1 <- 2. * cs * co + ss * so * cospsi
|
|
1216 t2 <- 0
|
|
1217 if (bt2 > 0) {
|
|
1218 t2 <- sin(bt2) * (2. * ds * doo + ss * so * cos(bt1) * cos(bt3))
|
|
1219 }
|
|
1220
|
|
1221 denom <- 2. * pi * pi
|
|
1222 frho <- ((pi - bt2) * t1 + t2) / denom
|
|
1223 ftau <- (-bt2 * t1 + t2) / denom
|
|
1224
|
|
1225 if (frho < 0) {
|
|
1226 frho <- 0
|
|
1227 }
|
|
1228 if (ftau < 0) {
|
|
1229 ftau <- 0
|
|
1230 }
|
|
1231 res <- list("chi_s" = chi_s, "chi_o" = chi_o, "frho" = frho, "ftau" = ftau)
|
|
1232 return(res)
|
|
1233 }
|