-
Notifications
You must be signed in to change notification settings - Fork 0
/
Exploration_presentation_test.R
139 lines (124 loc) · 5.33 KB
/
Exploration_presentation_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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#' ---
#' title: "Exploration"
#' author: All
#' output: github_document
#' ---
#+ r load, message = FALSE, warnings = FALSE, echo = FALSE
library(ggplot2)
library(dplyr)
library(tidyr)
library(rmarkdown)
veg<-read.csv('00_raw_data/vegetation_sampling/saddptqd.hh.data.csv', header = T)
#' We decided to use the the NIWOT saddle grid data to answer our questions.
#' This grid was sampled intermittently since 1989 and annually or biannually
#' since 2006.
#' How often is each year represented?
#' We first wanted to know how sampling differed across the years in the data set.
#' Based on this information we excluded 1996 from our analysis
#'
#+ r table, echo = FALSE
table(veg$year)
#'we also noticed sampling design seemed to vary. In some years middle hits were
#'taken and in others they were not. In some years top hits were recorded for every
#'plot but not in others.
#'
#'Ultimately we decide to isolate only the top hits and covert bottom only hits to top
#'hits in the early years
#+r , echo = FALSE
veg.proc = veg %>% filter(!year %in% 1996)
veg.proc = veg.proc %>% filter(!grepl('^2', USDA_code))
veg.proc %>%
filter(USDA_code %in% c('KOMY', 'GEROT', 'DECE')) %>%
group_by(year, hit_type, USDA_code) %>%
summarise(n = n()) %>%
ggplot() +
geom_line(aes(x = year, y = n, group = hit_type, colour = hit_type)) +
geom_point(aes(x = year, y = n, colour = hit_type)) +
theme(legend.position = 'bottom') +
facet_wrap(~ USDA_code)
#'Although intially we were primarily interested in the Deschapsia-Geum interaction
#'we noticed that kobresia myoseroides was the second most common hit at decided to
#'include it in our analysis.
#'
#'We next wanted to visualise how each species was changing throughout time. Was
#'Deschampsia becoming much more abundant across the whole site?
#+ r species trends, echo = FALSE
veg.clean<-read.csv('01_process_data/output/veg_all_predictors.csv')
veg.clean%>%
ggplot( aes(x = year, y = n.obs, color = species))+
geom_point(aes(x =year, y = n.obs, color = species))+
geom_smooth(method = 'lm')
#'Within certain habitat type it seems like relative abundance of certain species changes
#+ r veg class trends, echo = FALSE
veg.clean<-read.csv('01_process_data/output/veg_all_predictors.csv')
veg.clean%>%
drop_na(veg_class)%>%
ggplot( aes(x = year, y = n.obs, color = veg_class))+
geom_point(aes(x =year, y = n.obs, color = veg_class))+
geom_smooth(method = 'lm')+
facet_wrap(vars(species))
#'
#' Elevation 1 with veg classes
#+ r elev1, echo = FALSE
read.csv('01_process_data/output/veg_all_predictors.csv') %>%
distinct(plot, year, .keep_all = TRUE) %>%
filter(!is.na(veg_class)) %>%
ggplot() +
geom_point(aes(x = UTM_E, y = UTM_N, shape = veg_class,
colour = elev1),
position = position_jitter(height = 5, width = 5),
size = 6) +
theme(panel.background = element_rect(fill = 'black'),
panel.grid = element_blank(),
plot.background = element_rect(fill = 'black'),
axis.text = element_text(colour = 'white'),
axis.title = element_text(colour = 'white'),
legend.background = element_rect(fill = 'black'),
legend.text = element_text(colour = 'white'),
legend.title = element_text(colour = 'white'),
legend.key = element_rect(fill = 'gray55')) +
scale_color_gradient(low = 'gray22', high = 'white') +
scale_shape_manual(values = 1:7)
#'Snow depth: early models showed that that there was a lot of variation in space,
#'with somewhat less in time.The only explanatory variable we had that varied in
#'both space and time was snow depth. Snow depth was measured at stakes associated
#'with each plot starting in March on a biweekly (mostly) schedule, producing a
#'mountain of data.
#'
#'We felt the most biologically relavant data to be max depth and snow melt date.
#'Due to gaps in the data it was hard to estimate snow melt date reliably. Scott
#'tried valiantly to fit GAMs to the snow data to allow us to estimate snow melts.
#'But after several weeks and many discussions the GAMs grew too frustrating and we
#'opted to use mean snow depth in June - simpler metric that still included some of
#'that information.
#'
#'We were able use temperature data at the site level from the saddle data logger.
#'This provided us with a wealth of data at a fine temporal scale.
#+ r daily temp, echo = FALSE
daily = read.csv('01_process_data/output/daily_airtemp_all.csv')
daily = daily %>%
mutate(jd = paste('1970', month, day, sep = '-') %>% as.Date() %>% julian(),
wyear = year + as.numeric(month > 9),
wd = jd - 273 + ifelse(wyear == year, 365, 0))
head(daily)
daily %>%
ggplot() +
geom_line(aes(x = wd, y = avg_temp)) +
facet_wrap(~ wyear)
#'
#'However, we decided to use the summer mean daily maximum from June - August as our #'predictor
#'
#'We also included a spatial autocorrelation term to incorporate the effect of the
#'abuncance of each speces in the previous year weighted by distance to the plot.
#+ r spatial autocorrelation, echo = FALSE
veg.2010=veg.clean%>%
filter(year %in% 2010)%>%
filter(species %in% 'DECE')
coords<- cbind(veg.2010$UTM_E, veg.2010$UTM_N)
matrix <- as.matrix(veg.2010$UTM_E, veg.2010$UTM_N, veg.2010$n.obs)
colnames(coords) <- c("X", "Y")
rownames(coords)<-veg.2010$plot
distmat <- as.matrix(dist(coords))
idw<-1/distmat
idw[!is.finite(idw)]<-NA
idw