[R-sig-Geo] Area of polygons wrong after intersection using gDifference

[Roger Bivand]

On Thu, 12 Jun 2014, Matthew Slocum wrote:

> I originally posted this question to Stack Overflow, and I was told that
> this was the right place to go.

Yes (of course), effectively none of the developers here have time to read 
SO. Users are expected to respect the constraints on developers' time - we 
have day jobs (really).

>
> I performed an analysis on a set of X, Y coordinates producing a Voronoi
> diagram, and I then created a border. Here I intersect the diagram and the
> border, and attempt to get the area of the resulting polygons. The areas
> for the interior "hole" polys are correct, but the edge polygons appear to
> maintain their original exaggerated size. A link to my data is here:
>
> https://drive.google.com/a/ruths.ai/file/d/0B8QG4cbDqH0UaGM2VkkxZHZkZTA/edit?usp=sharing
>
> Code illustrating the problem is here:
>

You posted HTML, don't.

You have provided data files (note that real shapefiles should have a 
*.prj). However, from there on you are just guessing.

Never user maptools to read data, use rgdal (you do not load maptools in 
your code). If you use rgdal, you'll find that the essential comments for 
SpatialPolygons and Polygons objects are generated. If you don't, the 
exterior and interior rings may get mislaid (rgeos imposes comments on 
the fly internally, but there are less forgiving heuristics in maptools).

Never use @ to access slots unless you are a developer. Use subsetting and 
access methods.

Do not guess that plot order is anything other than its name. The hole 
slot is intended to be correct, as is ring direction. The comment in 
Polygons objects is authoritative.

No Polygons object can by definition consist of a single hole - it has to 
have at least one exterior ring.

All of this is documented in Bivand et al. (2013, 2nd ed.) - your email is 
from a commercial organisation in oil/gas consulting so you can buy the 
book - academics have libraries (you do not give an affiliation, but 
should).

> sapply(slot(slot(SPDF[339,], "polygons")[[1]], "Polygons"), slot, 
"hole")
[1] FALSE

339 is not a hole - why should it be?? SPDF are Voronoi polygons, so no 
holes by definition?

> gArea(SPDF[337,])
[1] 3932082

is the area of the polygon - not wrong!

Why gDifference?? Of the possibilities, gIntersection is the obvious one 
to use, if you want to clip SPDF to spP. You even use the word intersect, 
so how does your code represent your actual workflow? Your final sentence 
is equally unreflected, as edges have by definition no area. Do you intend 
to do:

D <- gIntersection(SPDF, SpP, byid=c(TRUE, FALSE))
plot(D, add=TRUE, col=sample(rainbow(339)))

summary(gArea(SPDF, byid=TRUE) - gArea(D, byid=TRUE))

Hope this clarifies,

Roger

>
> # Read in shapefiles.
> # Files are located at:
>
> require(sp)
> require(rgeos)
> SPDF <- readShapeSpatial("SPDF.shp")
> SpP <- readShapeSpatial("SpP.shp")
>
> # Examine plots
> plot(SPDF)
> SPDF at polygons[[337]]@area # Too large; want it cut off
> SPDF at polygons[[339]]@area # Hole poly; area correct
> gArea(SPDF[339,]) # Ditto
> gArea(SPDF[337,]) # Still wrong
>
>
> # Merge polys using gDifference
> D <- gDifference(SpP, SPDF, byid = TRUE)
> plot(D)
>
> # Seems to work, but areas now have a couple of problems.
> # I pick apart D using the plotOrder slot to separate
> # polys that are holes versus those that are not, allowing
> # me to get the correct area for "hole" polys, but the
> # edge polygons are still not correct, maintaining their
> # area estimates from the original SPDF data frame.
>
> areas <- vector()
> for (i in 337:339){ # 337 = exterior poly, 338 and 339 are holes
>  po <- D at polygons[[i]]@plotOrder
>
>  if (max(po) == 2) {
>    areas[i] <- D at polygons[[i]]@Polygons[[2]]@area
>  } else {
>    areas[i] <- D at polygons[[i]]@area
>
>  }
> }
>
> areas
>
> How does one get the right areas for the edges that should be cut
> off by the intersection?
>
> 	[[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