[R-sig-Geo] sp, area and SpatialPolygons

[Roger Bivand]

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