-
Notifications
You must be signed in to change notification settings - Fork 15
/
maps.R
92 lines (81 loc) · 2.49 KB
/
maps.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
# *******************************************
# This file contains code for making maps of
# the WHO super-regions and regions based on
# clustering by PM2.5.
# *******************************************
# Setup -------------------------------------------------------------------
library(sp)
library(spdep)
library(mapmisc)
library(maptools)
library(rworldmap)
library(scales)
library(RColorBrewer)
library(dplyr)
# load 'GM' SpatialPointsDataFrame
load("bayes-vis.RData")
# Map of WHO super-regions ----------------------------------------------
spdf <- joinCountryData2Map(GM@data, nameJoinColumn="iso3")
# code missings
spdf$super_region[is.na(spdf$super_region) & spdf$continent == "Eurasia"] <- 4
spdf$super_region[is.na(spdf$super_region) & spdf$continent == "South America"] <- 5
spdf$super_region[is.na(spdf$super_region) & spdf$continent == "Africa"] <- 7
spdf$super_region[spdf$NAME %in% c("Algeria", "Libya", "Yemen", "Syria", "Iraq")] <- 2
spdf$super_region[spdf$NAME %in% c("Taiwan", "Vietnam", "Cambodia", "Laos", "Papua New Guinea", "N. Korea", "Sri Lanka")] <- 6
spdf$super_region[spdf$NAME %in% "Nepal"] <- 3
spdf$super_region[spdf$NAME %in% c("Greenland", "W. Sahara")] <- NA
opar <- par()
png(filename = "plots/map-who.png", res = 500, width = 1000, height = 600)
par(mar=c(0,0,0,0))
mapCountryData(
spdf,
nameColumnToPlot = "super_region",
catMethod = "categorical",
colourPalette = rev(
scales::hue_pal(
h = c(0, 360) + 15,
c = 100,
l = 65,
h.start = 0,
direction = 1
)(7)
),
addLegend = FALSE,
mapTitle = "",
lwd = 0.25
)
dev.off()
par(opar)
# Map of regions obtained from clustering ---------------------------------
average <-
GM@data %>%
group_by(iso3) %>%
summarise(pm25 = mean(pm25))
d <- dist(average)
hh <- hclust(d)
clust <- cutree(hh,k = 6)
GM@data$cluster_region <-
sapply(GM@data$iso3, function(x) clust[which(average$iso3 == x)])
spdf2 <- joinCountryData2Map(GM@data, nameJoinColumn="iso3")
spdf2$cluster_region[is.na(spdf2$cluster_region) & !is.na(spdf2$ISO3) & spdf$ISO3 != "ATA" & spdf2$ISO3 != "GRL"] <- 7
png(filename = "plots/map-clusters.png", res = 500, width = 1000, height = 600)
par(mar=c(0,0,0,0))
maplist <- mapCountryData(
spdf2,
nameColumnToPlot = "cluster_region",
catMethod = "fixedWidth",
colourPalette = rev(
scales::hue_pal(
h = c(0, 360) + 15,
c = 100,
l = 65,
h.start = 0,
direction = 1
)(7)
),
addLegend = FALSE,
mapTitle = "",
lwd = 0.25
)
dev.off()
graphics.off()