在一个shapefile中创建网格

Question

我将尝试在一个shapefile中创建一个网格，类似于这样。然而，我无法生成这样的网格。我想知道是否有人知道如何实现这一目标。

这是我的代码 -

WWWL.Shape<- readOGR("E:/Juan Arango", "WWL_Commerce_OK")
WWWL.Shape
plot(WWWL.Shape)
proj4string(WWWL.Shape)

bb <- bbox(WWWL.Shape)
cs <- c(3.28084, 3.28084)*6000  # cell size 
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(WWWL.Shape)))
plot(sp_grd)

< p> WWL.Shape 的输出

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 334367, 334498.7, 4088915, 4089057  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : Id 
min values  :  0 
max values  :  0 

- Perro

---

“grid”是什么意思？是点阵还是多边形网格？ - loki

Answer 1

sf 版本

请参考下面的 rgdal 版本

首先，我们从一个 shapefile 开始。您可以使用 st_read 从任何地理空间文件中加载它。

library(sf)
library(raster)
library(ggplot2)

# load some spatial data. Administrative Boundary
shp <- getData('GADM', country = 'aut', level = 0)
shp <- st_as_sf(shp)
# ggplot() + 
  # geom_sf(data = shp)

现在你需要的是将st_make_grid和st_intersection组合起来：

grid <- shp %>% 
  st_make_grid(cellsize = 0.1, what = "centers") %>% # grid of points
  st_intersection(shp)                               # only within the polygon

# ggplot() + 
  # geom_sf(data = shp) + 
  # geom_sf(data = grid)

rgdal和版本

要创建像素网格，您可以使用sp :: makegrid函数。
让我们从一个可重复的示例开始：

library(raster)
shp <- getData(country = "FRA", level = 0)

现在我们有了一个（多）多边形。让我们将其转换为度量坐标系（因为您的数据和单元格大小也是度量单位）：

shp <- spTransform(shp, CRSobj = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(shp)

接下来，我们会在这个多边形内创建一个网格，根据您指定的单元格大小。

cs <- c(3.28084, 3.28084)*6000
grdpts <- makegrid(shp, cellsize = cs)

然后，我们将这个网格（基本上是中心点的矩阵）转换为一个SpatialPoints对象：

spgrd <- SpatialPoints(grdpts, proj4string = CRS(proj4string(shp)))

这样就可以将其转换为一个SpatialPixels对象。（注意：添加子集[shp，]只选择原始多边形内的点。）

spgrdWithin <- SpatialPixels(spgrd[shp,])
plot(spgrdWithin, add = T)

如果您需要将网格转换为多边形或网格，可以使用：

spgrdWithin <- as(spgrdWithin, "SpatialPolygons")
# or
spgrdWithin <- as(spgrdWithin, "SpatialGrid")

- loki