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

[Roger Bivand]

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