SpatialPolygons-根据坐标在R中创建一组面 [英] SpatialPolygons - Creating a set of polygons in R from coordinates

问题描述

我正在尝试从顶点位置创建一组多边形并以X、Y格式保存。

这里是我的数据的一个例子--每行代表一个多边形的顶点。这些多边形是正方形

square <- rbind(c(255842.4, 4111578, 255862.4, 4111578, 255862.4, 4111558, 
                  255842.4, 4111558, 255842.4, 4111578, 255842.4, 4111578),
                c(257397.0, 4111309, 257417.0, 4111309, 257417.0, 4111289, 
                  257397.0, 4111289, 257397.0, 4111309, 257397.0, 4111309))

ID <- c("SJER1", "SJER2")'

我使用的是SpatialPolygons，因此我的数据需要在列表中。因此，我创建了一个循环，试图将数据从矩阵转换为列表格式。

我在这个站点上的其他一些问题中找到的代码后面创建了一个循环。我打破了每一步，试图理解为什么我只得到一个多边形作为输出，即使我有两组点。

for (i in 1:2) {  
  pts <- rbind(c(square[i,1], square[i,2]), c(square[i,3], square[i,4]), 
               c(square[i,5],square[i,6]), c(square[i,7],square[i,8]), 
               c(square[i,9],square[i,10]))
  sp1 <- list(Polygon(pts))
  sp2 <- list(Polygons(sp1,i))
  sp = SpatialPolygons(sp2)  
}
plot(sp)

您能告诉我如何调整代码以写出两个多边形而不是只写出一个吗？还有，如果我使用一个矩阵(正方形)作为我的起始数据集，并且如果我分配一个字符ID，它会将我的所有数据转换为一个字符，那么我如何将ID分配给每个多边形。

我的最终目标是SpatialPolygons对象中的两个多边形，第一个具有IDSJER1，第二个具有存储在SpatialPolygons对象中的IDSJER2。

然后我将把它写出到一个shapefile中。

推荐答案

?'SpatialPolygons-class'中有一些信息，但您或多或少希望执行以下操作：

polys <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

plot(polys)

基本要点是您需要创建Polygon对象(例如，从两列矩阵创建，x坐标在第一列，y坐标在第二列)。它们组合在列表中以创建Polygons对象(每个对象都应该具有唯一的ID)。这些Polygons对象组合在一个列表中以创建SpatialPolygons对象。如果愿意，您可以向SpatialPolygons添加CRS，参数为proj4string(请参阅?SpatialPolygons)。

要将其写出到ESRI shapefile，您需要通过组合我们创建的polys对象和一些数据将其转换为SpatialPolygonsDataFrame对象。我们只是将ID作为数据添加，因为没有更有趣的东西。

polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

然后写出来...

writeOGR(polys.df, '.', 'fancysquares', 'ESRI Shapefile')

第二个参数('.')表示将其写出到当前工作目录。

---

编辑

若要在有多行描述面的情况下快速创建SpatialPolygonsDataFrame，可以使用以下命令：

# Example data
square <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square)))

# Create SP
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID))

# Create SPDF
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

plot(polys.df, col=rainbow(50, alpha=0.5))

这篇关于SpatialPolygons-根据坐标在R中创建一组面的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！