[R-sig-Geo] Unexpected behaviour of "area" slot in sp Polygons

[Roger Bivand]

On Mon, 26 Jan 2015, Gregory Duveiller wrote:

> Dear All,
>
> I noticed something unusual when trying to get the areas of some polygons I
> extracted from a raster.
>
> In the script below, I generate a random raster with 4 classes, which I
> then turn to a "SpatialPolygonsDataFrame". I then try to calculate areas
> per class in two different ways.
>
> ########
>
> require(sp)
> require(raster)
> require(rgeos)
>
> nrows=100; ncols=100
>
> rdum<-raster(nrows=nrows,ncols=ncols,xmn=0,ymn=0,ymx=100,xmx=100)

set.seed() # please!

> rdum[]<-round(runif(nrows*ncols,1,4))
>
> vdum <- rasterToPolygons(rdum,dissolve=T)
>
> # calculate areas
> df <- data.frame(cls=1:4)
> for(i in 1:4){
>  ddd <- sapply(vdum at polygons[[i]]@Polygons, function(x) x at area)
>  df$area1[i] <- sum(ddd)
>  df$area2[i] <- sum(as.vector(rdum==i))
> }
> #########
>
> When I try to calculate the area per class from the attributes of the
> Polygons, I do not always get the same number as when I simply count the
> number of pixels of each class in the raster. In one iteration, for
> example, I get the following number in which the areas in the Polygons is
> higher by 52 units:
>
>> df
>  cls area1 area2
> 1   1  1670  1670
> 2   2  3384  3352
> 3   3  3355  3335
> 4   4  1643  1643
>

Without a seed set, this can't be reproduced, but with:

rdum<-raster(nrows=nrows,ncols=ncols,xmn=0,ymn=0,ymx=100,xmx=100)
set.seed(1)
rdum[]<-round(runif(nrows*ncols,1,4))
vdum <- rasterToPolygons(rdum,dissolve=T)
df <- data.frame(cls=1:4)
for(i in 1:4){
   ddd <- sapply(vdum at polygons[[i]]@Polygons, function(x) x at area)
   df$area1[i] <- sum(ddd)
   df$area2[i] <- sum(as.vector(rdum==i))
}
df$area3 <- gArea(vdum, byid=TRUE)

we get

> df
   cls area1 area2 area3
1   1  1692  1692  1692
2   2  3389  3355  3355
3   3  3249  3219  3219
4   4  1734  1734  1734
> sapply(df, sum)
   cls area1 area2 area3
    10 10064 10000 10000

The purpose of the area slot is to order the polygons for plotting from 
largest to smallest, which was needed for colour fills before the graphics 
polygon plotting function knew how to handle holes (so the area is gross 
of holes), and is still needed for hatched filling. It isn't necessarily 
to give the area of a Polygons object, see ?"Polygons-class".

There are holes in vdum:

> table(unlist(sapply(slot(vdum, "polygons"), function(x) sapply(slot(x, 
"Polygons"), slot, "hole"))))

FALSE  TRUE
  4734    28

so this is a probable explanation.

Roger

> Is there a reason one would not expect these to be always equal?
>
> Thanks in advance...
>
> Gregory Duveiller
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no