Roger, Works for me now. Thanks for your patience with this.
I'll definitely look at sf for future projects, but for now, I just want to analyze my darn data! :) - Anne On 05/01/2017 05:23 AM, Roger Bivand wrote: > On Sun, 30 Apr 2017, Anne C. Hanna wrote: > >> Okay, I think I've tracked it down. This is... absurdly complicated and I >> don't know what the actual correct resolution is, but I'm pretty sure that at >> least I know what's going wrong. > > Good, and thanks for staying with this. I've pushed a fix to my github repo, > and pull request #29 to edzer/sp. > > Note that this also works out of the box: > > library(sf) > sf <- st_read(dsn = "2011 Voting District Boundary Shapefiles", layer = > "MCDS", stringsAsFactors = FALSE) > sfd <- st_cast(sf, "POLYGON") > > However, using sf::st_touches() is a bit less intuitive for now than the spdep > "nb" object generated by spdep::poly2nb(), which also provides Queen (touch at > one or more boundary points) and Rook (touch at more than one boundary point), > and sf::st_touches() is seven times slower than spdep::poly2nb() now. > sf::st_touches() could be speeded up using an STR tree (spdep::poly2nb() can > do this, but in any case only looks at candidate neighbours with intersecting > bounding boxes). So: > > sp_sfd <- as(sfd, "Spatial") > library(spdep) > sp_nb <- poly2nb(sp_sfd) > > seems a reasonable workflow to find the adjacencies under Queen/Rook and > snapping control. > > Please try out the disaggregate fix and report back - and consider shifting to > sf. > > Roger > >> >> Basically it comes down to the fact that createSPComment() does not actually >> check whether each individual polygon in the SpatialPolygons object has a >> comment (much less a correct comment). Instead, it looks at a top-level >> "comment" attribute. If that top-level comment attribute (i.e. sf@comment >> for >> my SpatialPolygonsDataFrame object stored as "sf") is "TRUE", then >> createSPComment() does nothing. >> >> And here's where the second layer of the problem comes in for my dataset --- >> my original data, pre-disaggregation, has complete and correct comments. >> Some >> of those comments (for objects where it's a single polygon with 0 or 1 holes) >> are being retained by the basic disaggregation algorithm, while others are >> not. Your patch attempts to fill in the holes by running createSPComment() >> on >> the final results. But... the thing you are running createSPComment() on is >> SpatialPolygons(p), where p is the list of disaggregated polygons. And the >> SpatialPolygons() constructor looks at the list of the polygons you feed it, >> and if *any* of those polygons have comments, it sets the top-level "comment" >> attribute of its output to "TRUE", even if other polygons are missing >> comments. (And it also doesn't check the comments it gets for correctness.) >> So the fact that some polygon comments are retained in p means that >> SpatialPolygons() lies to createSPComment() about whether it has any missing >> comments that need to be fixed. >> >> On the other hand, when I manually run createPolygonComment() on the >> individual Polygons objects in disaggregate()'s output, I don't even look at >> the top-level "comment" attribute, so it works just fine. >> >> My little toy test polygons, on the other hand, were just not complex enough >> to end up with partial comments in them, so they also didn't experience this >> error. >> >> I am not familiar enough with the design philosophy behind sp to suggest at >> what level this issue should be fixed, but I hope this is at least enough >> info >> to help those who do know what they're doing get it corrected. >> >> - Anne >> >> >> On 04/30/2017 12:18 PM, Anne C. Hanna wrote: >>> Roger, >>> >>> Unfortunately I have a use case for you of createSPComment() not working. >>> >>> I tried your new version of disaggregate() on the original test script I >>> sent >>> you, and it does seem to have fixed all the cases in there. However, when I >>> tried it on my actual data, it had the same failure mode as before. The >>> original dataset has what appear to be correct comments, but the >>> disaggregate() output is full of NULLs in place of comments for the >>> multi-part >>> and multi-hole polygons. >>> >>> createSPComment() appears to be what's failing, rather than just some logic >>> in >>> your disaggregate() function --- when I try to run createSPComment() on the >>> disaggregate() output, I still get the same set of NULLs. I can fix it if I >>> run createPolygonsComment() individually on every Polygons object in my >>> SpatialPolygonsDataFrame. >>> >>> I'm not sure what is different between my dataset and the test polygons I >>> was >>> using earlier, as createSPComment() seems to handle the test shapes just >>> fine >>> (and converting my SpatialPolygonsDataFrame into just a set of >>> SpatialPolygons >>> also doesn't help). But presumably it is better to have createSPComment() >>> fixed than to have to work around it in disaggregate(). So I guess I'll >>> look >>> at the code for that and see if I can find what's wrong. >>> >>> Just in case you want to see this in action, I've attached a little test >>> script. The data set I am working with is available at: >>> >>> http://aws.redistricting.state.pa.us/Redistricting/Resources/GISData/2011-Voting-District-Boundary-Shapefiles.zip >>> >>> >>> - Anne >>> >>> >>> On 04/29/2017 10:06 AM, Roger Bivand wrote: >>>> On Sat, 29 Apr 2017, Anne C. Hanna wrote: >>>> >>>>> Roger, >>>>> >>>>> This looks great, and I will try it out ASAP. I do have one reservation >>>>> though --- it seems you are using createSPComment() to reconstruct the >>>>> comments, and I have seen some discussion that that may not be reliable in >>>>> all cases (e.g. if the initial polygons are wonky in some way). I don't >>>>> know a lot about this, but is it possible that it would be preferable to >>>>> parse and appropriately disaggregate the original comments strings (if >>>>> they >>>>> exist), so as to deal slightly more smoothly with such cases (i.e., the >>>>> polygons would still be wonky, but at least the hole/polygon matching >>>>> would >>>>> track with whatever was in the original data)? >>>> >>>> Contributions very welcome - this was code moved from the raster package, >>>> so I >>>> really have little feeling for why it is necessary. >>>> >>>> Please provide references and use cases. GEOS is strict in its treatments >>>> of >>>> geometries, but does work in most cases. People just expressing opinions >>>> doesn't help, and may well be misleading. It is possible to clean >>>> geometries >>>> too, see the cleangeo package. createPolygonsComment probably isn't >>>> foolproof, >>>> but specific reports help, because then it is possible to do something >>>> about >>>> it. sp classes didn't use simple features when written because sf was not >>>> widely used then - introduced in 2004, when sp was being developed. Using >>>> sf >>>> helps, and rgeos::createPolygonsComment was a work-around from 2010 when >>>> rgeos >>>> was written. >>>> >>>> Even when you use sf in the sf package, you will still run into invalid >>>> geometries because the geometries are in fact invalid, the sf package also >>>> uses GEOS. >>>> >>>> Roger >>>> >>>>> >>>>> I also had no success using createSPComment to fix disaggregate()'s output >>>>> previously, even though my polygons are perfectly non-wonky, so I am >>>>> perhaps >>>>> a little more untrusting of it than I should be. But I'll let you know >>>>> how >>>>> this version works with my data. Thanks for addressing this so quickly! >>>>> >>>>> - Anne >>>>> >>>>> >>>>> On 04/28/2017 12:11 PM, Roger Bivand wrote: >>>>>> I've pushed the fix to my fork: >>>>>> >>>>>> https://github.com/rsbivand/sp >>>>>> >>>>>> and created a pull request: >>>>>> >>>>>> https://github.com/edzer/sp/pull/28 >>>>>> >>>>>> Only one part of a complicated set if nested if() in disaggregate() was >>>>>> adding >>>>>> comments, but in some settings the existing comments survived the >>>>>> disaggregation. Now the Polygons object comment attribute is re-created >>>>>> for >>>>>> all Polygons objects. This is version 1.2-6, also including code changes >>>>>> that >>>>>> internally affect rgdal and rgeos - you may like to re-install them from >>>>>> source after installing this sp version (shouldn't matter). >>>>>> >>>>>> Roger >>>>>> >>>>>> On Fri, 28 Apr 2017, Anne C. Hanna wrote: >>>>>> >>>>>>> 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 >>>>>>> >>>>>> >>>>> >>>>> >>>> >> >> >
signature.asc
Description: OpenPGP digital signature
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo