R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

setwd("~/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/R_files/")

lib_loc <- "/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library"
# Folders location
CODE_DIR <- "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/R_files/"
DATA_DIR <- "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Input_datasets/"
H100_DIR <- "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Site_indices_H100/"
RDI_DIR  <- "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Relative_density_index/"
PARAM_DIR<- "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Parameters/"

#Outputs
OUT_MOD_DIR<-"/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Outputs/"


#################################################################
#                                                                                       #
#                       Initialisation des données IFN                                  #
#                       Calcul de Dg, age_1, Ho, add regions                            #
#                                                                                       #
#                       @ Aude Valade, adapted from Patrick Vallet                      #
#                                                       first written March 2014        #
#                                                                                       #
#################################################################

source(paste(CODE_DIR,'F-function_init_IGN.R',sep=""))
#rgdal, nlme
lapply(c('gtools','sp','rgdal','zoo','lmtest','micEcon','frontier','nlme'), require,character.only=T,lib.loc=lib_loc)

Loading required package: gtools
Loading required package: sp
Loading required package: rgdal
Please note that rgdal will be retired during October 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-7, (SVN revision 1203)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.5.3, released 2022/10/21
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/rgdal/gdal
 GDAL does not use iconv for recoding strings.
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 9.1.0, September 1st, 2022, [PJ_VERSION: 910]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-1
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: lmtest
Loading required package: micEcon

If you have questions, suggestions, or comments regarding one of the 'micEcon' packages, please use a forum or 'tracker' at micEcon's R-Forge site:
https://r-forge.r-project.org/projects/micecon/
Loading required package: frontier

Please cite the 'frontier' package as:
Tim Coelli and Arne Henningsen (2013). frontier: Stochastic Frontier Analysis. R package version 1.1. http://CRAN.R-Project.org/package=frontier.

If you have questions, suggestions, or comments regarding the 'frontier' package, please use a forum or 'tracker' at frontier's R-Forge site:
https://r-forge.r-project.org/projects/frontier/
Loading required package: nlme

[[1]]
[1] TRUE

[[2]]
[1] TRUE

[[3]]
[1] TRUE

[[4]]
[1] TRUE

[[5]]
[1] TRUE

[[6]]
[1] TRUE

[[7]]
[1] TRUE

[[8]]
[1] TRUE

flag_get_H100   <-'calculate' # read or calculate #Read H100 values or recalculate them if updates in the script
flag_get_RDI    <-'calculate' # read or calculate #Read H100 values or recalculate them if updates in the script
#-------------------------------------------------------------------------------------------------
# Initialize definition dataframes and tools
#-------------------------------------------------------------------------------------------------
species         <- init_species()
species_all     <- init_species_all()
GRECO           <- read_GRECO()
serreg          <- read_ser(DATA_DIR)

Warning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among others

OGR data source with driver: ESRI Shapefile 
Source: "/Users/valade/EcoSols_Nextcloud/Foret_francaise/BiCaFF/code_modele_share_VB/Input_datasets/1-Donnees_IGN/fiches_GRECO/ser_l93", layer: "ser_l93"
with 473 features
It has 2 fields

ser0            <- data.frame(codes=unique(serreg$codeser),names=remove_accents(unique(serreg$NomSER)))
ser                     <- ser0[!is.na(ser0$names),]

#-------------------------------------------------------------------------------------------------
#       Read all years of data for trees and stands
#--------------------------------------------------------------------------------------------------
# xx VBL'intérêt du code est de pouvoir mettre à jour avec de nouvelles données IGN, mais là ça semble pas évident puisque tu utilises pas les fichiers bruts. Du coup, il faudrait mettre des commentaires pour expliquer comment tu as généré ces fichiers modifiés et/ou inclure les scripts ici. On pourra aussi simplifier en générant un fichier .RData à partir de ce code (moins transparent). Qu'en penses-tu ?
suffx                           <-'0812'
arbre                           <- read.table(paste(DATA_DIR,"1-Donnees_IGN/Arbres_foret_simplifReformateVolumeCorrige_0812.txt",sep=""), h=T, sep=";", colClasses=c(rep("factor", 2),"integer",rep("factor", 14), rep("numeric", 4), "factor", rep("numeric", 8)),na.strings=c("",".","NA"))
arbre                           <- arbre[order(arbre$idp, arbre$a),]
listidp_arbr            <- unique(arbre$idp)
arbres_morts_foret      <- read.table(paste(DATA_DIR,'1-Donnees_IGN/arbres_morts_foret_0812.txt',sep=''),h=T,sep=';')
placette2012_2          <- read.table(paste(DATA_DIR,"1-Donnees_IGN/placettes_foret_2012-2.csv",sep=""),h=T, sep=";", colClasses=c( "factor","numeric","numeric",rep("factor", 26)),na.strings=c("",".","NA"))
placette                        <- read.table(paste(DATA_DIR,"1-Donnees_IGN/Placettes_foret_0812.txt",sep=""), h=T, sep=";", colClasses=c(rep("factor", 2), rep("numeric", 2), rep("factor", 6), "numeric", rep("factor", 16),'character',rep("factor", 4)),na.strings=c("",".","NA"))
placette                        <- placette[order(placette$idp),]
placette                        <- merge(placette, placette2012_2[,c('idp','ess_age_1')],by='idp',all.x=T)
placette                        <- placette[placette$idp %in% listidp_arbr,]

#-------------------------------------------------------------------------------------------------
#       Adjust species description
#       (pool 'residual species together', select species from (esspre, ess_age_1), write in numerical format)
#--------------------------------------------------------------------------------------------------
placette$ess_age_1      <- ifelse(is.na(placette$ess_age_1.y)&!is.na(placette$ess_age_1.x),as.vector(as.factor(placette$ess_age_1.x)),as.vector(as.factor(placette$ess_age_1.y)))
placette$type           <- ifelse(as.numeric(gsub("[^0-9]","",placette$ess_age_1))>50,type<-'R',type<-'F')
placette[is.na(placette$ess_age_1),]$type<- ifelse(as.numeric(gsub("[^0-9]","",placette[is.na(placette$ess_age_1),]$esspre))>50,type<-'R',type<-'F')
placette$esspre_R       <- mapply(findresiduals,placette$esspre,placette$type)  #no esspre var in arbre db but been merged with placette db
placette$ess_age_R      <- mapply(findresiduals,placette$ess_age_1,placette$type)
placette$ess_age_R      <- ifelse(is.na(placette$ess_age_R),as.vector(as.factor(placette$esspre_R)),as.vector(as.factor(placette$ess_age_R)))
placette$es_num         <- as.vector(as.factor(as.numeric(gsub("[^0-9]","",placette$ess_age_R))))
placette                        <- placette[!is.na(placette$ess_age_R),]
#---- Small fixes
placette[is.na(placette$gest),]$gest    <-'999'
names(placette)[names(placette)=="tm2"] <-"tm_2"

#-------------------------------------------------------------------------------------------------
# Add geographical references (SER & GRECO)
#-------------------------------------------------------------------------------------------------
list_idp_ser            <- read.csv(paste(DATA_DIR,'1-Donnees_IGN/list_idp.csv',sep=''))
placette                        <- merge(placette,list_idp_ser,by="idp",all.x=T)
placette$GRECO          <- substr(placette$SER,1,1)


#--------------------------------------------------------------------------------------------------
# Make some adjustments and additions to the tree data: compute basal area (g), relative height expresse as a proportion of the height of the highest tree of the plot (hrel), dominant vs understorey (strate), xx? (gDeb), coniferous vs broadleaf (type), age only for trees of the dominant species (age_pur)
#--------------------------------------------------------------------------------------------------
# Basal area (G) of each tree in  m^2/ha  : G = pi*r^2 = pi*(c/2)^2 = c^2/(4*pi)
# Divide by 10000 to convert cm^2 to m^2 and multiply by w to have representativity on one hectare
arbre                           <- merge(arbre,list_idp_ser                             ,by="idp",all.x=TRUE)
arbre                           <- merge(arbre,placette [,c('idp','esspre_R','ess_age_R','cac','peupnr')],by='idp',all.x=T)

hmax                            <- tapply(arbre$htot, arbre$idp, max) # Identify dominant and understory from height relative to maximum height
arbre$g                         <- arbre$c13*arbre$c13*arbre$w/(40000*pi)
arbre$d13                       <- arbre$c13/pi
arbre$hrel                      <- arbre$htot/hmax[arbre$idp]
arbre$strate            <- "sousEtage" ; arbre$strate[arbre$hrel>=0.50] <- "etageDominant"
arbre$gDeb                      <- pi/10000*(arbre$c13/(2*pi)-arbre$ir5/10)^2*arbre$w
arbre$GRECO                     <- substr(arbre$SER,1,1)
arbre$type                      <- ifelse(as.numeric(gsub("[^0-9]","",arbre$espar))>50,type<-'R',type<-'F')
arbre$esparR            <- mapply(findresiduals,arbre$espar,arbre$type) #add flag 'compl' if take letters in the species code
arbre$esparR            <- as.numeric(gsub("[^0-9]","",arbre$esparR))
arbre$ess_age_R         <- as.numeric(gsub("[^0-9]","",arbre$ess_age_R))
arbre$esspre_R          <- as.numeric(gsub("[^0-9]","",arbre$esspre_R))
arbre$main_species      <- ifelse(as.numeric(gsub("[^0-9]","",arbre$ess_age_R))==as.numeric(gsub("[^0-9]","",arbre$esparR)),'Y','N')
arbre$age_pur           <- ifelse(arbre$main_species =='N',NA,arbre$age)
#Check for presence of trees of main species xx VB Delete what follows?

# ##### Problem some plots have no tree of species ess_age or esspre
# list_count                    <- data.frame(
# idp                   =unique(arbre$idp),
# count_main    =as.vector(tapply(   (arbre$main_species =='Y'), arbre$idp, sum)),
# count_main2   =as.vector(tapply(   (arbre$esparR ==arbre$ess_age_R), arbre$idp, sum)),
# count_all=as.vector(tapply(   arbre$main_species, arbre$idp, length)))
# arbre                         <- merge(arbre,list_count,by='idp',all.x=T,all.y=T)
# arbre$ess_age_R               <- ifelse(arbre$count_main == 0 & arbre$esspre_R!=arbre$ess_age_R,arbre$esspre_R,arbre$ess_age_R)
# arbre$main_species    <- ifelse(as.numeric(gsub("[^0-9]","",arbre$ess_age_R))==as.numeric(gsub("[^0-9]","",arbre$esparR)),'Y','N')
# list_count2                   <- data.frame(idp=unique(arbre$idp),count_main2=as.vector(tapply(   (arbre$main_species =='Y'), arbre$idp, sum)))
# arbre<                                - merge(arbre,list_count2,by='idp',all.x=T,all.y=T)
# arbre$main_count      <-arbre$main_count2
# arbre[arbre$idp %in% list_count2[list_count2$count_main==0,'idp'],c('cac')]
#-------------------------------------------------------------------------------------------------
# Read and add to dataset pre-calculated density index, tree volume increment, volume expansion factors
#-------------------------------------------------------------------------------------------------
list_iv5_arbre          <- read.csv(paste(DATA_DIR,'1-Donnees_IGN/list_iv5_arbres.H100.4.csv',sep=''))
list_VEF                        <- read.csv(paste(DATA_DIR,'1-Donnees_IGN/VEF_list.csv',sep=''))
arbre                   <- merge(arbre,list_iv5_arbre[c('idp','a','iv5')],all.x=T,by=c('idp','a'))



#h0<-aggregate(htot~idp+ess_age_R,data=arbre[!is.na(arbre$age),],mean)
#h0<-h0[order(h0$idp,h0$ess_age_R),]
#write.table(h0,'list_h0.2.txt',sep=';')

#--------------------------------------------------------------------------------------------------
# Make some adjustments and additions to the plot-level data: average plot age (age_1), average age of the dominant species (age_pur), basal area of dominant stratum and of the whole plot (Gdom and G), plot volume (V), real (* w) numver of trees and trees of the dominant stratum (N and Ndom), quadratic mean diameter of plot and dominant stratum (Dg and Dgdom), average plot diameter (d13), plot basal area five years ago (GDeb) xx Aude vérifie, increase in basal area over last five years (DeltaG5)
#--------------------------------------------------------------------------------------------------
placetteDendro                  <- data.frame(idp=unique(arbre$idp))
placetteDendro$age_1    <- as.vector(tapply(arbre$age , arbre$idp, mean,na.rm=TRUE)) #no age data at tree level before 2008
placetteDendro$age_pur  <- as.vector(tapply(arbre$age_pur,                           arbre$idp, mean,na.rm=TRUE)) ; placetteDendro$age_1[placetteDendro$age_pur==0]<-NA
placetteDendro$Gdom             <- as.vector(tapply(arbre$g*(arbre$strate=="etageDominant"), arbre$idp, sum)) # 50470 values
placetteDendro$V                <- as.vector(tapply(arbre$v*arbre$w,                         arbre$idp, sum))   ### ok simplif trees (corrected file P. Vallet)+ ok all species
placetteDendro$G                <- as.vector(tapply(arbre$g,                                 arbre$idp, sum))
placetteDendro$Ndom             <- as.vector(tapply(arbre$w*(arbre$strate=="etageDominant"), arbre$idp, sum))
placetteDendro$N                <- as.vector(tapply(arbre$w,                                 arbre$idp, sum))
placetteDendro$Dgdom    <- 100*sqrt(4*placetteDendro$Gdom/(placetteDendro$Ndom*pi))
placetteDendro$Dg               <- 100*sqrt(4*placetteDendro$G/(placetteDendro$N*pi))
placetteDendro$d13              <- as.vector(tapply(arbre$d13*arbre$w,                       arbre$idp, sum))/placetteDendro$N   ### ok simplif trees (corrected file P. Vallet)+ ok all species
placetteDendro$GDeb             <- as.vector(tapply(arbre$gDeb,                              arbre$idp, sum))                            ### CHECK Need to incorporate bark ratio!!
placetteDendro$DeltaG5  <- placetteDendro$G - placetteDendro$GDeb

# ------ conversion of stands into homogeneous stand xxVB: replace with following line? you are not converting stand or computing stand level characteristics here ...
# ------ compute average (per tree) characteristics of dominant species
# sum the cumulative variables of the stand for the dominant species
placetteDendro$G_ess    <- as.vector(tapply(arbre$g*         (arbre$main_species=='Y'), arbre$idp, sum))
placetteDendro$V_ess    <- as.vector(tapply(arbre$v* arbre$w* (arbre$main_species=='Y'), arbre$idp, sum))
placetteDendro$N_ess    <- as.vector(tapply(arbre$w*          (arbre$main_species=='Y'), arbre$idp, sum))

## PROBLEM when no tree of main species,  G_ess and N_ess=0 => Dg_ess=NA xxVB: you don't solve it so I guess you consider it's not a problem. Can you explain and rephrase?
# average the tree diameter and diameter increment measurements over the main species of each plot
placetteDendro$d13_ess  <- as.vector(tapply(arbre$d13*arbre$w*(arbre$main_species=='Y'), arbre$idp, sum))/placetteDendro$N_ess   ### ok simplif trees (corrected file P. Vallet)+ ok all species xxVB?? explanation needed or delete comment
placetteDendro$ir5_ess  <- as.vector(tapply(arbre$ir5*arbre$w*(arbre$main_species=='Y'), arbre$idp, sum))/placetteDendro$N_ess   ### ok simplif trees (corrected file P. xxVB?? explanation needed or delete comment

placetteDendro$Dg_ess   <- 100*sqrt(4*placetteDendro$G_ess/(placetteDendro$N_ess*pi))
placetteDendro$v_ind_ess<-placetteDendro$V_ess/placetteDendro$N_ess
placetteDendro$N_homog  <-placetteDendro$N_ess+(placetteDendro$V-placetteDendro$V_ess)/placetteDendro$v_ind_ess
#placetteDendro$mdr<- placetteDendro$d13_ess/placetteDendro$Dg # attempt to give a conversion factor between real-world diameters used for scenario description and modelled Dg diameters over all stand even when mixed
placetteDendro$IV5              <-as.vector(tapply(arbre$iv5*arbre$w, arbre$idp, sum,na.rm=T)) # NA for iv5 for trees for which NA for ir5. Put these ones in NR xxVB 1) should this be moved up? (not related to dominant species) 2) did you solve the NA problem? If so, explain how.
#------------------------------------------------------------------------------------------------
# Mortality data
#-------------------------------------------------------------------------------------------------
placetteMort                    <-data.frame(idp=unique(arbres_morts_foret$idp))
placetteMort$Vmort5             <-as.vector(tapply(arbres_morts_foret$v*(arbres_morts_foret$datemort==1)*arbres_morts_foret$w,arbres_morts_foret$idp,sum))
placette                                <-merge(placette, placetteMort, by='idp',all.x=T)


#-------------------------------------------------------------------------------------------------
# Merge everything together
#--------------------------------------------------------------------------------------------------
placette                                <-merge(placette, placetteDendro,by='idp')



# !!!! GO UNTIL HERE  AND REDO CALCULATION OF H100 and RDI with special attention to non recensables xxVB Unclear. Rephrase

#-------------------------------------------------------------------------------------------------
# Add RDI values : read in stored file or recalculate or recalculate
#-------------------------------------------------------------------------------------------------if(flag_get_RDI=='read'){
if(flag_get_RDI=='read'){
  list_RDI        <- read.csv(paste(RDI_DIR,'list_RDI_N_homog_Dg_ess_noshift.csv',sep=''))#shift suffx means that a parallel to the sfa curve has been calculated such that all points are below the self-thinning line
}else{
  source(paste(CODE_DIR,'get_density_indexes.2.R',sep=''))
}

[1] "============"
[1] " Density indices at scale national for dataset 0812"
[1] "============"
[1] "02  :  chêne pédonculé"
[1] "dim(bd_calc) 3968"
[1] "dim(bd_apply) 4993"
(Intercept) 
   12.98275 
log(Dg_ess) 
  -1.967139 
[1] "03  :  chêne sessile"
[1] "dim(bd_calc) 2983"
[1] "dim(bd_apply) 3747"
(Intercept) 
   13.29625 
log(Dg_ess) 
  -2.069658 
[1] "05  :  chêne pubescent"
[1] "dim(bd_calc) 1861"
[1] "dim(bd_apply) 2349"
(Intercept) 
   11.51596 
log(Dg_ess) 
  -1.456404 
[1] "06  :  chêne vert"
[1] "dim(bd_calc) 597"
[1] "dim(bd_apply) 757"
(Intercept) 
   9.198496 
log(Dg_ess) 
 -0.4613006 
[1] "09  :  hêtre"
[1] "dim(bd_calc) 2368"
[1] "dim(bd_apply) 2976"
(Intercept) 
   13.94508 
log(Dg_ess) 
    -2.1993 
[1] "10  :  châtaignier"
[1] "dim(bd_calc) 957"
[1] "dim(bd_apply) 1201"
(Intercept) 
   12.27637 
log(Dg_ess) 
  -1.625782 
[1] "11  :  charme"
[1] "dim(bd_calc) 511"
[1] "dim(bd_apply) 660"
(Intercept) 
   13.72552 
log(Dg_ess) 
  -2.199929 
[1] "17C  :  frêne commun"
[1] "dim(bd_calc) 1112"
[1] "dim(bd_apply) 1404"
(Intercept) 
   12.49372 
log(Dg_ess) 
  -1.801063 
[1] "98  :  residuels feuillus"
[1] "dim(bd_calc) 2919"
[1] "dim(bd_apply) 3681"
(Intercept) 
   11.02437 
log(Dg_ess) 
  -1.256723 
[1] "51  :  pin maritime"
[1] "dim(bd_calc) 1239"
[1] "dim(bd_apply) 1556"
(Intercept) 
    12.1024 
log(Dg_ess) 
  -1.682464 
[1] "52  :  pin sylvestre"
[1] "dim(bd_calc) 1663"
[1] "dim(bd_apply) 2090"
(Intercept) 
   12.48559 
log(Dg_ess) 
  -1.749981 
[1] "53CO  :  pin laricio"
[1] "dim(bd_calc) 309"
[1] "dim(bd_apply) 391"
(Intercept) 
   11.52072 
log(Dg_ess) 
  -1.443934 
[1] "54  :  pin noir d'Autriche"
[1] "dim(bd_calc) 288"
[1] "dim(bd_apply) 362"
(Intercept) 
   12.45827 
log(Dg_ess) 
  -1.743694 
[1] "57A  :  pin d'Alep"
[1] "dim(bd_calc) 336"
[1] "dim(bd_apply) 422"
(Intercept) 
   12.45088 
log(Dg_ess) 
  -1.832441 
[1] "61  :  sapin pectiné"
[1] "dim(bd_calc) 1051"
[1] "dim(bd_apply) 1320"
(Intercept) 
   12.98737 
log(Dg_ess) 
  -1.852882 
[1] "62  :  epicéa commun"
[1] "dim(bd_calc) 1095"
[1] "dim(bd_apply) 1378"
(Intercept) 
   13.22562 
log(Dg_ess) 
    -1.9095 
[1] "64  :  douglas"
[1] "dim(bd_calc) 717"
[1] "dim(bd_apply) 902"
(Intercept) 
   11.41146 
log(Dg_ess) 
  -1.421504 
[1] "99  :  residuels resineux"
[1] "dim(bd_calc) 599"
[1] "dim(bd_apply) 750"
(Intercept) 
   11.67514 
log(Dg_ess) 
  -1.482562 

placette        <- merge(placette,list_RDI[,c('idp','RDI')], all.x=T,by='idp')
avgRDIes        <-(tapply(placette$RDI,placette$es_num,mean,na.rm=T))
dfavgRDIes      <- data.frame(es_num=rownames(avgRDIes),avgRDI= as.vector(avgRDIes))
placette        <- merge(placette, dfavgRDIes,by='es_num')


#-------------------------------------------------------------------------------------------------
# Add H100 values : read in stored file or recalculate or recalculate
#-------------------------------------------------------------------------------------------------
if(flag_get_H100=='read'){
  list_H100       <- read.csv(paste(H100_DIR,'list_H100_age_1_sfo_diff.2.csv',sep=''), h=T) # list_H100_pur is when curves are fitted on (age_pur,h0)
}else{
  source(paste(CODE_DIR,'get_site_indexes.4.R',sep=''))
}

[1] "2"
[1] "2436"
[1] "try fitting iteration : 0"
[1] "try fitting iteration : 1"
[1] "try fitting iteration : 2"
[1] "2444"
[1] "try fitting iteration : 0"
[1] "try fitting iteration : 1"
[1] "try fitting iteration : 2"
[1] "try fitting iteration : 3"
[1] "3"
[1] "2064"
[1] "try fitting iteration : 0"
[1] "try fitting iteration : 1"
[1] "try fitting iteration : 2"
[1] "try fitting iteration : 3"
[1] "1633"
[1] "try fitting iteration : 0"
[1] "5"
[1] "748"
[1] "try fitting iteration : 0"
[1] "903"
[1] "try fitting iteration : 0"
[1] "652"
[1] "try fitting iteration : 0"
[1] "6"
[1] "741"
[1] "try fitting iteration : 0"
[1] "9"
[1] "2919"
[1] "try fitting iteration : 0"
[1] "10"
[1] "810"
[1] "try fitting iteration : 0"
[1] "355"
[1] "try fitting iteration : 0"
[1] "11"
[1] "277"
[1] "try fitting iteration : 0"
[1] "345"
[1] "try fitting iteration : 0"
[1] "17"
[1] "705"
[1] "try fitting iteration : 0"
[1] "658"
[1] "try fitting iteration : 0"
[1] "98"
[1] "1419"
[1] "try fitting iteration : 0"
[1] "2118"
[1] "try fitting iteration : 0"
[1] "51"
[1] "1523"
[1] "try fitting iteration : 0"
[1] "52"
[1] "2068"
[1] "try fitting iteration : 0"
[1] "53"
[1] "384"
[1] "try fitting iteration : 0"
[1] "54"
[1] "358"
[1] "try fitting iteration : 0"
[1] "57"
[1] "252"
[1] "try fitting iteration : 0"
[1] "165"
[1] "try fitting iteration : 0"
[1] "61"
[1] "1302"
[1] "try fitting iteration : 0"
[1] "62"
[1] "1359"
[1] "try fitting iteration : 0"
[1] "try fitting iteration : 1"
[1] "64"
[1] "891"
[1] "try fitting iteration : 0"
[1] "99"
[1] "740"
[1] "try fitting iteration : 0"

#xxVB J'ai commenté la ligne suivante car placette a déjà une colonne H100 (cf section "get site indices" plus haut). On peut le mettre dans le "if" de dessus, mais avant il faut vérifier car la colonne placette$H100 n'est pas égale à list_H100$H100. Tu vois et tu fais ce qu'il faut ?
## placette        <- merge(placette,list_H100[,c('idp','H100')],all.x=T,by='idp') # all.x because many plots missing in list_H100 due to lack of age_1 variables for before 2008
avgH100es       <-(tapply(placette$H100,placette$es_num,mean,na.rm=T))
dfavgH100es     <- data.frame(es_num=rownames(avgH100es),avgH100= as.vector(avgH100es))
placette        <- merge(placette, dfavgH100es,by='es_num')


#-------------------------------------------------------------------------------------------------
# Plots which are too young to be inventoried: attribute to them the species average site index and RDI
#--------------------------------------------------------------------------------------------------
NR_idp  <-      placette[placette$cac %in% c('NR','AA')& !is.na(placette$peupnr)&as.numeric(placette$peupnr)>0,'idp']

placette[is.na(placette$RDI) & placette$idp %in% NR_idp, c('RDI')]      <-placette[is.na(placette$RDI) & placette$idp %in% NR_idp,c('avgRDI')]
placette[is.na(placette$H100)& placette$idp %in% NR_idp,c('H100')]      <-placette[is.na(placette$H100)& placette$idp %in% NR_idp,c('avgH100')]


#-------------------------------------------------------------------------------------------------
# Add fertility class from quantile distribution per species
#-------------------------------------------------------------------------------------------------

placette$cac_H100<-NA
for(es in species_all$codes){
  data <- placette[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es)) ,]
  
  if(length(data$idp)>=400){
    H100_1<-quantile(data[!is.na(data$gest) & data$gest>1 & !is.na(data$H100),'H100'],0.25)
    H100_2<-quantile(data[!is.na(data$gest) & data$gest>1 & !is.na(data$H100),'H100'],0.5)
    H100_3<-quantile(data[!is.na(data$gest) & data$gest>1 & !is.na(data$H100),'H100'],0.75)
    
    placette$cac_H100[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es))& !is.na(placette$H100) & placette$H100>H100_2 & placette$H100<=H100_3]   <-3
    placette$cac_H100[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es))& !is.na(placette$H100) & placette$H100>H100_3 ]                                          <-4
  }else{
    H100_1<-quantile(data[!is.na(data$gest) & data$gest>1& !is.na(data$H100),'H100'],0.3)
    H100_2<-quantile(data[!is.na(data$gest) & data$gest>1& !is.na(data$H100),'H100'],0.6)
    placette$cac_H100[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es))& !is.na(placette$H100) & placette$H100>H100_2 ]                                          <-3
    
  }
  
  placette$cac_H100[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es)) & placette$H100<=H100_1]                                                 <-1
  placette$cac_H100[as.numeric(gsub("[^0-9]", "",placette$ess_age_R))==as.numeric(gsub("[^0-9]", "",es)) & placette$H100> H100_1 & placette$H100<=H100_2] <-2
}

placette$cac_age<-0
placette[placette$cac=="01",]$cac_age<-2.5
placette[placette$cac=="02",]$cac_age<-7.5
placette[placette$cac=="03",]$cac_age<-12.5
placette[placette$cac=="04",]$cac_age<-17.5
placette[placette$cac=="05",]$cac_age<-22.5
placette[placette$cac=="06",]$cac_age<-27.5
placette[placette$cac=="07",]$cac_age<-32.5
placette[placette$cac=="08",]$cac_age<-37.5
placette[placette$cac=="09",]$cac_age<-45.0
placette[placette$cac=="10",]$cac_age<-55.0
placette[placette$cac=="11",]$cac_age<-65.0
placette[placette$cac=="12",]$cac_age<-75.0
placette[placette$cac=="13",]$cac_age<-90.0
placette[placette$cac=="14",]$cac_age<- 110.0
placette[placette$cac=="15",]$cac_age<-130.0
placette[placette$cac=="16",]$cac_age<-150.0
placette[placette$cac=="17",]$cac_age<-170.0
placette[placette$cac=="18",]$cac_age<-190.0
placette[placette$cac=="19",]$cac_age<-220.0
placette[placette$cac=="20",]$cac_age<-250.0

placette_reset<-placette
## xx2 following line added to avoid re-running the script
save(file=paste(CODE_DIR,"inventory_data.Rdata",sep=""),list=c("placette_reset","species_all","GRECO"))
#-------------------------------------------------------------------------------------------------
#xxVB: I'm not sure what that is ...
for(i in 1:5){
  alarm()
  Sys.sleep(0.5)
}



ggplot(test,aes(x=age_1,y=Ndom,col=esspre))+
  geom_point()+
  theme_bw()
ggsave("age_1_N_NR.png")

Saving 7.29 x 4.51 in image

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.