[R-sig-Geo] sample.Polygons fails for shapes with more than four Polygons

[Roger Bivand]

On Mon, 7 Mar 2011, Roger Bivand wrote:

> On Mon, 7 Mar 2011, Aman Verma wrote:
>
>> Hi everybody,
>> 
>> I have discovered some strange behaviour in the function sample.Polygons in 
>> the sp package. Consider this code:
>> 
>> library(maptools)
>> library(sp)
>> # Download and unpack this shape file to your working directory:
>> # 
>> http://surveillance.mcgill.ca/countries_subdiv_Project_PP_project_region.zip
>> shapefile<-readShapePoly("countries_subdiv_Project_PP_project_region.shp")
>> 
>> # Try and sample some points from different shapes (with multiple 
>> polygons):
>> sample.Polygons(shapefile at polygons[[1930]],n=1,type='regular') # Has 2 
>> polygons
>> # Works just fine.
>> sample.Polygons(shapefile at polygons[[1933]],n=1,type='regular') # Has 10 
>> polygons
>> # Fails:
>> # Error: subscript out of bounds
>> 
>> I have done some debugging, and located the source of the error. In 
>> sample.Polygons, if there is more than one Polygon in the 'shape', this for 
>> loop is run:
>> 
>> for (i in seq(along = pls)) {
>> 	bbi <- .bbox2SPts(bbox(pls[[i]]), proj4string = proj4string) # 
>> Convert bbox to set of 4 points
>> 	# Check which other shapes overlap the bbox of this shape.
>> 	# For a two dimensional object, bb_in is a list of vectors with 
>> length 4:
>> 	# one for each point in the bbox.
>> 	bb_in <- lapply(pls[-i], function(x, pts) pointsInPolygon(pts, x), 
>> pts = bbi)
>> 	# For a two dimensional object, zzz must be a matrix with exactly 
>> four columns.
>> 	zzz <- do.call("rbind", bb_in)
>> 	# For a 2D object, zzz only has four columns, but i can exceed four!
>> 	# But 'i' represents the index of the polygon in the shape, which can 
>> exceed four.
>> 	# zzz[, i] will fail when i exceeds four.
>> 	if (holes[i] || (any(unlist(bb_in) == 1) && !(sum(zzz[, i])%%2) == 
>> 0)) smple[i] <- FALSE
>> }
>> 
>> The comments are my own. Basically, for two-dimensional shapes, the last 
>> line will fail when the number of Polygons in the shape exceeds four, and 
>> at least one of those shapes doesn't have a "hole". (If they do have holes, 
>> then zzz[,i] won't be evaluated.)
>
> Hi Aman,
>
> Thanks for a clear analysis and sample data! Bug fix committed to R-Forge 
> SVN, will be in tomorrow's R-Forge builds (revision 1039). If you need it 
> earlier, please make an anonymous checkout.

In revision 1040, the whole heuristic is removed - it is the user's 
responsibility to set the Polygon objects "hole" slots correctly - if in 
doubt, use checkPolygonsHoles() in the maptools package.

Roger

>
> pts <- spsample(shapefile, n=5000, type="random")
>
> now works.
>
>> 
>> I'm not sure exactly what this last line is supposed to be doing... which 
>> is why I can't make a suggestion for a correction. I suspect that zzz[,i] 
>> should be zzz[i,], but that is just a guess, really.
>> 
>
> The change adding this code on 16 July 2010 was made to try to improve hole 
> logic. Quite often, imported multipolygon objects (Polygons objects in sp) 
> have wrong hole status assigned to the geometries. If checkPolygonsHoles() in 
> the maptools package has not been run to clean them up, this may lead to 
> grief. The code is a guess at finding conditions that would identify a hole 
> even though it is declared as not a hole. It isn't bullet-proof, and did 
> contain a bug, which you identified.
>
> Best wishes,
>
> Roger
>
>> I have tested this on R version 2.12.2 with sp package 0.9-78 on both the 
>> Windows and UNIX (Ubuntu) platforms.
>> 
>> Thanks for any help.
>> 
>> aman
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>
>

-- 
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