Skip to content

Calculate each polygon percentage inside a circles

6 messages · Juan Tomas Sayago Gomez, Rolf Turner, ASANTOS

ASANTOS
# 
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
Juan Tomas Sayago Gomez
# 
you should check the intersect function. I used it like in this function in
the link:
http://gis.stackexchange.com/questions/140504/extracting-intersection-areas-in-r

Juan

On Tue, May 24, 2016 at 9:17 PM ASANTOS via R-sig-Geo <
r-sig-geo at r-project.org> wrote:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      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

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      ASANTOS
    

    

      
      #
    

  


  

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

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

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

        

            

      
      
      
    

        
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

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      ASANTOS
    

    

      
      #
    

  


  

      
Dear Rolf Turner,

         Amazing!!!!! I will never put R in doubt.

        My problem was solved with success, thanks very much again!!!

Best wishes,

Alexandre