-
Notifications
You must be signed in to change notification settings - Fork 1
/
dswe_hist.R
80 lines (67 loc) · 3.23 KB
/
dswe_hist.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
library(terra)
# vg extent
vg <-vect("/Users/jacktarricone/ch1_jemez_data/vector_data/valle_grande_aoi.geojson")
# cumlative swe change
cum_r <-rast("/Users/jacktarricone/ch1_jemez_data/gpr_rasters_ryan/final_swe_change/dswe_cum.tif")
cum <-crop(cum_r,vg)
plot(cum)
# pair 1
feb12_19_r <-rast("/Users/jacktarricone/ch1_jemez_data/gpr_rasters_ryan/final_swe_change/dswe_feb12-19_same_pixels.tif")
feb12_19 <-crop(feb12_19_r, vg)
plot(feb12_19)
# pair 2
feb19_26_r <-rast("/Users/jacktarricone/ch1_jemez_data/gpr_rasters_ryan/final_swe_change/dswe_feb19-26_same_pixels.tif")
feb19_26 <-crop(feb19_26_r, vg)
plot(feb19_26)
# convert to data frames for plotting
p1_df <-as.data.frame(feb12_19, xy = TRUE, cells = TRUE, na.rm = TRUE)
p2_df <-as.data.frame(feb19_26, xy = TRUE, cells = TRUE, na.rm = TRUE)
cum_df <-as.data.frame(cum, xy = TRUE, cells = TRUE, na.rm = TRUE)
# plot
# theme_set(theme_classic(base_size = 11))
ggplot()+
geom_vline(xintercept = 0, linetype=3, col = "black") +
geom_histogram(cum_df, mapping = aes(x=lyr.1, fill = "Feb 12-26th", color = "Feb 12-26th"), alpha=0.4, bins = 80) +
geom_histogram(p1_df, mapping = aes(x=lyr.1, fill = "Feb 12-19th", color = "Feb 12-19th"), alpha=0.4, bins = 80) +
geom_histogram(p2_df, mapping = aes(x=lyr.1, fill = "Feb 19-26th", color = "Feb 19-26th"), alpha=0.4, bins = 80) +
scale_colour_manual(name = "Date Range",
labels = c("Feb 12-19th","Feb 19-26th","Feb 12-26th"),
values=c("firebrick","goldenrod","black"))+
scale_fill_manual(name = "Date Range",
labels = c("Feb 12-19th","Feb 19-26th","Feb 12-26th"),
values = c("firebrick","goldenrod","black"))+
xlim(c(-2,1)) + xlab("dSWE [cm]")
# density plot
fancy_scientific <- function(l) {
# turn in to character string in scientific notation
l <- format(l, scientific = TRUE)
# quote the part before the exponent to keep all the digits
l <- gsub("^(.*)e", "'\\1'e", l)
# turn the 'e+' into plotmath format
l <- gsub("e", "%*%10^", l)
# return this as an expression
parse(text=l)
}
ggplot()+
geom_vline(xintercept = 0, linetype=3, col = "black") +
geom_density(cum_df, mapping = aes(x=lyr.1, y=stat(count),fill = "Feb 12-26th", color = "Feb 12-26th"), alpha=0.1) +
geom_density(p1_df, mapping = aes(x=lyr.1, y=stat(count),fill = "Feb 12-19th", color = "Feb 12-19th"), alpha=0.1) +
geom_density(p2_df, mapping = aes(x=lyr.1, y=stat(count),fill = "Feb 19-26th", color = "Feb 19-26th"), alpha=0.1) +
scale_colour_manual(name = "InSAR Dates",
labels = c("Feb 12-19th","Feb 19-26th","Feb 12-26th"),
values=c("firebrick","goldenrod","black"))+
scale_fill_manual(name = "InSAR Dates",
labels = c("Feb 12-19th","Feb 19-26th","Feb 12-26th"),
values = c("firebrick","goldenrod","black"))+
xlim(c(-2,1)) + xlab("dSWE [cm]")+ylab("Count")+
scale_y_continuous(labels=fancy_scientific) +
theme(legend.position = c(.22,.75))
setwd("/Users/jacktarricone/ch1_jemez_data/plots")
ggsave("dswe_hist.png",
width = 5,
height = 3,
units = "in",
dpi = 300)
# global(cum, mean, na.rm = TRUE)
# global(feb12_19, mean, na.rm = TRUE)
# global(feb19_26, mean, na.rm = TRUE)