-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfunction_calls.R
89 lines (62 loc) · 2.21 KB
/
function_calls.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# calls for cso precip..
setwd('/Users/grad/Desktop/umces/cso/precip')
test_gage = getGHCN(lat = 39.9, lon = -77.18, NSdeg = 3, EWdeg = 3, StartYear = 1984, EndYear = 2015, StnOpt = "Full")
test = read_ghcn(test_gage, cut_date = "2015-12-31", write_csv = FALSE)
pcp_l = test$pcp
tl = lapply(pcp_l, make_df, test)
require(pbapply)
tl = pblapply(tl, remove_nas)
# Have a look
require(pbapply)
all_gauges = pblapply(tl, gauge_spatial)
# look at a gauge dist
## !! call
plot_gauges(all_gauges, item = 2000)
# !! end call
## !! call
require(pbapply)
vor_l = pblapply(all_gauges, thiessen_poly)
# !! call, define outside to save time ...
direcotry = "/Users/grad/Desktop/umces/cso/css"
layer_name = "CSS_final_jan14_clean"
cso = readOGR(dsn = directory, layer = layer_name)
## !! end call
# !! call
cso = readOGR(dsn = directory, layer = layer_name)
require(pbapply)
out = mclapply(vor_l, extract_pcp, cso = cso, write_csv = FALSE)
# !! end call
# !! call
output = finish_pcp(out, test$dates, write_csv = FALSE)
# cso conv function.
# define the TT cso equation ....
# outside the other functions, because this makes it CBP specific ....
cso_est = function(x){
mm_out = x / 10
inches = mm_out*0.03937
cso_est = 1567*inches^2 - 46.955*inches + 1309.8
#cso_est[which(x <= threshold)] <- 0
cso_est = cso_est / 1e6
return(cso_est)
}
#calculate cso output on outp
output$cso = cso_est(output$pcp) * output$ac
out1 = aggregate.data.frame(output$pcp, by = list(output$id), FUN = sum)
out2 = aggregate.data.frame(output$ac, by = list(output$id), FUN = sum)
out = merge(out1, out2, by = 'Group.1')
colnames(out) = c('id', 'pcp', 'ac')
outp = out[order(out$id),]
# write it separte,
write.csv(outp, paste(getwd(), 'output_pcp.csv', sep = '/'), row.names = FALSE)
# end call
##plot the min and max number of stations, and min and max rainfall
# number of stations
lens = lapply(tl, nrow)
mins = which(lens == min(unlist(lens)))
plot.min = mins[[length(mins)]]
date.min = dates[[plot.min]]
maxs = which(lens == max(unlist(lens)))
plot.max = maxs[[length(maxs)]]
date.max = dates[[plot.max]]
theme_set(theme_gray(base_size = 22))
plot_extract_pcp(vor = vor_l, dates = dates, cso = cso, date_lab = dates[[which(dates == date.min)]])