[R-sig-Geo] sp::disaggregate() does not seem to add the necessary "comment" attribute to identify holes if Polygons does not contain exactly 1 polygon with 0 or 1 holes

[Anne C. Hanna]

Hello.  I first posted this issue report on the sp GitHub repo
(https://github.com/edzer/sp/issues/27) and it was suggested that I redirect
it to here.

I am working with a geographic dataset with complex borders. The data is
stored as a SpatialPolygonsDataFrame.  Each Polygons object in the data frame
may be composed of multiple disjoint polygons, and each polygon may have
multiple holes.  I want to disaggregate each of the Polygons objects into its
individual disjoint polygons and construct an adjacency matrix for all the
disjoint components, and I was using disaggregate() to do this.  However, when
I run gTouches() on the disaggregated data, in order to compute the
adjacencies, I get a number of warnings like this:

Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_touches") : Polygons
object missing comment attribute ignoring hole(s).  See function createSPComment.

Looking at the Polygons "comment" attributes in the SpatialPolygonsDataFrame
output by disaggregate(), I see that the only comment values are "0"
(indicating a single polygon with no holes), "0 1" (indicating a single
polygon with a single hole), and NULL (apparently no comment was written).
Since I know my dataset contains several Polygons objects which are composed
of multiple disjoint regions, and also several Polygons which contain more
than one hole, this is not the expected result.  In reading the disaggregate()
code in the sp GitHub repository (specifically, explodePolygons()), I also
can't see anywhere the comment is being added for the cases where a Polygons
object has more than two parts or more than two holes.  It actually seems like
it's getting carried along almost accidentally in the few cases that do get
comments, and neglected otherwise.

Assuming I'm not failing to understand the code and the desired behavior
(entirely possible, as I am new at working with this software!), this seems
suboptimal to me.  My dataset is pretty well-behaved (despite its complexity),
so I should be able to fix my issues with judicious application of
createPolygonsComment.  But I had a heck of a time figuring out what was going
wrong with gTouches, since Polygons comment management appears to be a pretty
obscure field (and createSPComment wasn't working for me, for whatever
reason).  So it seems like it might be better if disaggregate() just parses
and passes along the comments from its input correctly, or, if it's absolutely
necessary to not create comments, passes nothing and warns clearly in the
manual that comments and associated hole information are being lost.  Passing
along comments in some cases while silently dropping them in others seems like
kind of the worst of both worlds.

I've attached a set of tests I wrote to demonstrate desired/undesired
behavior: disaggregate_comment_tests.R.  My R version is 3.4.0, my sp version
is 1.2-4, my rgeos version is 0.3-23 (SVN revision 546), and my GEOS runtime
version is 3.5.1-CAPI-1.9.1 r4246.  I am using Debian Release 9.0 with kernel
version 4.9.0-2-amd64.  I hope this is useful; please let me know if you need
more info or if there is a better place to post this issue.

 - Anne
-------------- next part --------------
require(sp);
require(rgeos);

# Make two disjoint Polygon objects (p1 and p2) and two holes in p1 (p3 and p4), and one hole in p2 (p5)
p1 <- Polygon(cbind(c(0, 1, 0, 0), c(0, 0, 1, 0)), hole = FALSE);
p2 <- Polygon(cbind(c(1, 2, 1, 1), c(1, 1, 2, 1)), hole = FALSE);
p3 <- Polygon(cbind(c(0.25, 0.5, 0.25, 0.25), c(0.25, 0.25, 0.5, 0.25)), hole = TRUE);
p4 <- Polygon(cbind(c(0.0625, 0.125, 0.0625, 0.0625), c(0.0625, 0.0625, 0.125, 0.0625)), hole = TRUE);
p5 <- Polygon(cbind(c(1.25, 1.5, 1.25, 1.25), c(1.25, 1.25, 1.5, 1.25)), hole = TRUE);

# Combine Polygon objects in a few different ways 
P123 <- Polygons(list(p1, p2, p3), 1);
P12 <- Polygons(list(p1, p2), 1);
P13 <- Polygons(list(p1, p3), 1);
P134 <- Polygons(list(p1, p3, p4), 1);
P1 <- Polygons(list(p1), 1);
P1235 <- Polygons(list(p1, p2, p3, p5), 1);

sp123 <- SpatialPolygons(list(P123));
sp12 <- SpatialPolygons(list(P12));
sp13 <- SpatialPolygons(list(P13));
sp134 <- SpatialPolygons(list(P134));
sp1 <- SpatialPolygons(list(P1));
sp1235 <- SpatialPolygons(list(P1235));

# Add hole comments to the SpatialPolygons, so that everything is correct before disaggregation
attr(sp123 at polygons[[1]], "comment") <- createPolygonsComment(sp123 at polygons[[1]]);
attr(sp12 at polygons[[1]], "comment") <- createPolygonsComment(sp12 at polygons[[1]]);
attr(sp13 at polygons[[1]], "comment") <- createPolygonsComment(sp13 at polygons[[1]]);
attr(sp134 at polygons[[1]], "comment") <- createPolygonsComment(sp134 at polygons[[1]]);
attr(sp1 at polygons[[1]], "comment") <- createPolygonsComment(sp1 at polygons[[1]]);
attr(sp1235 at polygons[[1]], "comment") <- createPolygonsComment(sp1235 at polygons[[1]]);

# Disaggregate the SpatialPolygons
spd123 <- disaggregate(sp123, ignoreholes = FALSE);
spd12 <- disaggregate(sp12, ignoreholes = FALSE);
spd13 <- disaggregate(sp13, ignoreholes = FALSE);
spd134 <- disaggregate(sp134, ignoreholes = FALSE);
spd1 <- disaggregate(sp1, ignoreholes = FALSE);
spd1235 <- disaggregate(sp1235, ignoreholes = FALSE);

# Extract the comments from each result
ca123 <- unlist(lapply(lapply(spd123 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));
ca12 <- unlist(lapply(lapply(spd12 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));
ca13 <- unlist(lapply(lapply(spd13 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));
ca134 <- unlist(lapply(lapply(spd134 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));
ca1 <- unlist(lapply(lapply(spd1 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));
ca1235 <- unlist(lapply(lapply(spd1235 at polygons, function(p) {attr(p, "comment")}), function(x) {ifelse(is.null(x), NA, x)}));

# Compute what the comments *should* be for each case
cd123 <- unlist(lapply(spd123 at polygons, function(p) {createPolygonsComment(p)}));
cd12 <- unlist(lapply(spd12 at polygons, function(p) {createPolygonsComment(p)}));
cd13 <- unlist(lapply(spd13 at polygons, function(p) {createPolygonsComment(p)}));
cd134 <- unlist(lapply(spd134 at polygons, function(p) {createPolygonsComment(p)}));
cd1 <- unlist(lapply(spd1 at polygons, function(p) {createPolygonsComment(p)}));
cd1235 <- unlist(lapply(spd1235 at polygons, function(p) {createPolygonsComment(p)}));

# Shape consists of one polygon with a hole and one without
print(ca123);  # This is NA NA (wrong).
print(cd123);  # This is "0 1" "0" (correct).

# Shape consists of two polygons without holes
print(ca12);  # This is NA NA (wrong, but less important).
print(cd12);  # This is "0" "0" (correct).

# Shape consists of one polygon with one hole
print(ca13);  # This is "0 1" (correct).
print(cd13);  # This is "0 1" (correct).

# Shape consists of one polygon with two holes
print(ca134);  # This is NA (wrong).
print(cd134);  # This is "0 1 1" (correct).

# Shape consists of one polygon with no holes
print(ca1);  # This is "0" (correct).
print(cd1);  # This is "0" (correct).

# Shape consists of two polygons, each with one hole
print(ca1235);  # This is NA NA (wrong).
print(cd1235);  # This is "0 1" "0 1" (correct).
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 195 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20170428/842f5e93/attachment.sig>