-
Notifications
You must be signed in to change notification settings - Fork 7
/
Figure_2.r
72 lines (60 loc) · 2.19 KB
/
Figure_2.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
# Note: this relies on old libraries YMMV
# this is a simplified version of the figure
# in the final publication - for reference only
library(mapdata)
library(maptools)
library(maps)
library(raster)
library(RColorBrewer)
#---- set colours ----
red = colorRampPalette(brewer.pal(9,"YlOrRd"))(40)
blue = colorRampPalette(brewer.pal(9,"PuBu"))(20)
jet.col.arid = rev(c(rev(red),blue))
ylgn = colorRampPalette(brewer.pal(5,"YlGn"))(54)
ylorbr = colorRampPalette(brewer.pal(5,"YlOrBr"))(11)
ylgn = ylgn[c(-1:-4)]
ylorbr = ylorbr[-1]
jet.col.fcover = c(rev(ylorbr),ylgn)
#---- load data ----
fcover = raster("data/NCC_map_data/figure_2_fCover_change.tif")
aridity = raster("data/NCC_map_data/figure_2_aridity_index_change.tif")
#---- plot data ----
# margins etc
par(mfrow=c(2,1),oma=c(2,0,0,0),mar=c(4,4,4,4))
# plot
plot(fcover,box=F,axes=F,legend.only=F,zlim=c(-20,100),add=F,
ylim=c(22,52.875),
col=jet.col.fcover,
horiz = T,
lwd=1,
cex.lab=1.5,
legend.width=0.5, legend.shrink=0.9,
axis.args=list(cex.axis=1.5),
legend.args=list(
text='Fractional cover change (%)',
side=1,
font=1,
line=3,
cex=1.5)
)
# country outlines
map("worldHires","usa",resolution=0,add=T,lwd=1.1,xlim=c(-130,-70),ylim=c(22,52.875))
map("worldHires","canada",resolution=0,add=T,lwd=1.1,xlim=c(-130,-70),ylim=c(22,52.875))
map("worldHires","mexico",resolution=0,add=T,lwd=1.1,xlim=c(-130,-70),ylim=c(22,52.875))
legend('bottomright',legend='a',bty='n',cex=2.5)
plot(aridity,box=F,axes=F,legend.only=F,
zlim=c(-0.1,0.2),
add=F,
ylim=c(22,52.875),
col=jet.col.arid,
horiz = T,
lwd=1,
cex.lab=1.5,
legend.width=0.5, legend.shrink=0.9,
axis.args=list(cex.axis=1.5),
legend.args=list(text=expression(Delta*" Aridity Index"), side=1, font=1, line=3, cex=1.5))
# country outlines
map("worldHires","usa",resolution=0,add=T,lwd=1.3,xlim=c(-130,-70),ylim=c(25.125,52.875))
map("worldHires","canada",resolution=0,add=T,lwd=1.3,xlim=c(-130,-70),ylim=c(25.125,52.875))
map("worldHires","mexico",resolution=0,add=T,lwd=1.3,xlim=c(-130,-70),ylim=c(25.125,52.875))
legend('bottomright',legend='b',bty='n',cex=2.5)