我正在尝试将SpatialGrid类中存在 1 个多边形的子集SpatialPolygons。我怎样才能做到这一点？
我试过这样：
grd.clip <- grd[!is.na(over(grd, polygon))]
但我得到错误
Error in matrix(idx, gr@cells.dim[2], gr@cells.dim[1], byrow = TRUE)[rows, :
(subscript) logical subscript too long
问问题
3665 次
3 回答
4
我尝试了一种方法，结果证明不是理想的解决方案。我将其保留在这里以说明我的原始想法，但如有必要我愿意将其删除
我仍然不清楚最初问题的解决方案，但如有必要，我会做得更好
library(sp)
library(rgdal)
library(raster)
library(latticeExtra)
加载 shapefile
shp <- readOGR(dsn = "D:/Programacao/R/Stackoverflow/17962821", layer = "shp")
proj4string(shp)
创建网格拓扑
grid <- GridTopology(cellcentre.offset=c(731888.0,7457552.0),
cellsize=c(16,16),cells.dim=c(122,106))
grid <- SpatialGrid(grid, proj4string=CRS(proj4string(shp)))
将 SpatialGrid 转换为 rasterLayer
rgrid <- raster(extent(grid))
res(rgrid) <- c(16, 16)
给它一些数字
rgrid[] <- runif(ncell(grid), 1, 10)
proj4string(rgrid) <- CRS(proj4string(shp))
plot(rgrid)
使用 SPDF 遮罩栅格
rgrid_msk <- mask(rgrid,shp)
plot(rgrid_msk)
将其转换回保留属性值的网格
grid_ae <- as(rgrid_msk, 'SpatialPointsDataFrame')
grid_ae <- grid_ae[!is.na(grid_ae@data$layer), ]
gridded(grid_ae) <- TRUE
summary(grid_ae)
> summary(grid_ae)
Object of class SpatialPixelsDataFrame
Coordinates:
min max
x 731912 733816
y 7457560 7459224
Is projected: TRUE
proj4string :
[+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs]
Number of points: 7814
Grid attributes:
cellcentre.offset cellsize cells.dim
x 731920 16 119
y 7457568 16 104
Data attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.005 3.231 5.523 5.512 7.748 9.999
spplot(grid_ae) +
latticeExtra::layer(sp.polygons(shp, fill = NA, col = 'red'))
与常规区域相交后保留 SPDF 属性的解决方案
library(rgeos)
library(rgdal)
library(sp)
library(latticeExtra)
grid <- readOGR(dsn = 'S:/Temporarios', layer = 'grid')
proj4string(grid) <- CRS('+init=epsg:4326')
grid
class : SpatialPolygonsDataFrame
features : 110
extent : -9.6, -7.95, 36.45, 37.95 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : ID
min values : 652
max values : 761
summary(grid@data)
ID
Min. :652.0
1st Qu.:679.2
Median :706.5
Mean :706.5
3rd Qu.:733.8
Max. :761.0
polyg <- readOGR(dsn = 'S:/Temporarios', layer = 'polyg')
proj4string(polyg) <- CRS('+init=epsg:4326')
# plot it
spplot(grid, 'ID') +
latticeExtra::layer(sp.polygons(polyg, fill = NA, col = 'blue'))´
# clip
clipgrid <- gIntersection(grid, polyg, byid = T, id = as.character(grid@data$ID))
cells <- row.names(clipgrid)
cells <- split(cells, ' ')
clipspdf <- as(clipgrid, 'SpatialPolygonsDataFrame')
clipspdf@data$id <- as.numeric(row.names(clipspdf@data))
spplot(clipspdf, 'id')
summary(clipspdf@data)
dummy id
Min. :0 Min. :665.0
1st Qu.:0 1st Qu.:687.8
Median :0 Median :706.5
Mean :0 Mean :706.4
3rd Qu.:0 3rd Qu.:725.2
Max. :0 Max. :747.0
于 2014-03-12T23:37:43.627 回答
1
以下 R 代码有效（请注意,您的示例中缺少），运行此代码所需的一些对象在此答案的末尾给出：
meuse.grid[!is.na(over(meuse.grid, sr)),]
如果这不能解决您的问题，请提供一个可重现的示例来说明该问题。
需要的一些对象：
r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))
sr1=Polygons(list(Polygon(r1)),"r1")
sr2=Polygons(list(Polygon(r2)),"r2")
sr3=Polygons(list(Polygon(r3)),"r3")
sr4=Polygons(list(Polygon(r4)),"r4")
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), row.names=c("r1","r2","r3","r4")))
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y
于 2013-07-31T08:48:43.423 回答
0
这是我在赏金之后想出的：
require(sp)
require(rgdal)
grid = GridTopology(cellcentre.offset=c(731888.0,7457552.0),cellsize=c(16,16),cells.dim=c(122,106))
# convert object of class SpatialGrid to SpatialPixelsDataFrame.
grid = SpatialPixelsDataFrame(grid,
data=data.frame(id=1:prod(122,106)),
proj4string=CRS(proj4string(shp)))
plot(grid)
shp = readOGR(dsn = "...", layer = "shp")
bound = shp@polygons
bound = SpatialPolygons(bound, proj4string=CRS("+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs"))
plot(bound)
归功于 Paul Hiemstra的回答。
clip_grid = grid[!is.na(over(grid, bound)),]
plot(clip_grid)
于 2014-03-13T18:26:35.737 回答