Joining grid or polygon layers in R?

Question

I have two grids. Cells in each grid have an attribute. Size of grid cells differs between the two grids and the two grids do not overlay each other perfectly. I would like to create a new layer, or group of polygons, or grid, that results from joining the two existing grids. Each cell in the new grid should retain the attribute from each of the original grids. Area not occurring in both original grids can be discarded from the new grid. Such discarded areas would only have one attribute (and one missing observation for the second attribute).

How can I do this is R?

So far I have tried over in the sp package and unionSpatialPolygons. I have also considered the PBSmapping package, which allows for creating intersections of polygons. This post describes how to join polygons within a layer, but perhaps not how to join polygons in different layers: Joining polygons in R

I have written code to compute the (x,y) coordinates of the four corners of each cell that results from joining the two original grids. That code is at the very bottom of this post. Perhaps I can create a new polygon layer from the resulting data frame. However, I suspect that process has already been implemented in a package.

Below is code to create the two original grids, assign an attribute to each cell in each grid and plot the two grids.

There will be 16 cells in the new grid. Those cells generally will not be squares, but rather rectangles of various shapes and sizes.

library(rgdal)
library(maptools)
library(raster)
library(rgeos)
library(gridExtra)

setwd('c:/users/mmiller21/simple R programs')

set.seed(1234)

# create grid 1

grd1 <- GridTopology(c(2,2), c(1,1), c(4,4))
polys1 <- as(grd1, "SpatialPolygons")

# assign projection to grid

proj4string(polys1) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

# create fake atttribute data for each grid cell

poly.data1 = data.frame(f1=runif(length(row.names(polys1)), 1, 10))
row.names(poly.data1) <- paste0('g', 1:length(row.names(polys1)))

# convert grid to a SpatialPolygonsDataFrame:

poly.df1 = SpatialPolygonsDataFrame(polys1, poly.data1)

# create grid 2

grd2 <- GridTopology(c(3.8,3.8), c(2,2), c(2,2))
polys2 <- as(grd2, "SpatialPolygons")

# assign projection to grid

proj4string(polys2) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

# create fake atttribute data for each grid cell

poly.data2 = data.frame(f2=runif(length(row.names(polys2)), 100, 200))
row.names(poly.data2) <- paste0('g', 1:length(row.names(polys2)))

# convert grid to a SpatialPolygonsDataFrame:

poly.df2 = SpatialPolygonsDataFrame(polys2, poly.data2)

jpeg(filename = "joined.grid.jpeg")

plot(poly.df2, col='blue', xlim=c(1,8), ylim=c(1,8))
text(coordinates(poly.df2), labels=row.names(poly.df2))
axis(1, lty = 0, lwd = 0)
axis(2, lty = 0, lwd = 0)
plot(poly.df1, add=TRUE)
text(coordinates(poly.df1), labels=row.names(poly.df1))

dev.off()

Here is code to compute the (x,y) coordinates of each cell in the new grid:

x1 <- seq(1.5,5.5,1)
y1 <- seq(1.5,5.5,1)

x2 <- seq(2.8,6.8,2)
y2 <- seq(2.8,6.8,2)

x <- c(x1,x2)
x <- sort(x)
x
# [1] 1.0 2.0 2.8 3.0 3.8 4.0 4.8 5.0 5.8 6.8

new.x <- x[x>= min(x2) & x <= max(x1)]

y <- c(y1,y2)
y <- sort(y)
# [1] 1.0 2.0 2.8 3.0 3.8 4.0 4.8 5.0 5.8 6.8

new.y <- y[y>= min(y2) & y <= max(y1)]

z <- merge(new.x, new.y)

plot(z$x, z$y)

Answer 1

I used gintersection, as Jeffrey Evans suggested, and then I think I was able to assign the attribute from each original grid to the new grid with the following code:

both.polys <- gIntersection(poly.df1, poly.df2,  byid=TRUE)
both.polys
#class       : SpatialPolygons 
#features    : 16 
#extent      : 2.8, 5.5, 2.8, 5.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0

row.names(both.polys)

#  [1] "g2 g1"  "g2 g3"  "g3 g1"  "g3 g3"  "g4 g1"  "g4 g2"  "g4 g3"  "g4 g4"  "g6 g3"  "g7 g3"  "g8 g3"  "g8 g4" 
# [13] "g10 g3" "g11 g3" "g12 g3" "g12 g4"

row.names(poly.df1)
# [1] "g1"  "g2"  "g3"  "g4"  "g5"  "g6"  "g7"  "g8"  "g9"  "g10" "g11" "g12" "g13" "g14" "g15" "g16"

row.names(poly.df2)
# [1] "g1" "g2" "g3" "g4"

poly.df1$f1
#  [1] 2.023331 6.600695 6.483473 6.610415 8.748238 6.762795 1.085462 3.092955 6.994754 5.628260 7.242322 5.904774
# [13] 3.544602 9.310901 3.630843 8.535661

poly.df2$f2
# [1] 128.6223 126.6821 118.6723 123.2226

new.attribs <- data.frame(do.call(rbind, strsplit(row.names(both.polys), " ")))

poly.df1$X1 <- row.names(poly.df1)
poly.df2$X2 <- row.names(poly.df2)

new.attrib.data <- merge(new.attribs    , poly.df1, by='X1')
new.attrib.data <- merge(new.attrib.data, poly.df2, by='X2')

row.names(new.attrib.data) <- row.names(both.polys)

# convert new grid to a SpatialPolygonsDataFrame:

new.grid = SpatialPolygonsDataFrame(both.polys, new.attrib.data)

jpeg(filename = "new.grid.jpeg")
plot(new.grid)
dev.off()

> new.grid
class       : SpatialPolygonsDataFrame 
features    : 16 
extent      : 2.8, 5.5, 2.8, 5.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 4
names       : X2,  X1,       f1,       f2 
min values  : g1, g10, 1.085462, 118.6723 
max values  : g4,  g8, 7.242322, 128.6223

Answer 2

The gUnion or gIntersection function in "rgeos" should do the trick. The rgeos package is an interface to the GEOS Geometry Engine and is becoming the de facto R standard for handling topology operations.