-
Notifications
You must be signed in to change notification settings - Fork 0
/
example2.Rmd
162 lines (121 loc) · 6.87 KB
/
example2.Rmd
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
---
title: "2. Generating EU NUTS maps"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Generating EU NUTS maps}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## Introduction
For some projects we want to create a map of EU NUTS level results produced by GLOBIOM. This vignette provides several examples on how to do this. The EU NUTS map is part of the package and can be loaded using buildin commands (see below). The package also includes data on potentially disappeared fraction of global species (PDF/ha) for three scenarios (2010-2050) at the EU Nuts region (`eu_nuts_pdf`) from [Di Fulvio et al. (2019)](https://www.sciencedirect.com/science/article/pii/S0048969718333941?via%3Dihub), which will be used in the example.
## Preparing the data
We use the [sf package](https://cran.r-project.org/web/packages/sf/index.html) for most of the vector (e.g. polygon) processing so it needs to be loaded first. The package can be easily combined with [ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html) the core plotting package in R.
Before we can create the map we need to: (1) load the EU NUTS polygons (part of the package); (2) load the example data (part of the package) and (3) connect the the two.
```{r message=FALSE, warning=FALSE, include=TRUE}
library(globiomvis)
library(ggplot2)
library(sf)
library(dplyr)
# Load the EU NUTS polygons
eu_nuts_poly <- st_read(system.file("shp", "NUTS2_revised.shp", package = "globiomvis"))
# Load the test data from the into the global R environment (i.e. it will display in RStudio so you can see it)
data("eu_nuts_pdf", package="globiomvis")
# Join the data. The variable names in the polygon file and the data file are the same (NUTS2) so they are automatically linked. Note that the order is important because otherwise the polygon format will be lost.
eu_nuts <- left_join(eu_nuts_poly, eu_nuts_pdf)
```
## Create a map for one year and one scenario
```{r message=FALSE, warning=FALSE, include=TRUE}
# Prepare data. We select only the year 2050 and the BASE scenario.
# We remove NA values, if any.
df <- eu_nuts %>%
filter(year == 2050, scenario == "BASE", !is.na(value))
# Create the map. Note that we plot the map twice (2xgeom_sf). The first call ensures that the full map is depicted (with background color white). As we used a filter, it is possible we removed polygons for which data was not available, which would otherwise be dropped.
ggplot() +
geom_sf(data = eu_nuts_poly, colour = "grey30", fill = "white") +
geom_sf(data = df, colour = "grey30", aes(fill = value)) +
scale_fill_gradientn(colours = terrain.colors(10)) +
labs(fill = "", x = "", y = "", title = "PDF of global species for BASE scenario, 2050") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "transparent"),
plot.title = element_text(hjust = 0.5))
```
Like in the original figure 2, we add bin the data to better reveal the spatial pattern. I also use a similar color palette as in the original figure.
```{r message=FALSE, warning=FALSE, include=TRUE}
# Create bins
interval <- c(0, 1.5, 3, 6, 9, 12, 15, 20, 30, 60, 100, 110)
# Create legend labels
lab <- paste(interval, interval[-1], sep = "-")[c(1:11)]
# As the palette has only 9 colors and we need 12, we need to interpolate additional colors.
library(RColorBrewer)
n_col = length(interval)
getPalette = colorRampPalette(brewer.pal(9, "YlOrBr"))
# Add the binned data
df <- df %>%
mutate(value2 = cut(value, interval, include.lowest = T))
# Create the map
ggplot() +
geom_sf(data = eu_nuts_poly, colour = "grey30", fill = "transparent") +
geom_sf(data = df, colour = "grey30", aes(fill = value2)) +
scale_fill_manual(values = getPalette(n_col), labels = lab) +
labs(x = "", y = "", fill = "PDF/ha" , title = "EU nuts abandoned land for 2050") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "transparent"),
plot.title = element_text(hjust = 0.5),
legend.position = c(0.9, 0.6),
legend.background = element_rect(fill = "transparent"),
legend.text=element_text(size=6))
```
## Create a set of maps for multiple periods and one scenario
In this example we create multiple maps that show the change in PDF/ha over time and in space. As we are increasing the number of maps, we first 'simplify' the polygons to speed up the processing time. To simplify further (and speed up), increase the `dTolerance` parameter. The EU NUTS polygon is quite detailed (i.e. capturing complex border patterns), which consume a lot of memory when plotting. This information is not needed for exploratory mapping.
```{r, message=FALSE, warning=FALSE, include=TRUE}
# Simplify the map
eu_nuts_poly_s <- st_simplify(eu_nuts_poly, dTolerance = 5000)
# Select the BASE scenario
df2 <- left_join(eu_nuts_poly_s, eu_nuts_pdf) %>%
filter(scenario == "BASE" , !is.na(value))
# plot
ggplot() +
geom_sf(data = eu_nuts_poly_s, colour = "grey30", fill = "red") +
geom_sf(data = df2, colour = "grey30", aes(fill = value)) +
scale_fill_gradientn(colours = terrain.colors(10)) +
labs(fill = "", x = "", y = "", title = "PDF/ha of global species for 2010-2050") +
facet_wrap(~year) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "transparent"),
plot.title = element_text(hjust = 0.5),
strip.background = element_rect(fill = "transparent"),
legend.position = "bottom")
```
## Create a map for multiple periods and multiple scenarios
Now we show all maps covering both the time and scenario dimension. As the maps are very small, this is probably only useful if there are large differences over time and space that can be easily spotted.
```{r, message=FALSE, warning=FALSE, include=TRUE}
# Remove NA values
df2 <- left_join(eu_nuts_poly_s, eu_nuts_pdf) %>%
filter(!is.na(value))
# plot
ggplot() +
geom_sf(data = eu_nuts_poly_s, colour = "grey30", fill = "red") +
geom_sf(data = df2, colour = "grey30", aes(fill = value)) +
scale_fill_gradientn(colours = terrain.colors(10)) +
labs(fill = "", x = "", y = "", title = "PDF/ha of global species for 2010-2050") +
facet_grid(scenario~ year) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "transparent"),
plot.title = element_text(hjust = 0.5),
strip.background = element_rect(fill = "transparent"),
legend.position = "bottom")
```