-
Notifications
You must be signed in to change notification settings - Fork 0
/
max_mk_test.R
113 lines (97 loc) · 2.95 KB
/
max_mk_test.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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# max mk test
# nov 24 2020
library(BiocManager)
library(rhdf5)
library(tidyverse)
library(rgeos)
library(raster)
library(rgdal)
library(parallel)
library(data.table)
library(gtools)
library(terra)
library(spatialEco)
library(EnvStats)
library(Kendall)
library(remotes)
# read in scf rasters to stack
setwd("/Volumes/jt/projects/margulis/snow_metric_rasters/max_swe/rasters")
list <-list.files()
max_list <-lapply(list, function(x) raster(x))
max_stack <-stack(max_list)
crs(max_stack)<-"+proj=leac +ellps=clrk66"
plot(max_stack[[9]])
max_stack
# test mk code by first running it on .5 degrees lat near tahoe
#tahoe_crop <-crop(scf_stack, extent(-123.3, -117.6, 39, 39.5))
#tahoe_crop
#plot(tahoe_crop[[9]])
###### mk test, trend.slope is full stats, 2 is just p-val and slope stats
trend.slope <- function(y, p.value.pass = TRUE, z.pass = TRUE,
tau.pass = TRUE, confidence.pass = TRUE, intercept.pass = TRUE) {
options(warn = -1)
fit <- EnvStats::kendallTrendTest(y ~ 1)
fit.results <- fit$estimate[2]
if (tau.pass == TRUE) {
fit.results <- c(fit.results, fit$estimate[1])
}
if (intercept.pass == TRUE) {
fit.results <- c(fit.results, fit$estimate[3])
}
if (p.value.pass == TRUE) {
fit.results <- c(fit.results, fit$p.value)
}
if (z.pass == TRUE) {
fit.results <- c(fit.results, fit$statistic)
}
if (confidence.pass == TRUE) {
ci <- unlist(fit$interval["limits"])
if (length(ci) == 2) {
fit.results <- c(fit.results, ci)
}
else {
fit.results <- c(fit.results, c(NA, NA))
}
}
options(warn = 0)
return(fit.results)
}
trend.slope2 <- function(y, p.value.pass = TRUE, z.pass = FALSE,
tau.pass = FALSE, confidence.pass = FALSE, intercept.pass = FALSE) {
options(warn = -1)
fit <- EnvStats::kendallTrendTest(y ~ 1)
fit.results <- fit$estimate[2]
if (tau.pass == TRUE) {
fit.results <- c(fit.results, fit$estimate[1])
}
if (intercept.pass == TRUE) {
fit.results <- c(fit.results, fit$estimate[3])
}
if (p.value.pass == TRUE) {
fit.results <- c(fit.results, fit$p.value)
}
if (z.pass == TRUE) {
fit.results <- c(fit.results, fit$statistic)
}
if (confidence.pass == TRUE) {
ci <- unlist(fit$interval["limits"])
if (length(ci) == 2) {
fit.results <- c(fit.results, ci)
}
else {
fit.results <- c(fit.results, c(NA, NA))
}
}
options(warn = 0)
return(fit.results)
}
# run it in parallel to see if stripping is gone
beginCluster(n=7)
system.time(max_trends_full <- clusterR(max_stack, overlay, args=list(fun=trend.slope2)))
endCluster()
plot(max_trends_full[[1]])
plot(max_trends_full[[2]])
max_p_value_full <-max_trends_full[[2]]
max_slope_full <-max_trends_full[[1]]
writeRaster(max_p_value_full,"/Volumes/jt/projects/margulis/snow_metric_rasters/max_swe/mk_results/max_p_value_full.tif")
writeRaster(max_slope_full, "/Volumes/jt/projects/margulis/snow_metric_rasters/max_swe/mk_results/max_slope_full.tif")