SpatialPolygons - 在 R 中从坐标创建一组多边形

Posted 2023-03-17

【中文标题】SpatialPolygons - 在 R 中从坐标创建一组多边形

【英文标题】：SpatialPolygons - Creating a set of polygons in R from coordinates

【发布时间】：2014-12-24 13:38:30

【问题描述】：

我正在尝试从顶点位置创建一组多边形，以 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，因此我的数据需要在列表中。所以我创建了一个循环来尝试将我的数据从矩阵中转换为列表格式。

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

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 对象中的两个多边形，第一个具有 ID SJER1，第二个具有 ID SJER2 存储在 SpatialPolygons 对象中。

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

【参考方案1】：

?'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 坐标的2 列矩阵）。这些组合在列表中以创建 Polygons 对象（每个对象都应具有唯一的 ID）。这些Polygons 对象组合在一个列表中以创建SpatialPolygons 对象。您可以根据需要添加CRS，将proj4string 参数添加到SpatialPolygons（请参阅?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))

【讨论】：

非常感谢。我想我更了解这一点，但你能告诉我为什么：for (i in 1:2) polys <- SpatialPolygons(list( Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i]) )) 或 a <- vector('list', length(2)) for (i in 1:2) a[[i]]<-list( Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i]) ) SpatialPolygons(a) 不起作用吗？我需要在很多点上运行它，所以我需要某种循环。非常感谢您的时间和解释。我也确实查看了文档。

@LWasser 第一个每次循环都会覆盖polys，所以你最终会得到一个只包含最后一个多边形的SpatialPolygons 对象。如果您删除最外层的list()，第二个将起作用，即只需执行a[[i]] <- Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i])。这是因为a 已经是一个列表。

@LWasser，如果您有很多行，请参阅我的编辑以了解可用于自动化流程的方法。 （另外，很抱歉，我不知何故忘记包含创建 SPDF 的代码。）

@jmaums 非常感谢。我现在明白这是如何工作的。该代码完美运行。我有一个 github 页面，我一直在用它来记录这些东西，因为我弄清楚它们/获得帮助弄清楚它们或者当我从其他人那里获取代码时。我非常感谢你，因为我真诚地感谢你花时间解释/帮助我 -- Git hub URL to code

真可惜sp 没有内置函数来更简洁地处理这个问题