Search code examples
rplotscopingr-spr-sf

scoping issues while plotting a raster with an overlayed polygon in rasterVis

I'm working on a plotting function based on rasterVis::levelplot to which the user can pass either just a raster object, or a raster object and an sf polygon object.

The function is rather complex, but a minimum subset showing the problem reads as: 

library(sf)
library(raster)
library(rasterVis)

myplot <- function(in_rast, in_poly = NULL) {
  rastplot <- rasterVis::levelplot(in_rast, margin = FALSE)
  polyplot <- layer(sp::sp.polygons(in_poly))
  print(rastplot + polyplot)
}

The problem is that I see some strange (for me) results while testing it. Let's define some dummy data - a 1000x1000 raster and a sf POYGON oject with four polygons which split the raster -: 

in_rast  <- raster(matrix(nrow = 1000, ncol = 1000))
in_rast  <- setValues(in_rast, seq(1:1000000))

my_poly  <- structure(list(cell_id = 1:4, geometry = structure(list(structure(list(
    structure(c(0, 0.5, 0.5, 0, 0, 0, 0, 0.5, 0.5, 0), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0.5, 1, 1, 0.5, 0.5, 0, 0, 0.5, 0.5, 0), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0, 0.5, 0.5, 0, 0, 0.5, 0.5, 1, 1, 0.5), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0.5, 1, 1, 0.5, 0.5, 0.5, 0.5, 1, 1, 0.5), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, class = c("sfc_POLYGON", 
"sfc"), precision = 0, crs = structure(list(epsg = NA_integer_, 
    proj4string = NA_character_), .Names = c("epsg", "proj4string"
), class = "crs"), bbox = structure(c(0, 0, 1, 1), .Names = c("xmin", 
"ymin", "xmax", "ymax")))), .Names = c("cell_id", "geometry"), row.names = c(NA, 
4L), class = c("sf", "data.frame"), sf_column = "geometry", 
agr = structure(NA_integer_, class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = "cell_id"))

and test the function. In theory, I think this should work: 

my_poly  <- as(my_poly, "Spatial") # convert to spatial
myplot(in_rast, in_poly = my_poly)

but I get:

enter image description here

doing this: 

in_poly <- my_poly
in_poly <- as(in_poly, "Spatial")
myplot(in_rast, in_poly = in_poly)

still fails, but with a different outcome:

enter image description here

The only way I found to have it working is to give to the polygon object the same name that I use inside the function (i.e., in_poly) from the beginning :

in_poly <- structure(list(cell_id = 1:4, geometry = structure(list(structure(list(
    structure(c(0, 0.5, 0.5, 0, 0, 0, 0, 0.5, 0.5, 0), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0.5, 1, 1, 0.5, 0.5, 0, 0, 0.5, 0.5, 0), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0, 0.5, 0.5, 0, 0, 0.5, 0.5, 1, 1, 0.5), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(0.5, 1, 1, 0.5, 0.5, 0.5, 0.5, 1, 1, 0.5), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, class = c("sfc_POLYGON", 
"sfc"), precision = 0, crs = structure(list(epsg = NA_integer_, 
    proj4string = NA_character_), .Names = c("epsg", "proj4string"
), class = "crs"), bbox = structure(c(0, 0, 1, 1), .Names = c("xmin", 
"ymin", "xmax", "ymax")))), .Names = c("cell_id", "geometry"), row.names = c(NA, 
4L), class = c("sf", "data.frame"), sf_column = "geometry", 
agr = structure(NA_integer_, class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = "cell_id"))


in_poly  <- as(in_poly, "Spatial")
myplot(in_rast, in_poly = in_poly)

enter image description here

Can anyone explain what's happening here ? It's clearly (?) a scoping problem, but I really do not understand why the function behaves like this !

Thanks in advance !

Solution

  • The help page of latticeExtra::layer explains that:

    the evaluation used in layer is non-standard, and can be confusing at first: you typically refer to variables as if inside the panel function (x, y, etc); you can usually refer to objects which exist in the global environment (workspace), but it is safer to pass them in by name in the data argument to layer.

    When using layer inside a function, you can embed your object in a list and pass it in the data argument:

    myplot <- function(in_rast, in_poly = NULL) {
  rastplot <- levelplot(in_rast, margin = FALSE)
  polyplot <- layer(sp.polygons(x), 
                    data = list(x = in_poly))
  print(rastplot + polyplot)
}

    Now the function produces the desired result:

    myplot(in_rast, in_poly = my_poly)

    enter image description here