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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
