[R-sig-Geo] PBSmapping calcArea

[Roger Bivand]

On Thu, 25 Sep 2014, Kenny Bell wrote:

> Hi all,
>
> I am trying to use calcArea to get the areas of the following polygons at
> the following link:
>
> https://www.dropbox.com/s/v1fjsdkcyxw61qt/rsiggeoQ.RData?dl=0

Thanks for providing a reproducible example.

>
> I use:
>
> library(PBSmapping)
> areaY <- calcArea(smallYPolySet, rollup = 2)
>
> which gives the output:
>
>> areaY
>  PID SID        area
> 1   1   1  38.7704526
> 2   2   1 116.6810084
> 3   2   2   0.9844843
>
> The second polygon has a hole which can be seen with:
>
> plotLines(smallYPolySet)
>
> However, the rollup = 2 option should see the hole and subtract the area
> and return the difference, as I understand. The above shows they are
> specified separately.
>
> Can someone help me with any/all of:
>
> 1) working out whether the hole is identified in the smallYPolySet object
> 2) flagging an element of a polyset as a hole
> 2) automatically determining if a polygon is a hole
>

Do everything step-by-step under your own control, using the standard sp 
class toolset:

load("rsiggeoQ.RData")
tail(smallYPolySet)
library(PBSmapping)
library(maptools)
sp <- PolySet2SpatialPolygons(smallYPolySet)
# function does not try to detect holes, as they may be wrong
sapply(slot(sp, "polygons"), function(p) sapply(slot(p, "Polygons"), slot,
  "hole")) # no holes
sapply(slot(sp, "polygons"), comment)
sp <- createSPComment(sp)
sapply(slot(sp, "polygons"), comment) # still no holes,
# the comment encodes to which exterior ring "0" in each Polygons object
# each hole belongs
library(rgeos)
slot(sp, "polygons") <- lapply(slot(sp, "polygons"), checkPolygonsHoles)
sapply(slot(sp, "polygons"), function(p) sapply(slot(p, "Polygons"), slot,
  "hole")) # to show change
sapply(slot(sp, "polygons"), comment) # to show change 
# hole now correctly assigned
library(rgdal)
sp_utm <- spTransform(sp, CRS("+proj=utm +zone=24 +units=km"))
gArea(sp_utm, byid=TRUE)

smallYPolySet1 <- smallYPolySet[c(1:59, 65:60),]
calcArea(smallYPolySet1, rollup=2)

I think that the specific problem is that your PolySet PID 2 SID 2 does 
not have POS in descending order (see tail() above) - a change made in the 
S3 class to try to say whether a ring is a hole or not; the sp Polygon 
class has a "hole" slot which is set TRUE when we know that it is a hole, 
and which can be checked and corrected using maptools::checkPolygonsHoles 
on Polygons objects (those with both exterior and possibly interior 
rings). Using a convenience function like calcArea isn't convenient when 
the object isn't what you think it is.

Hope this clarifies,

Roger

> Any help would be very much appreciated!
>
> Thanks,
> Kenny
>
> 	[[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