Commit 7709800e authored by Facundo Muñoz's avatar Facundo Muñoz ®️
Browse files

bugfix: compute distance to point or lines

- Also, add tests for distance_map()

- fixes #14
parent 401eac59
Package: mapMCDA
Title: Produce an epidemiological risk map by weighting multiple risk factors
Version: 0.3.4
Date: 2018-06-27
Version: 0.3.5
Date: 2019-04-06
Authors@R: c(
person("Andrea", "Apolloni", email = "andrea.apolloni@cirad.fr", role = c("ctb"), comment = "Animal mobility algorithm"),
person("Elena", "Arsevska", email = "elena.arsevska@cirad.fr", role = c("ctb")),
......
......@@ -6,7 +6,7 @@
#'
#'
#' @param x Spatial* object with relevant spatial features.
#' @param boundaries SpatialPolygon* object with administrative
#' @param boundaries SpatialPolygons* object with administrative
#' borders, or Raster* object from where to pick the resoltion.
#' @param res numeric. Resolution for the outcome distance map. The default divides
#' the smallest dimension into 100 cells.
......@@ -21,9 +21,12 @@
distance_map <- function(x, boundaries, res = resolution(boundaries, min_ncells = 100)) {
if (inherits(boundaries, "Spatial")) {
max_res <- resolution(boundaries, min_ncells = 2)
stopifnot(res > 0, res < max_res)
ext_grid <- raster::raster(raster::extent(boundaries), resolution = res)
msk <- raster::rasterize(boundaries, ext_grid, field = 1, background = NA, fun = "mean")
} else if (inherits(boundaries, "Raster")) {
## Ignores res
msk <- boundaries
} else stop(paste(substitute(boundaries), "must be either a SpatialPolygon* or a Raster* object."))
......@@ -43,7 +46,8 @@ distance_map <- function(x, boundaries, res = resolution(boundaries, min_ncells
#' @param x Spatial* object
#' @param mask Raster* object
#'
#' @return Raster* object. Distance map in meters.
#' @return Raster* object. Distance map in meters (even for unprojected maps in
#' geographical coordinates).
#'
#' @import raster
#' @export
......@@ -54,7 +58,7 @@ distance_map <- function(x, boundaries, res = resolution(boundaries, min_ncells
#' raster::plot(dist_wb)
distance_to_vector <- function(x, mask) {
xr <- raster::rasterize(x, mask, field = 0, fun = "mean", background = NA)
dx <- raster::mask(raster::distance(xr), mask)
xr <- raster::rasterize(x, mask, field = 0, fun = "last", background = NA)
raster::mask(raster::distance(xr), mask)
}
......@@ -10,7 +10,7 @@ distance_map(x, boundaries, res = resolution(boundaries, min_ncells =
\arguments{
\item{x}{Spatial* object with relevant spatial features.}
\item{boundaries}{SpatialPolygon* object with administrative
\item{boundaries}{SpatialPolygons* object with administrative
borders, or Raster* object from where to pick the resoltion.}
\item{res}{numeric. Resolution for the outcome distance map. The default divides
......
......@@ -12,7 +12,8 @@ distance_to_vector(x, mask)
\item{mask}{Raster* object}
}
\value{
Raster* object. Distance map in meters.
Raster* object. Distance map in meters (even for unprojected maps in
geographical coordinates).
}
\description{
Compute a raster file with the distances to the nearest vector feature,
......
context("Distance maps")
bnd <-
SpatialPolygons(list(Polygons(list(Polygon(
cbind(c(0, 10, 10, 0, 0), c(0, 0, 10, 10, 0))
)),
"Boundary")),
pO = 1L,
proj4string = CRS("+proj=longlat"))
## Raster Mask
rmsk <- raster(
nrows = 21,
ncols = 51,
xmn = 0,
xmx = 10,
ymn = 0,
ymx = 10,
crs = crs("+proj=longlat"),
vals = 1
)
polyA <- Polygon(cbind(c(1, 2, 2, 1, 1), c(1, 1, 2, 2, 1)))
polyB <- Polygon(cbind(c(5, 6, 6, 5, 5), c(3, 3, 4, 4, 3)))
plg1 <-
SpatialPolygons(list(Polygons(list(polyA),
"A")),
1L, proj4string = CRS("+proj=longlat"))
plg2 <-
SpatialPolygons(list(Polygons(list(polyA), "A"),
Polygons(list(polyB), "B")),
1:2L, proj4string = CRS("+proj=longlat"))
pts1 <- SpatialPoints(cbind(c(3, 7), c(3, 3)), proj4string = CRS("+proj=longlat"))
ln1 <- SpatialLines(list(Lines(Line(cbind(c(1, 9), c(8, 2))), "Line")), proj4string = CRS("+proj=longlat"))
# plot(bnd)
# plot(plg1, col = "blue", add = TRUE)
# plot(plg2, col = "red", add = TRUE)
# plot(pts1, col = "black", add = TRUE)
# plot(ln1, col = "green", add = TRUE)
test_that("Distance from points", {
ans <- expect_error(
distance_to_vector(pts1, rmsk),
NA
)
# plot(ans)
expect_is(ans, "RasterLayer")
expect_true(all(values(ans) >= 0))
})
test_that("Distance from lines", {
ans <- expect_error(
distance_to_vector(ln1, rmsk),
NA
)
# plot(ans)
expect_is(ans, "RasterLayer")
expect_true(all(values(ans) >= 0))
})
test_that("Distance from polygons", {
## Single polygon
ans1 <- expect_error(
distance_to_vector(plg1, rmsk),
NA
)
# plot(ans1)
expect_is(ans1, "RasterLayer")
expect_true(all(values(ans1) >= 0))
## Multiple polygons
ans2 <- expect_error(
distance_to_vector(plg2, rmsk),
NA
)
# plot(ans2)
expect_is(ans2, "RasterLayer")
expect_true(all(values(ans2) >= 0))
expect_true(all(values(ans2) <= values(ans1)))
})
test_that("distance_map() takes either a RasterLayer or a SpatialPolygons*", {
expect_error(
distance_map(pts1, bnd),
NA
)
expect_error(
distance_map(pts1, rmsk),
NA
)
})
test_that("distance_map() handles resolution spec", {
expect_error(
distance_map(pts1, bnd, res = 1),
NA
)
expect_error(
distance_map(pts1, bnd, res = 0),
"res > 0 is not TRUE"
)
expect_error(
distance_map(pts1, bnd, res = 20),
"res < max_res"
)
})
......@@ -25,9 +25,25 @@ test_that("risk re-scaling of a raster", {
test_that("risk map from vector: compute distances", {
## Polygons
st <- c(0, 100)
dwb <- expect_error(
risk_layer(
cmr$water_bodies, boundaries = cmr$cmr_admin3, scale_target = st
),
NA # Expects no error
)
## Raster values must be between 0 and 100
expect_identical(range(values(dwb), na.rm = TRUE), c(0, 100))
## Points
set.seed(20190405)
pts <- spsample(cmr$cmr_admin3, 10, type = "random")
expect_error(
risk_layer(cmr$water_bodies, boundaries = cmr$cmr_admin3, scale_target = st),
risk_layer(
pts, boundaries = cmr$cmr_admin3, scale_target = rev(st)
),
NA # Expects no error
)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment