-
Notifications
You must be signed in to change notification settings - Fork 2
/
utils_spatial_misc.R
177 lines (164 loc) · 4.72 KB
/
utils_spatial_misc.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
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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#' @include internal.R
NULL
#' Shapefile field names
#'
#' Quickly extract the field names in a shapefile.
#'
#' @param x `character` path.
#'
#' @param inherits `character` class to return only column names that
#' contain a certain type of data.
#' Defaults to `NULL` such that all column names are returned.
#'
#' @return `character` vector of field names.
#'
#' @examples
#' shapefile_field_names(system.file("shape/nc.shp", package = "sf"))
#' @export
shapefile_field_names <- function(x, inherits = NULL) {
# assert arguments are valid
assertthat::assert_that(
assertthat::is.string(x),
endsWith(x, ".shp")
)
# find layer name
l <- sf::st_layers(x)$name[[1]]
# prepare query
qu <- paste0("SELECT * FROM \"", l, "\" WHERE FID <= 2")
# import shapefile
s <- sf::read_sf(dsn = x, query = qu)
# drop geometry data
s <- sf::st_drop_geometry(s)
# if inherits is specified, then subset columns with specified data type
if (!is.null(inherits)) {
s <- dplyr::select_if(s, function(x) inherits(x, inherits))
}
# return column names
names(s)
}
#' Repair spatial data
#'
#' Repair geometry in a spatial dataset.
#'
#' @param x [sf::st_sf()] object.
#'
#' @return [sf::st_sf()] object.
#'
#' @examples
#' path <- system.file("shape/nc.shp", package = "sf")
#' g <- sf::read_sf(path)
#' repair_spatial_data(g)
#' @export
repair_spatial_data <- function(x) {
# assert arguments are valid
assertthat::assert_that(
inherits(x, "sf")
)
# repair geometry
x <- sf::st_make_valid(x)
# return result
x
}
#' Reserve sizes
#'
#' Calculate the size of individual reserves in a solution.
#'
#' @param x `numeric` vector containing the solution values for each planning
#' unit.
#'
#' @param areas `numeric` vector containing the area of each planning unit.
#'
#' @param boundary_matrix [Matrix::sparseMatrix()] object with boundary length
#' data for the planning units. Or [`NA`] indicating that no boundary data is
#' available (see Details for more information)
#'
#' @return A `numeric` vector containing the total area of each reserve.
#'
#' @details For spatial uploads (shapefile) with many planning units, building
#' boundary data can result in a std::bad_alloc error. To avoid this, the
#' user can skip generating a boundary matrix on the `new_dataset_from_auto`
#' method. For these scenarios, reserve sizes can not be calculated when the
#' `boundary_matrix` is set to `NA`.
#'
#' @export
reserve_sizes <- function(x, areas, boundary_matrix) {
# assert the argument are valid
assertthat::assert_that(
## x
is.numeric(x),
assertthat::noNA(x),
## areas
is.numeric(areas),
assertthat::noNA(areas),
length(areas) == length(x),
## boundary matrix
inherits(boundary_matrix, c("dsCMatrix", "dgCMatrix", "logical")))
if (inherits(boundary_matrix, c("dsCMatrix", "dgCMatrix"))) {
assertthat::assert_that(
nrow(boundary_matrix) == ncol(boundary_matrix),
nrow(boundary_matrix) == length(x)
)
}
# return NA if boundary_matrix is NA
if (!inherits(boundary_matrix, c("dsCMatrix", "dgCMatrix"))) {
return(NA_real_)
}
# return NA if no planning units selected, then return NA
if (sum(x) < 0.5) {
return(NA_real_)
}
# create adjacency matrix to describe relationships among units
idx <- which(x > 0.5)
adj_matrix <- boundary_matrix[idx, idx]
adj_matrix@x <- rep(1, length(adj_matrix@x))
Matrix::diag(adj_matrix) <- 1
# subset areas to contain only selected planning units
areas <- areas[idx]
# create graph
g <- igraph::graph_from_adjacency_matrix(adj_matrix, mode = "undirected")
# identify components
clu <- igraph::components(g)
# calculate total size of each reserve
vapply(seq_len(clu$no), FUN.VALUE = numeric(1), function(i) {
sum(areas[clu$membership == i])
})
}
#' Calculate coverage
#'
#' Calculate how well a solution covers data.
#'
#' @param x `numeric` solution values.
#'
#' @param data [Matrix::sparseMatrix()] object.
#'
#' @return `numeric` vector.
#'
#' @examples
#' # load dependency
#' library(Matrix)
#'
#' # simulate solution values for 10 planning units
#' solution_values <- sample(c(0, 1), 10, replace = TRUE)
#'
#' # simulate data for 5 features in those 10 planning units
#' feature_data <- as(matrix(runif(10 * 5), ncol = 10), "dgCMatrix")
#'
#' # calculate coverage
#' calculate_coverage(solution_values, feature_data)
#' @export
calculate_coverage <- function(x, data) {
assertthat::assert_that(
is.numeric(x),
inherits(data, "dgCMatrix")
)
if (nrow(data) > 0) {
out <- Matrix::Matrix(
x, byrow = TRUE, ncol = 1, nrow = length(x), sparse = TRUE
)
out <- as.numeric(data %*% out) / Matrix::rowSums(data)
names(out) <- rownames(data)
} else {
out <- numeric(0)
}
out
}