-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPopulation_Trend_per_y.R
102 lines (88 loc) · 2.97 KB
/
Population_Trend_per_y.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
# Library
library(raster)
library(maptools)
library(gam)
# Load data
load("PopD_all_December.rdata")
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)
names(per.origin) <- origins@data[, 1]
# Function standardization
std <- function(x) {
b <- (x - min(x)) / (max(x) - min(x))
return(rev(b))
}
# Calculating mean and
global.means <- global.gams <- list()
for (j in 1:length(per.origin)) {
print(j)
originI <- Pop[per.origin[[j]][, 1], ]
time <- 21:4
originI <- na.exclude(originI)
b <- apply(originI, 1, std)
nJ <- nrow(originI)
predictions <- matrix(nrow = nJ, ncol = length(time))
for(i in 1:nJ) {
model <- gam(b[, i] ~ s(time, df = 15))
col <- sample(rainbow(100), 1)
predictions[i, ] <- predict(model)
}
global.means[[j]] <- apply(predictions, 2, mean)
global.gams[[j]] <- apply(predictions, 2, sd)
}
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Plot per region
tiff("all_origins_pop_trend.tif", width = 50, height = 40, res = 300,
units = 'cm')
par(mfrow = c(4, 5))
for (j in 1:length(global.means)) {
plot(seq(0, 1, length.out = length(time)) ~ time, col = "white", main = names(per.origin)[j],
xlim = c(21, 4), ylab = "Pop Density",
xlab = "Thousand of years ago")
x <- global.means[[j]]
y <- global.gams[[j]]
down <- x - y
up <- x + y
lines(y = down, x = time, lty = 3, col = "gray40", lwd = 1)
lines(y = up, x = time, lty = 3, col = "gray40", lwd = 1)
lines(y = x, x = time, lwd = 2)
if (origin.time.region[j] == 1) {
polygon(cbind(c(12, 8.2, 8.2, 12, 12), c(-1, -1, 2, 2, -1)),
col = rgb(0, 1, 0, alpha = .2), border = F)
}
if (origin.time.region[j] == 2) {
polygon(cbind(c(8.2, 4.2, 4.2, 8.2, 8.2), c(-1, -1, 2, 2, -1)),
col = rgb(.28, 0, .28, alpha = .2), border = F)
}
}
dev.off()
# Save it in pdf
pdf("all_origins_pop_trend.pdf", width = 40, height = 50)
par(mfrow = c(5, 4), mar = c(5, 7, 7, 5))
for (j in 1:length(global.means)) {
plot(seq(0, 1, length.out = length(time)) ~ time, col = "white", main = names(per.origin)[j],
xlim = c(21, 4), ylab = "Pop Density",
xlab = "Thousand of years ago", cex.lab = 3, cex.main = 4, cex.axis = 2)
x <- global.means[[j]]
y <- global.gams[[j]]
down <- x - y
up <- x + y
lines(y = down, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = up, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = x, x = time, lwd = 4)
if (origin.time.region[j] == 1) {
polygon(cbind(c(12, 8.2, 8.2, 12, 12), c(-1, -1, 2, 2, -1)),
col = rgb(0, 1, 0, alpha = .2), border = F)
}
if (origin.time.region[j] == 2) {
polygon(cbind(c(8.2, 4.2, 4.2, 8.2, 8.2), c(-1, -1, 2, 2, -1)),
col = rgb(.28, 0, .28, alpha = .2), border = F)
}
}
dev.off()