如何 cut/crop/mask autoKrige 结果与边界多边形?
How to cut/crop/mask autoKrige result with boundary polygon?
我正忙于开发一个模块,该模块在单元格大小为 20x20m 的土壤样本数据上实现 autoKrige(自动地图库)。克里金法完成后,我想 crop/mask 克里金法结果与字段的边界。
cropping/masking 的问题在于(由于 20x20 单元格)字段的边框产生 "step" 效果。我正在寻找的是一个平滑的边界(穿过单元格)。
下面是生成上述两个场景的代码:
library(sp)
library(rgeos)
#create polygon
r1 <- cbind(c(641777, 642290, 642276, 641794), c(7036885, 7036743, 7036154, 7036146))
r2 <- cbind(c(642320, 642478, 642494, 642314), c(7036732, 7036699, 7036296, 7036290))
sr1 <- Polygons(list(Polygon(r1)),"r1")
sr2 <- Polygons(list(Polygon(r2)),"r2")
boundary.sp <- SpatialPolygons(list(sr1,sr2))
boundary.sp@proj4string <- CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs')
#create bounding box grid
bbox <- bbox(boundary.sp)
boundary.grid <- expand.grid(x = seq(from = bbox[1], to = bbox[3], by = 20), y = seq(from = bbox[2], to = bbox[4], by = 20))
coordinates(boundary.grid) <- ~x + y
gridded(boundary.grid) <- TRUE
boundary.grid@proj4string <- boundary.sp@proj4string
#create SpatialPixels grid
boundary.grid.stepped <- boundary.grid[!is.na(over(boundary.grid, boundary.sp)),]
plot(boundary.grid.stepped)
#cut grid with polygon to create SpatialPolygons grid
boundary.poly.grid <- as.SpatialPolygons.GridTopology(getGridTopology(boundary.grid), proj4string = CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs'))
boundary.grid.smooth <- gIntersection(boundary.poly.grid, boundary.sp, byid=TRUE)
plot(boundary.grid.smooth)
当前调用 autoKrige 时,上述网格 (boundary.grid.stepped) 作为 'new_data' 参数传递。
哪种方法更好,我该如何实施:
1) 事先准备好目标网格并将其用作'new_data' or,
2) 在边界框网格上进行克里金法然后切割/crop/mask?
我会说直接在目标网格上进行网格化,通过将其馈送到 new_data
。克里金结果仅取决于输入数据,所以我怀疑 cropping/calculating 和 calculating/cropping.
之间几乎没有区别
我正忙于开发一个模块,该模块在单元格大小为 20x20m 的土壤样本数据上实现 autoKrige(自动地图库)。克里金法完成后,我想 crop/mask 克里金法结果与字段的边界。
cropping/masking 的问题在于(由于 20x20 单元格)字段的边框产生 "step" 效果。我正在寻找的是一个平滑的边界(穿过单元格)。
下面是生成上述两个场景的代码:
library(sp)
library(rgeos)
#create polygon
r1 <- cbind(c(641777, 642290, 642276, 641794), c(7036885, 7036743, 7036154, 7036146))
r2 <- cbind(c(642320, 642478, 642494, 642314), c(7036732, 7036699, 7036296, 7036290))
sr1 <- Polygons(list(Polygon(r1)),"r1")
sr2 <- Polygons(list(Polygon(r2)),"r2")
boundary.sp <- SpatialPolygons(list(sr1,sr2))
boundary.sp@proj4string <- CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs')
#create bounding box grid
bbox <- bbox(boundary.sp)
boundary.grid <- expand.grid(x = seq(from = bbox[1], to = bbox[3], by = 20), y = seq(from = bbox[2], to = bbox[4], by = 20))
coordinates(boundary.grid) <- ~x + y
gridded(boundary.grid) <- TRUE
boundary.grid@proj4string <- boundary.sp@proj4string
#create SpatialPixels grid
boundary.grid.stepped <- boundary.grid[!is.na(over(boundary.grid, boundary.sp)),]
plot(boundary.grid.stepped)
#cut grid with polygon to create SpatialPolygons grid
boundary.poly.grid <- as.SpatialPolygons.GridTopology(getGridTopology(boundary.grid), proj4string = CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs'))
boundary.grid.smooth <- gIntersection(boundary.poly.grid, boundary.sp, byid=TRUE)
plot(boundary.grid.smooth)
当前调用 autoKrige 时,上述网格 (boundary.grid.stepped) 作为 'new_data' 参数传递。
哪种方法更好,我该如何实施:
1) 事先准备好目标网格并将其用作'new_data' or,
2) 在边界框网格上进行克里金法然后切割/crop/mask?
我会说直接在目标网格上进行网格化,通过将其馈送到 new_data
。克里金结果仅取决于输入数据,所以我怀疑 cropping/calculating 和 calculating/cropping.