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