Skip to content

overlapping polygons in shapefile

3 messages · Roger Bivand, Mihai Valcu

Original
Roger Bivand
#
On Wed, 12 Mar 2014, Alice C. Hughes wrote:
Please do not take threads off-list - I am not offering private advice, 
and what is said on the list is intended to be public so that others can 
contribute.
The whole point of using R is that users are encouraged to create and modify 
tools from parts (aka functions or methods). Most research problems are coerced 
into fixed subsets of problems that can be tackled with "other people's" tools. 
If researchers are to take responsibility for actually tackling research 
questions with the best alignment of data collection, data representation, 
workflow in data handling, analysis, it is their respopnsibility to ensure that 
the tools are as well aligned with the problem as possible. Otherwise, your 
ability to infer about the actual processes will have been amputated by the 
limitations imposed by using tools not matching your problem.

Change of support is particularly important here, so your ourput map of species 
incidence (don't call it a heatmap, it isn't measuring heat) depends crucially 
on steps in the workflow (for example, what is the error involved in 
delineating the range polygons?).

I was mislead by your phrasing to think that you also needed the areas of each 
species x species overlap (counts are a bit trivial), and you haven't said what 
you mean by area (is it a measure, a polygon, or a raster cell?). What does: 
"number o[f] species sugested to be in an area" mean? Do you mean a raster 
cell? Or the polygonal intersection output for each unique count? Or what?

gOverlaps(..., byid=TRUE) should give counts (it returns a matrix, and the 
rows/columns should be the counts if the polygons are species). However, if you 
have many polygons, do use the STR tree approach only to search for overlaps 
among polygons that actually have intersecting bounding boxes. You can count 
the number of overlaps in this way. If someone can provide random overlapping 
polygons in a SpatialPolygons object, it would be easier to show. All R lists 
do ask for reproducible cases with built-in or simulated data.

Roger

  
    

      
      
    

  


  


    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
Initial example below, please scroll down:

        
On Wed, 12 Mar 2014, Roger Bivand wrote:

        

            

      
      
      
    

        
A first cut at an example:

library(sp)
library(rgeos)
box <- readWKT("POLYGON((0 0, 0 1000, 1000 1000, 1000 0, 0 0))")
plot(box)
set.seed(1)
pts <- spsample(box, n=2000, type="random")
pols <- gBuffer(pts, byid=TRUE, width=50)
plot(pols, add=TRUE)
STR <- gBinarySTRtreeQuery(pols, pols)
length(STR)
summary(sapply(STR, length))
res <- vector(mode="list", length=length(STR))
for(i in seq(along=STR)) res[[i]] <- gOverlaps(pols[i], pols[STR[[i]]],
   byid=TRUE)
summary(sapply(res, sum)-1) # to remove self-counting, i overlaps i

So then we have # overlap counts per input polygon - but what are the 
output reporting units? Should we actually be counting the number of 
polygons that a fine grid of points over box belong to? If we focus on 
this, then maybe:

GT <- GridTopology(c(0.5, 0.5), c(1, 1), c(1000, 1000))
SG <- SpatialGrid(GT)
o <- over(SG, pols, returnList=TRUE)
ct <- sapply(o, length)
summary(ct)
SGDF <- SpatialGridDataFrame(SG, data=data.frame(ct=ct))
spplot(SGDF, "ct", col.regions=bpy.colors(20))

is closer? I haven't checked why the overlaps and the grid over counts 
differ - homework for someone? How does the output resolution affect the 
detected number of hits?

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Mihai Valcu
    

    

      
      #
    

  


  

      
Perhaps, rangeMapper can be of use. It is using spsample(, type = 
"regular") and over()  so it is similar with Roger's example below.

# Species richness of American wrens -----
require(rangeMapper); require(RSQLite); require(rgdal)
dbcon = rangeMap.start(file = "test.sqlite", overwrite = TRUE, dir = 
tempdir() )
f = system.file(package = "rangeMapper", "extdata", "wrens", 
"vector_combined")
r = readOGR(f, "wrens")
global.bbox.save(con = dbcon, bbox = f)
gridSize.save(dbcon, gridSize = 2) ; canvas.save(dbcon)
processRanges(spdf = r, con =  dbcon, ID = "sci_name" )
rangeMap.save(dbcon) # FUN is missing so species richness is computed by 
default. see ?rangeMap.save for other args
m = rangeMap.fetch(dbcon)
plot(m)
####

Mihai

        
On 3/12/2014 11:00 AM, Roger Bivand wrote: