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

[Rolf Turner]

Dear Alexandre,

Your tasks can all be done very simply in spatstat.  The code follows.
Note that I had to reverse the point order for sr2 and sr3 because
spatstat requires that the vertices of (exterior) boundaries of 
polygonal windows be given in anticlockwise order.

I commented out the plotting of the discs centred at each point of the 
simulated pattern since I found the plot to be unenlightening and messy 
looking.

The resulting percentages that you are after are in "pct".

If you want more precision for the disc areas, set npoly equal to a 
large number than the default 128 (e.g.1024 or 2048) in the calls to disc().

cheers,

Rolf Turner

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

#==================================================================================
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
}
#==================================================================================

On 25/05/16 13:17, ASANTOS via R-sig-Geo wrote:
> Dear members,
>
>         I will try to calculate each polygon percentage inside a circles
> given an arbitrary radius in a shapefile object with the code below and
> my output needs to be (Two first columns with center os circle
> coordinates and values of each polygon percentage):
>
> "pts$x"   "pts$y" "ID1" "ID2" "ID3" "ID4"
> 180486.1  330358.8  16.3   0.2  NA   17.3
> 179884.4  331420.8  88.3   NA   96.3 NA
> 179799.6  333062.5  25.3   22.3 0.5  NA
> ...
>
>       For this I try to do:
>
> #Packages
> require(shapefiles)
> require(rgdal)
> require(maptools)
> require(spatstat)
> require(berryFunctions)
>
>
> #Create 4 polygons (I create polygons because is more simple that given
> a shapefile in a post)
>
> sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477,
> 181294, 181007, 180409,
>    180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
>    332618, 332413, 332349)))),'1')
> sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314,
> 179955, 179142, 179437,
>    179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
>    331133, 331623, 332152, 332357, 332373)))),'2')
> sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712,
> 180752, 180329, 179875,
>    179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
>    c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
>    329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
> sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
>    c(332791, 333204, 333635, 333058, 332791)))),'4')
>
> #Spatial Polygons
> sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
> srdf=SpatialPolygonsDataFrame(sr,
> data.frame(row.names=c('1','2','3','4'), PIDS=1:4))
>
> #Create shapefile
> writeOGR(srdf, getwd(), 'POLY', 'ESRI Shapefile')
>
> #Read shapefile
> contorno_line_X <- readShapeSpatial("POLY.shp")
> plot(contorno_line_X)
>
> #Create ~28 random points in my window
> contorno_line_spat <- as(contorno_line_X,"owin")
> pts <- rpoispp(lambda=0.000005, win=contorno_line_spat)
> points(pts, col="blue")
>
> #Set the radius for the plots
> radius <- 600  # radius in meters
>
> #Draw the circles using berryFunctions package
> for(i in 1:length(pts$x)) circle(pts$x[i],pts$y[i],r=radius,
> col=rgb(1,0.5,0,alpha=0.4), border=NA)
> #
>
> Now I'd like to calculate de polygon percentage around each point in a
> 600 meters radius. This is possible? There are any function for this?
>
> Thanks,
>
> Alexandre
>