prioritization_scripts/3-num_features.R at main · matthewharris02/prioritization_scripts · GitHub

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
# 0. Load libraries ====
print("Loading libraries")
library(data.table)
library(arrow)
library(tidyverse)
library(terra)
# library(matrixStats)
library(viridisLite)

dir_wd <- "/mnt/sda/MH_restoration"
dir_src <- dir_wd
runid <- ""           # Additional ID to distinguish runs
source(file.path(dir_src, "script_tools/1.1-OPTIONS_20.R"))
dir_output <- file.path(dir_out, "output", ifelse(runid == "", "default", runid))

rast_template <- rast(
    crs = crs(EPSG),
    res = c(1000 * RES, 1000 * RES),
    ext = ext(EXT)
)


grid_cell <- open_dataset(file.path(dir_proc, "global_cells"),
                          partitioning = c("ISONUM"))  |>
    select(-id) |>
    collect()

ft <- grid_cell |>
    select(-lulc_converted, -lulc_other, -pu, -ISONUM) |>
    mutate(
        across(
            .cols = !c( x, y),
            .fns = ~if_else(is.na(.x), 0, 1)
        )
    )


# Create a raster where each cell is the number of features that exist for that cell

num_all <- ft[, `:=` (SUM = rowSums(.SD)),
                        .SDcols = c("ecoregions","ft_carbon",
                                    "ft_coastal", "ft_iucnrichness", "ft_kba",
                                    "ft_mangroves", "ft_pollination", "ft_ramsar", "ft_saltmarshes",
                                    "ft_usefulplants", "ft_waterquality")]

num_noEco <- ft[, `:=` (SUM = rowSums(.SD)),
                        .SDcols = c("ft_carbon",
                                    "ft_coastal", "ft_iucnrichness", "ft_kba",
                                    "ft_mangroves", "ft_pollination", "ft_ramsar", "ft_saltmarshes",
                                    "ft_usefulplants", "ft_waterquality")]

num_noC <- ft[, `:=` (SUM = rowSums(.SD)),
                       .SDcols = c(
                                   "ft_coastal", "ft_iucnrichness", "ft_kba",
                                   "ft_mangroves", "ft_pollination", "ft_ramsar", "ft_saltmarshes",
                                   "ft_usefulplants", "ft_waterquality")]

num_noCoastal <- ft[, `:=` (SUM = rowSums(.SD)),
                       .SDcols = c("ft_carbon",
                                   "ft_iucnrichness", "ft_kba",
                                   "ft_pollination", "ft_ramsar",
                                   "ft_usefulplants", "ft_waterquality")]


r_num_features <- rast(num_noCoastal[, c("x", "y", "SUM")],
                       crs = crs(EPSG),
                       extent = ext(rast_template))

writeRaster(r_num_features, file.path(dir_output, "num_features_noCoastal.tif"))

r_num_features_all <- rast(num_all[, c("x", "y", "SUM")],
                       crs = crs(EPSG),
                       extent = ext(rast_template))

writeRaster(r_num_features_all, file.path(dir_output, "num_features_all.tif"))

# Analyze solution ranking and number of features
solution <- read_csv(file.path(dir_output, "solution_full_lp_5km_0.01g_1t_1b.csv"))

all <- left_join(num_all, solution, join_by(x == x, y == y)) |>
    mutate(SUM = as.factor(SUM))


colours <- viridis(n = 12)


ggplot(all, aes(x = final, fill = SUM)) +
    scale_x_continuous(breaks = seq(1,10)) +
    geom_bar(just = 0.5) +
    scale_fill_manual(values = colours) +
    labs(
        x = "Solution Ranking",
        y = "Number of cells",
        fill = "Number of features"
    )