Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
merge.data.SWISS.R 5.65 KiB
1 
### MERGE Swiss DATA
2 
rm(list = ls())
3 
source("./R/format.function.R")
4 
source("./R/FUN.TRY.R")
5 
library(reshape)
6 
library(foreign)
7
8 
######################### READ DATA read individuals tree data
9 
data.swiss1 <- read.csv("./data/raw/DataSwiss/LFI12.csv", header = TRUE, stringsAsFactors = FALSE)
10 
data.swiss2 <- read.csv("./data/raw/DataSwiss/LFI23.csv", header = TRUE, stringsAsFactors = FALSE)
11 
data.swiss3 <- read.csv("./data/raw/DataSwiss/LFI34.csv", header = TRUE, stringsAsFactors = FALSE)
12 
data.swiss <- rbind(data.swiss1, data.swiss2, data.swiss3)
13 
rm(data.swiss1, data.swiss2, data.swiss3)
14 
data.swiss <- data.swiss[order(data.swiss$BANR), ]
15 
data.swiss <- data.swiss[, c("INVNR", "CLNR", "BANR", "X", "Y", "BAUMART", "TEXT",
16 
    "VEGPER", "BHD1", "BHD2", "BA1", "BA2", "BAI", "BHD_DIFF", "RPSTZ1", "RPSTZ2",
17 
    "HOEHE1", "HOEHE2")]
18
19 
colnames(data.swiss) <- c("Inventid", "siteid", "tree.id", "x", "y", "sp", "sp.name",
20 
    "year", "dbh1", "dbh2", "ba1", "ba2", "ba_diff", "dbh_diff", "repfactor1", "repfactor2",
21 
    "ht1", "ht2")
22
23 
#################
24 
##### change the projection of xy to wgs84 lat long.
25
26 
#### change coordinates system of x y to be in lat long WGS84
27 
library(sp)
28 
library(dismo)
29 
library(rgdal)
30
31 
data.sp <- data.swiss[, c("tree.id", "x", "y")]
32 
coordinates(data.sp) <- c("x", "y")  # define x y
33 
proj4string(data.sp) <- CRS("+init=epsg:21781")  # define projection system of our data ## EPSG:21781 CH1903 / LV03 CHECK WITH NICK STRANGE
34 
summary(data.sp)
35
36 
library(maptools)
37 
bioregion <- readShapePoly("./data/raw/DataSwiss/bio_regions/biogregions.shp")
38 
plot(bioregion)
39 
plot(data.sp,add=TRUE)
40 
test <-  over(data.sp,bioregion) ## need to check proj with nick
41
42
43 
detach(package:rgdal)
44 
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326"))  ## change projection in WGS84 lat lon
45 
data.swiss$Lon <- coordinates(data.sp2)[, "x"]
46 
data.swiss$Lat <- coordinates(data.sp2)[, "y"]
47 
## ## ## plot on world map
48 
library(rworldmap)
49 
newmap <- getMap(resolution = 'coarse')
50 
## # different resolutions available
51 
plot(newmap,xlim=c(5,11),ylim=c(45,48))
52 
points(data.sp2,cex=0.2,col='red')
53 
points(data.swiss$x,data.swiss$y)
54
55
56 
rm(data.sp, data.sp2)
57
58 
## Do not need to read in spp list as it is already available in data.swiss
59
60 
###################################### MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
61 
###################################### height to add variables to the traits data base Because we have two heights,
62 
###################################### then take the max of the two heights and then bootstrap
63 
res.quant.boot <- t(sapply(levels(factor(data.swiss[["sp"]])), FUN = f.quantile.boot,
64 
    R = 1000, x = log10(apply(data.swiss[, c("ht1", "ht2")], 1, max, na.rm = T)),
65 
    fac = factor(data.swiss[["spcode"]])))
66
67 
## create data base
68 
data.max.height <- data.frame(code = rownames(res.quant.boot), Max.height.mean = res.quant.boot[,
69 
    1], Max.height.sd = res.quant.boot[, 2], Max.height.nobs = res.quant.boot[, 3], stringsAsFactors =FALSE)
70 
rm(res.quant.boot)
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125
write.csv(data.max.height,file='./data/process/data.max.height.swiss.csv') ########################################## FORMAT INDIVIDUAL TREE DATA data.swiss$G <- 10 * (data.swiss$dbh_diff)/data.swiss$year ## diameter growth in mm per year data.swiss$D <- data.swiss[["dbh1"]] data.swiss$D[data.swiss$D == 0] <- NA ## diameter in cm data.swiss$dead <- rep(NA, nrow(data.swiss)) ## Mortality - MISSING data.swiss$plot.id <- data.swiss$siteid; data.swiss$siteid <- NULL ## plot code data.swiss$htot <- data.swiss$ht1 ## height of tree in m data.swiss$census <- data.swiss$Inventid data.swiss$obs.id <- apply(data.swiss[,c("tree.id","Inventid")],1,paste,collapse="_") data.swiss$weights <- data.swiss$repfactor1/10000 data.swiss$sp <- paste("sp.",data.swiss$sp,sep="") ###################### ECOREGION Ecoregion not available for swiss data data.swiss$ecocode <- rep("A", nrow(data.swiss)) ###################### PERCENT DEAD variable percent dead/cannot do with since dead variable is ###################### missing compute numer of dead per plot to remove plot with disturbance perc.dead <- tapply(data.swiss[["dead"]], INDEX = data.swiss[["plot"]], FUN = function.perc.dead2) data.swiss <- merge(data.swiss, data.frame(plot = names(perc.dead), perc.dead = perc.dead, stringsAsFactors =FALSE), by = "plot", sort = FALSE) ########################################################### PLOT SELECTION FOR THE ANALYSIS data.climate <- read.dbf(file = "./data/raw/DataSwiss/LFI14_climate.dbf") data.climate <- data.climate[, c("CLNR","swb_100","tave_68","prec_122","prec_35","prec_68","prec_911")] data.climate$MAP <- apply(data.climate[, c("prec_122","prec_35","prec_68","prec_911")], 1, sum) data.swiss <- merge(data.swiss, data.frame(siteid = data.climate$CLNR, swb = data.climate$swb_100, MAT = data.climate$tave_68, MAP = data.climate$MAP, stringsAsFactors =FALSE), sort = F, all.x = T) rm(data.climate) ############################################## #################### GENERATE ONE OBJECT PER ECOREGION # vector of ecoregion name ecoregion.unique <- unique(data.swiss[["ecocode"]]) ### read TRY data TRY.DATA.FORMATED <- readRDS("./data/process/TRY.DATA.FORMATED.rds") ## create lookup table for swiss species.lookup <- data.frame(sp=data.swiss[!duplicated(data.swiss[["sp"]]),"sp"], Latin_name=data.swiss[!duplicated(data.swiss[["sp"]]),"sp.name"], Latin_name_syn=data.swiss[!duplicated(data.swiss[["sp"]]),"sp.name"], stringsAsFactors =FALSE) #### lapply function to generate data for each ecoregion system.time(lapply(ecoregion.unique, FUN = fun.data.per.ecoregion, data.tot = data.swiss, plot.name = "plot", weight.full.plot = NA, name.country = "Swiss", data.TRY = TRY.DATA.FORMATED, species.lookup = species.lookup))