R/internalDynamics.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

.checkIDs<-function(treeData, plotData) {
  ids = unique(treeData$ID)
  mis = which(!(ids %in% rownames(plotData)))
  nmis = length(mis)
  if(nmis>0) {
    stop(paste0(nmis, " plot(s) mentioned in 'treeData' do not have environmental data as rows of 'plotData': ", paste(ids[mis], collapse = ",")))
  }
}
.checkSpeciesCodes<-function(treeData) {
  sp = unique(treeData$Species)
  spnames = speciesNamesModels(sp)
  mis = which(spnames=="")
  nmis = length(mis)
  if(nmis>0) {
    stop(paste0(nmis, " species code(s) mentioned in 'treeData' do not have models defined: ", paste(sp[mis], collapse = ",")))
  }
}
.IFNgrowthInternal<-function(dataInput, numYears, params= "defaultGrowthParams",
                             response = "DI", useProvince = TRUE, useRP = TRUE,
                             stochastic = FALSE, verbose=FALSE) {
  #Get params
  if(class(params)=="character") params = get(params)

  sp_vec = as.character(params$species)
  s = strsplit(sp_vec, split=",")

  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }
  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  if(verbose) {
    cat("\n Growth...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  INC = rep(NA, nrow(dataInput))
  vars = .getVarsGrowth()
  RP_vars = .getVarsRP()
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    p <- getParamList(spp[i])
    sel = (dataInput$Species==spp[i])

    if(!is.null(p)) {
      intercept = rep(as.numeric(p["Intercept"]), sum(sel))
      prsum = rep(0, sum(sel))
      rpsum = rep(0, sum(sel))
      if(useProvince) {
        v = dataInput$Province[sel]
        prstr = paste0("P", as.numeric(v))
        selpr = (!is.na(v))
        if(sum(selpr)>0) {
          selpr[selpr][!(prstr[selpr] %in% names(p))] = FALSE #Avoids Provinces that are not in p
          prsum[selpr] = as.numeric(p[prstr[selpr]])
          prsum[is.na(prsum)] = 0
        }
      }
      if(useRP) {
        for(j in 1:length(RP_vars)) {
          v = dataInput[[RP_vars[j]]][sel]
          rpstrj = paste0(RP_vars[j], v)
          selrp = (!is.na(v))
          if(sum(selrp)>0) {
            selrp[selrp][!(rpstrj[selrp] %in% names(p))] = FALSE #Avoids RP that are not in p
            rpsumj = as.numeric(p[rpstrj[selrp]])
            rpsumj[is.na(rpsumj)] = 0
            rpsum[selrp]  = rpsum[selrp] + rpsumj
          }
        }
      }
      intercept = intercept+prsum+rpsum
      # print(intercept)
      if(sum(sel)>1) {
        ss = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
      } else {
        ss = sum(as.numeric(dataInput[sel,vars])*as.numeric(p[vars]), na.rm=T)
      }
      INC[sel] = exp(intercept+ss)*(numYears[sel]/10) #Original model calibrated for 10 yr
      if(stochastic) {
        pred = INC[sel]
        gd = as.numeric(p["gamma_disp"])
        for(j in 1:length(pred)) {
          pred[j] = rgamma(1, shape = 1/gd, scale = pred[j]*gd)
        }
        INC[sel] = pred
      }
    } else {
      INC[sel] = NA
    }
  }
  if(sum(is.na(INC))>0){
    stop("NA INC")
  }

  # BAI model
  # if(response=="DI") res = 2*sqrt((INC+pi*(dataInput$d/2)^2)/pi) - dataInput$d
  # else if(response=="DBH") res = 2*sqrt((INC+pi*(dataInput$d/2)^2)/pi)
  # else res = INC
  # DI model
  if(response=="BAI") res = (pi*((dataInput$d+INC)/2)^2)-(pi*(dataInput$d/2)^2)
  else if(response=="DBH") res = dataInput$d+INC
  else res = INC
  return(res)
}

.IFNsurvivalInternal<-function(dataInput, numYears, params = "defaultSurvivalParams",
                               killVeryLargeTrees = TRUE, verbose = FALSE) {
  #Get params
  if(class(params)=="character") params = get(params)
  sp_vec = as.character(params$species)
  s = strsplit(sp_vec, split=",")

  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }

  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  species_selection = get("species_selection")
  getMaxDiameter<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(species_selection$max_diameter[i])
    return(1000000.0)
  }


  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  nr = nrow(dataInput)
  if(verbose) {
    cat("\n Survival...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  Psurv = rep(NA, nr)
  vars = .getVarsSurvival()
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    p <- getParamList(spp[i])
    sel = (dataInput$Species==spp[i])
    if(!is.null(p)) {
      if(sum(sel)>1) {
        ss = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
      } else {
        ss = sum(as.numeric(dataInput[sel,vars])*as.numeric(p[vars]), na.rm=T)
      }
      Psurv[sel] = 1/(1+exp(-1*(as.numeric(p["Intercept"])+ss)))
      if(killVeryLargeTrees) {
        Psurv[dataInput$d[sel] > getMaxDiameter(spp[i])] = 0
      }
    } else {
      Psurv[sel] = NA
    }
  }
  Psurv = Psurv^(numYears/10)
  return(Psurv)
}

.IFNsurvivalPGInternal<-function(dataInput, numYears, params = "defaultSurvivalPGParams",
                                 killVeryLargeTrees = TRUE, verbose = FALSE) {
  #Get params
  if(class(params)=="character") params = get(params)
  sp_vec = as.character(params$species)
  s = strsplit(sp_vec, split=",")

  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }

  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  species_selection = get("species_selection")
  getMaxDiameter<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(species_selection$max_diameter[i])
    return(1000000.0)
  }


  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  nr = nrow(dataInput)
  if(verbose) {
    cat("\n Survival (past growth)...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  Psurv = rep(NA, nr)
  vars = .getVarsSurvivalPG()
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    sel = (dataInput$Species==spp[i])
    p <- getParamList(spp[i])
    if(!is.null(p)) {
      if(sum(sel)>1) {
        ss = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars, drop=FALSE]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
      } else {
        ss = sum(as.numeric(dataInput[sel,vars, drop=FALSE])*as.numeric(p[vars]), na.rm=T)
      }
      Psurv[sel] = 1/(1+exp(-1*(as.numeric(p["Intercept"])+ss)))
      if(killVeryLargeTrees) {
        Psurv[dataInput$d[sel] > getMaxDiameter(spp[i])] = 0
      }
    } else { # Species not found
      Psurv[sel] = NA
    }
  }
  Psurv = Psurv^(numYears/10)
  return(Psurv)
}

.IFNingrowthInternal<-function(dataInput, numYears,
                               paramsB= "defaultIngrowthBParams",
                               paramsN= "defaultIngrowthNParams",
                               paramsDBH  = "defaultIngrowthDBHParams",
                               thresholdIngrowth = FALSE,
                               stochastic = FALSE, verbose=FALSE) {
  #Get params
  if(class(paramsB)=="character") paramsB = get(paramsB)
  if(class(paramsN)=="character") paramsN = get(paramsN)
  if(class(paramsDBH)=="character") paramsDBH = get(paramsDBH)


  B_sp_vec = as.character(paramsB$species)
  B_s = strsplit(B_sp_vec, split=",")
  getBParamList<-function(sp) {
    for(i in 1:length(B_s)) if(as.character(sp) %in% B_s[[i]]) return(as.list(paramsB[i,-c(1,2)]))
    return(NULL)
  }

  N_sp_vec = as.character(paramsN$species)
  N_s = strsplit(N_sp_vec, split=",")
  getNParamList<-function(sp) {
    for(i in 1:length(N_s)) if(as.character(sp) %in% N_s[[i]]) return(as.list(paramsN[i,-c(1,2)]))
    return(NULL)
  }
  DBH_sp_vec = as.character(paramsDBH$species)
  DBH_s = strsplit(DBH_sp_vec, split=",")
  getDBHParamList<-function(sp) {
    for(i in 1:length(DBH_s)) if(as.character(sp) %in% DBH_s[[i]]) return(as.list(paramsDBH[i,-c(1,2)]))
    return(NULL)
  }
  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  nr = nrow(dataInput)
  if(verbose) {
    cat("\n Recruiment density and diameter...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  vars.B = .getVarsIngrowthB()
  vars.N = .getVarsIngrowthN()
  vars.DBH = .getVarsIngrowthDBH()
  B = rep(NA, nr)
  P = rep(NA, nr)
  N = rep(NA, nr)
  DBH = rep(NA, nr)
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    B_p <- getBParamList(spp[i])
    N_p <- getNParamList(spp[i])
    DBH_p <- getDBHParamList(spp[i])
    sel = (dataInput$Species==spp[i])
    if(!is.null(N_p)) {
      if(sum(sel)>1) {
        ssB = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.B]),2,as.numeric(B_p[vars.B]), FUN="*"), na.rm=T))
        ssN = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.N]),2,as.numeric(N_p[vars.N]), FUN="*"), na.rm=T))
        ssDBH = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.DBH]),2,as.numeric(DBH_p[vars.DBH]), FUN="*"), na.rm=T))
      } else {
        ssB = sum(as.numeric(dataInput[sel,vars.B])*as.numeric(B_p[vars.B]), na.rm=T)
        ssN = sum(as.numeric(dataInput[sel,vars.N])*as.numeric(N_p[vars.N]), na.rm=T)
        ssDBH = sum(as.numeric(dataInput[sel,vars.DBH])*as.numeric(DBH_p[vars.DBH]), na.rm=T)
      }
      P[sel] = 1/(1+exp(-1*(as.numeric(B_p["Intercept"])+ssB)))
      B[sel] = ifelse(P[sel]>B_p["CutOff"],1,0)
      N[sel] = exp(as.numeric(N_p["Intercept"])+ssN)
      DBH[sel] = exp(as.numeric(DBH_p["Intercept"])+ssDBH)
      if(stochastic) {
        P[sel] = rbinom(sum(sel),1,prob = P[sel])
        predN = N[sel]
        predDBH = DBH[sel]
        gdN = as.numeric(N_p["gamma_disp"])
        gdDBH = as.numeric(DBH_p["gamma_disp"])
        for(j in 1:length(predN)) {
          predN[j] = rgamma(1, shape = 1/gdN, scale = predN[j]*gdN)
          predDBH[j] = rgamma(1, shape = 1/gdDBH, scale = predDBH[j]*gdDBH)
        }
        DBH[sel] = predDBH
        N[sel] = predN
      }
    } else { # Species not found
      DBH[sel] = NA
      N[sel] = NA
    }
  }
  N = N*(numYears/10)
  if(thresholdIngrowth) {N = B*N} # B is either 0 or 1, so that potential stochasticity is only included in N
  else {N = P*N}
  df = data.frame(ID = dataInput$ID, Species = dataInput$Species, N = N, DBH = DBH)
  df = df[df$N>0,]
  return(df)
}

.IFNingrowthdispInternal<-function(dataInput, numYears,
                               paramsB= "defaultIngrowthBParams_disp",
                               paramsN= "defaultIngrowthNParams_disp",
                               paramsDBH  = "defaultIngrowthDBHParams_disp",
                               thresholdIngrowth = FALSE,
                               stochastic = FALSE, verbose=FALSE) {
  #Get params
  if(class(paramsB)=="character") paramsB = get(paramsB)
  if(class(paramsN)=="character") paramsN = get(paramsN)
  if(class(paramsDBH)=="character") paramsDBH = get(paramsDBH)


  B_sp_vec = as.character(paramsB$species)
  B_s = strsplit(B_sp_vec, split=",")
  getBParamList<-function(sp) {
    for(i in 1:length(B_s)) if(as.character(sp) %in% B_s[[i]]) return(as.list(paramsB[i,-c(1,2)]))
    return(NULL)
  }

  N_sp_vec = as.character(paramsN$species)
  N_s = strsplit(N_sp_vec, split=",")
  getNParamList<-function(sp) {
    for(i in 1:length(N_s)) if(as.character(sp) %in% N_s[[i]]) return(as.list(paramsN[i,-c(1,2)]))
    return(NULL)
  }
  DBH_sp_vec = as.character(paramsDBH$species)
  DBH_s = strsplit(DBH_sp_vec, split=",")
  getDBHParamList<-function(sp) {
    for(i in 1:length(DBH_s)) if(as.character(sp) %in% DBH_s[[i]]) return(as.list(paramsDBH[i,-c(1,2)]))
    return(NULL)
  }
  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  nr = nrow(dataInput)
  if(verbose) {
    cat("\n Ingrowth density and diameter (including dispersal) ...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  vars.B = .getVarsIngrowthB_disp()
  vars.N = .getVarsIngrowthN_disp()
  vars.DBH = .getVarsIngrowthDBH_disp()
  B = rep(NA, nr)
  P = rep(NA, nr)
  N = rep(NA, nr)
  DBH = rep(NA, nr)
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    B_p <- getBParamList(spp[i])
    N_p <- getNParamList(spp[i])
    DBH_p <- getDBHParamList(spp[i])
    sel = (dataInput$Species==spp[i])
    if(!is.null(N_p)) {
      if(sum(sel)>1) {
        ssB = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.B]),2,as.numeric(B_p[vars.B]), FUN="*"), na.rm=T))
        ssN = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.N]),2,as.numeric(N_p[vars.N]), FUN="*"), na.rm=T))
        ssDBH = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars.DBH]),2,as.numeric(DBH_p[vars.DBH]), FUN="*"), na.rm=T))
      } else {
        ssB = sum(as.numeric(dataInput[sel,vars.B])*as.numeric(B_p[vars.B]), na.rm=T)
        ssN = sum(as.numeric(dataInput[sel,vars.N])*as.numeric(N_p[vars.N]), na.rm=T)
        ssDBH = sum(as.numeric(dataInput[sel,vars.DBH])*as.numeric(DBH_p[vars.DBH]), na.rm=T)
      }
      P[sel] = 1/(1+exp(-1*(as.numeric(B_p["Intercept"])+ssB)))
      B[sel] = ifelse(P[sel]>B_p["CutOff"],1,0)
      N[sel] = exp(as.numeric(N_p["Intercept"])+ssN)
      DBH[sel] = exp(as.numeric(DBH_p["Intercept"])+ssDBH)
      if(stochastic) {
        P[sel] = rbinom(sum(sel),1,prob = P[sel])
        predN = N[sel]
        predDBH = DBH[sel]
        gdN = as.numeric(N_p["gamma_disp"])
        gdDBH = as.numeric(DBH_p["gamma_disp"])
        for(j in 1:length(predN)) {
          predN[j] = rgamma(1, shape = 1/gdN, scale = predN[j]*gdN)
          predDBH[j] = rgamma(1, shape = 1/gdDBH, scale = predDBH[j]*gdDBH)
        }
        DBH[sel] = predDBH
        N[sel] = predN
      }
    } else { # Species not found
      DBH[sel] = NA
      N[sel] = NA
    }
  }
  N = N*(numYears/10)
  if(thresholdIngrowth) {N = B*N} # B is either 0 or 1, so that potential stochasticity is only included in N
  else {N = P*N} # potential stochasticity is only included in N
  df = data.frame(ID = dataInput$ID, Species = dataInput$Species, N = N, DBH = DBH)
  df = df[df$N>0,]
  return(df)
}

.IFNheightInternal<-function(dataInput, params = "defaultHeightParams",
                             useProvince = TRUE, useRP = TRUE, verbose = FALSE) {

  #Get params
  if(class(params)=="character") params = get(params)

  sp_vec = as.character(params$species)
  s = strsplit(sp_vec, split=",")
  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }
  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  nr = nrow(dataInput)
  if(verbose) {
    cat("\n Calculating height...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }

  RP_vars = .getVarsRP()

  H = rep(NA, nr)
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    p <- getParamList(spp[i])
    if(!is.null(p)) {
      sel = (dataInput$Species==spp[i])
      prsum = rep(0, sum(sel))
      rpsum = rep(0, sum(sel))
      numerator = p[["A0"]] + p[["B1"]]*dataInput$PPET[sel]+p[["B2"]]*dataInput$SWHC[sel] + p[["B3"]]*dataInput$Rad[sel]+p[["B4"]]*dataInput$elevation[sel]
      if(useProvince) {
        v = dataInput$Province[sel]
        prstr = paste0("E", as.numeric(v))
        selpr = (!is.na(v))
        if(sum(selpr)>0) {
          selpr[selpr][!(prstr[selpr] %in% names(p))] = FALSE #Avoids Provinces that are not in p
          prsum[selpr] = as.numeric(p[prstr[selpr]])
          prsum[is.na(prsum)] = 0
        }
      }
      if(useRP) {
        for(j in 1:length(RP_vars)) {
          v = dataInput[[RP_vars[j]]][sel]
          rpstrj = paste0("F",RP_vars[j], v)
          selrp = (!is.na(v))
          if(sum(selrp)>0) {
            selrp[selrp][!(rpstrj[selrp] %in% names(p))] = FALSE #Avoids RP that are not in p
            rpsumj = as.numeric(p[rpstrj[selrp]])
            rpsum[selrp]  = rpsum[selrp] + rpsumj
          }
        }
      }
      numerator = numerator + prsum + rpsum
      H[sel] = (numerator)/(1+(p[["C1"]]/dataInput$DBH[sel])+(p[["C2"]]/(dataInput$DBH[sel]^2)))
      if(sum(is.na(H[sel]))>0) {
        warning(paste("Warning: NA values in height estimation for",sum(is.na(H[sel])), "records."))
      } else {
        if(sum(H[sel]<0)>0) {
          warning(paste0(sum(H[sel]<0)," negative height values truncated to zero."))
          h = H[sel]
          h[h<0] = 0
          H[sel] = h
        }
      }
    } else {#Species not found
      H[sel] = NA
    }
  }
  return(H)
}

.IFNheightgrowthInternal<-function(dataInput, numYears, params= "defaultHeightGrowthParams",
                                   useProvince = TRUE, useRP = TRUE,
                                   stochastic = FALSE, verbose=FALSE) {
  #Get params
  if(class(params)=="character") params = get(params)
  HI_byProvince = get("HI_byProvince")


  sp_vec = as.character(params$species)
  s = strsplit(sp_vec, split=",")

  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }
  getHI_by_ProvinceList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(HI_byProvince[i,-c(1,2)]))
    return(NULL)
  }
  spp = unique(as.character(dataInput$Species))
  ns = length(spp)
  if(length(numYears)==1) numYears = rep(numYears, nrow(dataInput))

  if(verbose) {
    cat("\n Height growth...\n")
    if(ns>1) pb = txtProgressBar(1, ns, style=3)
  }
  INC = rep(NA, nrow(dataInput))
  vars = .getVarsHeightGrowth()
  RP_vars = .getVarsRP()
  for(i in 1:ns) {
    if(verbose && ns>1) setTxtProgressBar(pb,i)
    p <- getParamList(spp[i])
    sel = (dataInput$Species==spp[i])

    if(!is.null(p)) {
      intercept = rep(as.numeric(p["Intercept"]), sum(sel))
      prsum = rep(0, sum(sel))
      rpsum = rep(0, sum(sel))
      if(useProvince) {
        v = dataInput$Province[sel]
        prstr = paste0("P", as.numeric(v))
        selpr = (!is.na(v))
        if(sum(selpr)>0) {
          selpr[selpr][!(prstr[selpr] %in% names(p))] = FALSE #Avoids Provinces that are not in p
          prsum[selpr] = as.numeric(p[prstr[selpr]])
          prsum[is.na(prsum)] = 0
        }
      }
      if(useRP) {
        for(j in 1:length(RP_vars)) {
          v = dataInput[[RP_vars[j]]][sel]
          rpstrj = paste0(RP_vars[j], v)
          selrp = (!is.na(v))
          if(sum(selrp)>0) {
            selrp[selrp][!(rpstrj[selrp] %in% names(p))] = FALSE #Avoids RP that are not in p
            rpsumj = as.numeric(p[rpstrj[selrp]])
            rpsumj[is.na(rpsumj)] = 0
            rpsum[selrp]  = rpsum[selrp] + rpsumj
          }
        }
      }
      intercept = intercept+prsum+rpsum
      # print(intercept)
      if(sum(sel)>1) {
        ss = as.numeric(rowSums(sweep(as.matrix(dataInput[sel,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
      } else {
        ss = sum(as.numeric(dataInput[sel,vars])*as.numeric(p[vars]), na.rm=T)
      }
      INC[sel] = exp(intercept+ss)*(numYears[sel]/10) #Original model calibrated for 10 yr
      if(stochastic) {
        pred = INC[sel]
        gd = as.numeric(p["gamma_disp"])
        for(j in 1:length(pred)) {
          pred[j] = rgamma(1, shape = 1/gd, scale = pred[j]*gd)
        }
        INC[sel] = pred
      }
    } else {
      INC[sel] = NA
    }

    hi_byProv = getHI_by_ProvinceList(spp[i])
    prstr = paste0("P", as.numeric(as.numeric(dataInput$Province[sel])))
    allowed = as.logical(hi_byProv[prstr[selpr]])
    INC[sel][!allowed] = NA #Set not allowed combination of species-province to NA
  }
  res = INC

  return(res)
}