将两个多边形区域合并为R中的单个多边形区域 [英] Merge two polygon regions into a single polygon region in R

问题描述

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

I'm new to working with spatial data and polygons in R.

我有两个分别从Google Earth中提取的两个多边形的形状文件.因此，基本上第一个形状文件是一个位置(例如购物中心等)，第二个形状文件是第一个位置周围三公里的半径.我将两个形状文件都读取为R作为SpatialPolygonsDataFrames.我使用以下代码:

I have two separate shape files of two polygons that I extract from Google Earth. So basically the first shape file is a location (such as a shopping mall etc.) and the second shape file is a three kilometre radius around the first location. I read both shape files into R as SpatialPolygonsDataFrames. I use the following code:

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)

但是，我没有看到两个多边形的组合视图(一个多边形在另一个多边形内的视图).

However, I don't get a combined view of the two polygons, (view of the one polygon within the other).

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

If anybody have some advice on how I can code this the merge of two polygon regions into a single polygon region in R, I'd appreciate it very much!

推荐答案

如果您承诺使用spatstat包，这可能不会很有帮助.如果没有，请继续阅读...

If you're committed to using thespatstat package, this probably will not be very helpful. If not, read on...

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

If all you want do to is plot the polygons as separate layers, use 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()

如果您需要将它们合并为一个spacearyPolygonDataFrame，请使用此方法.这里的细微差别是，如果两层具有相同的多边形ID，则不能使用spRbind(...).因此，调用spChFIDs(...)会将poly2(圆圈)中一个多边形的ID更改为"R.3km".

If you need to combine them into one spatialPolygonDataFrame, use this. The nuance here is that you can't use spRbind(...) if the two layers have common polygon IDs. So the call to spChFIDs(...) changes the ID of the one polygon in poly2 (the circle) to "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)

您会注意到，在这些图中，圆圈"不是.这是因为我们使用的是long/lat，而您位于赤道以南.数据需要重新投影.事实证明，适用于南非的CRS是 utm-33 .

You'll notice that, on these plots, the "circle" isn't. This is because we're using long/lat and you're pretty far south of the equator. The data needs to be reprojected. Turns out the appropriate CRS for South Africa is 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()

现在是圆圈.

这篇关于将两个多边形区域合并为R中的单个多边形区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！