Skip to content
GitLab
Select Git revision
  • dev default
  • master
  • 136-decreased-performance-of-calibration-execution-time
  • 115-add-pcict-date-format-managment
  • 118-plot-outputsmodel-remove-case-sensitivity-to-the-which-argument
  • 96-parallelize-the-grid-screening-step-during-calibration
  • 60-migrate-transfoparam-functions-and-createcaliboptions-to-s3-methods
  • airGRplus
  • cleanAirGRplus
  • CRAN_v1.7.0
  • CRAN_v1.6.12
  • CRAN_v1.6.11
  • CRAN_v1.6.10.4
  • CRAN_v1.6.9.27
  • CRAN_v1.6.9.24
  • CRAN_v1.6.9.23
  • CRAN_v1.6.9.21
  • CRAN_v1.4.3.65
  • CRAN_v1.4.3.60
  • CRAN_v1.4.3.52
  • CRAN_v1.3.2.42
  • CRAN_v1.3.2.41
  • CRAN_v1.3.2.40
  • CRAN_v1.3.2.39
  • CRAN_v1.3.2.23
  • CRAN_v1.3.2.17
  • CRAN_v1.2.13.16
  • CRAN_v1.2.13.13
  • CRAN_v1.0.14.1
Find file BlameHistoryPermalink
Forked from HYCAR-Hydro / airGR
Source project has a limited visibility.
merge.data.SWISS.R 5.69 KiB
1 
### MERGE Swiss DATA Edited by FH
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 
detach(package:rgdal)
37 
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326"))  ## change projection in WGS84 lat lon
38 
data.swiss$Lon <- coordinates(data.sp2)[, "x"]
39 
data.swiss$Lat <- coordinates(data.sp2)[, "y"]
40 
## ## ## plot on world map
41 
## library(rworldmap)
42 
## newmap <- getMap(resolution = 'coarse')
43 
## ## # different resolutions available
44 
## plot(newmap,xlim=c(5,11),ylim=c(45,48))
45 
## points(data.sp2,cex=0.2,col='red')
46 
## points(data.swiss$x,data.swiss$y)
47 
## rm(data.sp, data.sp2)
48
49
50
51 
## Do not need to read in spp list as it is already available in data.swiss
52
53 
###################################### MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
54 
###################################### height to add variables to the traits data base Because we have two heights,
55 
###################################### then take the max of the two heights and then bootstrap
56 
res.quant.boot <- t(sapply(levels(factor(data.swiss[["spcode"]])), FUN = f.quantile.boot,
57 
    R = 1000, x = log10(apply(data.swiss[, c("ht1", "ht2")], 1, max, na.rm = T)),
58 
    fac = factor(data.swiss[["spcode"]])))
59
60 
## create data base
61 
data.max.height <- data.frame(code = rownames(res.quant.boot), Max.height.mean = res.quant.boot[,
62 
    1], Max.height.sd = res.quant.boot[, 2], Max.height.nobs = res.quant.boot[, 3], stringsAsFactors =FALSE)
63 
rm(res.quant.boot)
64 
write.csv(data.max.height,file='./data/process/data.max.height.swiss.csv')
65
66 
########################################## FORMAT INDIVIDUAL TREE DATA change unit and names of variables to be the same
67 
########################################## in all data for the tree
68 
data.swiss$G <- 10 * (data.swiss$dbh_diff)/data.swiss$year  ## diameter growth in mm per year - SOME EXTREMELY NEGATIVE HERE!
69 
data.swiss$D <- data.swiss[["dbh1"]]
70 
data.swiss$D[data.swiss$D == 0] <- NA
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123
## diameter in cm data.swiss$dead <- rep(NA, length(data.swiss[["dbh1"]])) ## Mortality - MISSING data.swiss$plot <- data.swiss$siteid ## plot code data.swiss$htot <- data.swiss$ht1 ## height of tree in m data.swiss$census <- data.swiss$Inventid data.swiss$obs.id <- as.character(1:nrow(data.swiss)) 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.dead) # ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OTHER VARIABLES IF # AVAILABLE (disturbance record) 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 Remove data with dead == 1 ## Nothing to remove data.climate <- read.dbf(file = "./data/raw/DataSwiss/LFI14_climate.dbf") data.climate <- data.climate[, c(1, 7, 15:19)] data.climate$MAP <- apply(data.climate[, 4:7], 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))