Newer
Older
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
#' Plot global distribution of lpjml simulated biomes
#'
#' Plots a map with the biome distribution as derived from a lpjml run based
#' on the "classify_biomes" function
#'
#' @param biome_data output (list) from classify_biomes()
#'
#' @param file_name directory for saving the plot (character string)
#' @param display_area boolean, adding occupied area per biome (default FALSE)
#' @param to_robinson logical to define if robinson projection should be used
#' for plotting
#' @param order_legend in which order the biomes should be displayed
#' default: c(1,2,9,10,11,3,4,5,12,13,14,6,7,8,15,16,17,18,19)
#'
#' @param bg_col character, specify background possible (`NA` for transparent)
#'
#' @examples
#' \dontrun{
#' plot_biomes(biome_data = biomes,
#' file_name ="/p/projects/open/Johanna/R/biomes.pfd")
#' }
#'
#' @md
#' @export
plot_biomes <- function(biome_data,
file_name = NULL,
display_area = FALSE,
to_robinson = TRUE,
cellarea = NULL,
order_legend = c(1, 2, 9, 10, 11, 3, 4, 5, 12, 13, 14,
6, 7, 8, 15, 16, 17, 18, 19),
bg_col = "white") {
# load required data: bbox, countries
lpjml_extent <- c(-180, 180, -60, 85)
bounding_box <- system.file("extdata", "ne_110m_wgs84_bounding_box.shp",
package = "biospheremetrics") %>%
rgdal::readOGR(layer = "ne_110m_wgs84_bounding_box", verbose = FALSE) %>%
{ if(to_robinson) sp::spTransform(., sp::CRS("+proj=robin")) else . } # nolint
countries <- system.file("extdata", "ne_110m_admin_0_countries.shp",
package = "biospheremetrics") %>%
rgdal::readOGR(layer = "ne_110m_admin_0_countries", verbose = FALSE) %>%
raster::crop(., lpjml_extent) %>%
{ if(to_robinson) sp::spTransform(., sp::CRS("+proj=robin")) else . } # nolint
biome_cols <- c("#993404", "#D95F0E", "#004529", "#238443",
"#D9F0A3", "#4EB3D3", "#2B8CBE", "#c4e2f4",
"#FE9929", "#FEC44F", "#FEE391", "#A8DDB5",
"#E0F3DB", "#F7FCF0", "#c79999", "#0868AC",
"#FFFFD4", "white", "#dad4d4")
biome_mapping <- system.file("extdata", "biomes.csv",
package = "biospheremetrics") %>%
readr::read_delim(delim = ";", col_types = readr::cols())
names(biome_cols) <- biome_mapping$short_name
biome_cols_legend <- biome_cols[order_legend]
biome_names_legend <- biome_mapping$short_name[order_legend]
biomes_lpjml <- to_raster(lpjml_array = biome_data$biome_id,
boundary_box = bounding_box,
ext = lpjml_extent,
to_robinson = to_robinson)
if (!is.null(file_name)) {
file_extension <- strsplit(file_name, split = "\\.")[[1]][-1]
switch(file_extension,
`png` = {
grDevices::png(file_name,
width = 8 * 1.8,
height = 4 * 2,
units = "cm",
res = 600,
pointsize = 7
)
},
`pdf` = {
grDevices::pdf(file_name,
width = 8 * 1.8 / 2.54,
height = (4 * 2) / 2.54,
pointsize = 7
)
}, {
stop("File extension ", dQuote(file_extension), " not supported.")
}
)
}
brk <- seq(min(biome_mapping$id) - 0.5,
max(biome_mapping$id, na.rm = TRUE) + 0.5, 1)
graphics::par(mar = c(4, 0, 0, 0), xpd = TRUE, bg = bg_col)
Fabian Stenzel
committed
raster::image(biomes_lpjml, asp = 1, xaxt = "n", yaxt = "n",
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
xlab = "", ylab = "", col = biome_cols, breaks = brk,
lwd = 0.1, bty = "n")
raster::plot(countries, add = TRUE, lwd = 0.3,
border = "#5c565667", usePolypath = FALSE)
if (to_robinson == TRUE) {
ypoint <- (-6736039)
} else {
ypoint <- (-67)
}
legend_text <- biome_names_legend[1:19]
if (display_area) {
if (is.null(cellarea)) stop("Cellarea needs to be supplied for displaying.")
biome_area <- rep(0, length(order_legend))
names(biome_area) <- biome_mapping$short_name
for (i in seq_along(order_legend)) {
biome_area[i] <- sum(cellarea[which(biome_data$biome_id == i)])
}
biome_area <- round(biome_area / sum(cellarea) * 100, 3)
legend_text <- paste(legend_text, paste0("(", biome_area, " %)"))
}
graphics::legend(0, y = ypoint, xjust = 0.45, yjust = 1, cex = 0.8,
legend_text,
fill = biome_cols_legend[1:19],
horiz = FALSE, border = NULL, bty = "o", box.col = "white",
bg = bg_col, ncol = 4)
if (!is.null(file_name)) grDevices::dev.off()
}
# convert lpjml vector to raster and change projection to robinson
to_raster <- function(lpjml_array,
boundary_box,
ext,
to_robinson) {
crs_init <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
lpj_ras <- raster::raster(res = 0.5, crs = crs_init)
lpj_ras[raster::cellFromXY(lpj_ras, cbind(lon, lat))] <- lpjml_array
if (to_robinson) {
ras_to <- raster::raster(xmn = -18000000,
xmx = 18000000,
ymn = -9000000,
ymx = 9000000,
crs = "+proj=robin",
nrows = 2 * 360,
ncols = 2 * 720)
out_ras <- raster::crop(lpj_ras, ext) %>%
raster::projectRaster(to = ras_to, crs = "+proj=robin",
na.rm = TRUE, method = "ngb") %>%
raster::mask(boundary_box) %>%
suppressWarnings()
} else {
out_ras <- raster::crop(lpj_ras, ext)
}
return(out_ras)
}