Nothing
#' R implementation of the foursail2 model with 2 canopy layers.
#'
#' The foursail2 model is a two layer implementation of the
#' foursail model described in Verhoef and Bach (2007).
#' Layers are assumed identical in particle inclination and
#' hotspot, but may differ in the relative density and types of
#' particles (see foursail2b for a layer specific inclination angle).
#' In comparison to foursail, the background (soil),
#' can now be non-Lambertain, having it own 4-stream
#' BDRF (not implemented here but may be input by the user).
#' There are two types of particles, generalized
#' to primary and secondary (originally termed "green"
#' and "brown" particles). The realtive abundance of
#' the secondary particle in the top canopy is regulated by
#' the dissociation paramerter.The model 4SAIL2 combines with
#' prospect, libery or procosine for the reflectance
#' and transmittance of the particles, and with the the foursail
#' or Hapke elements for the background reflectance.
#' If run alone, these require direct inputs which could be
#' measured leaf reflectance.
#'
#' @param rhoA primary particle reflectance from 400-2500nm (can be measured or modeled)
#' @param tauA primary particle transmittance from 400-2500nm (can be measured or modeled)
#' @param rhoB secondary particle reflectance from 400-2500nm (can be measured or modeled)
#' @param tauB secondary particle reflectance from 400-2500nm (can be measured or modeled)
#' @param rsobgr : background bidirectional reflectance (rso)
#' @param rdobgr : background directional hemispherical reflectance (rdo)
#' @param rsdbgr : background hemispherical directional reflectance (rsd)
#' @param rddbgr : background bi-hemispherical diffuse reflectance (rdd)
#' @param bgr background reflectance. Usual input is
#' soil reflectance spectra from 400-2500nm (can be measured or modeled)
#' @param param A named vector of 4SAIL2 parameter values (note:
#' program ignores case):
#' \itemize{
#' \item [1] = Leaf angle distribution function parameter a (LIDFa)
#' \item [2] = Leaf angle distribution function parameter b (LIDFb)
#' \item [3] = Leaf angle distribution function type (TypeLidf, see ?lidfFun)
#' \item [4] = Total Leaf Area Index (LAI), including primary and secondary
#' particles (brown and green leafs).
#' \item [5] = fraction secondary particles ("brown leaf fraction", fb)
#' \item [6] = Canopy dissociation factor for secondary particles ("diss")
#' \item [7] = Hot spot effect parameter (hspot). Often defined as the
#' ratio of mean leaf width and canopy height.
#' \item [7] = vertical crown coverage fraction (Cv), models clumping in combination
#' with parameter zeta.
#' \item [7] = tree shape factor (zeta), defined as the ratio of crown diameter and height.
#' \item [6] = Solar zenith angle (tts)
#' \item [7] = Observer zenith angle (tto)
#' \item [8] = Sun-sensor azimuth angle (psi)
#' }
#'
#' @return spectra matrixwith 4 reflectance factors and canopy transmission
#' for wavelengths 400 to 2500nm:
#' \itemize{
#' \item [1] = bi-hemispherical reflectance (rddt). White-sky albedo: the reflectance of the canopy
#' under diffuse illumination. The BRDF integrated over all viewing and illumination directions.
#' Diffuse reflectance for diffuse incidence.
#' \item [2] = hemispherical directional reflectance (rsdt). Black-sky albedpo: reflectance of a surface
#' under direct light without a diffuse component. It is the integral of the BRDF over all viewing
#' directions. Diffuse reflectance for direct solar incidence.
#' \item [3] = directional hemispherical reflectance (rdot). Diffuse reflectance in the vieweing
#' direction.
#' \item [4] = bi-directional reflectance (rsot). The ratio of reflected radiance in the viewing direction
#' to the incoming radiant flux in the solar direction.
#' }
#'
#' @examples
#' ## see ?foursail for lower-level implementations
#' fRTM(rho~prospect5+foursail2)
#'
#' @references Verhoef, W., Bach, H. (2007). Coupled soil-leaf-canopy and atmosphere radiative transfer
#' modeling to simulate hyperspectral multi-angular surface reflectance and TOA radiance data.
#' Remote Sens. Environ. 109, 166-182.
#' @export
foursail2 <- function(rhoA,tauA,rhoB=NULL,tauB=NULL,bgr,
rsobgr=NULL,rdobgr=NULL,rsdbgr=NULL,rddbgr=NULL,
param){
if(!is.null(rsobgr)) stop("SOIL BDRF not implemented currently")
if(is.null(rhoB)) { ## standard if B is empty
rhoB <- rhoA
tauB <- tauA
}
names(param) <- tolower(names(param))
## input paramters
LIDFa <- as.numeric(param["lidfa"])
LIDFb <- as.numeric(param["lidfb"])
TypeLIDF <- as.numeric(param["typelidf"])
lai <- as.numeric(param["lai"])
q <- as.numeric(param["hspot"])
Cv <- as.numeric(param["cv"])
Zeta0 <- as.numeric(param["zeta"])
fb <- as.numeric(param["fb"])
diss <- as.numeric(param["diss"])
psoil <- as.numeric(param["psoil"]) ## change to background
#skyl <- as.numeric(param["skyl"])
tts <- as.numeric(param["tts"])
tto <- as.numeric(param["tto"])
psi <- as.numeric(param["psi"])
## Leaf angle distribution function
## TO DO: build generic N angle LIA models
na <- 18 # number of angle
rd <- pi/180
lna <- seq(0,17,length.out=18)
thetal <- 2.5+5*lna
## reventing back to 4sail standard
## leaf angles
lna <- c(5,15,25,35,45,55,65,75,81,83,85,87,89)
## number of angles
na <- length(lna)
lidfpar <- c(na,LIDFa, LIDFb)
names(lidfpar) <- c("na","a","b")
class(lidfpar) <- paste0("lidf.",TypeLIDF)
lidFun <- lidf(lidfpar)
## special case when lai = 0
if(lai<=0){
## gap fractions (direct transmittance)
tss <- 1
too <- 1
tsstoo <- 1 # (bi-directional direct transmittance)
## reflectance and transmission
## hemispherical diffuse
rdd <- 0
tdd <- 1
## directional hemispheric for direct solar (s) and in viewing angle/observer (o)
rsd <- 0
tsd <- 0
rdo <- 0
tdo <- 0
## bi-directional in viewing angle
rso <- 0
rsos <- 0
rsod <- 0
rddt <- bgr
rsdt <- bgr
rdot <- bgr
rsodt <- 0
rsost <- bgr
rsot <- bgr
} else {
## Sensor geometry: Compute geometric quantities
## output = list(cts,cto,ctscto,dso)
sensGeom <- sail_sensgeom(tts,tto,psi,rd)
cts <- sensGeom[[1]]
cto <- sensGeom[[2]]
ctscto <- sensGeom[[3]]
dso <- sensGeom[[4]]
## Crown and vegatation clumping effects
## similar to flim implementation
Cs <- 1 # initialize
Co <- 1 # initialize
if(Cv<=1){
Cs <- 1-(1-Cv)^(1/cts)
Co <- 1-(1-Cv)^(1/cto)
}
Overlap <- 0
if(Zeta0>0){
Overlap <- min(Cs*(1-Co),Co*(1-Cs))*exp(-dso/Zeta0)
}
Fcd <- Cs*Co+Overlap
Fcs <- (1-Cs)*Co-Overlap
Fod <- Cs*(1-Co)-Overlap
Fos <- (1-Cs)*(1-Co)+Overlap
Fcdc <- 1-(1-Fcd)^(.5/cts+.5/cto)
## Part depending on diss, fb, and leaf optical properties
## First save the input fb as the old fb, as the following change is only artificial
## Better define an fb that is actually used: fbu, so that the input is not modified!
fbu <- fb
if(fb==0){ # if only green leaves
fbu <- 0.5
rhoB <- rhoA
tauB <- tauA
}
if(fb==1){ # only brown
fbu <- 0.5
rhoA <- rhoB
tauA <- tauB
}
s <- (1-diss)*fbu*(1-fbu)
## mixing of primary and secondary particles
## rho1 & tau1 : green foliage
## rho2 & tau2 : brown foliage (bottom layer)
rho1 <- ((1-fbu-s)*rhoA+s*rhoB)/(1-fbu)
tau1 <- ((1-fbu-s)*tauA+s*tauB)/(1-fbu)
rho2 <- (s*rhoA+(fbu-s)*rhoB)/fbu
tau2 <- (s*tauA+(fbu-s)*tauB)/fbu
## Angular distance - shadow length compensation
## output = list(ks,ko,sob,sof,sdb,sdf,dob,dof,ddb,ddf)
suitRes <- SUITS(na,lna,lidFun,tts,tto,cts,cto,psi,ctscto)
ks <- suitRes[[1]]
ko <- suitRes[[2]]
sob <- suitRes[[3]]
sof <- suitRes[[4]]
sdb <- suitRes[[5]]
sdf <- suitRes[[6]]
dob <- suitRes[[7]]
dof <- suitRes[[8]]
ddb <- suitRes[[9]]
ddf <- suitRes[[10]]
## devide the leaf area index in two layers
lai1 <- (1-fbu)*lai
lai2 <- fbu*lai
## 2 layer hot spot effect
hspotRes <- HotSpot2(lai,lai1,fbu,q,tss,ks,ko,dso)
tsstoo <- hspotRes[[1]]
sumint1 <- hspotRes[[2]]
sumint2 <- hspotRes[[3]]
## Calculate reflectances and transmittances for 2 layers
## Input LAI and leaf rho and tau
## and geometric factors
refTransRes <- ReflTrans2(rho1,rho2,tau1,tau2,ks,ko,
lai1,lai2,sdf,sdb,dof,dob,
sob,sof,ddb,ddf,
sumint1,sumint2)
dn <- refTransRes[[1]]
tup <- refTransRes[[2]]
tdn <- refTransRes[[3]]
rsdt <- refTransRes[[4]]
rdot <- refTransRes[[5]]
rsodt <- refTransRes[[6]]
rsost <- refTransRes[[7]]
rsot <- refTransRes[[8]]
rddt_t <- refTransRes[[9]]
rddt_b <- refTransRes[[10]]
tsst <- refTransRes[[11]]
toot <- refTransRes[[12]]
tsdt <- refTransRes[[13]]
tdot <- refTransRes[[14]]
tddt <- refTransRes[[15]]
## Apply clumping effects to vegetation layer
rddcb <- Cv*rddt_b
rddct <- Cv*rddt_t
tddc <- 1-Cv+Cv*tddt
rsdc <- Cs*rsdt
tsdc <- Cs*tsdt
rdoc <- Co*rdot
tdoc <- Co*tdot
tssc <- 1-Cs+Cs*tsst
tooc <- 1-Co+Co*toot
## New weight function Fcdc for crown contribution
## (W. Verhoef, 22-05-08)
rsoc <- Fcdc*rsot
tssooc <- Fcd*tsstoo+Fcs*toot+Fod*tsst+Fos
## Add the soil background
dn <- 1-rddcb*bgr
tup <- (tssc*bgr+tsdc*bgr)/dn
tdn <- (tsdc+tssc*bgr*rddcb)/dn
rddt <- rddct+tddc*bgr*tddc/dn
rsdt <- rsdc+tup*tddc
rdot <- rdoc+tddc*(bgr*tdoc+bgr*tooc)/dn
rsot <- rsoc+tssooc*bgr+tdn*bgr*tooc+tup*tdoc
}
taus <- (tsdc+tssc) ## direct and diffuse transmission of light
## rddt Bi-hemispherical reflectance
## rsdt Directional-hemispherical reflectance for solar incident flux
## rdot Hemispherical-directional reflectance in viewing direction
## rsot Bi-directional reflectance factor
finalOut <- list(rddt,rsdt,rdot,rsot,taus)
reflMat <- as.matrix(do.call(cbind,finalOut))
colnames(reflMat) <- c('rddt','rsdt','rdot','rsot','tau')
return(reflMat)
}
## 2 layer Hotspot function
HotSpot2 <- function(lai,lai1,fbu,q,tss,ks,ko,dso){
## Treatment of the hotspot-effect for 2 layers
tss <- exp(-ks*lai) ## full canopy
ck <- exp(-ks*lai1) ## only top
alf <- 1e6
## Apply correction 2/(K+k) suggested by F.M. Breon
if(q>0) alf <- (dso/q)*2/(ks+ko)
if(alf>200) alf <- 200 ## inserted H. Bach 1/3/04
if(alf==0){
## The pure hotspot - no shadow
tsstoo <- tss
sumint1 <- (1-ck)/(ks*lai)
sumint2 <- (ck-tss)/(ks*lai)
} else {
## Outside the hotspot
fhot <- lai*sqrt(ko*ks)
## Integrate by exponential Simpson method in 20 steps
## the steps are arranged according to equal partitioning
## of the slope of the joint probability function
x1 <- 0
y1 <- 0
f1 <- 1
ca <- exp(alf*(fbu-1))
fint <- (1-ca)*.05
sumint1 <- 0
for(i in 1:20){
if(i<20) {
x2 <- -log(1-i*fint)/alf
} else {
x2 <- 1-fbu
}
y2 <- -(ko+ks)*lai*x2+fhot*(1-exp(-alf*x2))/alf
f2 <- exp(y2)
sumint1 <- sumint1+(f2-f1)*(x2-x1)/(y2-y1)
x1 <- x2
y1 <- y2
f1 <- f2
}
fint <- (ca-exp(-alf))*.05
sumint2 <- 0
for(i in 1:20){
if(i<20) {
x2 <- -log(ca-i*fint)/alf
} else {
x2 <- 1
}
y2 <- -(ko+ks)*lai*x2+fhot*(1-exp(-alf*x2))/alf
f2 <- exp(y2)
sumint2 <- sumint2+(f2-f1)*(x2-x1)/(y2-y1)
x1 <- x2
y1 <- y2
f1 <- f2
}
tsstoo <- f1
}
return(list(tsstoo,sumint1,sumint2))
} # HotSpot2
ReflTrans2 <- function(rho1,rho2,tau1,tau2,ks,ko,
lai1,lai2,sdf,sdb,dof,dob,
sob,sof,ddb,ddf,
sumint1,sumint2){
## Bottom layer (rho2,tau2,lai2)
bottom <- cReflTransSingleLayer(rho2,tau2,lai2,
ks,ko,
sdf,sdb,dof,dob,
sob,sof,ddb,ddf)
rddb <- bottom[[1]]
rsdb <- bottom[[2]]
rdob <- bottom[[3]]
rsodb <- bottom[[4]]
tddb <- bottom[[5]]
tsdb <- bottom[[6]]
tdob <- bottom[[7]]
toob <- bottom[[8]]
tssb <- bottom[[9]]
w2 <- bottom[[10]]
## Top layer (rho1,tau1,lai1)
top <- cReflTransSingleLayer(rho1,tau1,lai1,
ks,ko,
sdf,sdb,dof,dob,
sob,sof,ddb,ddf)
rdd <- top[[1]]
rsd <- top[[2]]
rdo <- top[[3]]
rsod <- top[[4]]
tdd <- top[[5]]
tsd <- top[[6]]
tdo <- top[[7]]
too <- top[[8]]
tss <- top[[9]]
w1 <- top[[10]]
## Combine with bottom layer reflectances and
## transmittances (adding method)
lai <- lai1+lai2 ## total LAI
dn <- 1-rdd*rddb
tup<-(tss*rsdb+tsd*rddb)/dn
tdn<-(tsd+tss*rsdb*rdd)/dn
rsdt<-rsd+tup*tdd
rdot<-rdo+tdd*(rddb*tdo+rdob*too)/dn
rsodt<-rsod+(tss*rsodb+tdn*rdob)*too+tup*tdo
rsost<-(w1*sumint1+w2*sumint2)*lai
rsot<-rsost+rsodt
## Diffuse reflectances at the top and the bottom are now different
rddt_t<-rdd+tdd*rddb*tdd/dn
rddt_b<-rddb+tddb*rdd*tddb/dn
## Transmittances of the combined canopy layers
tsst<-tss*tssb
toot<-too*toob
tsdt<-tss*tsdb+tdn*tddb
tdot<-tdob*too+tddb*(tdo+rdd*rdob*too)/dn
tddt<-tdd*tddb/dn
return(list(dn,tup,tdn,rsdt,rdot,rsodt,rsost,rsot,rddt_t,
rddt_b,tsst,toot,tsdt,tdot,tddt))
}
## calculate transmittance and reflectance of a single layer
ReflTransSingleLayer <- function(rho,tau,lai,
ks,ko,
sdf,sdb,dof,dob,
sob,sof,ddb,ddf){
## if 3 or leaf refl-- Geometric leaf reflectance
## Input rho and tau from PROSPECT/ other particle model
## Correct using SUITS factors
## Here the reflectance and transmission come in
## with the suits coefficients
RTgeomRes <- RTgeom(rho,tau,ddb,ddf,sdb,sdf,dob,dof,sob,sof)
sigb <- RTgeomRes[[1]]
att <- RTgeomRes[[2]]
m <- RTgeomRes[[3]]
sb <- RTgeomRes[[4]]
sf <- RTgeomRes[[5]]
vb <- RTgeomRes[[6]]
vf <- RTgeomRes[[7]]
w <- RTgeomRes[[8]]
tss <- exp(-ks*lai);
too <- exp(-ko*lai);
J1ks <- Jfunc1(ks,m,lai)
J2ks <- Jfunc2(ks,m,lai)
J1ko <- Jfunc1(ko,m,lai)
J2ko <- Jfunc2(ko,m,lai)
z <- Jfunc3(ks,ko,lai)
J4 <- Jfunc4(m,lai)
m_val1 <- which(m>0.01)
m_val2 <- which(m<0.01)
## Normal case
rinf<-rinf2<-e1<-e2<-denom<-re<-rdd<-tdd <- 0*m
e1 <- exp(-m[m_val1]*lai)
e2[m_val1] <- e1[m_val1]*e1[m_val1]
rinf[m_val1] <- ((att[m_val1]-m[m_val1])/sigb[m_val1])
rinf2[m_val1] <- rinf[m_val1]*rinf[m_val1]
re[m_val1] <- rinf[m_val1]*e1[m_val1]
denom[m_val1] <- 1-rinf2[m_val1]*e2[m_val1]
Ps<-Qs<-Pv<-Qv<-rdo<-tdo<-rds<-tds<-tsd<-rsd<- 0*m
Ps[m_val1] <- (sf[m_val1]+sb[m_val1]*rinf[m_val1])*J1ks[m_val1]
Qs[m_val1] <- (sf[m_val1]*rinf[m_val1]+sb[m_val1])*J2ks[m_val1]
Pv[m_val1] <- (vf[m_val1]+vb[m_val1]*rinf[m_val1])*J1ko[m_val1]
Qv[m_val1] <- (vf[m_val1]*rinf[m_val1]+vb[m_val1])*J2ko[m_val1]
rdd[m_val1] <- rinf[m_val1]*(1-e2[m_val1])/denom[m_val1]
tdd[m_val1] <- (1-rinf2[m_val1])*e1[m_val1]/denom[m_val1]
tsd[m_val1] <- (Ps[m_val1]-re[m_val1]*Qs[m_val1])/denom[m_val1]
rsd[m_val1] <- (Qs[m_val1]-re[m_val1]*Ps[m_val1])/denom[m_val1]
tdo[m_val1] <- (Pv[m_val1]-re[m_val1]*Qv[m_val1])/denom[m_val1]
rdo[m_val1] <- (Qv[m_val1]-re[m_val1]*Pv[m_val1])/denom[m_val1]
g1 <- g2 <- 0*m
g1[m_val1] <- (z-J1ks[m_val1]*too)/(ko+m[m_val1])
g2[m_val1] <- (z-J1ko[m_val1]*tss)/(ks+m[m_val1])
Tv1 <- Tv2 <- T1 <- T2 <- T3 <- 0*m
Tv1[m_val1] <- (vf[m_val1]*rinf[m_val1]+vb[m_val1])*g1[m_val1]
Tv2[m_val1] <- (vf[m_val1]+vb[m_val1]*rinf[m_val1])*g2[m_val1]
T1[m_val1] <- Tv1[m_val1]*(sf[m_val1]+sb[m_val1]*rinf[m_val1])
T2[m_val1] <- Tv2[m_val1]*(sf[m_val1]*rinf[m_val1]+sb[m_val1])
T3[m_val1] <- (rdo[m_val1]*Qs[m_val1]+
tdo[m_val1]*Ps[m_val1])*rinf[m_val1]
## Multiple scattering contribution to
## bidirectional canopy reflectance
rsod <- 0*m
rsod[m_val1] <- (T1[m_val1]+T2[m_val1]-T3[m_val1])/(1-rinf2[m_val1])
## Near or complete conservative scattering
amsig <- apsig <- rtp <- rtm <- 0*m
amsig[m_val2] <- att[m_val2]-sigb[m_val2]
apsig[m_val2] <- att[m_val2]+sigb[m_val2]
rtp[m_val2] <- (1-amsig[m_val2]*J4[m_val2])/
(1+amsig[m_val2]*J4[m_val2])
rtm[m_val2]<-(-1+apsig[m_val2]*J4[m_val2])/
(1+apsig[m_val2]*J4[m_val2])
rdd[m_val2]<-5*(rtp[m_val2]+rtm[m_val2])
tdd[m_val2]<-5*(rtp[m_val2]-rtm[m_val2])
dns<-dno<-cks<-cko<-dks<-dko<-ho<-0*m
dns[m_val2]<-ks*ks-m[m_val2]*m[m_val2]
dno[m_val2]<-ko*ko-m[m_val2]*m[m_val2]
cks[m_val2]<-(sb[m_val2]*(ks-att[m_val2])-
sf[m_val2]*sigb[m_val2])/dns[m_val2]
cko[m_val2]<-(vb[m_val2]*(ko-att[m_val2])-
vf[m_val2]*sigb[m_val2])/dno[m_val2]
dks[m_val2]<-(-sf[m_val2]*(ks+att[m_val2])-
sb[m_val2]*sigb[m_val2])/dns[m_val2]
dko[m_val2]<-(-vf[m_val2]*(ko+att[m_val2])-
vb[m_val2]*sigb[m_val2])/dno[m_val2]
ho[m_val2]<-(sf[m_val2]*cko[m_val2]+sb[m_val2]*dko[m_val2])/(ko+ks)
rsd[m_val2]<-cks[m_val2]*(1-tss*tdd[m_val2])-dks[m_val2]*rdd[m_val2]
rdo[m_val2]<-cko[m_val2]*(1-too*tdd[m_val2])-dko[m_val2]*rdd[m_val2]
tsd[m_val2]<-dks[m_val2]*(tss-tdd[m_val2])-
cks[m_val2]*tss *rdd[m_val2]
tdo[m_val2]<-dko[m_val2]*(too -tdd[m_val2])-
cko[m_val2]*too*rdd[m_val2]
## Multiple scattering contribution to
## bidirectional canopy reflectance
rsod[m_val2]<-ho[m_val2]*(1-tss*too)-
cko[m_val2]*tsd[m_val2]*too-dko[m_val2]*rsd[m_val2]
return(list(rdd,rsd,rdo,rsod,tdd,tsd,tdo,too,tss,w))
} # ReflTransSingleLayer
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.