Initiation function for even harvest #13

Added initiation distrib from issue #13 on github

Also edited the tests because new def_init is causing test fails

Author: gowachin

DESCRIPTION

@@ -45,4 +45,4 @@ LazyData: true
 LazyDataCompression: bzip2
 URL: https://github.com/gowachin/matreex
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1

NAMESPACE

@@ -42,6 +42,7 @@ export(def_init)
 export(def_initBA)
 export(def_init_even)
 export(def_init_k)
+export(def_init_planting)
 export(delay)
 export(disturb_fun)
 export(disturb_fun_mixt)

R/class_species.R

@@ -428,6 +428,37 @@ def_init_k <- function(x){
     return(fun)
 }
 
+#' Init population from plantation, after even clear cut
+#'
+#' @param species name of the species. This allows to extract the parameters
+#' from matreex::distrib_planting dataset, that have been computed by Georges in
+#' issue #13.
+#'
+#' @details
+#' This function may replace def_init_even() in the long term.
+#'
+#'
+#' @export
+def_init_planting <- function(species){
+    # dev
+    # species <- "Picea_abies"
+    distrib <- matreex::distrib_planting
+    pms <- distrib[distrib$species == species,]
+    if(nrow(pms) != 1){
+        stop("This species does not retrieve any row from matreex::distrib_planting dataset. Check the format of the species column.")
+    }
+    force(params)
+    fun <- function(mesh, SurfEch = 0.03) {
+
+        x <- dnorm(log(mesh), mean = pms$mean_DBH_log, sd = pms$sd_DBH_log)/
+            (sum(dnorm(log(mesh), mean = pms$mean_DBH_log, sd = pms$sd_DBH_log)))
+        x[x < 1e-11] <- 0 # this teshold is the same as the IPM minimal non null values
+
+        return(x * SurfEch)
+    }
+
+    return(fun)
+}
 
 #' Default population harvest
 #'

R/distrib_planting.R

@@ -0,0 +1,7 @@
+#' Species distribution parameter for young stands.
+#'
+#' @name distrib_planting
+#'
+#' @source Georges Kunstler
+#'
+"distrib_planting"

R/matreex.R

@@ -7,4 +7,4 @@
 #'
 #' @name matreex
 #' @docType package
-NULL
+"_PACKAGE"

data/distrib_planting.rda

@@ -0,0 +1,46 @@
+"","species","mean_DBH_log","sd_DBH_log","Nha","BAha"
+"1","Abies alba",4.70218549577077,0.0495809812130126,517.661497058824,5.41172158515756
+"2","Acer campestre",4.69021478810264,0.0229491520261313,206.183766666667,1.9017154702442
+"3","Acer pseudoplatanus",4.7110533387838,0.037044737558841,442.281173333333,4.481553224364
+"4","Alnus glutinosa",4.72228044928083,0.0627169270786061,369.777573707849,4.03433579604816
+"5","Alnus sp",4.71078238563064,0.0339154927124127,842.02096,8.79999945002538
+"6","Betula sp",4.71352082053515,0.0463938807874255,282.29810822709,2.93144266175903
+"7","Carpinus betulus",4.70842261081656,0.0581175873300147,331.342175357143,3.53293192490761
+"8","Castanea sativa",4.71322723646938,0.0576733957110027,625.385746924825,6.5825770483407
+"9","Eucalyptus camaldulensis",4.72113227559557,0.0836156442799526,288.196760443228,3.04503871031746
+"10","Eucalyptus globulus",4.73635502053714,0.0985885082207965,447.105139686778,4.91044956388888
+"11","Fagus sylvatica",4.71471148673444,0.0625589907059087,604.240923948785,6.4944952337039
+"12","Fraxinus angustifolia",4.67689121819777,0.197738049873788,163.716623610466,3.51404610058426
+"13","Fraxinus excelsior",4.6956048830764,0.0297308087373162,314.217673577297,3.13943794613797
+"14","Juniperus oxycedrus",4.70322258115659,0.0314976991059382,193.277400741428,1.91874865612117
+"15","Juniperus thurifera",4.71800125274837,0.102156348647792,272.445577438099,2.95473192131547
+"16","Larix decidua",4.70118933007119,0.0663804484729467,400.1715,3.87032614001525
+"17","Larix kaempferi",4.74655021548501,0.0523477139968134,1272.01855,13.2804700521601
+"18","Olea europaea",4.72329762114183,0.0842188126502913,287.657823069796,3.07421186342592
+"19","Picea abies",4.71391698732338,0.0549258868794927,554.30328116423,5.81616466566444
+"20","Picea sitchensis",4.68447641748224,0.0277434397471149,352.89115,3.45317349617457
+"21","Pinus contorta",4.72205093735784,0.075883711426374,392.582192960009,4.15032430555556
+"22","Pinus halepensis",4.73065495555649,0.116309449280215,344.17175202085,3.78384832192317
+"23","Pinus nigra",4.71770137770318,0.0874038286224909,422.0333604511,4.53370277221821
+"24","Pinus pinaster",4.70908191963463,0.0597032021328669,387.877221793085,4.04038045893054
+"25","Pinus pinea",4.74932247954807,0.0687127671922933,324.132515576554,3.51209450569149
+"26","Pinus radiata",4.76502411150768,0.142227689390037,477.464829275685,5.36547499999999
+"27","Pinus sylvestris",4.72355624381171,0.0611146710645294,287.184169902517,3.05368625086921
+"28","Pinus uncinata",4.71558595746972,0.0568791860360685,223.594503705949,2.26539015861161
+"29","Populus nigra",4.69794613568356,0.0271635322521335,284.3002,2.58986559174726
+"30","Populus sp",4.70646107072577,0.0543639300819205,470.006666666667,4.55245740649722
+"31","Populus tremula",4.71721661636695,0.0388965259535102,283.873404787988,2.91486038731533
+"32","Prunus padus",4.65554525240519,0.0020745190845443,113.04,1.00437331580267
+"33","Pseudotsuga menziesii",4.72362157464496,0.0559649421456623,551.975063432745,5.69872716068797
+"34","Quercus faginea",4.71653756548537,0.0731074159796038,385.107342938483,4.10933972785039
+"35","Quercus ilex",4.71797541980459,0.0769973139597457,409.234076051567,4.31240975661893
+"36","Quercus petraea",4.72275057112809,0.0576060013731273,587.732822466617,6.14638891217109
+"37","Quercus pubescens",4.72620358488652,0.0659217683414891,440.365722298375,4.75538259486895
+"38","Quercus pyrenaica",4.72764545638401,0.0967788182623929,471.27761374457,5.12762126577823
+"39","Quercus robur",4.71854679648392,0.0542881150244748,535.488522923456,5.61847598446725
+"40","Quercus rubra",4.7025161285923,0.0391011708167626,543.47855,5.47228416031546
+"41","Quercus suber",4.72519630147892,0.142016971915139,172.97195000328,2.10987500462963
+"42","Robinia pseudacacia",4.71226747547154,0.0470528425569732,424.696153846153,4.4769811823802
+"43","Salix caprea",4.72083877933273,0.0413963352956474,257.442778048919,2.69437893753388
+"44","Salix sp",4.71678579103847,0.0333416179745275,533.486066666667,5.33333257535536
+"45","Tilia cordata",4.62359195460859,0,88.4100000000001,0.720430709203283

dev/ParamSizeDistYoungStand.csv

@@ -0,0 +1,145 @@
+---
+title: "Debug cases from Jasper"
+output:
+    github_document
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+    collapse = TRUE,
+    comment = "#>"
+)
+options(cli.progress_show_after = 600)
+# 10 minutes before showing a cli progress bar
+
+library(tidyr)
+```
+
+# Common
+
+```{r make_IPM_Fagus}
+# library(matreex)
+load_all()
+library(dplyr)
+library(ggplot2)
+
+# Load fitted model for a species
+# fit_species # list of all species in dataset
+data("fit_Fagus_sylvatica")
+
+# Load associated climate
+data("climate_species")
+climate <- subset(climate_species, N == 2 & sp == "Fagus_sylvatica", select = -sp)
+# see ?climate_species to understand the filtering of N.
+climate
+
+Fagus_ipm <- make_IPM(
+    species = "Fagus_sylvatica",
+    climate = climate,
+    fit = fit_Fagus_sylvatica,
+    clim_lab = "optimum clim",
+    mesh = c(m = 700, L = 90, U = get_maxdbh(fit_Fagus_sylvatica) * 1.1),
+    BA = 0:80, # Default values are 0:200, smaller values speed up this vignette.
+    verbose = TRUE
+)
+```
+
+```{r make_IPM_sylvatica}
+data("fit_Abies_alba")
+Abies_ipm <- make_IPM(
+    species = "Abies_alba",
+    climate = climate, # this variable is defined at the top of the doc.
+    fit = fit_Abies_alba,
+    clim_lab = "optimum clim",
+    mesh = c(m = 700, L = 90, U = get_maxdbh(fit_Abies_alba) * 1.1),
+    BA = 0:80, # Default values are 0:200, smaller values speed up this vignette.
+    verbose = TRUE
+)
+```
+
+```{r hide warn, echo=FALSE}
+options(W_matreex_edist = FALSE)
+```
+
+# Uneven Harvest and Disturbance
+
+```{r uneven_forest}
+
+U_Forest <- forest(
+    species = list(
+        Fagus = species(IPM = Fagus_ipm, init_pop = def_initBA(30),
+                        harvest_fun = Uneven_harv, disturb_fun = disturb_fun,
+                        harv_lim = c(dth = 175, dha = 575, hmax = 1), 
+                        disturb_coef = "sp")
+    )
+)
+
+disturb <- data.frame(type = "storm", intensity = 0.8, IsSurv = FALSE, t = 100)
+
+U_Forest$species[[1]]$disturb_coef
+
+Uneven_sim <- sim_deter_forest(
+    U_Forest,
+    tlim = 260, equil_time = 260, equil_dist = 10, equil_diff = 1,
+    harvest = "Uneven", targetBA = 20,
+    SurfEch = 0.03, disturbance  = disturb,
+    verbose = TRUE
+)
+
+Uneven_sim  %>%
+    dplyr::filter(var %in% c("BAsp", "BAstand", "N"), ! equil) %>%
+    ggplot(aes(x = time, y = value)) +
+    facet_wrap(~ var, scales = "free_y") +
+    geom_line(linewidth = .2) + geom_point(size = 0.4)
+
+```
+
+# Harvest and Disturbance competition
+
+
+```{r uneven_forest}
+load_all()
+E_Forest <- forest(
+    species = list(
+        Fagus = species(IPM = Fagus_ipm, 
+                        init_pop = def_init_planting("Fagus_sylvatica"),
+                        harvest_fun = Even_harv, disturb_fun = disturb_fun,
+                        disturb_coef = "sp", rdi_coef = "sp") #,
+        # Abies = species(IPM = Abies_ipm, init_pop = def_init_even,
+                        # harvest_fun = Even_harv, disturb_fun = disturb_fun,
+                        # disturb_coef = "sp", rdi_coef = "sp")
+    ), 
+    harv_rules = c(Pmax = 0.25, dBAmin = 3, freq = 5, alpha = 1)
+)
+
+disturb <- data.frame(type = "storm", intensity = .5, IsSurv = FALSE, t = 150)
+
+E_Forest$species[[1]]$IPM$info["delay"] <- 48
+as.numeric(E_Forest$species[[1]]$IPM$info["delay"])
+
+Even_sim <- sim_deter_forest(
+    E_Forest,
+    tlim = 300, equil_time = 300, equil_dist = 10, equil_diff = 1,
+    harvest = "Even",  targetRDI = 0.6, targetKg = 0.7,
+    final_harv = 120,
+    disturbance  = disturb,
+    SurfEch = 0.03,
+    verbose = TRUE
+)
+
+Even_sim  %>%
+    # dplyr::filter(var %in% c("H", "N", "BAsp"), ! equil, time >90, time < 120) %>%
+    dplyr::filter(var %in% c("BAsp", "H", "N"), ! equil) %>%
+    ggplot(aes(x = time, y = value, color = species)) +
+    facet_wrap(~ var, scales = "free_y") +
+    geom_line(linewidth = .2) + geom_point(size = 0.4)
+
+Even_sim  %>%
+    # dplyr::filter(var %in% c("H", "N", "BAsp"), ! equil, time >90, time < 120) %>%
+    dplyr::filter(var %in% c("n", "h"), time %in% c(118:121, 165:175)) %>%
+    ggplot(aes(x = size, y = value, color = species)) +
+    facet_wrap(time ~ var, scales = "free_y") +
+    geom_line(linewidth = .2) + geom_point(size = 0.4)
+
+```
+

dev/dev_jasper_b.Rmd

@@ -0,0 +1,89 @@
+# Compute mean and sd of species dist for young stand for Jasper
+
+#FUNDIV <- read.csv("data/NewFundiv/FUNDIV_tree.csv")
+FUNDIV <- read.csv("data/NewFundiv/FUNDIV_data.csv", sep = " ")
+
+hist(FUNDIV$dbh2, breaks = 100)
+abline(v = 100, col = "red")
+
+FUNDIV <- FUNDIV[FUNDIV$dbh2 >100, , ]
+FUNDIV <- FUNDIV[FUNDIV$treestatus %in% c(1,2), ]
+ # 1 2 in growth alive
+library(dplyr)
+library(tidyr)
+library(Hmisc)
+
+# Compute dominant species in ba
+fun_dom_sp <- function (sp,ba, Nha){
+  res <- data.frame(sp = sp, ba_ha2 = ba, Nha = Nha) %>% group_by(sp) %>% summarise(BA = sum(ba_ha2), N = sum(Nha)) %>%
+    arrange(desc(BA)) %>%ungroup() %>% mutate(BA_per = BA/sum(BA))
+  return((res[1,]))
+}
+
+df <- FUNDIV %>% group_by(plotcode) %>%
+  summarise(country = unique(country),
+            MQD_log = sum(log(dbh2)*Nha2)/sum(Nha2),
+            MQD_log_sd = sqrt(wtd.var(log(dbh2), Nha2)),
+            domsp = list(fun_dom_sp(sp=species, ba = Nha2 *(dbh2/1000)^2/4*pi, Nha = Nha2)))%>%
+  unnest_wider(domsp)
+
+sp_sel <- names(sort(table(df$sp), decreasing = TRUE)[1:45])
+c("Pinus sylvestris", "Picea abies", "Fagus sylvatica", "Quercus ilex",
+  "Pinus pinaster", "Pinus halepensis", "Quercus robur", "Quercus petraea",
+  "Pinus nigra", "Betula sp")
+res <- df %>% filter(sp %in% sp_sel ,
+                     BA_per >0.9,
+                     MQD_log < log(120))%>%
+  group_by(sp) %>% summarise(MQD_log = mean(MQD_log),
+                             MQD_log_sd = mean(MQD_log_sd),
+                             Nha = mean(N),
+                             BA = mean(BA))
+
+res <- recrut[, -1]
+names(res) <- c("sp", "MQD_log", "MQD_log_sd", "Nha", "BAha")
+
+# Plot of the simulated distribution
+par(mfrow = c(2,5))
+for (i in 1:10){
+  hist(exp(rnorm(10000, mean = res$MQD_log[i], sd = res$MQD_log_sd[i])),
+       main = res$sp[i], breaks = 40, probability= TRUE, ylim = c(0, 1))
+  # size mesh for tree size
+  mesh <- 100:150
+  # compute relative probability density for each mesh based on lognormal proba
+  prob_dens <- dnorm(log(mesh), mean = res$MQD_log[i], sd = res$MQD_log_sd[i])/
+    (sum(dnorm(log(mesh), mean = res$MQD_log[i], sd = res$MQD_log_sd[i])))
+  lines(mesh, prob_dens, col = "blue")
+  abline(v = 99, col = "red")
+  # compute basal area from size dist and mean density and compare to mean basal area in data
+  print(res$sp[i])
+  print(sum(mesh^2/(4*1000*1000) *pi *  prob_dens *res$Nha[i]))
+  print(res$BA[i])
+}
+names(res) <- c("species", "mean_DBH_log", "sd_DBH_log", "Nha", "BAha")
+
+write.csv(res, file = "ParamSizeDistYoungStand.csv")
+#How to introduce this new cohort ? Start after the time lag of the species? As plantation shorter than the full time lag ? By how much ?10 years?
+# A given percentage of teh time lag ? 50 % ?
+
+
+#Not possible to get the structure in the lag Y years before by inverting the IPM because the IPM is non invertible (its determinant is zero)
+
+# writing a function to use these parameters
+matreex::distrib_planting
+load_all()
+data("fit_Picea_abies")
+data("climate_species")
+climate <- subset(climate_species, N == 2 & sp == "Picea_abies", select = -sp)
+
+Picea_ipm <- make_IPM(
+    species = "Picea_abies",
+    climate = climate,
+    fit = fit_Picea_abies,
+    clim_lab = "optimum clim",
+    mesh = c(m = 700, L = 90, U = get_maxdbh(fit_Picea_abies) * 1.1),
+    BA = 0:60, # Default values are 0:200, smaller values speed up this vignette.
+    verbose = TRUE
+)
+
+fun <- def_init_planting("Picea_abies")
+fun(Picea_ipm$mesh, 0.03)