我正在开发一个绘图函数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)
谁能解释这里发生了什么?这显然是(?)一个范围界定问题,但我真的不明白为什么函数会这样!
提前致谢 !