我在R中有一个SpatialPointsDataFrame,看起来像这样:
coordinates id order hole piece group box_id
326 (-94.4, 27.6586) 47 1 FALSE 1 47.1 1
327 (-93.64, 27.6232) 47 2 FALSE 1 47.1 1
328 (-92.04, 27.7649) 47 3 FALSE 1 47.1 1
329 (-90.36, 27.0903) 47 4 FALSE 1 47.1 1
330 (-91.12, 25.6929) 47 5 FALSE 1 47.1 1
331 (-92.92, 25.6569) 47 6 FALSE 1 47.1 1
332 (-93.44, 26.0169) 47 7 FALSE 1 47.1 1
333 (-94.4, 25.9809) 47 8 FALSE 1 47.1 1
334 (-94.4, 27.6586) 47 9 FALSE 1 47.1 1
335 (-92.04, 27.7649) 48 1 FALSE 1 48.1 2
336 (-93.64, 27.6232) 48 2 FALSE 1 48.1 2
337 (-94.4, 27.6586) 48 3 FALSE 1 48.1 2
338 (-94.4, 27.8356) 48 4 FALSE 1 48.1 2
339 (-93.64, 27.7649) 48 5 FALSE 1 48.1 2
340 (-90.28, 28.1182) 48 6 FALSE 1 48.1 2
341 (-90.56, 27.9417) 48 7 FALSE 1 48.1 2
342 (-92.04, 27.7649) 48 8 FALSE 1 48.1 2
100 (-94.4, 27.8356) 20 1 FALSE 1 20.1 3
101 (-94.4, 28.0829) 20 2 FALSE 1 20.1 3
102 (-90.28, 28.1182) 20 3 FALSE 1 20.1 3
103 (-93.64, 27.7649) 20 4 FALSE 1 20.1 3
104 (-94.4, 27.8356) 20 5 FALSE 1 20.1 3
(行名/编号不正确,因为我按box_id列进行了排序)
这些点是多边形的节点(由box_id标识)。我想将其编写为多边形形状文件以读入GIS程序,但是在将其转换为SpatialPolygonDataFrame(然后使用writeOGR)或将其直接写入shp文件时遇到了麻烦。任何帮助将不胜感激。我是使用R的GIS的新手,如果我遗漏了一些明显的内容,我们深表歉意。
最佳答案
因此,这是另一种解决方案,在概念上类似于@ JoshO'Brien's。这将容纳子多边形(给定ID的多个组),并创建一个SpatialPolygonsDataFrame类的对象,其中包含数据部分(例如,属性表)。这可以完全导出到ESRI shapefile。
首先,作为示例,我们创建一个SpatialPointsDataFrame,其格式与您的格式相同。 id
列标识多边形ID
,而group
列标识给定多边形的子多边形。在此示例中,我们有3个多边形,其中两个具有1个子多边形,第三个具有3个子多边形。我们还创建了一个数据框data
(属性表),在这种情况下,它将多边形ID与box_id相关联。
# this code just creates a sample SpatialPointsDataFrame
x <- c(-10,-10,10,10,-10)
y <- c(-10,10,10,-10,-10)
df.1 <- data.frame(x,y,id=47, order=1:5,hole=F,piece=1,group=47.1,box_id=1)
coordinates(df.1)=c("x","y")
x <- c(+15+3*cos(2*pi/5*(0:5)))
y <- c(-15+3*sin(2*pi/5*(0:5)))
df.2 <- data.frame(x,y,id=48, order=1:6,hole=F,piece=1,group=48.1,box_id=2)
coordinates(df.2)=c("x","y")
x <- c(-15+2*cos(2*pi/8*(0:8)))
y <- c(+15+2*sin(2*pi/8*(0:8)))
df.3.1 <- data.frame(x,y,id=20, order=1:9,hole=F,piece=1,group=20.1,box_id=3)
coordinates(df.3.1)=c("x","y")
x <- c(0+2*cos(2*pi/8*(0:8)))
y <- c(+15+2*sin(2*pi/8*(0:8)))
df.3.2 <- data.frame(x,y,id=20, order=1:9,hole=F,piece=1,group=20.2,box_id=3)
coordinates(df.3.2)=c("x","y")
x <- c(+15+2*cos(2*pi/8*(0:8)))
y <- c(+15+2*sin(2*pi/8*(0:8)))
df.3.3 <- data.frame(x,y,id=20, order=1:9,hole=F,piece=1,group=20.3,box_id=3)
coordinates(df.3.3)=c("x","y")
df <- rbind(df.1,df.2,df.3.1,df.3.2,df.3.3)
df
# coordinates id order hole piece group box_id
# 1 (-10, -10) 47 1 FALSE 1 47.1 1
# 2 (-10, 10) 47 2 FALSE 1 47.1 1
# 3 (10, 10) 47 3 FALSE 1 47.1 1
# 4 (10, -10) 47 4 FALSE 1 47.1 1
# 5 (-10, -10) 47 5 FALSE 1 47.1 1
# 6 (18, -15) 48 1 FALSE 1 48.1 2
# 7 (15.92705, -12.14683) 48 2 FALSE 1 48.1 2
# 8 (12.57295, -13.23664) 48 3 FALSE 1 48.1 2
# ...
data <- data.frame(box_id=unique(df$box_id),row.names=unique(df$id))
现在,这是
points2polygons(...)
的代码。使用现有的SpatialPointsDataFrame,这就是您所需要的。points2polygons <- function(df,data) {
get.grpPoly <- function(group,ID,df) {
Polygon(coordinates(df[df$id==ID & df$group==group,]))
}
get.spPoly <- function(ID,df) {
Polygons(lapply(unique(df[df$id==ID,]$group),get.grpPoly,ID,df),ID)
}
spPolygons <- SpatialPolygons(lapply(unique(df$id),get.spPoly,df))
SpatialPolygonsDataFrame(spPolygons,match.ID=T,data=data)
}
spDF <- points2polygons(df,data)
plot(spDF,col=spDF$box_id+1)
library(rgdal)
writeOGR(spDF,dsn=".",layer="myShapefile", driver="ESRI Shapefile")
关于r - 在R中将SpatialPointsDataframe转换为SpatialPolygons,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/21759134/