在 rasterVis 中绘制具有重叠多边形的栅格时的范围问题

Question

我正在开发一个基于 rasterVis::levelplot 的绘图函数，用户可以将光栅对象或光栅对象和 sf 多边形对象传递给该函数。

该函数相当复杂，但显示问题的最小子集如下：

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)
}

问题是我在测试时看到了一些奇怪的（对我来说）结果。让我们定义一些虚拟数据 - 一个 1000x1000 的光栅和一个 sf POYGON 对象，它有四个分割光栅的多边形 -:

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"))

并测试功能。理论上，我认为这应该有效：

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

但我得到：

这样做：

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

仍然失败，但结果不同：

我发现让它工作的唯一方法是为多边形对象赋予我在函数中使用的相同名称（即，in_poly）从一开始：

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)

谁能解释一下这里发生了什么？这显然（？）是一个范围问题，但我真的不明白为什么函数会这样！

提前致谢！

Answer 1

latticeExtra::layer 的帮助页面解释说：

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.

在函数内部使用 layer 时，您可以将对象嵌入列表中并将其传递到 data 参数中：

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)
}

现在函数产生了想要的结果：

myplot(in_rast, in_poly = my_poly)

在 rasterVis 中绘制具有重叠多边形的栅格时的范围问题

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

r

plot

scoping

sp

sf