3

我是在 R 中处理空间数据和多边形的新手。

我有两个从谷歌地球提取的两个多边形的单独形状文件。所以基本上第一个形状文件是一个位置（例如购物中心等），第二个形状文件是第一个位置周围三公里的半径。我将两个形状文件作为 SpatialPolygonsDataFrames 读入 R。我使用以下代码：

library(maptools)
library(sp)
library(spatstat)
options(digits=10) 

# Read polygon a

a <- readShapeSpatial(file.choose())
class(a)

spatstat.options(checkpolygons=FALSE)

r <- slot(a,"polygons")
r <- lapply(r, function(a) { SpatialPolygons(list(a)) })
windows <- lapply(r, as.owin)
Ploy_One <- tess(tiles=windows)

# Read polygon b

b <- readShapeSpatial(file.choose())
class(b)

spatstat.options(checkpolygons=FALSE)

s <- slot(b,"polygons")
s <- lapply(s, function(b) { SpatialPolygons(list(b)) })

windows <- lapply(s, as.owin)
Poly_Two <- tess(tiles=windows)

# Read polygon b

Combined_Region <- intersect.tess(Poly_One, Poly_Two)
plot(Combined_Region)

但是，我没有得到两个多边形的组合视图（一个多边形内的一个多边形的视图）。

如果有人对我如何将两个多边形区域合并为 R 中的单个多边形区域进行编码有一些建议，我将非常感激！

4

1 回答 1

4

如果您致力于使用该spatstat软件包，这可能不会很有帮助。如果没有，请继续阅读...

如果您只想将多边形绘制为单独的图层，请使用ggplot

library(ggplot2)
library(sp)
library(maptools)

setwd("<directory with all your files...>")

poly1 <- readShapeSpatial("Polygon_One")
poly2 <- readShapeSpatial("Polygon_Two")
# plot polygons as separate layers,,,
poly1.df <- fortify(poly1)
poly2.df <- fortify(poly2)
ggplot() +
  geom_path(data=poly1, aes(x=long,y=lat, group=group))+
  geom_path(data=poly2, aes(x=long,y=lat, group=group), colour="red")+
  coord_fixed()

如果您需要将它们组合成一个 spatialPolygonDataFrame，请使用它。spRbind(...)这里的细微差别是，如果两个图层具有共同的多边形 ID ，则不能使用。因此，调用将（圆）spChFIDs(...)中的一个多边形的 ID 更改poly2为“R.3km”。

# combine polygons into a single shapefile
poly.combined <- spRbind(poly1,spChFIDs(poly2,"R.3km"))
# plot polygons using ggplot aesthetic mapping
poly.df <- fortify(poly.combined)
ggplot(poly.df) + 
  geom_path(aes(x=long, y=lat, group=group, color=group)) + 
  scale_color_discrete("",labels=c("Center", "3km Radius")) +
  coord_fixed()
# plot using plot(...) method for spatialObjects
plot(poly.combined)

您会注意到，在这些图中，“圆圈”不是。这是因为我们使用的是经度/纬度，而您位于赤道以南很远的地方。数据需要重新投影。事实证明，适合南非的 CRS 是utm-33。

wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
utm.33 <- "+proj=utm +zone=33 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
proj4string(poly.combined) <- CRS(wgs.84)
poly.utm33 <- spTransform(poly.combined,CRS(utm.33))
poly.df    <- fortify(poly.utm33)
ggplot(poly.df) + 
  geom_path(aes(x=long, y=lat, group=group, color=group)) + 
  scale_color_discrete("",labels=c("Center", "3km Radius")) +
  theme(axis.text=element_blank()) + labs(x=NULL,y=NULL) +
  coord_fixed()

现在是圈子。

于 2014-01-25T23:17:37.147 回答