Merge two polygon regions into a single polygon region in R

Question

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

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).

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!

Answer 1

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

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()

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)

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()

Now the circle is.