sample.Polygons fails for shapes with more than four Polygons

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.

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