R/Landscape_manipulation_functions.R

In tcarruth/SSES: Spatial Social Ecological Systems

# ==============================================================================================================================================================
# === Landscape manipulation functions =========================================================================================================================================
# ==============================================================================================================================================================





#' Subset the landscape object by lake
#'
#' @param obj an object of class 'Landscape'
#' @param ind integer vector: the lakes to subset by in the position they are listed in the lakenam slot
#' @author T. Carruthers
#' @export Lsub
Lsub<-function(obj,ind){
  obj@nl<-length(ind)
  obj@lakenam<-obj@lakenam[ind]
  obj@longnam<-obj@longnam[ind]
  obj@lakex<-obj@lakex[ind]
  obj@lakey<-obj@lakey[ind]
  obj@lakearea<-obj@lakearea[ind]
  dim0=dim(obj@lxslev)
  dim0[2]<-length(ind)
  obj@lxslev<-array(obj@lxslev[,ind,],dim0)
  obj@GDD<-t(matrix(obj@GDD[,ind]))
  obj@TDS<-t(matrix(obj@TDS[,ind]))
  obj@BagLim<-t(matrix(obj@BagLim[,ind]))
  obj@stockable<-obj@stockable[ind]
  dim0=dim(obj@pcxl)
  dim0[3]<-length(ind)
  obj@pcxl<-array(obj@pcxl[,,ind],dim0)
  dim0=dim(obj@lxattr)
  dim0[2]<-length(ind)
  obj@lxattr<-array(obj@lxattr[,ind,],dim0)
  obj@Scosts<-array(obj@Scosts[,ind,],dim0)
  dim0=dim(obj@exeff)
  dim0[2]<-length(ind)
  obj@exeff<-array(obj@exeff[,ind],dim0)
  obj@acc<-array(obj@acc[,ind],c(dim(obj@acc)[1],obj@nl))
  dim0<-dim(obj@costs)
  obj@costs<-array(obj@costs[,ind,,],c(dim0[1],obj@nl,dim0[3:4]))
  dim0<-dim(obj@eff)
  dim0[4]<-length(ind)
  if(length(dim(obj@eff))>2)obj@eff<-array(obj@eff[,,,ind,],dim0)
  obj
}



#' Add stochasticity to the landscape object
#'
#' For all landscape characteristics that are uncertain, these parameters are sampled.
#'
#' @param obj an object of class 'Landscape'
#' @param nsim integer: the number of simulations (stochastic samples of parameters)
#' @author T. Carruthers
#' @export Stoch
Stoch<-function(obj,nsim){

  # Dimensions
  nl<-obj@nl
  obj@nsim<-nsim
  npc<-obj@npc
  nattr<-obj@nattr
  na<-obj@na
  nage<-obj@nage
  nst<-obj@nst
  ncat<-obj@ncat

  # Create an nsim dimension to each attribute (the first dimension allowing for easy parallel computing)
  obj@pcsize<-array(geterr(nsim*npc,fac2num(obj@errs$pcsize[2]))*rep(obj@pcsize,each=nsim),dim=c(nsim,npc))
  #obj@pcxl<-array(geterr(nsim*npc*nl,fac2num(obj@errs$pcxl[2]))*rep(obj@pcxl,each=nsim),dim=c(nsim,npc,nl))
  errvec<-rep(geterr(nsim,fac2num(obj@errs$GDDamong[2])),nl)*rep(geterr(nl,fac2num(obj@errs$GDDacross[2])),each=nsim)
  obj@GDD<-array(errvec*rep(obj@GDD,each=nsim),c(nsim,nl))
  errvec<-rep(geterr(nsim,fac2num(obj@errs$TDSamong[2])),nl)*rep(geterr(nl,fac2num(obj@errs$TDSacross[2])),each=nsim)
  obj@TDS<-array(errvec*rep(obj@TDS,each=nsim),c(nsim,nl))

  obj@pcxa<-getmvlogerr(obj@pcxa,fac2num(obj@errs$pcxa[2]),nsim)
  obj@costs<-array(geterr(nsim*nl*(nattr-3)*max(fac2num(obj@attr[1,4:obj@nattr])),fac2num(obj@errs$cost[2]))*
                     rep(obj@costs,each=nsim),dim=c(nsim,nl,nattr-1,ncat))
  obj@Scosts<-array(geterr(nsim*nst*nl,fac2num(obj@errs$cost[2]))*
                      rep(obj@Scosts,each=nsim),dim=c(nsim,nl,nst))
  obj@effval<-array(geterr(nsim*na,fac2num(obj@errs$effval[2]))*rep(obj@effval,each=nsim),dim=c(nsim,na))
  obj@licval<-array(geterr(nsim,fac2num(obj@errs$licval[2]))*obj@licval,dim=c(nsim,1))
  obj@aq<-array(geterr(nsim*na,fac2num(obj@errs$aq[2]))*rep(obj@aq,each=nsim),dim=c(nsim,na))

  nofish<-1-obj@apr
  aprerr<-rnorm(nsim*na,0,fac2num(obj@errs$apr[2]))+log(rep(obj@apr,each=nsim))
  obj@apr<-array(exp(aprerr)/(exp(aprerr)+rep(nofish,each=nsim)),dim=c(nsim,na))

  # error in the way in which anglers are attracted to lake characteristics
  lnind<-(1:nattr)[obj@attr[2,]=="lognorm"]
  mvind<-(1:nattr)[obj@attr[2,]=="mvlogistic"]
  nln<-length(lnind)
  nmv<-length(mvind)
  obj@axattr<-array(rep(obj@axattr,each=nsim),dim=c(nsim,dim(obj@axattr)[2:4]))
  obj@axattr[,,lnind,]<-obj@axattr[,,lnind,]*geterr(nsim*na*nln*max(fac2num(obj@attr[1,])),rep(fac2num(obj@attr[3,lnind]),each=nsim*na))
  obj@axattr[,,mvind,]<-obj@axattr[,,mvind,]+rnorm(nsim*na*nmv*max(fac2num(obj@attr[1,])),0,rep(fac2num(obj@attr[3,mvind]),each=nsim*na))

  # error in the population variables
  npe<-ncol(obj@poperr)
  nind<-(1:npe)[obj@poperr[1,]=="norm"]
  lnind<-(1:npe)[obj@poperr[1,]=="lognorm"]
  nnm<-length(nind)
  nlog<-length(lnind)
  popvals<-array(NA,dim=c(nsim,ncol(obj@poperr),nst))
  popvals[,nind,]<-rnorm(nsim*nnm*nst,1,rep(fac2num(obj@poperr[2,nind]),each=nsim*nst))*rep(obj@popval[,nind,],each=nsim)
  popvals[,lnind,]<-getbias(nsim*nlog*nst,rep(fac2num(obj@poperr[2,lnind]),each=nsim*nst))*rep(obj@popval[,lnind,],each=nsim)
  #popvals<-as.data.frame(popvals)
  #names(popvals)<-names(obj@popval)
  obj@popval<-popvals

  obj@Mage<-array(getbias(nsim*nage*nst,fac2num(obj@errs$Mage[2]))*rep(obj@Mage,each=nsim),dim=c(nsim,nage,nst))
  obj

}


#' Landscape simplification (to ensure identifiability in IFD effort - a distinct lake x population centre distance matrix)
#'
#' For all landscape characteristics that are uncertain, these parameters are sampled.
#'
#' @param obj an object of class 'Landscape'
#' @param h postive real number, the height of the dendrogram for cutting (population centre simplification)
#' @author T. Carruthers
#' @export pc_simp
pc_simp<-function(obj,h=1000){

  npc<-obj@npc
  nl<-obj@nl
  mat<-array(NA,c(npc,npc))
  for(i in 1:npc){
    for(j in 1:npc){
      dif<-((obj@pcxl[1,i,]-obj@pcxl[1,j,])^2)^0.5
      mat[i,j]<-sum(dif)/nl
    }
  }
  rownames(mat)<-obj@pcnam
  mat<-as.dist(mat)
  #mat<-cbind(obj@pcx,obj@pcy)
  #rownames(mat)<-obj@pcnam
  out<-hclust(dist(mat),method="ward.D")
  plot(out)
  rect.hclust(out,h=h)
  groups<-cutree(out, h=h)

  npc<-max(groups)

  pcsize<-matrix(aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x,nrow=1)
  pcx<-aggregate(obj@pcx*obj@pcsize[1,],by=list(as.vector(groups)),sum)$x/aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x
  pcy<-aggregate(obj@pcy*obj@pcsize[1,],by=list(as.vector(groups)),sum)$x/aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x

  pcnam<-rep(NA,npc)
  szg<-aggregate(rep(1,length(groups)),by=list(as.vector(groups)),sum)$x
  for(i in 1:npc){
    if(szg[i]==1)pcnam[i]<-obj@pcnam[match(i, groups)]
    if(szg[i]>1)pcnam[i]<-paste(obj@pcnam[match(i, groups)],"_etal",sep="")
    #pcnam[i]<-obj@pcnam[match(i, groups)]
  }

  apr<-array(NA,c(1,npc,4))
  for(i in 1:4){
    apr[1,,i]<-aggregate(obj@apr[1,,i]*obj@pcsize[1,],by=list(as.vector(groups)),sum)$x/
      aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x
  }

  pcxl<-array(NA,c(1,npc,nl))
  for(i in 1:nl){
    pcxl[1,,i]<-aggregate(obj@pcxl[1,,i]*obj@pcsize[1,],by=list(as.vector(groups)),sum)$x/
      aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x
  }


  totno<-obj@pcxa[1,,]*obj@pcsize[1,]
  pcxa<-array(NA,c(1,npc,4))
  for(i in 1:max(groups)){
    ind<-(1:length(groups))[i==groups]
    if(length(ind)>1){
      temp<-apply(totno[ind,],2,sum)
    }else{
      temp<-totno[ind,]
    }
    pcxa[1,i,]<-temp/sum(temp)
  }
  obj@FTroutAng<-aggregate(obj@FTroutAng*obj@pcsize[1,],by=list(as.vector(groups)),sum)$x/
    aggregate(obj@pcsize[1,],by=list(as.vector(groups)),sum)$x
  obj@npc<-npc
  obj@pcsize<-pcsize
  obj@pcx<-pcx
  obj@pcy<-pcy
  obj@pcnam<-pcnam
  obj@apr<-apr
  obj@pcxl<-pcxl
  obj@pcxa<-pcxa
  obj

}


#' Calculate couch effect (the attractivity of other activities given a fitted landscape object)
#'
#' @param obj an object of class 'Landscape'
#' @author T. Carruthers
#' @export calccouch2
calccouch2<-function(obj){
  mspa<-as.matrix(expand.grid(1:obj@nmanage,1:obj@nsim,1:obj@npc,1:obj@na))
  U<-apply(obj@U,c(1,2,3,5),sum) # sum of utility across lakes
  maxdays<-as.matrix(obj@maxdays)
  couch<-array(NA,c(obj@nmanage,obj@nsim,obj@npc,obj@na))
  couch[mspa]<-((U[mspa]*unlist(maxdays[mspa[,c(2,4)]]))/as.numeric(obj@apr))-U[mspa]
  obj@couch=couch
  obj
}



#' Expand a Lansdscape object to mulitple management regimes
#'
#' @param obj an object of class 'Landscape'
#' @param nmanage positive integer greater than 1, the number of replicate management options
#' @param nsim positive integer greater than 1, the number of replicate simulations
#' @author T. Carruthers
#' @export replic8
replic8<-function(obj,nmanage=NA,nsim=NA){ # replicates the object over multiple management options

  if(!is.na(nmanage)){
    nl<-obj@nl
    npc<-obj@npc
    nst<-obj@nst
    obj@nmanage<-nmanage
    obj@lxattr<-array(rep(obj@lxattr,each=nmanage),c(nmanage,dim(obj@lxattr)[2:3]))
    obj@lxslev<-array(rep(obj@lxslev,each=nmanage),c(nmanage,dim(obj@lxslev)[2:3]))
    obj@pcxl<-array(rep(obj@pcxl,each=nmanage),c(nmanage,dim(obj@pcxl)[2:3]))
    obj@BagLim<-array(rep(obj@BagLim,each=nmanage),dim=c(nmanage,nl))
    obj@couch<-array(rep(obj@couch,each=nmanage),dim=c(nmanage,dim(obj@couch)[2:4]))
    if(length(dim(obj@eff))>2)obj@eff<-array(rep(obj@eff,each=nmanage),dim=c(nmanage,dim(obj@eff)[2:5]))

  }else if(!is.na(nsim)){

    obj@nsim<-nsim
    nl<-obj@nl
    npc<-obj@npc
    nattr<-obj@nattr
    na<-obj@na
    nage<-obj@nage
    nst<-obj@nst
    ncat<-obj@ncat

    obj@pcsize<-array(rep(obj@pcsize,each=nsim),dim=c(nsim,npc))
    obj@GDD<-array(rep(obj@GDD,each=nsim),dim=c(nsim,nl))
    obj@TDS<-array(rep(obj@TDS,each=nsim),dim=c(nsim,nl))
    obj@pcxa<-array(rep(obj@pcxa,each=nsim),dim=c(nsim,dim(obj@pcxa)[2:3]))
    #obj@costs<-array(rep(obj@costs,each=nsim),dim=c(nsim,dim(obj@costs)[2:4]))
    obj@Scosts<-array(rep(obj@Scosts,each=nsim),dim=c(nsim,dim(obj@Scosts)[2:3]))
    obj@effval<-array(rep(obj@effval,each=nsim),dim=c(nsim,na))
    obj@licval<-array(obj@licval,dim=c(nsim,1))
    obj@aq<-array(rep(obj@aq,each=nsim),dim=c(nsim,na))
    obj@apr<-array(rep(obj@apr,each=nsim),dim=c(nsim,dim(obj@apr)[2:3]))
    obj@axattr<-array(rep( obj@axattr,each=nsim),dim=c(nsim,dim( obj@axattr)[2:4]))
    obj@popval<-array(rep( obj@popval,each=nsim),dim=c(nsim,dim( obj@popval)[2:3]))
    obj@Mage<-array(rep(  obj@Mage,each=nsim),dim=c(nsim,dim(  obj@Mage)[2:3]))
    obj@acc<-array(rep(obj@acc,each=nsim),dim=c(nsim,nl))
    obj@DR<-array(rep(obj@DR,each=nsim),dim=c(nsim,na))
    obj@sel<-array(rep(obj@sel,each=nsim),dim=c(nsim,nage))
    obj@DD<-rep(obj@DD,nsim)
    obj@maxdays<-array(rep(obj@maxdays,each=nsim),dim=c(nsim,na))



  }else{
    message("You did not specify a number of replicate management options or replicate simulations!")
  }

  obj

}


#' Expand a Landscape object by copying an existing management scenario
#'
#' @param obj an object of class 'Landscape'
#' @param tocopy positive integer greater than 1, the number of replicate management options
#' @author T. Carruthers
#' @export AddMan
AddMan<-function(obj,tocopy=1){

  obj2<-obj
  nmanage<-obj@nmanage+1
  obj2@nmanage<-nmanage
  nl<-obj@nl
  npc<-obj@npc
  nst<-obj@nst

  tc_lxattr<-array(obj@lxattr[tocopy,,],c(1,dim(obj@lxattr)[2:3]))
  tc_lxslev<-array(obj@lxslev[tocopy,,],c(1,dim(obj@lxslev)[2:3]))
  tc_pcxl<-array(obj@pcxl[tocopy,,],c(1,dim(obj@pcxl)[2:3]))
  tc_BagLim<-array(obj@BagLim[tocopy,],c(1,nl))
  tc_couch<-array(obj@couch[tocopy,,,],c(1,dim(obj@couch)[2:4]))
  tc_eff<-array(obj@eff[tocopy,,,,],c(1,dim(obj@eff)[2:5]))

  obj2@lxattr<-abind(obj@lxattr,tc_lxattr,along=1)
  obj2@lxslev<-abind(obj@lxslev,tc_lxslev,along=1)
  obj2@pcxl<-  abind(obj@pcxl,  tc_pcxl,  along=1)
  obj2@BagLim<-abind(obj@BagLim,tc_BagLim,along=1)
  obj2@couch<- abind(obj@couch, tc_couch, along=1)
  obj2@eff<-   abind(obj@eff,   tc_eff,   along=1)

  obj2

}

#' Delete a management option from a Landscape object
#'
#' @param obj an object of class 'Landscape'
#' @param tocopy positive integer greater than 1, the number of replicate management options
#' @author T. Carruthers
#' @export DelMan
DelMan<-function(obj,todel=1){
  if(obj@nmanage>1){
    keepind<-((1:obj@nmanage)!=todel)

    obj2<-obj

    obj2@lxattr<-obj@lxattr[keepind,,,drop=F]
    obj2@lxslev<-obj@lxslev[keepind,,,drop=F]
    obj2@pcxl<-  obj@pcxl[keepind,,,drop=F]
    obj2@BagLim<-obj@BagLim[keepind,,drop=F]
    obj2@couch<- obj@couch[keepind,,,,drop=F]
    obj2@eff<-   obj@eff[keepind,,,,,drop=F]

    obj2@nmanage<-obj@nmanage-1

    obj2
  }

}



#' Returns the names of slots that have a dimension 'simulation'
#'
#' @param obj an object of class 'Landscape' that has only 1 simulation (nsim=1)
#' @author T. Carruthers
#' @export get_sim_slots
get_sim_slots<-function(obj){

  obj2<-replic8(obj,nsim=3)
  replist<-NULL
  slots<-slotNames(obj)
  for(i in 1:length(slots)){

    if(class(slot(obj,slots[i]))!="character"){ # not a character

      if(!is.null(dim(slot(obj,slots[i]))[1])){        # at least 2 dimensions

        if(!all(dim(slot(obj,slots[i]))==dim(slot(obj2,slots[i])))) replist<-c(replist,slots[i])

      }

    }

  }

  replist

}


#' Returns the names of slots that have a dimension 'management'
#'
#' @param obj an object of class 'Landscape' that has only 1 management policy (nmanage=1)
#' @author T. Carruthers
#' @export get_manage_slots
get_manage_slots<-function(obj){

  obj2<-replic8(obj,nmanage=3)
  replist<-NULL
  slots<-slotNames(obj)
  for(i in 1:length(slots)){

    if(class(slot(obj,slots[i]))!="character"){ # not a character

      if(!is.null(dim(slot(obj,slots[i]))[1])){        # at least 2 dimensions

        if(!all(dim(slot(obj,slots[i]))==dim(slot(obj2,slots[i])))) replist<-c(replist,slots[i])

      }

    }

  }

  replist

}