[R-sig-Geo] [FORGED] Calculate each polygon percentage inside a circles

[Rolf Turner]

On 26/05/16 00:23, ASANTOS wrote:

> Dear Rolf Turner,
>
>            It's much better a clean code with a minimum packages, thank
> you very much for your answer. But "pct" object give me a total polygon
> percentage around each point and I need too an identification (in
> columns) of individual contribution of each polygon. In my simulation, I
> find 50.38001% for the point 1, but this is a total percentage of
> polygons and I'd like to know a percentage contribution for each polygon
> (e.g. ID1 = 0.00000 + ID1 = 10.00000 + ID3 = 0.00000 + ID4 = 40.38001 =
> 50.38001 total), this is possible?


Of course it's possible!  This is R. :-)

You just need to look at the component polygons as separate window 
objects and do the same thing to them that I did to the overall "W".
There is a small gotcha that has to be coped with when the intersection
of the disc with a polygon is empty (as often occurs).  (See below.)

There is a number of different ways to write the code for this.

One way is as follows:

==========================================================================
library(spatstat)

sr1 <- owin(poly=cbind(c(180114, 180553, 181127, 181477, 181294,
                      181007, 180409, 180162, 180114),
                    c(332349, 332057, 332342, 333250, 333558,
                      333676, 332618, 332413, 332349)))

sr2 <- owin(poly=cbind(rev(c(180042, 180545, 180553, 180314, 179955,
                          179142, 179437, 179524, 179979, 180042)),
                    rev(c(332373, 332026, 331426, 330889, 330683,
                          331133, 331623, 332152, 332357, 332373))))

sr3 <- owin(poly=cbind(rev(c(179110, 179907, 180433, 180712, 180752,
                          180329, 179875, 179668, 179572, 179269,
                          178879, 178600, 178544, 179046, 179110)),
                    rev(c(331086, 330620, 330494, 330265, 330075,
                          330233, 330336, 330004, 329783, 329665,
                          329720, 329933, 330478, 331062, 331086))))

sr4 <- owin(poly=cbind(c(180304, 180403,179632,179420,180304),
                    c(332791, 333204, 333635, 333058, 332791)))

wins <- solist(sr1,sr2,sr3,sr4)

W <- union.owin(wins)

set.seed(42)
X <- rpoispp(28/area.owin(W),win=W)
N <- npoints(X)
plot(X,cols="blue")
AD <- area.owin(disc(radius=600))

pct <- matrix(nrow=N,ncol=4)
rownames(pct) <- paste("point",1:N,sep=".")
colnames(pct) <- paste("sr",1:4,sep=".")
for(i in 1:npoints(X)) {
     Di <- disc(radius=600,centre=c(X$x[i],X$y[i]))
     for(j in 1:4) {
         Aij <- intersect.owin(Di,wins[[j]],fatal=FALSE)
         pct[i,j] <- if(is.null(Aij)) 0 else 100*area.owin(Aij)/AD
     }
}
==========================================================================

HTH

cheers,

Rolf Turner

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276