[R-sig-Geo] How does sf do rgeos::gUnaryUnion or maptools::unionSpatialPolygons?

[Edzer Pebesma]

On 09/05/17 15:32, Roger Bivand wrote:
> While https://github.com/edzer/sfr/wiki/migrating is very helpful,
> rgeos::gUnaryUnion is not the typical use case of sf::st_union. The
> typical use case, from ASDAR 1st edition, chapter 5 and the maptools
> "combine_maptools" vignette (the shapefiles are shipped with maptools),
> is grouping features that should belong to the same statistical entity:
> 
> library(maptools)
> vignette("combine_maptools")
> 
> ###################################################
> ### chunk number 16:
> ###################################################
> 
> library(sf)
> nc90 <- st_read(system.file("shapes/co37_d90.shp", package = "maptools"))
> st_crs(nc90) <- "+proj=longlat +datum=NAD27"
> table(table(paste(nc90$ST, nc90$CO, sep="")))
> 
> The point is that two counties are represented by multiple features of
> "POLYGON" objects, rather than single "MULTIPOLYGON" objects, in the
> input shapefile. I've tried doing:
> 
> ids <- factor(paste(nc90$ST, nc90$CO, sep=""))
> nc90a <- st_cast(nc90, to="MULTIPOLYGON", ids=as.integer(ids))
> 
> but:
> 
>> dim(nc90a)
> [1] 104   9
> 
> all turned into "MULTIPOLYGON", but not grouped, even though I think
> ids= are as they should be:
> 
>> table(table(as.integer(ids)))
> 
>  1  2  4
> 98  1  1

this is at least documented: nc90 is of class `sf`, and st_cast.sf
has no ids argument; ... is ignored. If it would merge, it'd need
some guidance what to do with non-geometry feature attributes.

> 
> This may be to avoid dropping data.frame rows. It looks as though I can
> get there using:
> 
> nc90a <- st_cast(st_geometry(nc90), to="MULTIPOLYGON", ids=as.integer(ids))
> 
>> length(nc90a)
> [1] 100
> 
> but all are "MULTIPOLYGON", not a mixture of "POLYGON" and
> "MULTIPOLYGON" features, and I've no idea which is which - that is how
> the order of nc90a relates to the counties of nc90. How do I associate
> the features in nc90a with their county ids? Where do i find the ids I
> gave to st_cast in nc90a?

st_cast uses base::split, which converts ids to a factor,
and returns a list with the order of the levels of that factor
(alphabetically?). Neither convenient, nor documented. I guess more
intuitive would be to squeeze, but keep original order, right?

I'd suggest to use instead

nc90a = aggregate(nc90, list(ids = ids), head, n = 1)

which returns a GEOMETRY sf, with 2 MULTIPOLYGON and 98 POLYGON. You
could st_cast that to MULTIPOLYGON.

Alternatively:

library(dplyr)
bind_cols(nc90, ids=ids) %>%
	group_by(ids) %>%
	summarise(ST=head(ST,1), CO=head(CO,1), do_union=FALSE)

converts straight into MULTIPOLYGON.

I'll update the sp -> sf migration wiki.

> 
> What am I missing? I'm writing here rather than raising an issue on GH
> because others may want to know too, and Edzer needs others to reply for
> him if they know the answer.
> 
> Puzzled,
> 
> Roger
> 

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20170509/aef97118/attachment.sig>