-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlatitude_analysis_2
151 lines (112 loc) · 5.42 KB
/
latitude_analysis_2
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
library(raster)
library(xts)
library(lubridate)
library(ggstatsplot)
library(rsq)
library(extrafont)
library(car)
library(bfast)
library(remotes)
library(ggpmisc)
library(dplyr)
library(lubridate)
library(mice)
library(VIM)
setwd("F:\\Global_water_sink_source\\Figures\\A_B")
multiband_tif <- stack("diff_indirect_direct_ratio.tif")
##############
multi_frame <-
as.data.frame(multiband_tif, xy = TRUE) %>%
na.omit() %>% mutate(across(c(x, y), round, digits = 0))
multi_long = multi_frame %>% reshape2::melt(id=c("x", "y"))
statis_plot = multi_long %>%
group_by(y,variable) %>%
summarise(median = median(value, na.rm = TRUE),
Q75 = quantile(value,probs=0.75, na.rm = TRUE),
Q25 = quantile(value,probs=0.25, na.rm = TRUE) )
ggplot(statis_plot, aes(x=y, y=median))+
geom_ribbon(aes(ymin=Q25, ymax=Q75),alpha=0.3)+
coord_flip()+geom_line(size = 0.8)+geom_hline(yintercept=0, linetype="dashed", color = "black")+
labs(x="Lontitude", y = "Sensitivity") +scale_y_continuous(breaks = seq(-0.8,0.8,0.4))+
theme(axis.title=element_text(size=20,face="bold",color="black"),
axis.text = element_text(size=18,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
text=element_text(size=18),
legend.position="bottom",
legend.text = element_text(size=18
),legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#########################################################################################compare the sensitivity and the trend due to KNDVI change from 1982~2000
#first for the sensitivity
setwd("F:\\Global_water_sink_source\\Data\\scale_10000")
sensi_1982 <- stack("coeff_model1_kndvi_EI.tif")
sensi_2000 <- stack("coeff_model2_kndvi_EI.tif")
##############
multi_frame <-
as.data.frame(sensi_1982 , xy = TRUE) %>%
na.omit() %>% mutate(across(c(x, y), round, digits = 3))
multi_frame$sensi2000= extract(sensi_2000 ,multi_frame[,1:2])
colnames(multi_frame) = c("x","y","sensi1982","sensi2000")
multi_long = multi_frame %>%mutate(across(c(x, y), round, digits = 0)) %>% reshape2::melt(id=c("x", "y"))
statis_plot = multi_long %>%
group_by(variable) %>%
summarise(median = median(value, na.rm = TRUE),
Q75 = quantile(value,probs=0.75, na.rm = TRUE),
Q25 = quantile(value,probs=0.25, na.rm = TRUE) )
ggplot(statis_plot, aes(x=y, y=median,color = variable))+
geom_line(size = 0.8)+geom_hline(yintercept=0, linetype="dashed", color = "black")+
labs(x="Lontitude", y = "Sensitivity") +scale_y_continuous(breaks = seq(-0.8,0.8,0.1))+
theme(axis.title=element_text(size=20,face="bold",color="black"),
axis.text = element_text(size=18,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
text=element_text(size=18),
legend.position="bottom",
legend.text = element_text(size=18
),legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#first for the trend change
setwd("F:\\Global_water_sink_source\\Data\\ratio_change_value")
sensi_1982 <- stack("trend_ndvi_EI_1982_2000.tif")
sensi_2000 <- stack("trend_ndvi_EI_2001_2019.tif")
##############
multi_frame <-
as.data.frame(sensi_1982 , xy = TRUE) %>%
na.omit() %>% mutate(across(c(x, y), round, digits = 3))
multi_frame$sensi2000= extract(sensi_2000 ,multi_frame[,1:2])
colnames(multi_frame) = c("x","y","sensi1982","sensi2000")
multi_long = multi_frame %>%mutate(across(c(x, y), round, digits = 0)) %>% reshape2::melt(id=c("x", "y"))
statis_plot = multi_long %>%
group_by(variable) %>%
summarise(median = median(value, na.rm = TRUE),
Q75 = quantile(value,probs=0.75, na.rm = TRUE),
Q25 = quantile(value,probs=0.25, na.rm = TRUE) )
ggplot(statis_plot, aes(x=y, y=median,color = variable))+
geom_line(size = 0.8)+geom_hline(yintercept=0, linetype="dashed", color = "black")+
labs(x="Lontitude", y = "Sensitivity") +scale_y_continuous(breaks = seq(-2,2,0.5))+scale_x_continuous(breaks = seq(-2,2,0.5))+
theme(axis.title=element_text(size=20,face="bold",color="black"),
axis.text = element_text(size=18,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
text=element_text(size=18),
legend.position="bottom",
legend.text = element_text(size=18
),legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
###################################################################################################
setwd("F:\\Global_water_sink_source\\Data\\scale_10000")
sensi_1982 <- stack("trend_gimms_ndvi_1982_2000.tif")
sensi_2000 <- stack("trend_gimms_ndvi_2000_2019.tif")
##############
multi_frame <-
as.data.frame(sensi_1982 , xy = TRUE) %>%
na.omit() %>% mutate(across(c(x, y), round, digits = 3))
multi_frame$sensi2000= extract(sensi_2000 ,multi_frame[,1:2])
colnames(multi_frame) = c("x","y","sensi1982","sensi2000")
multi_long = multi_frame %>%mutate(across(c(x, y), round, digits = 0)) %>% reshape2::melt(id=c("x", "y"))
statis_plot = multi_long %>%
group_by(variable) %>%
summarise(median = median(value, na.rm = TRUE),
Q75 = quantile(value,probs=0.75, na.rm = TRUE),
Q25 = quantile(value,probs=0.25, na.rm = TRUE) )