Commit cc6aa5d3 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

strat to add germany

parent 57d53071
......@@ -47,7 +47,7 @@ $(D3)/Done.txt: scripts/process.data/remove.out.R scripts/process.data/summarise
R CMD BATCH scripts/process.data/summarise.data.R
R CMD BATCH scripts/process.data/plot.data.all.R
$(D3)/Done.t.txt: scripts/process.data/merge.all.processed.data.R
$(D3)/Done.t.txt: scripts/process.data/merge.all.processed.data.R $(sites)
R CMD BATCH $<
#-------------------------------------------------------
......
......@@ -175,12 +175,23 @@ load.data.for.lmer <- function(trait, cat.TF,
return( res)
}
fun.merge.biomes <- function(data){
data$biomes <- factor(data$biomes)
levels.name <- levels(data$biomes)
levels.name[levels.name == 'Tundra'] <- 'Boreal forest'
levels.name[levels.name == 'Subtropical desert'] <- 'Temperate grassland desert'
levels(data$biomes) <- levels.name
return(data)
}
load.and.save.data.for.lmer <- function(trait,
min.obs= 10,
type.filling = species,
fname = 'data.all.no.log.rds',
base.dir = "output/processed"){
base.dir = "output/processed",
merge.biomes.TF = TRUE){
data.tree.tot <- readRDS(file.path(base.dir, fname))
if(merge.biomes.TF) data.tree.tot <- fun.merge.biomes(data.tree.tot)
df <- fun.select.data.for.lmer(data.tree.tot, trait,
type.filling = type.filling)
if (fname == 'data.all.no.log.rds') type <- 'no.log'
......
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.biomes",
pars = c('intercept' , 'mean_logD', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_logD_species', 'sigma_sumBn_species', 'sigma',
'sigma_sumBn_biomes', 'sigma_sumTfBn_biomes',
'sigma_sumTnBn_biomes', 'sigma_sumTnTfBn_abs_biomes',
'param_Tf_biomes', 'param_sumBn_biomes',
'param_sumTfBn_biomes', 'param_sumTnBn_biomes',
'param_sumTnTfBn_abs_biomes'),## parameters to save
init.fun = function(chain_id= 1, list){
step <- (chain_id-1)/10
init <-
list(intercept = 0+step,
intercept_species = rnorm(list$N_species, 0, 0.3+step),
param_sumBn_species = rnorm(list$N_species, 0, 0.1+step),
param_logD_species = rnorm(list$N_species, 0, 0.3+step),
intercept_plot = rnorm(list$N_plot, 0, 0.3 +step),
intercept_set = rnorm(list$N_set, 0, 0.1 +step),
mean_logD = (2/3 - 0.1) +step,
param_Tf_biomes= rep(-0.1+step, list$N_biomes),
param_sumBn_biomes= rep(-0.1+step, list$N_biomes),
param_sumTfBn_biomes= rep(-0.1+step, list$N_biomes),
param_sumTnBn_biomes= rep(-0.1+step, list$N_biomes),
param_sumTnTfBn_abs_biomes= rep(-0.1+step, list$N_biomes),
sigma_inter_species = 0.3 +step,
sigma_inter_plot = 0.3 + +step,
sigma_inter_set = 0.1 + +step,
sigma = 0.5 +step,
sigma_sumBn_species = 0.1 +step,
sigma_sumBn_biomes = 0.1 +step,
sigma_sumTfBn_biomes = 0.1 +step,
sigma_sumTnBn_biomes = 0.1 +step,
sigma_sumTnTfBn_abs_biomes = 0.1 +step,
sigma_logD_species = 0.3 +step)
return(init)
},
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_biomes;
int<lower=0> N_set;
int species_id[N_indiv];
int plot_id[N_indiv];
int biomes_id[N_indiv];
int set_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
real Tf[N_indiv];
real sumTfBn[N_indiv];
real sumTnBn[N_indiv];
real sumTnTfBn_abs[N_indiv];
}
parameters{
real<lower=-5,upper=5> intercept;
real<lower=-5,upper=5> mean_logD;
real<lower=0.0001,upper=5> sigma_inter_species;
real<lower=0.0001,upper=5> sigma_inter_set;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_logD_species;
real<lower=0.0001,upper=5> sigma_sumBn_species;
real<lower=0.0001,upper=5> sigma_sumBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTfBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_biomes;
real<lower=0.0001,upper=5> sigma;
real<lower=-5,upper=5> param_Tf_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_species[N_species];
real<lower=-5,upper=5> param_sumTfBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnTfBn_abs_biomes[N_biomes];
real<lower=-5,upper=5> param_logD_species[N_species];
real<lower=-5,upper=5> intercept_species[N_species];
real<lower=-5,upper=5> intercept_plot[N_plot];
real<lower=-5,upper=5> intercept_set[N_set];
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_species[species_id[i]]
+ intercept_plot[plot_id[i]]
+ intercept_set[set_id[i]]
+ (mean_logD + param_logD_species[species_id[i]])*logD[i]
+ (param_Tf_biomes[biomes_id[i]])*Tf[i]
+ ( param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
+ (param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
+ (param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
;
}
}
model {
# constants for prior
real sigma0;
sigma0 <- 10;
###############################################
########### Hierarchical parameters ########
### species random param
param_logD_species ~ normal(0,sigma_logD_species);
param_sumBn_species ~ normal(0,sigma_sumBn_species);
intercept_species ~ normal(0, sigma_inter_species);
### plot random param
intercept_plot ~ normal(0,sigma_inter_plot);
### set random effect
intercept_set ~ normal(0, sigma_inter_set);
### biomes random effect
param_sumBn_biomes ~ normal(0, sigma_sumBn_biomes);
param_sumTfBn_biomes ~ normal(0, sigma_sumTfBn_biomes);
param_sumTnBn_biomes ~ normal(0, sigma_sumTnBn_biomes);
param_sumTnTfBn_abs_biomes ~ normal(0, sigma_sumTnTfBn_abs_biomes);
###############################################
########### Non-Hierarchical parameters ########
# slope and intercept
intercept ~ normal(0,sigma0);
mean_logD ~ normal(0,sigma0);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
}
......@@ -10,6 +10,9 @@ model.files.stan.Tf.1 <-
c("R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.fixed.R",
"R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.fixed.biomes.R")
model.files.stan.Tf.2 <-
c("R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.biomes.R")
### NEED TO MODIFY THIS FUNCTION work with all parameters in stan
......@@ -52,7 +55,7 @@ fun.call.stan.and.save <- function(stan.model, list.stan, path.out,
var.sample){
start <- Sys.time()
set_cppo(mode = "fast")
inits <- fun.generate.init(stan.model$init.fun,
inits <- fun.generate.init(stan.model$init.fun,
list.stan, chains)
stan.output <- stan(model_code = stan.model$stan,
data = list.stan,
......@@ -122,7 +125,7 @@ run.stan <- function (model.file, trait, init.TF,
var.sample = var.sample,
chain_id = chain_id,
...)
}
}
}
gc()
}
......
......@@ -34,7 +34,6 @@ GetClimate <-function(lats,lons) {
#get the set of tiles corresponding to the lons/lats arguments
tiles <- extract(NumRast,cbind(lons,lats))
unique.tiles <- unique(tiles)
plot.mat <- numeric(length(tiles))
......@@ -42,7 +41,7 @@ GetClimate <-function(lats,lons) {
#download each of the tiles
for (i.tile in unique.tiles) {
## browser()
temp.tile <- getData('worldclim',
var = 'bio',
res = 0.5,
......@@ -58,7 +57,14 @@ GetClimate <-function(lats,lons) {
cbind(lons[tiles==i.tile],
lats[tiles==i.tile]),
layer=12,nl=1) #retrieve MAP
## plot.map <- extract(rast,
## cbind(lons,
## lats),
## nl=1) #retrieve MAP
## plot.mat <- extract(rast1,
## cbind(lons,
## lats),
## nl=1)/10 #retrieve MAP
nas <- which(tiles==i.tile & is.na(plot.mat))
#assign climate of nearest plot to all nas
......
......@@ -96,6 +96,12 @@ plot.ellips.on.biome.map(MAP = clim$MAP/10,
MAT = clim$MAT,
set = clim$set, cols = cols.vec)
##+ Fig1b, fig.cap="**Number of tree per Whittaker biomes**", fig.width=7, fig.height=7, echo = FALSE, results = 'hide', warning=FALSE, message=FALSE
data.all <- readRDS.root('output/processed/data.all.no.log.rds')
par(mar = c(3, 13, 3,3))
barplot(table(data.all$biomes), horiz = TRUE, las = 2)
rm(data.all)
## \newpage
## #Results
......@@ -104,11 +110,15 @@ plot.ellips.on.biome.map(MAP = clim$MAP/10,
## re-sampling with weights equal to inverse focal species abundance
list.all.results.species.id <-
readRDS.root('output/list.lmer.out.all.species.id.no.log.old2.rds')
readRDS.root('output/list.lmer.out.all.sp.name.no.log.rds')
## resampling with weights equal to inverse focal ecocode abundance
list.all.results.ecocode.id <-
readRDS.root('output/list.lmer.out.all.ecocode.id.no.log.old2.rds')
readRDS.root('output/list.lmer.out.all.ecocode.no.log.rds')
## No resampling
list.all.results <-
readRDS.root('output/list.lmer.out.all.NA.no.log.rds')
##+ check.convergence, echo = FALSE, results = 'hide'
......@@ -126,7 +136,7 @@ names.biomes <- c('Subtrop', 'Temp grassland', 'Medit', 'Temp forest',
##+ Fig2, fig.cap='**Effect size parameters with data re-sampled per ecocode**', fig.width=10, fig.height=5, echo = FALSE
plot.param(list.all.results.ecocode.id ,
model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species',
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.species',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
......@@ -189,7 +199,27 @@ pandoc.table(mat.param, caption = "Full effect size")
##+ Fig3, fig.cap='**effect size parameters with data re-sampled per species.**', fig.width=10, fig.height=5, echo = FALSE
plot.param(list.all.results.species.id ,
model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species',
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
param.names = c(expression(paste("Direct trait effect ",
(m[1] %*% t[f]))),
expression(paste('Comp mean ', (c[0]))),
expression(paste('Compet response ',
(c[r] %*% t[f]))),
expression(paste('Compet effect ',
(c[e] %*% t[n]))),
expression(paste('Limit simil ',
(c[l] %*% abs(t[n] - t[f]))))) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.42, 0.18))
##+ Fig3b, fig.cap='**effect size parameters with no subsampling.**', fig.width=10, fig.height=5, echo = FALSE
plot.param(list.all.results ,
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
......@@ -207,6 +237,7 @@ plot.param(list.all.results.species.id ,
names.bio = names.biomes ,
xlim = c(-0.42, 0.18))
## The type of re-sampling of the data (ecocode or species) doesn't change a lot the results. I have done this analysis 4 times (so each time with a different re-sampling) and the results are quite stable.
## \newpage
......
......@@ -24,7 +24,7 @@ data.canada <- subset(data.canada,
subset = (data.canada$SubPlotSize>=0.02 &
data.canada$InitDBH>=10) )
sum(is.na(data.canada$species))
sum(is.na(data.canada$species_FIACode))
data.canada$sp <- data.canada$species_FIACode
data.canada$sp[data.canada$sp==8349] <- 762
#MCV replace incorrect sp code for Prunus serotina
......@@ -48,15 +48,15 @@ data.canada$sp.name[is.na(data.canada$sp.name)] <- data.canada$sp[is.na(data.can
## FORMAT INDIVIDUAL TREE DATA
data.canada$FinalDBH[data.canada$FinalDBH < 0] <- NA
data.canada$G <- 10 * (data.canada$FinalDBH - data.canada$InitDBH)/
data.canada$Interval ## diameter growth in mm per year
data.canada$IntervalYears ## diameter growth in mm per year
data.canada$BA.G <- (pi*(data.canada$FinalDBH/2)^2 -
pi*(data.canada$InitDBH/2)^2)/data.canada$Interval
pi*(data.canada$InitDBH/2)^2)/data.canada$IntervalYears
## diameter growth in mm per year
data.canada$G[which(data.canada$InitDBH == 0 |
data.canada$FinalDBH == -999)] <- NA
data.canada$BA.G[which(data.canada$InitDBH == 0 |
data.canada$FinalDBH == -999)] <- NA
data.canada$year <- data.canada$Interval; data.canada$IntervalYears <- NULL
data.canada$year <- data.canada$IntervalYears; data.canada$IntervalYears <- NULL
## number of year between measurement
data.canada$D <- data.canada[["InitDBH"]]
data.canada$D[data.canada$D == 0] <- NA ## diameter in cm
......@@ -105,7 +105,9 @@ data.canada <- merge(data.canada, data.frame(plot = names(perc.dead),
source("R/utils/climate.R")
clim <- GetClimate(data.canada$Lat, data.canada$Lon)
data.canada$MAT <- clim$MAT
data.canada$MAP2 <- clim$MAP
data.canada$MAP <- clim$MAP
######################## BIOMES
......
#!/usr/bin/env Rscript
### Download and merge Germaby
source("R/format.data/format-fun.R")
library(reshape, quietly = TRUE)
dir.create("output/formatted/Germany", recursive = TRUE, showWarnings = FALSE)
#############
## DOWNLOAD DATA FROM WEBSITE
data_path <- "data/raw/Germany"
gdi.web1 <- "https://gdi.ti.bund.de/geoserver/bwi_2009dt/ows?typeName=bwi_2009dt:"
gdi.web2 <- "&service=WFS&version=1.0.0&request=GetFeature&outputFormat=SHAPE-ZIP"
files <- c("bwi2002dt_wzp4", "bwi2002dt_trakte",
"bwi2002dt_ecken", "bwi2002dt_waldecken",
"bwi1987dt_wzp4", "bwi1987dt_trakte",
"bwi1987dt_ecken", "bwi1987dt_waldecken")
download.file2 <- function(x, web1, web2, path){
download.file(paste0(web1, x, web2),
destfile = file.path(path, paste0(x, ".zip")),
method = "curl",
quiet = TRUE)
}
## download all data
lapply(files, FUN = download.file2 ,
web1 = gdi.web1,
web2 = gdi.web2,
path = data_path)
#unzip data
fun.unzip <- function(file, data_path){
unzip(file.path(data_path, paste0(file, ".zip")), exdir = file.path(data_path, file))
}
lapply(files, FUN = fun.unzip ,
data_path = data_path)
## load data
library(rgdal)
library(sp)
## library(GISTools)
## wzp_02<- readShapePoly(file.path(data_path, "bwi2002dt_wzp4",
## "bwi2002dt_wzp4"),
## proj4string=CRS("+init=epsg:5676"))
wzp_02 <- readOGR(file.path(data_path, "bwi2002dt_wzp4"), "bwi2002dt_wzp4",
p4s = "+init=epsg:31466")
wzp_87 <- readOGR(file.path(data_path, "bwi1987dt_wzp4"), "bwi1987dt_wzp4",
p4s = "+init=epsg:31466")
class(wzp_02)
## save as rds to load faster later on
saveRDS(wzp_87, file.path(data_path,'wzp_87.rds'))
saveRDS(wzp_02, file.path(data_path,'wzp_02.rds'))
centroids_02<- getSpPPolygonsLabptSlots(wzp_02)
centroids_87<- getSpPPolygonsLabptSlots(wzp_87)
colnames(centroids_87) <- colnames(centroids_02) <- c('x', 'y')
wzp_02_df<- cbind(as.data.frame(wzp_02), centroids_02)
wzp_87_df<- cbind(as.data.frame(wzp_87), centroids_87)
####### change coordinates system of x y to be in lat long WGS84
library(sp)
library(dismo)
library(rgdal)
data.sp <- wzp_87_df[, c("x", "y")]
coordinates(data.sp) <- c("x", "y")
proj4string(data.sp) <- CRS("+init=epsg:2166")
# define projection system of our data
summary(data.sp)
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326"))
## change projection in WGS84 lat lon
wzp_87_df$Lon <- coordinates(data.sp2)[, "x"]
wzp_87_df$Lat <- coordinates(data.sp2)[, "y"]
## plot on world map
library(rworldmap)
newmap <- getMap(resolution = 'high')
## different resolutions available
plot(newmap, xlim = c(7, 16), ylim = c(47, 55))
points(data.sp2, cex = 0.2, col = 'red')
plot(wzp_87_df$x, wzp_87_df$y, cex = 0.2, col = 'red')
## Need to ask Bjoern if good EPSG seems ok
### VARIABLES SELECTION
vars.select1 <- c('bhd_130',# diamter at 130 cm corrected in mm
'm_hoe',# measured height in dcm
'hori',# horizontal distance in cm
'al_ba1987', # age in year
'n_ha',
'tree.id',
'tnr',
'enr',
'bnr',
'ba_gattung',
'ba_art') #weight to have number of tree per ha
vars.select2 <- c('bhd_130',# diamter at 130 cm corrected in mm
'm_hoe',# measured height in dcm
'hori',# horizontal distance in cm
'tot', # dead
'al_ba2002', # age in year
'n_ha',
'tree.id',
'tnr',
'enr',
'bnr',
'ba_gattung',
'ba_art') #weight to have number of tree per ha
#### NEED TO FIND VARIABLES ON DISTURBANCE HARVESTING AND PLANTATION
wzp_87_df$tree.id <- paste(wzp_87_df$tnr, wzp_87_df$enr, wzp_87_df$bnr,sep='_')
wzp_02_df$tree.id <- paste(wzp_02_df$tnr, wzp_02_df$enr, wzp_02_df$bnr,sep='_')
data_germany<- merge(wzp_87_df[, vars.select1], wzp_02_df[, vars.select2],
by = c("tree.id"), all.x = TRUE, all.y = FALSE)
## check growth
plot(data_germany$bhd_130.x, data_germany$bhd_130.y - data_germany$bhd_130.x)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment