[R-sig-Geo] unionSpatialPolygons - single polygon instead of multipart polygon

[Roger Bivand]

On Sun, 14 Sep 2008, hadley wickham wrote:

> Hi All,
>
> I'm trying to figure out exactly what output unionSpatialPolygons
> gives me - it looks like it merges my singleton polygons into
> multi-polygons (see output of str below), but what I really want is
> just a single singleton polygon, containing the boundary of the
> unioned areas.  How can I get that?

By defining the IDs= argument (as you do) for the output Polygons objects. 
Could you make your test file available - there may be slivers there 
making it look as though county boundaries are identical, but are 
separated by a river? See:

> nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
  proj4string=CRS("+proj=longlat +datum=NAD27"))
lps <- coordinates(nc1)
ID <- cut(lps[,1], quantile(lps[,1]), include.lowest=TRUE)
reg4 <- unionSpatialPolygons(nc1, ID)
length(slot(reg4, "polygons"))
# [1] 4
sapply(slot(reg4, "polygons"), function(x) length(slot(x, "Polygons")))
# [1] 1 1 1 6

where the islands end up as separate Polygon objects in the 4th Polygons 
object, because they don't share a boundary. If your dataset has slivers, 
that might be avoided using the threshold argument, but I wouldn't 
describe that as convenient, really.

Hope this helps,

Roger

>
>> str(unioned[4]@polygons[[1]])
> Formal class 'Polygons' [package "sp"] with 5 slots
>  ..@ Polygons :List of 3
>  .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
>  .. .. .. ..@ labpt  : num [1:2] -99.6  28.6
>  .. .. .. ..@ area   : num 0.063
>  .. .. .. ..@ hole   : logi FALSE
>  .. .. .. ..@ ringDir: int 1
>  .. .. .. ..@ coords : num [1:151, 1:2] -99.8 -99.7 -99.7 -99.6 -99.6 ...
>  .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
>  .. .. .. ..@ labpt  : num [1:2] -98.2  28.5
>  .. .. .. ..@ area   : num 0.112
>  .. .. .. ..@ hole   : logi FALSE
>  .. .. .. ..@ ringDir: int 1
>  .. .. .. ..@ coords : num [1:131, 1:2] -98 -98 -98 -98 -98 ...
>  .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
>  .. .. .. ..@ labpt  : num [1:2] -96.7  28.5
>  .. .. .. ..@ area   : num 0.0558
>  .. .. .. ..@ hole   : logi FALSE
>  .. .. .. ..@ ringDir: int 1
>  .. .. .. ..@ coords : num [1:190, 1:2] -96.8 -96.8 -96.8 -96.8 -96.8 ...
>  ..@ plotOrder: int [1:3] 2 1 3
>  ..@ labpt    : num [1:2] -98.2  28.5
>  ..@ ID       : chr "007"
>  ..@ area     : num 0.231
>
>
> I've included my full code below to give more of a flavour of what I'm
> trying to do - simplify the boundaries of TX counties and then get the
> data into a regular R data frame.  This doesn't work at the moment,
> because I'm only getting the first part of the multipart polygons.
>
> counties <- readShapeSpatial("tx-counties")
>
> attr <- as.data.frame(counties)
> names(attr) <- tolower(names(attr))
> attr <- attr[c("fips", "name", "area", "perimeter")]
>
> polys <- split(row.names(attr), attr$fips)
>
> cp <- polygons(counties)
> unioned <- unionSpatialPolygons(cp, invert(polys))
>
> coords <- function(x) x at polygons[[1]]@Polygons[[1]]@coords
> ccoords <- lapply(seq_len(254), function(i) coords(unioned[i]))
>
> id <- function(x) x at polygons[[1]]@ID
> names(ccoords) <-  sapply(seq_len(254), function(i) id(unioned[i]))
>
> cdf <- do.call("rbind", lapply(seq_along(ccoords), function(i) {
>  df <- as.data.frame(ccoords[[i]])
>  names(df) <- c("x", "y")
>  df <- as.data.frame(dp(df, 0.01))
>  df$id <- names(ccoords)[i]
>  df
> }))
>
>
> Hadley
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no