#' Area-based clear and partial cuts
#'
#' Selects cells for clear and partial cuts on an area-based
#'
#' @param land A \code{landscape} data frame with forest stand records in rows
#' @param params A list of default parameters generated by the function \code{default.params()} or
#' a customized list of model parameters
#' @param cc.area A data frame with the number of cells to clear cut per mgmt unit returned by the function \code{timber.area()}
#' @param pc.area A data frame with the number of cells to partial cut per mgmt unit returned by the function \code{timber.partial.area()}
#' @param km2.pixel Area in km^2^ of grid cells
#'
#' @return A list of four items:
#' \itemize{
#' \item{\code{cc.cells}: A vector with the \code{cell.id} of the clear cut cells}
#' \item{\code{pc.cells}: A vector with the \code{cell.id} of the partial cut cells}
#' \item{\code{track.cut}: A data frame with managed areas by different prescriptions}
#' \item{\code{spp.track}: A data frame with the area and volume clear-cut and partial cut per species}
#' }
#'
#' @export
#'
#' @examples
#' data(landscape)
#' params = default.params()
#' cc.area = timber.area(landscape, params)
#' pc.area = timber.partial.area(landscape, params, pc.step=5)
#' harvest = harvest.area(landscape, params, cc.area, pc.area, km2.pixel=4)
#'
harvest.area <- function(land, params, cc.area, pc.area, km2.pixel){
cat(" Area-based selection of clearcut and partial cut cells", "\n" )
land2 <- land[!is.na(land$mgmt.unit) & !land$spp=="NonFor",]
land2$vol <- stand.volume(land2)*km2.pixel*100
units <- as.character(sort(unique(land2$mgmt.unit[!is.na(land2$mgmt.unit)])))
units <- units[-which(units == c("AGNA"))]
s.inc <- filter(land2, !is.na(mgmt.unit) & is.na(exclus)) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
s.inc.mat <- filter(land2, !is.na(mgmt.unit) & is.na(exclus) & age>age.matu) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
## For those locations that can be harvested (included), differentiate those that have been burnt or killed
## by an outbreak, and then count the young (cannot be salvaged) vs the mature (can be salvaged)
s.inc.burnt <- filter(land2, !is.na(mgmt.unit) & is.na(exclus) & tsfire==0) %>%
group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
s.inc.mat.burnt <- filter(land2, !is.na(mgmt.unit) & is.na(exclus) & tsfire==0 & age>=age.matu) %>%
group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
s.inc.kill <- filter(land2, !is.na(mgmt.unit) & is.na(exclus) & tssbw%in%c(0,5)) %>%
group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
s.inc.mat.kill <- filter(land2, !is.na(mgmt.unit) & is.na(exclus) & tssbw%in%c(0,5) & age>=age.matu) %>%
group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
## Also, look for zones at defforestation risk, both included and excluded
reg.fail.ex <- filter(land2, !is.na(mgmt.unit) & spp %in% c("EPN", "SAB", "OTH.RES.N"),
tsfire==0, age<=params$age.seed, temp < 1.5, runif(length(land2$temp))<params$p.failure, !is.na(exclus)) %>%
group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
reg.fail.inc <- filter(land2, !is.na(mgmt.unit) & spp %in% c("EPN", "SAB", "OTH.RES.N"),
tsfire==0, age<=params$age.seed, temp < 1.5, runif(length(land2$temp))<params$p.failure, is.na(exclus)) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
#land2 <- mutate(land2, rndm=runif(nrow(land2)))
#even <- land2$spp %in% c("EPN", "PET", "SAB", "OthCB", "OthCT", "OthDB") & is.na(land2$exclus) & land2$rndm<=0.95
#sum(even)
#even[land2$spp %in% c("BOJ", "ERS", "OthDT")& is.na(land2$exclus) & land2$rndm>0.95] <- 1
land2$even[land2$tsfire==0 & land2$spp!="NonFor"] <- 1 # NO change for NonFor cells
land.coniferes <- land2[land2$even==1,]
land.feuillu.tol <- land2[land2$even==0,]
land.ea <- land.coniferes
s.ea <- group_by(land.ea, mgmt.unit) %>% summarise(x=length(mgmt.unit))
### possibilité
poss.init <- cc.area
## Subset the mature even-aged cells from those that are harvestable
land.rec <- filter(land.ea, age>=age.matu)
# initialisation des variables
cc.cells.salv.tot <- cc.cells.unaff.tot <- cc.cells <- numeric(0)
# selection
for(unit in units){
harv.level.u <- poss.init[poss.init$mgmt.unit == unit,2]
land.ea.u <- land.ea[land.ea$mgmt.unit==unit,]
s.ea.u <- length(land.ea.u$cell.id)
# Subset of harvestable (mature even-aged) cells
land.ea.mat.u <- land.ea.u[land.ea.u$age >= land.ea.u$age.matu,]
subland.salv.mature.burn <- land.ea.u[(land.ea.u$age >= (land.ea.u$age.matu-params$diff.prematurite)) &
land.ea.u$tsfire ==0, ]
# table(land.ea.u$tsfire)
liste.recuperables.FMU <- sample(subland.salv.mature.burn$cell.id, round(params$salvage.rate.FMU*nrow(subland.salv.mature.burn)))
subland.salv.mature.burn2 <- subland.salv.mature.burn[subland.salv.mature.burn$cell.id %in% liste.recuperables.FMU,]
if(length(subland.salv.mature.burn2$cell.id) > round(params$salvage.rate.event* harv.level.u)) {
liste.2 <- sample(subland.salv.mature.burn2$cell.id,round((params$salvage.rate.event* harv.level.u)), replace=FALSE)
subland.salv.mature.burn3 <- subland.salv.mature.burn2[subland.salv.mature.burn2$cell.id %in% liste.2,]
} else {
subland.salv.mature.burn3 <- subland.salv.mature.burn2
}
subland.salv.mature.sbw <- land.ea.u[(land.ea.u$age >= (land.ea.u$age.matu-params$diff.prematurite)) &
land.ea.u$tssbw %in% c(0,5), ]
subland.salv.mature <- rbind(subland.salv.mature.burn3,subland.salv.mature.sbw)
# sélection de cellules récupérables en tenant compte
# des contraintes a priori (maximum salvage rate, etc.)
cell.salv.available <- subland.salv.mature$cell.id #sample(subland.salv.mature$cell.id, round(salvage.rate.event*nrow(subland.salv.mature)), replace=FALSE)
#############################################
# Randomly select cells among the even-aged mature cells present in non-protected areas
# Prioritize clear cuts in disturbed areas (salvage logging)
x <- length(cell.salv.available )# x = disponible
cc.cells.salv <- numeric(0)
xx <-0 # récolté
while (x > 0 & length(cell.salv.available)>0 & xx < harv.level.u) {
paquet <- ifelse((x <= harv.level.u) ,x,harv.level.u ) # pour accélérer le calcul, paquets de 5 cellules
cc.cells.salv.x <- sample(cell.salv.available,paquet)#cell.salv.available[1:paquet] #
cc.cells.salv <- c(cc.cells.salv,cc.cells.salv.x)
cell.salv.available <- cell.salv.available[-which(cell.salv.available%in%cc.cells.salv.x)]
# fonction neutralisée pour le monent car fonctionne mal (vol=NA pour certaines cellules)
x <- 1 #sum(subland.salv.mature[subland.salv.mature$cell.id %in%cell.salv.available,]$vol )
xx <- length(cc.cells.salv )
}
# When salvaged cells were not enough to satisfy sustained yield level, then harvest some mature
# forests unaffected by disturbances (cc.cells.unaff).
subland.non.pertu <- land.ea.mat.u[land.ea.mat.u$tsfire!=0 & land.ea.mat.u$tssbw>5, ]
x <- length(subland.non.pertu$cell.id ) # x = disponible
cc.cells.unaff <- numeric(0)
# arrête la récolte lorsqu'on est rendu à < 40000m3 du but, il y aura parfois des dépassements
while(x > 0 & length(subland.non.pertu$cell.id)>0 & xx < (harv.level.u)) {
paquet <- ifelse(harv.level.u-xx >= 0 ,min(x,harv.level.u-xx) ,0 )
cc.cells.unaff.x <- sample( subland.non.pertu$cell.id, paquet)
cc.cells.unaff <- c(cc.cells.unaff,cc.cells.unaff.x)
subland.non.pertu <- subland.non.pertu[-which(subland.non.pertu$cell.id%in%cc.cells.unaff.x),]
xx <- length(c(cc.cells.unaff,cc.cells.salv ))
}
# Combine all types of harvested cells with those already harvested in other FMUs during the same period
cc.cells <- c(cc.cells, cc.cells.salv, cc.cells.unaff)
cc.cells.salv.tot <- c(cc.cells.salv.tot, cc.cells.salv)
cc.cells.unaff.tot <- c(cc.cells.unaff.tot, cc.cells.unaff)
}
################## partial cuts ############################
land.uea <- land.feuillu.tol
s.uea <- group_by(land.uea, mgmt.unit) %>% summarise(x=length(mgmt.unit))
# volume par cellule. Moitié du volume accessible
land.uea$vol <- land.uea$vol/2
## The maturity age for partial cuts is half the maturity age for a clear cut
land.uea$age.matu.pc <- round(land.uea$age.matu,-1)/2
## Subset of harvestable (i.e. mature uneven-aged, ot recently partial cut) cells
land.rec.pc <- filter(land.uea, tspcut >=age.matu.pc )
## Get the number of cells to be managed under a partial-cut regime
s.uea <- group_by(land.uea, mgmt.unit) %>% summarise(x=length(mgmt.unit))
s.mat.pc <- group_by(land.rec.pc, mgmt.unit) %>% summarise(x=length(mgmt.unit))
harv.level.pc <- pc.area # read.table("InitialVolumePC.txt", header=T)
pc.cells <- numeric()
for(unit in unique(land.uea$mgmt.unit)){ #unit=9351
poss.cp.ua <- harv.level.pc$x[harv.level.pc$mgmt.unit==unit]
poss.cp.ua <- ifelse(length(poss.cp.ua) == 0,0,poss.cp.ua)
cell.dispo.ua <- land.rec.pc[land.rec.pc$mgmt.unit==unit, ]
x <- length(cell.dispo.ua$cell.id ) #disponible
pc.cells.ua <- numeric(0)
xx <- 0
# arrête la récolte lorsqu'on est rendu à < 20000m3 du but, il y aura parfois des dépassements
while(x > 0 & length(cell.dispo.ua$cell.id)>0 & xx < (poss.cp.ua)) {
paquet <- ifelse(poss.cp.ua-xx >= 0 ,min(x,poss.cp.ua-xx) ,0 )
pc.cells.ua.x <- sample( cell.dispo.ua$cell.id, paquet)
pc.cells.ua <- c(pc.cells.ua,pc.cells.ua.x)
cell.dispo.ua <- cell.dispo.ua[-which(cell.dispo.ua$cell.id%in%pc.cells.ua.x),]
#x <- sum(cell.dispo.ua$vol,na.rm=T ) # volume restant
xx <- length(c(pc.cells.ua )) # aire récoltée
}
## Add the cells partially cut in this UA:
pc.cells <- c(pc.cells, pc.cells.ua)
}
################################################# TRACKING #################################################
## Area salvaged and non-salvaged, clearcut
## areas are in cells, volumes are in m3
a.salv <- filter(land2, cell.id %in% cc.cells.salv.tot) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
a.unaff <- filter(land2, cell.id %in% cc.cells.unaff.tot) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
### Volume logged, clearcut
v.salv <- filter(land.ea, cell.id %in% cc.cells.salv.tot) %>% group_by(mgmt.unit) %>% summarise(x=sum(vol))
v.unaff <- filter(land.ea, cell.id %in% cc.cells.unaff.tot) %>% group_by(mgmt.unit) %>% summarise(x=sum(vol))
## Area clearcut per species per management unit (clearcut)
spp.ccut <- filter(land2, cell.id %in% c(cc.cells.salv.tot, cc.cells.unaff.tot)) %>%
group_by(mgmt.unit, spp) %>% summarise(x=length(mgmt.unit))
## Volume clearcut per species per management unit (clearcut)
spp.ccut.vol <- filter(land2, cell.id %in% c(cc.cells.salv.tot, cc.cells.unaff.tot)) %>%
group_by(mgmt.unit, spp) %>% summarise(x=sum(vol))
## Area partial cut per management unit
a.pcut <- filter(land2, cell.id %in% pc.cells) %>% group_by(mgmt.unit) %>% summarise(x=length(mgmt.unit))
## Volume partial cut per management unit
v.pcut <- filter(land2, cell.id %in% pc.cells) %>% group_by(mgmt.unit) %>% summarise(x=sum(vol))
## Area partialcut per species per management unit (clearcut)
spp.pcut <- filter(land2, cell.id %in% pc.cells) %>% group_by(mgmt.unit, spp) %>% summarise(x=length(mgmt.unit))
## Volume partialcut per species per management unit (clearcut)
spp.pcut.vol <- filter(land2, cell.id %in% pc.cells) %>% group_by(mgmt.unit, spp) %>% summarise(x=sum(vol))
## Merge all the info, FMU level
## @@@@@@@@ MATHIEU, here the names of the variables area "area" but they are in fact NUMBER OF CELLS!
## I don't change anything, but it is confusing ;-)
track <- left_join(s.inc, s.ea, by="mgmt.unit") %>% left_join(s.mat.pc, by="mgmt.unit") %>%
left_join(s.inc.burnt, by="mgmt.unit") %>% left_join(s.inc.mat.burnt, by="mgmt.unit") %>%
left_join(s.inc.kill, by="mgmt.unit") %>% left_join(s.inc.mat.kill, by="mgmt.unit") %>%
left_join(reg.fail.ex, by="mgmt.unit") %>% left_join(reg.fail.inc, by="mgmt.unit") %>%
left_join(a.salv, by="mgmt.unit") %>% left_join(a.unaff, by="mgmt.unit") %>%
left_join(v.salv, by="mgmt.unit") %>% left_join(v.unaff, by="mgmt.unit") %>%
left_join(a.pcut, by="mgmt.unit") %>% left_join(v.pcut, by="mgmt.unit")
names(track)[2:ncol(track)] <- c("a.inc", "a.even.age", "a.mat.pc", "a.inc.burnt", "a.inc.mat.burnt",
"a.inc.kill", "a.inc.mat.kill", "a.reg.fail.ex", "a.reg.fail.in", "a.salvaged", "a.unaff","v.salv",
"v.unaff","a.pcut","v.pcut")
track[is.na(track)] <- 0
track[,c(2:12,15)] = track[,c(2:12,15)]*km2.pixel
#### merge, species level
spp.track <- left_join(spp.ccut, spp.ccut.vol, by=c("mgmt.unit", "spp")) %>%
left_join(spp.pcut, by=c("mgmt.unit", "spp")) %>% left_join(spp.pcut.vol, by=c("mgmt.unit", "spp"))
names(spp.track)[3:ncol(spp.track)] <- c("area.ccut","vol.ccut", "area.pcut","vol.pcut")
spp.track[is.na(spp.track)] <- 0
spp.track[,c(3,5)] = spp.track[,c(3,5)]*km2.pixel
## Return the cell.id of the cut locations and the tracking info
return(list(cc.cells=cc.cells, pc.cells=pc.cells, track.cut=track, spp.track=spp.track))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.