R/getPlotData.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

getRegionesProcedencia<-function(sp, RPfile="RegionesProcedencia/All_species.rds") {
  RPs = readRDS(RPfile)
  RPs_sp = as_Spatial(RPs)
  sp_2 = spTransform(sp, RPs_sp@proj4string)
  RPs_plots = over(sp_2,RPs_sp)
  for(i in 2:ncol(RPs_plots)) {
    RPs_plots[,i] = as.character(RPs_plots[,i])
  }
  return(RPs_plots)
}
getDynamicDataIFN<-function(plotDataIFN, ncfile, yearIni = 2001, yearFin = 2100,
                            includeInput = FALSE) {
  library(meteoland)
  IDs = row.names(plotDataIFN)


  numSteps = (yearFin-yearIni+1)/10
  if(includeInput) numSteps = numSteps+1

  dates = seq(as.Date(paste0(yearIni,"-01-01")), as.Date(paste0(yearFin, "-12-01")), by="month")

  cat(paste("A. Reading ncfile..."))
  sp = readNetCDFpoints(ncfile)
  dates_file = readNetCDFdates(ncfile)
  dates = dates[dates %in% dates_file]
  IDsfound = IDs[IDs %in% row.names(sp@coords)]
  if(length(IDsfound)<length(IDs)) cat(paste0(" [", length(IDsfound), "/", length(IDs)-length(IDsfound)," plots found / not found.]"))
  spm = readmeteorologypoints(ncfile, dates = dates, stations = IDsfound)
  cat("\n")

  cat(paste("B. Processing climate data...\n"))
  pb=txtProgressBar(0, 3, style=3)
  PETmat = extractvars(spm, "PET")@data
  setTxtProgressBar(pb, 1)
  Precmat = extractvars(spm, "Precipitation")@data
  setTxtProgressBar(pb, 2)
  Tempmat = extractvars(spm, "MeanTemperature")@data
  setTxtProgressBar(pb, 3)
  cat("\n")

  cat(paste("C. Calculating soil water deficit...\n"))
  #Piedallu, C., Gégout, J.-C., Perez, V., & Lebourgeois, F. 2012. Soil water balance performs better than climatic water variables in tree species distribution modelling. Global Ecology and Biogeography. doi: 10.1111/geb.12012
  nplots = length(IDsfound)
  nmonths = ncol(PETmat)
  SWHC =  plotDataIFN[IDsfound, "SWHC"]
  SWD = PETmat
  SWD[] = 0
  CWD = PETmat
  CWD[] = 0
  SWC = SWHC #Init SWC to SWHC
  AET = rep(NA, nplots)
  pb=txtProgressBar(1, nmonths, style=3)
  for(i in 1:nmonths) {
    setTxtProgressBar(pb, i)
    CD = Precmat[,i] - PETmat[,i]
    pos = (CD>= 0)
    pos[is.na(pos)] = TRUE
    SWCprev = SWC
    SWC[pos] = pmin(SWCprev[pos]+CD[pos], SWHC[pos]) #Fill soils if positive balance
    AET[pos] = PETmat[pos,i]
    SWC[!pos] = SWCprev[!pos]*exp(CD[!pos]/SWHC[!pos]) #Decrease soils if negative balance
    AET[!pos] = SWCprev[!pos] + (Precmat[!pos,i]-SWC[!pos]) #Water that actually evaporated including drawing from soil
    CWD[,i] =  pmax(0,- CD)
    SWD[,i] =  PETmat[,i]-AET
    SWD[is.na(SWD[,i]),i] = 0 #Set NA to zero
  }
  cat("\n")

  cat(paste("D. Arranging output...\n"))
  cutyears = cut(as.Date(colnames(PETmat)), breaks="years")
  Prec_year = t(apply(Precmat, 1, function(x) {tapply(x, cutyears, sum, na.rm=T)}))
  Temp_year = t(apply(Tempmat, 1, function(x) {tapply(x, cutyears, mean, na.rm=T)}))
  PET_year = t(apply(PETmat, 1, function(x) {tapply(x, cutyears, sum, na.rm=T)}))
  CSWD_year = t(apply(SWD, 1, function(x) {tapply(x, cutyears, sum, na.rm=T)}))

  nplots = length(IDs)

  plotDynamicData = data.frame(ID = rep(IDs, numSteps), Step = as.numeric(gl(numSteps, nplots)),
                               CSWD= NA, Prec = NA, Temp = NA)
  pb=txtProgressBar(1, numSteps, style=3)
  iniStep = 1
  if(includeInput) {
    iniStep = 2
    for(i in 1:length(IDs)) {
      id = IDs[i]
      if(id %in% IDsfound) {
        plotDynamicData[i, "CSWD"] = plotDataIFN[id, "CSWD"]
        plotDynamicData[i, "Prec"] = plotDataIFN[id, "Prec"]
        plotDynamicData[i, "Temp"] = plotDataIFN[id, "Temp"]
        plotDynamicData[i, "PET"] = plotDataIFN[id, "PET"]
      }
    }
  }
  for(s in iniStep:numSteps) {
    setTxtProgressBar(pb, s)
    if(includeInput) {
      yearIniStep = yearIni + 10*(s-2)
    } else {
      yearIniStep = yearIni + 10*(s-1)
    }
    yearFinStep = yearIniStep + 10

    yearSel = paste0(seq(yearIniStep, yearFinStep-1, by=1),"-01-01")
    yearSel = yearSel[yearSel%in% levels(cutyears)]
    for(i in 1:length(IDs)) {
      id = IDs[i]
      r = i + nplots*(s-1)
      if(id %in% IDsfound) {
        plotDynamicData[r, "CSWD"] = mean(CSWD_year[id,yearSel])
        plotDynamicData[r, "Prec"] = mean(Prec_year[id,yearSel])
        plotDynamicData[r, "Temp"] = mean(Temp_year[id,yearSel])
        plotDynamicData[r, "PET"] = mean(PET_year[id,yearSel])
      }
    }
  }
  return(plotDynamicData)
}