progress on Germany

scripts/format.data/Canada.R

@@ -46,10 +46,10 @@ data.canada$sp[data.canada$sp==8349] <- 762
 ### read species names and merge with data.canada
 species.clean <- read.csv("data/raw/Canada/FIA_REF_SPECIES.csv",
                           stringsAsFactors = FALSE)
-species.clean <- mutate(species.clean,
+species.clean <- dplyr::mutate(species.clean,
                         Latin_name = paste(GENUS,
                                            SPECIES, sep = " "))
-species.clean <- mutate(species.clean,
+species.clean <- dplyr::mutate(species.clean,
                         Latin_name_syn = paste(GENUS,
                                                SPECIES, sep = " "))
 
@@ -59,11 +59,11 @@ data.canada <- left_join(data.canada,
                                   stringsAsFactors = F),
                        by = "sp")
 
-data.canada <- mutate(data.canada, sp =paste("sp", sp, sep = "."))
-data.canada <- select(data.canada,
+data.canada <- dplyr::mutate(data.canada, sp =paste("sp", sp, sep = "."))
+data.canada <- dplyr::select(data.canada,
                       sp.name = Latin_name,
                       matches("."))
-data.canada <- mutate(data.canada,
+data.canada <- dplyr::mutate(data.canada,
                       sp.name = ifelse(is.na(sp.name), sp, sp.name))
 
 ## MASSAGE TRAIT DATA HEIGHT DATA FOR TREE MISSING BRING US DATA FOR HEIGHT OVER
@@ -117,9 +117,9 @@ data.canada <- subset(data.canada, subset = (data.canada$ecocode != "-251") &
 ##### compute numer of dead per plot to remove plot with disturbance
 perc.dead <- tapply(data.canada[["dead"]], INDEX = data.canada[["plot"]],
                    FUN = function.perc.dead2)
-data.canada <- merge(data.canada, data.frame(plot = names(perc.dead),
+data.canada <- left_join((data.canada, data.frame(plot = names(perc.dead),
                                              perc.dead = perc.dead),
-                      by = "plot", sort = FALSE)
+                      by = "plot")
 
 
 #### get wc climate AND OVERRIDE THE DATA PROVIDED BECAUSE UNIT ERROR
@@ -150,7 +150,7 @@ source("R/utils/ecoregions.R")
 data.p <- data.canada[!duplicated(data.canada$plot), ]
 data.p$koppen <- GetKoppen(data.p$Lon, data.p$Lat)
 data.p$wwf <- GetEcoregions(data.p$Lon, data.p$Lat)
-data.p$wwf <- as.character(data.p$wwf) 
+data.p$wwf <- as.character(data.p$wwf)
 data.canada <- merge(data.canada, data.p[, c('plot', 'koppen', 'wwf')], by = 'plot')
 
 

scripts/format.data/Germany.R

@@ -38,7 +38,7 @@ lapply(files, FUN = download.file2 ,
 fun.unzip <- function(file, data_path){
   unzip(file.path(data_path, paste0(file, ".zip")), exdir = file.path(data_path, file))
   obj <- readOGR(file.path(data_path, file), file, encoding = "latin1")
-  saveRDS(obj, file.path(data_path, paste0(file, ".rds")))  
+  saveRDS(obj, file.path(data_path, paste0(file, ".rds")))
   unlink(file.path(data_path, file), recursive = TRUE)
 }
 
@@ -50,7 +50,7 @@ library(rgdal)
 library(sp)
 ## library(GISTools)
 ## wzp_02<- readShapePoly(file.path(data_path, "bwi2002dt_wzp4",
-##                                  "bwi2002dt_wzp4"), 
+##                                  "bwi2002dt_wzp4"),
 ##                        proj4string=CRS("+init=epsg:5676"))
 
 wzp_02 <- readOGR(file.path(data_path, "bwi2002dt_wzp4"), "bwi2002dt_wzp4")
@@ -64,7 +64,7 @@ saveRDS(wzp_02, file.path(data_path,'wzp_02.rds'))
 }else{
 wzp_87 <- readRDS(file.path(data_path,'wzp_87.rds'))
 wzp_02 <- readRDS(file.path(data_path,'wzp_02.rds'))
-}    
+}
 # get proj
 p4.bwi <- proj4string(wzp_87)
 # get centroids
@@ -74,19 +74,32 @@ 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")  
+coordinates(data.sp) <- c("x",  "y")
 proj4string(data.sp) <- p4.bwi
-# define projection system of our data 
+# 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"]
+
+#02
+data.sp <- wzp_02_df[,  c("x",  "y")]
+coordinates(data.sp) <- c("x",  "y")
+proj4string(data.sp) <- p4.bwi
+# 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_02_df$Lon <- coordinates(data.sp2)[,  "x"]
+wzp_02_df$Lat <- coordinates(data.sp2)[,  "y"]
+
 ## plot on world map
 library(rworldmap)
 newmap <- getMap(resolution = 'high')
@@ -126,6 +139,22 @@ vars.select2 <-  c('bhd_130',# diamter at 130 cm corrected in mm
 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='_')
 
+## add species name
+wzp_87_df$sp.name <- paste(wzp_87_df$ba_gattung, wzp_87_df$ba_art)
+wzp_02_df$sp.name <- paste(wzp_02_df$ba_gattung, wzp_02_df$ba_art)
+
+sp_vec <- sort(table(wzp_02_df$sp.name))
+sp_table <- data.frame(Species = names(sp_vec), Nobs = as.vector(sp_vec))
+write.csv(sp_table,
+          file.path("output", "formatted", "Germany", "species_Germany.csv"),
+          row.names = FALSE)
+
+write.csv(wzp_02_df[ , c("gid", "cell_id", "tnr",
+                         "enr", "bnr", "netz",
+                         "Lon", "Lat")],
+          file.path("output", "formatted", "Germany", "coord_Germany.csv"),
+          row.names = FALSE)
+
 data_germany<- merge(wzp_87_df[, vars.select1], wzp_02_df[, vars.select2],
                      by = c("tree.id"), all.x = TRUE, all.y = FALSE)
 

scripts/process.data/explore.processed.data.report.R

@@ -16,9 +16,9 @@ source("R/analysis/lmer.run.R")
 
 
 ### read all data
-data.all <-  readRDS(file.path('output', 'processed', "data.all.no.log.rds"))
+data.all <-  readRDS(file.path('output', 'processed', "data.all.no.log.t.rds"))
 
-data.glob <-  readRDS(file.path('output', 'processed', "data.all.global.rds"))
+data.glob <-  readRDS(file.path('output', 'processed', "data.all.global.t.rds"))
 
 
 library(dplyr)
@@ -178,7 +178,7 @@ data.all.t$biomes <-  factor(data.all.t$biomes, c("Boreal forest",
                                   "Tropical forest savanna",
                                   "Tropical rain forest"))
 
-pdf('figs/range_traits_per_biomes.pdf', height = 12, width = 12)
+pdf('figs/range_traits_per_biomes2.pdf', height = 6, width = 6)
 m <- matrix(c(1:4), 2, 2)
 layout(m, widths=c(4,4),
         heights=c(1,1))
@@ -204,6 +204,17 @@ legend("center",legend = levels(data.all.t$biomes),
        bty = 'n', cex = 1.2)
 dev.off()
 
+
+library(dplyr)
+
+df_biomes <- group_by(data.all.t, biomes) %>% summarise(SLAmean = mean(SLA.focal, na.rm = TRUE),
+                                WDmean = mean(Wood.density.focal, na.rm = TRUE),
+                                MHmean = mean(Max.height.focal, na.rm = TRUE),
+                                SLAmax = max(SLA.focal, na.rm = TRUE),
+                                WDmax = max(Wood.density.focal, na.rm = TRUE),
+                                MHmax = max(Max.height.focal, na.rm = TRUE)
+                                )
+
 rm(data.all.t)
 
 ## plot of correlation between BA.G and D the best is log(BA.G) vs log(D)
@@ -321,6 +332,37 @@ names.biomes.t[5] <-  'Trop forest'
 names.biomes.t[6] <-  'Trop rain forest'
 
 
+
+
+
+pdf('figs/CWMrange_traits_per_biomes2.pdf', height = 12, width = 12)
+m <- matrix(c(1:4), 2, 2)
+layout(m, widths=c(4,4),
+        heights=c(1,1))
+par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
+boxplot(mean.SLA ~ biomes, data.temp,
+        col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
+        ylab = expression(paste('SLA ' ,(mm^2/mg))),
+        outline = TRUE, xaxt = 'n',  cex.y = 1.5,cex.lab=1.2)
+par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
+boxplot(mean.Wood.density ~ biomes, data.temp,
+        col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
+        ylab = expression(paste('Wood density ' ,( mg/mm^3))),
+        outline = TRUE, xaxt = 'n',  cex.y = 1.5,cex.lab=1.2)
+par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
+boxplot(mean.Max.height ~ biomes, data.temp,
+        col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
+        ylab = expression(paste('Max height ' ,(m))),
+        outline = TRUE, xaxt = 'n',  cex.y = 1.5,cex.lab=1.2)
+par(mar=c(0, 0, 0, 0))
+plot(1,1,type="n", axes=F, xlab="", ylab="")
+legend("center",legend = levels(data.temp$biomes),
+       fill = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
+       bty = 'n', cex = 1.2)
+dev.off()
+
+
+
 ## need to look at what is Temperate grassland desert
 pdf('figs/biomes_div_sd_traits.pdf',
     width = 12,