sp, area and SpatialPolygons
On Sat, 17 Apr 2010, Patrick Giraudoux wrote:
I did not find a function in sp to extract area from SpatialPolygons objects. Thus I was thinking about writing some function to do the job similar to the functions written in order to extract lengths from SpatialLines objects (see ?SpatialLines). I have incidentally discovered a strange thing. Let us define a list of 3 polygons: polylist <- list(structure(c(180016, 180034, 180452, 180588, 180615, 180533, 180225, 180016, 332182, 331756, 331774, 332074, 332418, 332518, 332319, 332182), .Dim = c(8L, 2L)), structure(c(179907, 180325, 180397, 180152, 179781, 179672, 179735, 179907, 331611, 331611, 331266, 330931, 330967, 331266, 331466, 331611), .Dim = c(8L, 2L)), structure(c(179499, 179971, 180343, 180161, 179753, 179418, 179499, 330577, 330768, 330468, 330096, 330078, 330369, 330577 ), .Dim = c(7L, 2L))) And make a simple SpatialPolygons object with them: SP<-SpatialPolygons(list(Polygons(list(Polygon(polylist[[1]])),ID="P1"),Polygons(list(Polygon(polylist[[2]])),ID="P2"),Polygons(list(Polygon(polylist[[3]])),ID="P3")),pO=as.integer(c(1,2,3))) It looks like if the area of Polygons #1 was not computed:
Patrick, There is a bug, not not a simple one. The ring direction of your first Polygon object signals that it is a hole, but in Polygons(), singleton holes are made into islands, and the largest hole in an all-hole list of Polygon objects is also made into an island. The bug was that the hole vector used for assigning the area had not been updated, so the area slot in the Polygons object was assigned 0 (because it was still treated as a hole). Had you said hole=FALSE for P1, it would have been OK. The updating is now in place. Thanks, Roger
SP at polygons[[1]]@area [1] 0 Whilst the area of the Polygons are: SP at polygons[[1]]@Polygons[[1]]@area [1] 322111 However, no trouble with the next Polygons:
SP at polygons[[2]]@area
[1] 384740
SP at polygons[[2]]@Polygons[[1]]@area
[1] 384740
SP at polygons[[3]]@area
[1] 423937
SP at polygons[[3]]@Polygons[[1]]@area
[1] 423937 Did I do something wrong, miss something or is there a bug in the SpatialPolygons function ? Patrick
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