Japan.R 9.6 KB
Newer Older
fhui28's avatar
fhui28 committed
1
2
3
4
#!/usr/bin/env Rscript

### MERGE japan DATA
rm(list = ls())
Georges Kunstler's avatar
Georges Kunstler committed
5
source("R/format.data/format-fun.R")
fhui28's avatar
fhui28 committed
6
7
8
9
10
11
12
13
14
15
16
17
library(reshape,quietly=TRUE)

dir.create("output/formatted/Japan", recursive=TRUE,showWarnings=FALSE)

######################### READ DATA read individuals tree data 
filedir <- "data/raw/Japan/TreeData/"
filenames <- c("AI-BC1","AU-DB1","GR-DB1","KK-DB1","OG-DB1","SD-EB1","TM-DB4","AM-EB1","AU-EC1",
	"HY-EC1","KM-DB1","OS-EC1","SI-DB1","TN-EB1","AO-BC1","AY-EB1","IC-BC1","KS-DB1","OT-EC1",
	"TB-DB1","UR-BC1","AS-DB1","CC-DB1","KA-EB1","KY-DB1","OW-EB1","TM-DB1","WK-EC1","AS-DB2",
	"CC-DB2","KG-EC1","NB-EC1","OY-DB1","TM-DB2","YK-EB1","AS-DB3","CC-DB3","KJ-EB1","NP-DB1",
	"OZ-DB1","TM-DB3","YN-EB1")

Georges Kunstler's avatar
Georges Kunstler committed
18
19
data.japan <- data.frame(x = 1, y = 1, id = 1, plot = 1, sp = 1, sp.name = 1, dbh1 = 1, dbh2 = 1,
                         date1 = 1, date2 = 1, status1 = 1, status2 = 1,er1=NA,er2=NA)
fhui28's avatar
fhui28 committed
20
21
22
23
24
for(i in 1:length(filenames)) {
	make.filename <- paste(filenames[i],"-tree-ver1.csv",sep="")
	cat("Working on dataset", make.filename, "\n")

	if(i == 6) next;
Georges Kunstler's avatar
Georges Kunstler committed
25
26
        if(filenames[i] %in% c('CC-DB1', 'CC-DB2', 'CC-DB3', 'OG-DB1', 'OY-DB1', 'SD-EB1', 'TM-DB4')) next;
        
fhui28's avatar
fhui28 committed
27
28
29
30
31
32
33
34
35
36
37
	tmp.dat <- read.csv(paste(filedir,make.filename,sep=""), header = T, stringsAsFactors = FALSE, fill = T)
	how.many.measures <- length(grep("gbh",colnames(tmp.dat)))
	if(how.many.measures == 1) next; 
	
	format.tmp.dat <- tmp.dat[,1:8] ## Relies on the first 8 columns being the same in each dataset
	format.tmp.dat <- format.tmp.dat[rep(1:nrow(format.tmp.dat),each=how.many.measures-1),]
	rownames(format.tmp.dat) <- 1:nrow(format.tmp.dat)

	get.gbh.cols <- grep("gbh",colnames(tmp.dat))
	get.dl.cols <- grep("dl",colnames(tmp.dat))
	get.date.cols <- grep("s_date",colnames(tmp.dat))
Georges Kunstler's avatar
Georges Kunstler committed
38
39
40
	get.error.cols <- grep("error",colnames(tmp.dat))
	format.tmp.dat <- data.frame(format.tmp.dat, dbh1 = as.vector(unlist(t(tmp.dat[,
                                                         get.gbh.cols[-length(get.gbh.cols)]]))),
fhui28's avatar
fhui28 committed
41
42
43
44
		dbh2 = as.vector(unlist(t(tmp.dat[,get.gbh.cols[-1]]))),
		status1 = as.vector(unlist(t(tmp.dat[,get.dl.cols[-length(get.dl.cols)]]))),
		status2 = as.vector(unlist(t(tmp.dat[,get.dl.cols[-1]]))),
		date1 = as.vector(unlist(t(tmp.dat[,get.date.cols[-length(get.date.cols)]]))),
Georges Kunstler's avatar
Georges Kunstler committed
45
46
47
48
49
50
51
52
53
54
55
56
57
58
		date2 = as.vector(unlist(t(tmp.dat[,get.date.cols[-1]]))),
                er1 = as.vector(unlist(t(tmp.dat[,get.error.cols[-length(get.error.cols)]]))),
		er2 = as.vector(unlist(t(tmp.dat[,get.error.cols[-1]])))                                     )

#Plots were subdivided into 10×10-m grid cells (Fig. 2).
#Exceptions were HY-EC1 which was subdivided into 10×10 or 10×5-m gird cells (see Appendix.pdf), and CC-DB1 which was subdivided into sixteen 25×25-m grid cells. Should not lead to problems

## IN CC-DB1, CC-DB2, CC-DB3, OG-DB1, OY-DB1, SD-EB1, TM-DB4 where stem locations were not measured        
        
	format.tmp.dat$x <- format.tmp.dat$grid_xcord+format.tmp.dat$stem_xcord; 
	format.tmp.dat$y <- format.tmp.dat$grid_ycord+format.tmp.dat$stem_ycord;
## plot(format.tmp.dat$x,format.tmp.dat$y)
        format.tmp.dat$tag_no <- NULL;
	format.tmp.dat$id <- paste(format.tmp.dat$indv_no,format.tmp.dat$tag_no,sep=""); format.tmp.dat$indv_no <- NULL
fhui28's avatar
fhui28 committed
59
60
61
62
63
	format.tmp.dat$stem_xcord <- NULL; 
	format.tmp.dat$stem_ycord <- NULL
	format.tmp.dat$plot <- filenames[i]
	format.tmp.dat$sp_jpn <- NULL; 
	format.tmp.dat$sp.name <- format.tmp.dat$sp; 
Georges Kunstler's avatar
Georges Kunstler committed
64
65
66
67
	format.tmp.dat$dbh1 <-  format.tmp.dat$dbh1/pi
	format.tmp.dat$dbh2 <-  format.tmp.dat$dbh2/pi

        format.tmp.dat <- format.tmp.dat[,colnames(data.japan)]
fhui28's avatar
fhui28 committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
	
	#data.japan <- merge(data.japan, format.tmp.dat, by = "sp")
	data.japan <- rbind(data.japan, format.tmp.dat)
	cat("Combined dataset now has", length(unique(data.japan$id)), "individuals\n")
	cat("Combined dataset now has", nrow(data.japan), "rows\n\n")
	rm(tmp.dat)
	}
	
data.japan <- data.japan[-1,]
data.japan$date1 <- as.Date(as.character(data.japan$date1), format = "%Y%m%d")
data.japan$date2 <- as.Date(as.character(data.japan$date2), format = "%Y%m%d")
	
site.data <- read.csv(file = "data/raw/Japan/ERDP_2011_01.5.1-SiteList.csv", header = T)
data.japan <- merge(data.japan, data.frame(plot = site.data$PlotID, Lat = site.data$Latitude, 
Georges Kunstler's avatar
Georges Kunstler committed
82
	Lon = site.data$Longitude,mat = site.data$Temp30yr, map = site.data$Rain30yr), by = "plot")
Georges Kunstler's avatar
Georges Kunstler committed
83
84
85
86
87
88
## library(rworldmap,quietly=TRUE)
## newmap <- getMap(resolution = 'coarse')
## plot(newmap)
## points(site.data$Longitude,site.data$Latitude,col='red')
## x11()
## plot(site.data$Temp30yr,site.data$Rain30yr)
fhui28's avatar
fhui28 committed
89

Georges Kunstler's avatar
Georges Kunstler committed
90
91


fhui28's avatar
fhui28 committed
92
################################## FORMAT INDIVIDUAL TREE DATA
Georges Kunstler's avatar
Georges Kunstler committed
93
94
95
96
97
98
99
100
data.japan$year <- as.numeric(difftime(data.japan$date2, data.japan$date1, units = "weeks")/52) ## number of year between measurement
data.japan$G <- 10*(data.japan[["dbh2"]] - data.japan[["dbh1"]])/data.japan$year ## diameter growth in mm per year
data.japan$BA.G <-  (pi*(data.japan[["dbh2"]]/2)^2 -pi*(data.japan[["dbh1"]]/2)^2)/data.japan[["year"]] ## BA growth in cm2
data.japan$G[data.japan$status1 > 0 | data.japan$status2>0] <- NA
data.japan$BA.G[data.japan$status1 > 0 | data.japan$status2>0] <- NA
data.japan$G[data.japan$er1 > 0 | data.japan$er2>0] <- NA
data.japan$BA.G[data.japan$er1 > 0 | data.japan$er2>0] <- NA

fhui28's avatar
fhui28 committed
101
102
103
104
data.japan$D <- data.japan[["dbh1"]]  ## diameter in mm convert to cm
data.japan$D[data.japan$status1 > 0] <- NA
data.japan$dead1 <- as.numeric(data.japan$status1 > 0)  ## dummy variable for dead tree 0 alive 1 dead
data.japan$dead2 <- as.numeric(data.japan$status2 > 0)  ## dummy variable for dead tree 0 alive 1 dead
Georges Kunstler's avatar
Georges Kunstler committed
105
data.japan$census <- substr(data.japan$date1,1,4)
Georges Kunstler's avatar
Georges Kunstler committed
106
107
108
data.japan$cluster <- data.japan$plot ## cluster number
data.japan$plot <- NA ## no plot

fhui28's avatar
fhui28 committed
109
data.japan$htot <- rep(NA,nrow(data.japan)) ## height of tree in m - missing
110
data.japan$tree.id <- paste(data.japan$cluster,data.japan$id) ## tree unique id
fhui28's avatar
fhui28 committed
111
112
113
114
115
116
data.japan$obs.id <- 1:nrow(data.japan) ## obs uniquue id

data.japan <- subset(data.japan,subset= (data.japan$dead1 == 0))
data.japan$dead1 <- NULL
data.japan$dead <- data.japan$dead2
data.japan$dead2 <- NULL
Georges Kunstler's avatar
Georges Kunstler committed
117

118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
### TODO CHANGE PLOT SIZE AS other

#remove tree with no x y
data.japan <- data.japan[!(is.na(data.japan$x) | is.na(data.japan$y)),]

## create small plot
square.s.t <- 20
make.quad <- do.call("rbind",lapply(unique(data.japan$cluster),
                           FUN=fun.quadrat.cluster,
                           cluster=data.japan$cluster,
                           tree.id=data.japan$obs.id,
                           x=data.japan$x,
                           y=data.japan$y,
                           square.s=square.s.t))
data.japan <- merge(data.japan, make.quad,by.x='obs.id',by.y="tree.id")


Georges Kunstler's avatar
Georges Kunstler committed
135
## ## #plot each plot
136
137
138
139
## pdf("figs/plots.japan.pdf")
## lapply(unique(data.japan[["cluster"]]),FUN=fun.circles.plot,
## data.japan[['x']],data.japan[['y']],
##        data.japan[["cluster"]],data.japan[["D"]],
Georges Kunstler's avatar
Georges Kunstler committed
140
## inches=0.1,fg.l=data.japan$make.quad)
141
## dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
142
## #
143
144

data.japan$plot <- data.japan$make.quad
Georges Kunstler's avatar
Georges Kunstler committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185

## LOAD WORLDCLIM
#### get wc climate 
source("R/utils/climate.R")
clim <- GetClimate(data.japan$Lat,data.japan$Lon)
data.japan$MAT <- clim$MAT
data.japan$MAP <- clim$MAP

## ## LOAD ecoregion WWF
## source("R/utils/ecoregions.R")
## data.japan$ecocode <- GetEcoregions(data.japan$Lon,data.japan$Lat)
## data.japan$ecocode <- as.character(data.japan$ecocode)
## ### NEED TO VERIFY THAT ECOREGIONS ARE ASSIGNED PROPERLY
## table(as.character(data.japan$ecocode),data.japan$cluster)
## data.japan$ecocode <- NULL

filenames <- c("AI-BC1","AU-DB1","GR-DB1","KK-DB1","OG-DB1","SD-EB1","TM-DB4","AM-EB1","AU-EC1",
	"HY-EC1","KM-DB1","OS-EC1","SI-DB1","TN-EB1","AO-BC1","AY-EB1","IC-BC1","KS-DB1","OT-EC1",
	"TB-DB1","UR-BC1","AS-DB1","CC-DB1","KA-EB1","KY-DB1","OW-EB1","TM-DB1","WK-EC1","AS-DB2",
	"CC-DB2","KG-EC1","NB-EC1","OY-DB1","TM-DB2","YK-EB1","AS-DB3","CC-DB3","KJ-EB1","NP-DB1",
	"OZ-DB1","TM-DB3","YN-EB1")
biomes <- rep(NA,length(filenames))
## biomes from the data publication
sa <- c('HY','OT')
ct <- c('UR','AS','TM','NP','AO','OG','CC','OY','NB','OZ','KM','KK','KY','OS','KS','AU','GR','TB','SI')
wt <- c('KG','AI','KA','WK','IC','SD','KJ','AY','TN','YK')
st <- c('AM','YN','OW')
biomes[substr(filenames, 1, 2) %in% sa] <- 'sa'
biomes[substr(filenames, 1, 2) %in% ct] <- 'ct'
biomes[substr(filenames, 1, 2) %in% wt] <- 'wt'
biomes[substr(filenames, 1, 2) %in% st] <- 'st'
data.biomes <- data.frame(cluster=filenames,ecocode=biomes, stringsAsFactors = FALSE)

data.japan <- merge(data.japan,data.biomes,by='cluster')
table(data.japan$ecocode)

## to start rm sa
data.japan <- subset(data.japan,subset=data.japan$ecocode!='sa')


data.japan$perc.dead <- NA
fhui28's avatar
fhui28 committed
186
		
fhui28's avatar
fhui28 committed
187

Georges Kunstler's avatar
Georges Kunstler committed
188
189
data.japan <- subset(data.japan,subset=data.japan[["D"]]>=10 & !is.na(data.japan[["D"]])
                                       & !is.na(data.japan[["x"]]) & !is.na(data.japan[["y"]])) ### select only tree above 10cm of dbh at census 1
Georges Kunstler's avatar
Georges Kunstler committed
190

Georges Kunstler's avatar
Georges Kunstler committed
191
192
193
194
195
196
197
#### get wc climate 
source("R/utils/climate.R")
clim <- GetClimate(data.japan$Lat,data.japan$Lon)
data.japan$MAT <- clim$MAT
data.japan$MAP <- clim$MAP


Georges Kunstler's avatar
Georges Kunstler committed
198
###################### PLOT SELECTION FOR THE ANALYSIS
199
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
Georges Kunstler's avatar
Georges Kunstler committed
200
    'Lon','Lat',"x", "y", "census",'MAT','MAP')
Georges Kunstler's avatar
Georges Kunstler committed
201
202
203
204
data.tree <- subset(data.japan, select = c(vec.basic.var))

data.tree <- fun.convert.type.B(data.tree)
write.csv(data.tree,file="output/formatted/Japan/tree.csv",row.names = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
205

fhui28's avatar
fhui28 committed
206
### write data plot with variables only at the plot level. I should delete them from the tree table but so far I keep them to not distroy later code
Georges Kunstler's avatar
Georges Kunstler committed
207
208
209

vec.basic.var.p <- c("plot", "cluster", "Lon","Lat","ecocode")
data.plot <- subset(data.tree, subset=!duplicated(data.tree$cluster),select = c(vec.basic.var.p))
fhui28's avatar
fhui28 committed
210
write.csv(data.plot,file="output/formatted/Japan/plot.csv",row.names = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
211