Skip to content
Prev 24434 / 29559 Next

[FORGED] Calculate each polygon percentage inside a circles

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?

Thanks again,

Best wishes,

Alexandre


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

bdry <- list(cbind(c(180114, 180553, 181127, 181477, 181294,
                      181007, 180409, 180162, 180114),
                    c(332349, 332057, 332342, 333250, 333558,
                      333676, 332618, 332413, 332349)),

              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))),

              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))),

              cbind(c(180304, 180403,179632,179420,180304),
                    c(332791, 333204, 333635, 333058, 332791)))

W <- owin(poly=bdry)

set.seed(42)
X  <- rpoispp(28/area.owin(W),win=W)
plot(X,cols="blue")

for(i in 1:npoints(X)) plot(disc(radius=600,centre=c(X$x[i],X$y[i])),
      add=TRUE,col=rgb(1,0.5,0,alpha=0.4),border=NA)

# To be fair, calculate the area of the discs using area.owin()
# rather than as pi*600^2.

AD <- area.owin(disc(radius=600))

pct <- numeric(npoints(X))
for(i in 1:npoints(X)) {
     Di <- disc(radius=600,centre=c(X$x[i],X$y[i]))
     pct[i] <- 100*area.owin(intersect.owin(Di,W))/AD
}
#==================================================================================