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
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
#' 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 to_robinson logical to define if robinson projection should be used
#' for plotting
#'
#' @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,
to_robinson = TRUE,
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
order_legend <- c(1, 2, 9, 10, 11, 3, 4, 5, 12, 13, 14, 6, 7, 8, 15, 16,
17, 18, 19)
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` = {
png(file_name,
width = 8 * 1.8,
height = 4 * 2,
units = "cm",
res = 600,
pointsize = 7)
},
`pdf` = {
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)
par(mar = c(4, 0, 0, 0), xpd = T, bg = bg_col)
image(biomes_lpjml, asp = 1, xaxt = "n", yaxt = "n",
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(0, y = ypoint, xjust = 0.45, yjust = 1, cex = 0.8,
biome_names_legend[1:19],
fill = biome_cols_legend[1:19],
horiz = F, border = NULL, bty = "o", box.col = "white",
bg = bg_col, ncol = 4)
if (!is.null(file_name)) 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)
#TODO replace lpjmliotools
lpj_ras[raster::cellFromXY(lpj_ras, cbind(lon, lat))] <-
lpjml_array
if (to_robinson) {
ras_to <- 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)
}