6

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

enter image description here

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

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.

  • Thank you. I obtained the intersection using: both.polys <- gIntersection(poly.df1, poly.df2, byid=TRUE). Unfortunately, the attribute data have been lost. Ideally, each of the new cells will retain the attribute assigned to that area in each of the two original grids. – Mark Miller Mar 23 '14 at 18:33
4

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

enter image description here

> 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
  • I do not understand what was wrong with my original solution. The code ran and appeared to return the correct answer. The code returned errors after being edited by a different user. I believe I have edited the edited answer so that the code again returns the correct answer. The answer provided now seems to be the same as the answer provided originally, but the code now uses a couple of names provided by the other user. – Mark Miller Oct 17 '14 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.