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

[Patrick Giraudoux]

Except this inexpected 'bug' discovery point the fact that I draw the 
hole inadvertently ! I wanted to draw standards polygons but did the 
first counterclockwise (so a hole...)... Chaos added to chaos = bug fix. 
Let us call that the French wine method... hips !

Anyway, thanks a lot Roger,

Cheers,

Patrick






Roger Bivand a écrit :
> 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
>>
>>
>>
>>
>>
>>
>>
>>
>