Rasterise network: filter national sub-network if necessary.

Fix #44

DESCRIPTION

 Package: mapMCDA
 Title: Produce an epidemiological risk map by weighting multiple risk
         factors
-Version: 0.4.51
-Date: 2020-04-29
+Version: 0.4.52
+Date: 2020-05-23
 Authors@R: c( person("Andrea", "Apolloni", email =
         "andrea.apolloni@cirad.fr", role = c("ctb"), comment = "Animal
         mobility algorithm"), person("Elena", "Arsevska", email =

R/network.R

@@ -133,7 +133,7 @@ must be respected:",
     } else {
       ## Aggregate volume with warning
       warning(
-        "Aggregating volumes for multiple links.",
+        "Aggregating volumes for multiple links. ",
         "If this is not what you want, stop now and aggregate manually."
       )
       dat <- 
@@ -239,7 +239,9 @@ setOldClass("geonetwork")
 #' @param ... Other arguments passed on to
 #'   \code{\link[raster]{rasterize}}.
 #' @importMethodsFrom raster rasterize
-#' @importFrom igraph vertex.attributes
+#' @importFrom raster extract
+#' @importFrom igraph vertex.attributes induced_subgraph
+#' @importFrom sf st_transform st_crs as_Spatial
 #' @importFrom stats setNames
 #' @export
 #' @name rasterize_geonetwork
@@ -251,7 +253,27 @@ setMethod(
     function(
       x, y, field = "importance", ...) {
       
-      etx <- epidemic_threshold(x, beta = 1)
+      ## National subnetwork extraction
+      geo_nodes <- sf::st_transform(
+        attr(x, "geom_node"),
+        crs = sf::st_crs(y)
+      )
+      raster_value_at_nodes <- raster::extract(y, sf::as_Spatial(geo_nodes))
+      
+      if (any(idx <- is.na(raster_value_at_nodes))) {
+        sub_x <- igraph::induced_subgraph(
+          graph = x,
+          vids = which(!idx),
+          impl = "copy_and_delete"
+        )
+        
+        class(sub_x) <- class(x)
+        attr(sub_x, "geom_node") <- attr(x, "geom_node")[!idx,]
+      } else {
+        sub_x <- x
+      }
+      
+      etx <- epidemic_threshold(sub_x, beta = 1)
       
       ## If the graph is weighted, use weights
       if (is.null(epiR0 <- etx$weighted)) {
@@ -274,9 +296,9 @@ setMethod(
       )
       
       nodes <- cbind(
-        as.data.frame(igraph::vertex.attributes(x)),
+        as.data.frame(igraph::vertex.attributes(sub_x)),
         setNames(
-          data.frame(sf::st_coordinates(attr(x, "geom_node"))),
+          data.frame(sf::st_coordinates(attr(sub_x, "geom_node"))),
           c("Lon", "Lat")
         )
       )

tests/testthat/test-network.R

@@ -98,3 +98,38 @@ test_that("Network rasterisation", {
   expect_equal(vals_at_B, 100)
   expect_true(all(vapply(vals_not_B, all.equal, TRUE, 0)))
 })
+
+
+test_that("Network rasterisation takes into account the sub-network within boundaries", {
+  
+  r <- raster(nrow = 10, ncol = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10, vals = 1)
+  
+  ## A geonetwork with a single source and several destinations:
+  ## A -> B; A -> C; A -> D
+  ## Node C is beyond raster boundaries
+  edges <- data.frame(
+    O = rep("A", 3),
+    D = LETTERS[2:4],
+    stringsAsFactors = FALSE
+  )
+  nodes <- data.frame(
+    node = LETTERS[1:4],
+    x = rep(c(2, 8), times = 2),
+    y = c(2, 2, 11, 8)
+  )
+  gn <- geonetwork::geonetwork(edges, nodes)
+  
+  # plot(r)
+  # plot(gn, add = TRUE)
+  
+  expect_error(gnr <- rasterize(gn, r), NA)
+
+  ## Remove edges related with node C off-boundaries
+  edges2 <- edges[-2, ]
+  nodes2 <- nodes[-3, ]
+  gn2 <- geonetwork::geonetwork(edges2, nodes2)
+  expect_error(gnr2 <- rasterize(gn2, r), NA)
+
+  ## The results are identical
+  expect_identical(gnr, gnr2)  
+})