Handle epidemiological units without CRS

fixes #27

DESCRIPTION

 Package: mapMCDA
-Title: Produce an epidemiological risk map by weighting multiple risk factors
-Version: 0.4.2
-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")),
-  person("Françoise", "Chirara", email = "francoise.chirara@cirad.fr", role = c("ctb"), comment = "Logo designer"),
-  person("Sylvain", "Falala", email = "sylvain.falala@cirad.fr", role = c("aut"), comment = "Shiny interface"),
-  person("Jean Marc", "Feussom", email = "mfeussom@gmail.com", role = c("ctb")),
-  person("Renaud", "Lancelot", email = "renaud.lancelot@cirad.fr", role = c("ctb")),
-  person("Raphaëlle", "Metras", email = "raphaelle.metras@cirad.fr", role = c("ctb")),
-  person("Facundo", "Muñoz", email = "facundo.munoz@cirad.fr", role = c("aut", "cre"), comment = c("Package developer", ORCID = "0000-0002-5061-4241"))
-  )
-Description: Use a Shiny interface to help users import multiple layers of risk
-  factors, scale and combine them using a Multi-Criteria Decision Analysis
-  approach to produce a risk map as an outcome.
+Title: Produce an epidemiological risk map by weighting multiple risk
+        factors
+Version: 0.4.3
+Date: 2019-04-11
+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")),
+        person("Françoise", "Chirara", email =
+        "francoise.chirara@cirad.fr", role = c("ctb"), comment = "Logo
+        designer"), person("Sylvain", "Falala", email =
+        "sylvain.falala@cirad.fr", role = c("aut"), comment = "Shiny
+        interface"), person("Jean Marc", "Feussom", email =
+        "mfeussom@gmail.com", role = c("ctb")), person("Renaud",
+        "Lancelot", email = "renaud.lancelot@cirad.fr", role =
+        c("ctb")), person("Raphaëlle", "Metras", email =
+        "raphaelle.metras@cirad.fr", role = c("ctb")),
+        person("Facundo", "Muñoz", email = "facundo.munoz@cirad.fr",
+        role = c("aut", "cre"), comment = c("Package developer", ORCID
+        = "0000-0002-5061-4241")) )
+Description: Use a Shiny interface to help users import multiple layers
+        of risk factors, scale and combine them using a Multi-Criteria
+        Decision Analysis approach to produce a risk map as an outcome.
 Depends: R (>= 3.2)
 License: GPL-3 + file LICENSE
 Encoding: UTF-8
 LazyData: true
-Imports:
-  classInt,
-  deldir,
-  geojsonio,
-  geojsonlint,
-  geonetwork,
-  igraph,
-  methods,
-  maps,
-  maptools,
-  plyr,
-  raster,
-  rasterVis,
-  RColorBrewer,
-  rgdal,
-  rgeos,
-  rhandsontable,
-  shiny,
-  shinydashboard,
-  shinyFiles,
-  sp,
-  stringr,
-  utils
-Suggests:
-    devtools,
-    knitr,
-    mapview,
-    rmarkdown,
-    roxygen2,
-    testthat
+Imports: classInt, deldir, geojsonio, geojsonlint, geonetwork, igraph,
+        methods, maps, maptools, plyr, raster, rasterVis, RColorBrewer,
+        rgdal, rgeos, rhandsontable, shiny, shinydashboard, shinyFiles,
+        sp, stringr, utils
+Suggests: devtools, knitr, mapview, rmarkdown, roxygen2, testthat
 RoxygenNote: 6.1.1
 URL: https://Cirad-ASTRE.github.io/mapMCDA
 BugReports: https://github.com/Cirad-ASTRE/mapMCDA/issues

R/risk_unit.R

@@ -22,6 +22,28 @@
 #'   sp::spplot(cmr$cmr_admin3[, "rv"], cuts = 3)
 risk_unit <- function(r, eu, fun = mean) {
   rgrid <- methods::as(r, "SpatialGridDataFrame")  # needed for overlay methods
+  
+  ## Possible differences in CRS
+  if (!identicalCRS(r, eu)) {
+    ## r should have some CRS defined at this point
+    stopifnot(!is.na(proj4string(r)))
+    
+    ## If the epidemiological units don't have CRS, assume geographical
+    ## whenever possible and reproject to match r's CRS
+    if (is.na(proj4string(eu))) {
+      if (couldBeLonLat(eu)) {
+        proj4string(eu) <- CRS("+proj=longlat +datum=WGS84")
+      } else {
+        stop(
+          "Missing Coordinate Reference System (CRS) in the layer of epidemiological units.\n",
+          "Please load a layer with the CRS information."
+        )
+      }
+    }
+    
+    eu <- spTransform(eu, CRS(proj4string(r)))
+  }
+  
   funrisk_poly <- over(eu, rgrid, fn = fun)[[1]]
   
   ## Small-polygon correction

inst/interface/server.R

@@ -435,7 +435,7 @@ server <- function(input, output, session) {
     
     if(is.null(epidUnitLayer)) return(NULL)
     
-    if( !is.projected(epidUnitLayer)) {
+    if( !isTRUE(is.projected(epidUnitLayer))) {
       warning("This map is not projected. This can lead to very
               inaccurate computations of distances and areas, depending
               on the location and size of the region of interest.

tests/testthat/test-risk_unit.R

 context("risk_unit")
 
 ## 10 x 10 raster with sequential values
-r <- raster(nrow = 10, ncol = 10)
-r[] <- seq.int(ncell(r))
+r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, nrow = 10, ncol = 10)
+r <- setValues(r, seq.int(ncell(r)))
 
 ## polygons covering the region by squares of 2x2 pixels
 r0 <- r
 res(r0) <- 2*res(r)
 epiunits <- as(r0, "SpatialPolygons")
 
-## index for the polygons (for each pixel, what block it belongs to)
-block <- 
-  do.call(
-    paste0,
-    expand.grid(x = rep(1:5, each = 2), y = rep(1:5, each = 2))
-  )
-block <- factor(block, levels = unique(block))
-
 test_that("risk_unit returns a vector of length n_polygons", {
 
   riskvalues <- risk_unit(r, epiunits)  
@@ -27,7 +19,16 @@ test_that("risk_unit returns a vector of length n_polygons", {
 })
 
 test_that("alternative risk summaries", {
+
+  ## index for the polygons (for each pixel, what block it belongs to)
+  block <- 
+    do.call(
+      paste0,
+      expand.grid(x = rep(1:5, each = 2), y = rep(1:5, each = 2))
+    )
+  block <- factor(block, levels = unique(block))
   
+
   expect_summary <- function(fun) {
     expect_identical(
       risk_unit(r, epiunits, fun = fun),
@@ -42,3 +43,32 @@ test_that("alternative risk summaries", {
   expect_summary(sum)
 
 })
+
+
+test_that("Missing CRS in epidemiological units", {
+  
+  eu_na <- epiunits
+  proj4string(eu_na) <- CRS()
+
+  ## If no CRS but possibly geographical coordinates, 
+  ## assume them, and yield results with a Warning
+  res1 <- expect_warning(risk_unit(r, eu_na), "CRS is NA")
+  expect_true(!any(is.na(res1)))
+  
+  ## Also, reproject if necessary to match CRS of r
+  # proj.4 projection description
+  newproj <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
+  rt <- suppressWarnings(projectRaster(r, crs = CRS(newproj)))
+  
+  res2 <- expect_warning(risk_unit(rt, eu_na), "CRS is NA")
+  expect_true(!any(is.na(res2)))
+  
+  ## However, if the CRS is missing but coordinates are projected,
+  ## fail.
+  eu_pr <- epiunits
+  eu_pr <- spTransform(eu_pr, CRS(newproj))
+  proj4string(eu_pr) <- CRS()
+  
+  expect_error(risk_unit(rt, eu_pr), "Missing Coordinate Reference System")
+})
+