我正在研究基于rasterVis::levelplot
的绘图函数,用户可以仅通过栅格对象或栅格对象以及多边形对象。在栅格中绘制具有覆盖多边形的栅格时确定问题的范围
功能是相当复杂的,但显示该问题的最小子集读作:
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)
}
的问题是,我看到一些奇怪的(对我来说)的结果,同时测试它。让我们来定义一些虚拟的数据 - 一个1000×1000的光栅和sf
POYGON oject四个多边形其分裂光栅 - :
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)
任何人都可以解释这里发生了什么?很明显(?)是一个范围问题,但我真的不明白为什么函数的行为像这样!
在此先感谢!