-
Notifications
You must be signed in to change notification settings - Fork 0
/
fig4_script.R
145 lines (129 loc) · 7.94 KB
/
fig4_script.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
139
140
141
142
143
144
145
#This script is for making the four panel choropleth maps shwoing country level PA effectiveness as well as other stats at the country level
#-------prepare the data by joing the Atotal AGB values and the ancillary information by country ---------------------
totalAGB_anc <- read.csv("/gpfs/data1/duncansongp/GEDI_global_PA/WDPA_L4b_output/l4bAGB_countryxBiome_byCountry_mar08_ancillary_v2.csv")
totalAGB <- read.csv("/gpfs/data1/duncansongp/GEDI_global_PA/WDPA_L4b_output/l4bAGB_countryxBiome_byCountry_mar08_usa.csv") %>%
dplyr::select(totalAGB , AGB_stderr , AGB_type, ISO3, continent, iso3Status) %>%
pivot_wider(names_from=AGB_type, values_from=c(totalAGB,AGB_stderr)) %>%
left_join(totalAGB_anc, by=c("ISO3"="iso3"))
#-----plotting the choropleth maps-----
library(rnaturalearth)
library(BAMMtools)
library(sf)
# all2_mod <- all2 %>% mutate(extra_AGB_in_PA=as.numeric(extra_AGB_in_PA)) %>% mutate(totalPAArea_all=totalPAArea_all/1000000)
all2_mod <- totalAGB
sf_use_s2(FALSE)
world <-ne_countries(scale = "medium", returnclass = "sf")
class(world)
dim(world)
world2=world %>% left_join(all2_mod, by=c("iso_a3"="ISO3")) %>%
dplyr::mutate(totalAGB_allPAExtraAGB= totalAGB_allPAExtraAGB*0.49) %>%
dplyr::mutate(totalAGB_allPAAGB= totalAGB_allPAAGB*0.49)
# %>%
# mutate(meanAGBD_PA= meanagbd_1,meanAGBD_control=meanagbd_0, sdAGBD_PA=sdagbd_1, sdAGBD_control=sdagbd_0)%>%
# dplyr::mutate(absolute_diff_AGBD=meanAGBD_PA-meanAGBD_control) %>%
# dplyr::mutate(percent_diff_AGBD=100*(meanAGBD_PA-meanAGBD_control)/meanAGBD_control) %>%
# dplyr::mutate(absolute_diff_rh98=(meanrh98_1-meanrh98_0)/100) %>%
# dplyr::mutate(percent_diff_rh98=100*(meanrh98_1-meanrh98_0)/meanrh98_0) %>%
# dplyr::mutate(absolute_diff_cover=meancov_1-meancov_0) %>%
# dplyr::mutate(percent_diff_cover=100*(meancov_1-meancov_0)/meancov_0) %>%
# dplyr::mutate(absolute_diff_pai=meanpai_1-meanpai_0) %>%
# dplyr::mutate(percent_diff_pai=100*(meanpai_1-meanpai_0)/meanpai_0)
# br <- append(getJenksBreaks(world2$totalAGB_allPAExtraAGB,7), 1.5*2) %>% unique() %>% sort() %>% round(digits = 1)
br <- append(getJenksBreaks(world2$totalAGB_allPAExtraAGB,7), 1.5*0.49) %>% unique() %>% sort() %>% round(digits = 1) %>% append(0.05) %>% sort() %>% unique()
t1=ggplot(data = world2) +
geom_sf(aes(fill = cut(totalAGB_allPAExtraAGB, breaks=br, include.lowest=TRUE, right=FALSE)), size=0.4, colour="#494a4a") +
coord_sf(crs = "+proj=cea +lon_0=Central Meridian +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84")+
# coord_sf(crs = "+proj=moll")+
scale_fill_manual(values=c("#dfc27d", "#eef5dc","#d9f0a3","#addd8e", "#41ab5d", "#1b7038", "#004529"), na.translate = F) +
# scale_colour_manual(values=NA) +
# guides(colour=guide_legend("No PAs/No matches", override.aes=list(colour="black", fill="white"))) +
labs(title = "Additional total preserved C in PAs",
fill = "Preserved C (Gt) ")+
theme_bw() +
theme(text=element_text(family="Helvetica", face="bold", size=12),
legend.title=element_text(size=12),
legend.text=element_text( size=12), axis.text=element_text(size=12)) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())+
theme(axis.ticks.x = element_blank()) +
theme(legend.key.size = unit(0.5, 'cm'))+
theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))+
labs(tag = "A")
t11=t1+coord_sf(xlim = c(-179, 180), ylim = c(-55, 55), expand = FALSE)
t11
# br2 <- getJenksBreaks(world2$totalAGB_allPAAGB,8) %>% unique() %>% sort() %>% round(digits = 1)
br2 <- getJenksBreaks(world2$totalAGB_allPAAGB,8) %>% unique() %>% sort() %>% round(digits = 1)
br2[8] <- 19.4
t2=ggplot(data = world2) +
geom_sf(aes(fill = cut(totalAGB_allPAAGB, breaks=br2, include.lowest=TRUE, right=FALSE)), size=0.4, colour="#494a4a") +
coord_sf(crs = "+proj=cea +lon_0=Central Meridian +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84")+
# coord_sf(crs = "+proj=moll")+
scale_fill_manual(values=c("#ffffcc", "#c7e9b4","#7fcdbb", "#41b6c4", "#1d91c0",'#225ea8','#0c2c84'),na.translate = F) +
# scale_colour_manual(values=NA) +
# guides(colour=guide_legend("No PAs/No matches", override.aes=list(colour="black", fill="white"))) +
labs(title = "Total aboveground C in PAs",
fill = str_wrap(" Total aboveground C in PAs (Gt)", width = 20))+
theme_bw() +
theme(text=element_text(family="Helvetica", face="bold", size=12),
legend.title=element_text(size=12),
legend.text=element_text( size=12), axis.text=element_text(size=12)) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())+
theme(axis.ticks.x = element_blank()) +
theme(legend.key.size = unit(0.5, 'cm'))+
theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))+
labs(tag = "C")
t22=t2+coord_sf(xlim = c(-179, 180), ylim = c(-55, 55), expand = FALSE)
t22
br3 <-append(seq(-65, 75, len=7), 0) %>% sort() %>% round(digits = 1) #append(getJenksBreaks(world2$weighted_agbd_diff,7), 0) %>% unique() %>% sort() %>% round(digits = 2)
t3=ggplot(data = world2) +
geom_sf(aes(fill = cut(weighted_agbd_diff, breaks=br3, include.lowest=TRUE, right=FALSE)), size=0.4, colour="#494a4a") +
coord_sf(crs = "+proj=cea +lon_0=Central Meridian +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84")+
scale_fill_manual(values=c("#bf812d", "#dfc27d", "#f6e8c3","#c7eae5", "#80cdc1", "#35978f", "#0a4a45"),na.translate = F) +
# scale_colour_manual(values=NA) +
# guides(colour=guide_legend("No PAs/No matches", override.aes=list(colour="black", fill="white"))) +
labs(title = "Mean additional preserved AGBD in PAs",
fill = str_wrap("Mean additional AGBD (Mg/ha)", width = 20))+
theme_bw() +
theme(text=element_text(family="Helvetica", face="bold", size=12),
legend.title=element_text(size=12),
legend.text=element_text( size=12),
axis.text=element_text(size=12)) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())+
theme(axis.ticks.x = element_blank()) +
theme(legend.key.size = unit(0.5, 'cm'))+
theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))+
labs(tag = "B")
t33=t3+coord_sf(xlim = c(-179, 180), ylim = c(-55, 55), expand = FALSE)
t33
world2[world2$adm0_a3=="BRA",]$total_pa_area <- 2.52*1000000
br4 <- getJenksBreaks(world2$total_pa_area/1000000,8) %>%round(digits = 2)
br4[8] <- 2.52
t4=ggplot(data = world2) +
geom_sf(aes(fill = cut(total_pa_area/1000000, breaks=br4, include.lowest=TRUE, right=FALSE), colour=""), size=0.4)+
coord_sf(crs = "+proj=cea +lon_0=Central Meridian +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84")+
# coord_sf(crs = "+proj=moll")+
scale_fill_manual(values=c("#fffad6", "#fff7bc","#fee391", "#fec44f", "#fe9929", "#ec7014", "#cc4c02", "#993404","#662506"), na.translate = F) +
scale_colour_manual(values=NA) +
guides(colour=guide_legend("No PAs/No matches", override.aes=list(colour="black", fill="#f0f0f0"))) +
labs(title = "Total PA Area",
fill = "PA area (Million ha) ")+
theme_bw() +
theme(text=element_text(family="Helvetica", face="bold", size=12),
legend.title=element_text(size=12),
legend.text=element_text( size=12),
axis.text=element_text(size=12)) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())+
theme(axis.ticks.x = element_blank()) +
theme(legend.key.size = unit(0.5, 'cm'))+
theme(legend.spacing.y = unit(0.1, "cm"))+
theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))+
labs(tag = "D")
t44=t4+coord_sf(xlim = c(-179, 180), ylim = c(-55, 55), expand = FALSE)
# scale_y_continuous(breaks = c(-50, -25, 0, 25,50))
t44
ggplot2::ggsave(filename=paste('/gpfs/data1/duncansongp/GEDI_global_PA/figures/JAN21_FIGS/iso_level_biomass_results_fig4_July29_C_v3.png', sep=''),
plot= grid.arrange(t11, t33,t22, t44, ncol=1),
width=8, height=10, units = "in", bg = "transparent")