Sunday, 8 April 2018

Merge a list of spatial polygon objects in R


I have a list of spatial buffers (30000 buffers) that I built with the function lapply:


buff.pts <- lapply(1:nrow(pts.prj), FUN=function(l){
buff <- gBuffer(pts.prj[l,], width=1000) ## 1km
return(buff)
}))

> head(buff.pts)

[[1]]
class : SpatialPolygons
features : 1
extent : 307941.8, 311941.8, 4994518, 4998518 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0

[[2]]
class : SpatialPolygons
features : 1
extent : 307226, 311226, 4991153, 4995153 (xmin, xmax, ymin, ymax)

coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0

From this list, how can I merge all spatial buffers to obtain a shapefile with the 30000 buffers (or features) ? (This shapefile will then be used in the function aggregate to aggregate spatial polygons by attributes.)


I tested this code but I obtain this error message:


test <- as.data.frame(do.call("rbind", buff.pts))
Error in as.data.frame(do.call("rbind", buff.pts)) :
error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': Error in validObject(res) :
invalid class “SpatialPolygons” object: non-unique Polygons ID slot values

Answer



Given a list of SpatialPolygons objects, here's how to construct a spatial polygons data frame with one feature per original SpatialPolygons feature.



Sample data: spl is a list of 12 SpatialPolygons objects - make sure your object gives the same results as this, and test on a small sample before running on 30,000:


> length(spl)
[1] 12
> class(spl)
[1] "list"
> class(spl[[1]])
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"


You want to create a single Spatial Polygons object with all the features in it in order to then make a Spatial Polygons Data Frame:


> joined = SpatialPolygons(lapply(spl, function(x){x@polygons[[1]]}))
> plot(joined)

This takes the first polygons slot from the object (and there should be only one, since each list element is currently a single feature) and then constructs a list of Polygons objects which is what you feed to SpatialPolygons to make a multi-feature SpatialPolygons. Plot this, and you should see all your features. Next, if you want to save as a shapefile, you need to add some data. In the absence of anything else, I create a simple 1 to 12 ID column:


> jdata = SpatialPolygonsDataFrame(Sr=joined, data=data.frame(i=1:12),FALSE)

The FALSE flag just stops R trying to rearrange the spatial and non-spatial data to match up. You might want to put the buffer sizes in the data frame or something.


Job done.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...