R/internals.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

.getSpainProv<-function() {
  SpainProv = c("01","02","03","04","05","06","07", "08","09","10",
                "11","12","13","14","15","16","17", "18","19","20",
                "21","22","23","24","25","26","27", "28","29","30",
                "31","32","33","34","35","36","37", "38","39","40",
                "41","42","43","44","45","46","47", "48","49","50")
  return(SpainProv)
}
.getProvinceFromID<-function(x) {
  return(substr(x,1, nchar(x)-4))
}
.getVarsRP<-function() {
  vars =c("RP_Abiesalba", "RP_Abiespinsapo", "RP_Fagussylvatica",
          "RP_Pinushalepensis", "RP_Pinusnigra", "RP_Pinuspinaster",
          "RP_Pinuspinea", "RP_Pinussylvestris", "RP_Pinusuncinata",
          "RP_Quercuscanariensis", "RP_Quercusfaginea", "RP_Quercushumilis",
          "RP_Quercusilex","RP_Quercuspetraea", "RP_Quercuspyrenaica",
          "RP_Quercusrobur", "RP_Quercussuber")
  return(vars)
}
.speciesModels<-function(x) {
  species_selection = get("species_selection")
  sp_list = strsplit(as.character(species_selection$species), split=",")
  ms = rep("", length(x))
  for(i in 1:length(sp_list)) {
    ms[as.character(x) %in% sp_list[[i]]] = i
  }
  return(ms)
}
.getVarsGrowth<-function() {
  vars = c("d","inv.d", "ln.d", "sq.d", "sqrt.d","BAL","BAL.ln.d", "BAL.ext", "ln.G", "Temp15sq","PPET", "SWHC", "Rad", "slope")
  return(vars)
}
.getVarsSurvival<-function(){
  return(c("d", "ln.d", "sq.d", "sqrt.d", "inv.d", "h", "BAL","BAL.ext", "ln.G","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsSurvivalPG<-function(){
  return(c("d", "ln.d", "sq.d", "sqrt.d", "inv.d", "h", "BAL","BAL.ext", "ln.d.inc.prev.ln.d", "ln.d.inc.prev","ln.G","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthB<-function(){
  return(c("Nsp.N","ln.N", "G", "ln.G","sd.DBH","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthN<-function(){
  return(c("Nsp.N","ln.N", "G", "ln.G", "sd.DBH", "Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthDBH<-function(){
  return(c("G", "ln.G", "sd.DBH","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthB_disp<-function(){
  return(c("Nsp.N_s0.5","Nsp.N_s1","Nsp.N_s2","Nsp.N_s4","ln.N", "G", "ln.G","sd.DBH","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthN_disp<-function(){
  return(c("Nsp.N_s0.5","Nsp.N_s1","Nsp.N_s2","Nsp.N_s4","ln.N", "G", "ln.G", "sd.DBH", "Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsIngrowthDBH_disp<-function(){
  return(c("G", "ln.G", "sd.DBH","Temp15sq","PPET", "SWHC", "Rad", "slope"))
}
.getVarsHeightGrowth<-function() {
  vars = c("inv.h","ln.h", "h.d","BAL.ln.d", "BAL.ext", "G", "Temp15sq","PPET", "SWHC", "Rad", "slope")
  return(vars)
}
.densityFactor<-function(d) {
  factor<-c(127.3239546, 31.83098865,14.14710607, 5.092958185)
  FACTOREXP = rep(NA, length(d))
  FACTOREXP[which(d<12.5)] = 1
  FACTOREXP[which(d>=12.5 & d<22.5)] = 2
  FACTOREXP[which(d>=22.5 & d<42.5)] = 3
  FACTOREXP[which(d>=42.5)] = 4
  return(factor[FACTOREXP])
}
.biomassModel1<-function(DBH, H, alpha, a, b, c) {
  return(as.numeric(alpha + a*(DBH^b)*(H^c)))
}
.biomassModel2<-function(DBH, H, alpha, t, u, v, w) {
  return(as.numeric(alpha + t*(DBH^u)+v*(H^w)))
}
.biomassModel3<-function(DBH, H, Zth, g,i,j,k,l,m,n) {
  Z = ifelse(DBH>Zth, 1,0)
  return(as.numeric(Z*(g*((i+DBH)^j)+k*((l+DBH)^m)*(H^n))))
}
.biomassModel4<-function(DBH, H, alpha, p, q) {
  return(as.numeric(alpha + p*(DBH^2)*H+ q*(DBH*H)))
}
.biomassModel5<-function(DBH, H, e, f) {
  return(as.numeric(e*(((DBH^2)*H)^f)))
}
.biomassModel6<-function(DBH, H, alpha, r, s, x) {
  return(as.numeric(alpha+ (r*(DBH^2))+ (s*H)+ ((x*(DBH^2))*H)))
}
.biomassModel7<-function(DBH, beta, gamma) {
  return(as.numeric(beta*exp(gamma*DBH)))
}

.biomass<-function(DBH, H, par) {
  model = par[["Model"]]
  if(model==1) {
    return(.biomassModel1(DBH, H, par[["alpha"]], par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==2) {
    return(.biomassModel2(DBH, H, par[["alpha"]], par[["t"]], par[["u"]], par[["v"]], par[["w"]]))
  }
  else if(model==3) {
    return(.biomassModel3(DBH, H, par[["Zth"]], par[["g"]], par[["i"]], par[["j"]], par[["k"]], par[["l"]], par[["m"]], par[["n"]]))
  }
  else if(model==4) {
    return(.biomassModel4(DBH, H, par[["alpha"]], par[["p"]], par[["q"]]))
  }
  else if(model==5) {
    return(.biomassModel5(DBH, H, par[["e"]], par[["f"]]))
  }
  else if(model==6) {
    return(.biomassModel6(DBH, H, par[["alpha"]], par[["r"]], par[["s"]], par[["x"]]))
  }
  else if(model==7) {
    return(.biomassModel7(DBH, par[["beta"]], par[["gamma"]]))
  }
  stop("Wrong parameter 'model' (must be an integer between 1 and 7)")
}

.getBiomassParams<-function(sp, compartment, area = NA) {
  sel = biomass_parameters$ID_TAXON==sp & biomass_parameters$Compartment==compartment
  if(!is.na(area)) {
    selArea = biomass_parameters$Area==area
    selArea[is.na(selArea)] = TRUE
    sel = sel & selArea
  }
  df = biomass_parameters[sel,]
  return(df)
}



.volumeModel1<-function(DBH, HT, a, b, c) {
  return(as.numeric(a + b*(DBH^2)*HT))
}
.volumeModel7<-function(VCC, a, b, c) {
  return(as.numeric(a + b*VCC + c*(VCC^2)))
}
.volumeModel8<-function(VCC, a, b, c) {
  return(as.numeric(a + b*VCC + c*(VCC^2)))
}
.volumeModel9<-function(DBH, q) {
  return(as.numeric(DBH^q))
}
.volumeModel10<-function(VCC, a, b, c) {
  return(as.numeric(a + b*VCC + c*(VCC^2)))
}
.volumeModel11<-function(DBH, HT, p, q, r) {
  return(as.numeric(p*(DBH^q)*(HT^r)))
}
.volumeModel12<-function(DBH, p, q) {
  return(as.numeric(p*(DBH^q)))
}
.volumeModel13<-function(DBH, a, b, dnm) {
  return(as.numeric(a+b*(DBH-dnm)))
}
.volumeModel14<-function(DBH, p, q) {
  return(as.numeric(p*(DBH^q)))
}
.volumeModel15<-function(DBH, a, b, dnm) { # TO DO!!!!
  return(as.numeric(NA))
}
.volumeModel16<-function(DBH, a, b) {
  return(as.numeric(a+b*(DBH^2)))
}
.volumeModel17<-function(DBH, a, b,c) {
  return(as.numeric(a+b*DBH+c*(DBH^2)))
}
.volumeModel18<-function(DBH, p, q) {
  return(as.numeric(p*exp(q*DBH)))
}
.volumeModel19<-function(DBH, a, b, c, d) {
  return(as.numeric(a + (b*DBH)+(c*(DBH^2))+ (d*(DBH^3))))
}
.volumeModel20<-function(DBH, a, b, d) {
  return(as.numeric(a + (b*DBH)+ (d*(DBH^3))))
}
.volumeModel21<-function(DBH, c, d) {
  return(as.numeric((c*(DBH^2))+ (d*(DBH^3))))
}
.volumeModel25<-function(DBH, HT, p, q, r) {
  return(as.numeric(p*(DBH^q)*(HT^r)))
}

.volume<-function(DBH, HT, VCC, par) {
  model = par[["Model"]]
  if(model==1) {
    return(.volumeModel1(DBH, HT, par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==7) {
    return(.volumeModel7(VCC, par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==8) {
    return(.volumeModel8(VCC, par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==9) { #No formula available!
    return(.volumeModel9(DBH, par[["r"]]))
  }
  else if(model==10) {
    return(.volumeModel10(VCC, par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==11) {
    return(.volumeModel11(DBH, HT, par[["p"]], par[["q"]], par[["r"]]))
  }
  else if(model==12) {
    return(.volumeModel12(DBH, par[["p"]], par[["q"]]))
  }
  else if(model==13) {
    return(.volumeModel13(DBH, par[["a"]], par[["b"]], par[["dnm"]]))
  }
  else if(model==14) {
    return(.volumeModel14(DBH, par[["p"]], par[["q"]]))
  }
  else if(model==15) {
    return(.volumeModel15(DBH, par[["a"]], par[["b"]], par[["dnm"]]))
  }
  else if(model==16) {
    return(.volumeModel16(DBH, par[["a"]], par[["b"]]))
  }
  else if(model==17) {
    return(.volumeModel17(DBH, par[["a"]], par[["b"]], par[["c"]]))
  }
  else if(model==18) {
    return(.volumeModel18(DBH, par[["p"]], par[["q"]]))
  }
  else if(model==19) {
    return(.volumeModel19(DBH, par[["a"]], par[["b"]], par[["c"]], par[["d"]]))
  }
  else if(model==20) {
    return(.volumeModel20(DBH, par[["a"]], par[["b"]], par[["d"]]))
  }
  else if(model==21) {
    return(.volumeModel21(DBH, par[["c"]], par[["d"]]))
  }
  else if(model==25) {
    return(.volumeModel25(DBH, HT, par[["p"]], par[["q"]], par[["r"]]))
  }
  stop(paste("Wrong parameter 'model' ", model))
}
.getVolumeParams<-function(ifn, prov, sp, param, fc, code_missing = "99") {
  sel = (volume_parameters$IFN %in% ifn) & (volume_parameters$Province==prov) & (volume_parameters$PARAM==param) & (volume_parameters$FC==fc)
  sel[is.na(sel)] = FALSE
  df = volume_parameters[sel,]
  if(nrow(df)>0) { # select species
    sel = rep(FALSE, nrow(df))
    spp_s = strsplit(as.character(df$ID_TAXON), split=",")
    for(i in 1:length(spp_s)) if(as.character(sp) %in% spp_s[[i]]) sel[i] = TRUE
    if(sum(sel)==0) sel[df$ID_TAXON==code_missing] = TRUE #Set otras frondosas (code = "99") to true
    return(df[sel,])
  }
  return(df)
}

.domheight<-function(h, n) {
  o <-order(h, decreasing=TRUE)
  h = h[o]
  n = n[o]
  ncum = 0
  for(i in 1:length(h)) {
    ncum = ncum + n[i]
    if(ncum>100) return(sum(h[1:i]*n[1:i])/sum(n[1:i]))
  }
  return(sum(h*n)/sum(n))
}
.domdiameter<-function(dbh, n) {
  o <-order(dbh, decreasing=TRUE)
  dbh = dbh[o]
  n = n[o]
  ncum = 0
  for(i in 1:length(dbh)) {
    ncum = ncum + n[i]
    if(ncum>100) return(sum(dbh[1:i]*n[1:i])/sum(n[1:i]))
  }
  return(sum(dbh*n)/sum(n))
}
.hb<-function(h, n) {
  return((100/.domheight(h,n))*sqrt(10000/sum(n)))
}
.diameterSD<-function(DBH, N) {
  if(length(N)==0) return(NA)
  sx2 = sum(N*(DBH^2), na.rm=T)
  sx = sum(DBH*N, na.rm=T)
  n = sum(N, na.rm=T)
  var = (sx2 - ((sx^2)/n))/(n-1)
  if(!is.na(var)) if(var>=0) return(sqrt(var))
  return(0)
}

.mergeTrees<-function(x) {
  mergeTreesSizeID<-function(x) {
    ntree = nrow(x)
    if(ntree>0) {
      BA = x$N*pi*(x$DBH/200)^2
      BAsp = tapply(BA, x$Species, FUN = sum)
      Nsp = as.numeric(tapply(x$N, x$Species, FUN = sum))
      DBHsp =  2*sqrt(10000*as.numeric(BAsp)/(pi*Nsp))
      y = data.frame(ID = rep(x$ID[1], length(BAsp)),
                     Species = names(BAsp),
                     DBH = DBHsp, N = Nsp,
                     row.names = 1:length(BAsp),
                     stringsAsFactors = FALSE)
      if("H" %in% names(x)) y$H = as.numeric(tapply(x$H*BA, x$Species, FUN = sum)/BAsp)
      if("DBHincprev" %in% names(x)) y$DBHincprev = as.numeric(tapply(x$DBHincprev*BA, x$Species, FUN = sum)/BAsp)
      if("Run" %in% names(x)) y <- y %>% add_column(Run = x$Run[1], .before = "ID")
      if("Step" %in% names(x)) y <- y %>% add_column(Step = x$Step[1], .before = "ID")
      if("Name" %in% names(x)) y <- y %>% add_column(Name = speciesNames(y$Species), .after = "Species")
      if("Ingrowth" %in% names(x)) y <- y %>% add_column(Ingrowth = FALSE, .after = "N")
      if("Colonization" %in% names(x)) y <- y %>% add_column(Colonization = FALSE, .after = "Ingrowth")
      return(y)
    }
    return(x)
  }
  mergeTreesID<-function(x) {
    sel0a = (x$DBH >= 47.5) & (x$DBH < 52.5) &(x$N < 5.1)
    sel0b = (x$DBH >= 42.5) & (x$DBH < 47.5) &(x$N < 5.1)
    sel1a = (x$DBH >= 37.5) & (x$DBH < 42.5) & (x$N < 14.1)
    sel1b = (x$DBH >= 32.5) & (x$DBH < 37.5) & (x$N < 14.1)
    sel1c = (x$DBH >= 27.5) & (x$DBH < 32.5) & (x$N < 14.1)
    sel1d = (x$DBH >= 22.5) & (x$DBH < 27.5) & (x$N < 14.1)
    sel2a = (x$DBH >= 17.5) & (x$DBH < 22.5) & (x$N < 31.8)
    sel2b = (x$DBH >= 12.5) & (x$DBH < 17.5) & (x$N < 31.8)
    sel3 = (x$DBH < 12.5) & (x$N < 127.3)
    if("Ingrowth" %in% names(x)) sel3 = sel3 & (x$Ingrowth == FALSE)
    nosel = !(sel0a | sel0b | sel1a | sel1b| sel1c| sel1d| sel2a| sel2b | sel3)
    y = x[nosel,, drop = FALSE]
    if(sum(sel0a)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel0a,, drop = FALSE]))
    if(sum(sel0b)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel0b,, drop = FALSE]))
    if(sum(sel1a)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel1a,, drop = FALSE]))
    if(sum(sel1b)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel1b,, drop = FALSE]))
    if(sum(sel1c)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel1c,, drop = FALSE]))
    if(sum(sel1d)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel1d,, drop = FALSE]))
    if(sum(sel2a)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel2a,, drop = FALSE]))
    if(sum(sel2b)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel2b,, drop = FALSE]))
    if(sum(sel3)>0) y <- bind_rows(y, mergeTreesSizeID(x[sel3,, drop = FALSE]))
    return(y)
  }
  IDs = unique(x$ID)
  y = NULL
  for(id in IDs) {
     if(is.null(y)) y = mergeTreesID(x[x$ID == id,,drop=FALSE])
     else y = bind_rows(y, mergeTreesID(x[x$ID == id,,drop=FALSE]))
  }
  return(y)
}