Testing for non-random overlap of polygons

Question

I have two polygon shapefiles and I want to see to what extent the observed amount of overlap is due to chance.

I'm thinking of some kind of permutation test, but not sure of the best way to proceed.

One idea is randomly relocating the polygons in one shapefile to new locations, measuring the total area of intersection, then repeating this a thousand times or so to get a distribution of areas of intersection. Then I can measure how likely the observed area is compared to the distribution of random areas. I have a recent GIS question relates to this approach, it was a comment on that question motivated this question.

Another idea is to combine the two polgyon shapefiles into one shapefile, then randomly reassign the labels of the two types of polygon, measure the total area of intersection under the new labelling scheme, then repeat, etc. as above.

The specific use-case is something like this: one polygon shapefile contains polygons of human activity areas (a few larger polygons). The other polygon shapefile contains polygons of rock outlines (many smaller polygons). The question is whether the observed arrangement of rocks is determined by the activity areas or not: did people put rocks on their activity areas or are we just looking at a natural distribution of rocks?

Are either of the two ideas to test this on the right track, or should I be thinking of something completely different? Any insights would be most welcome!

Answer 1

In this case you should use spatial pattern analysis to identify relations between human activity and distribution of rocks. Rocks have to be represented as a point shp-file. Adding another level of abstract to distribution of rocks (making them polygons instead of points) is highly questionable approach that doesn't seem to provide any help for analysis. You will need a spatstat package (with the great tutorial) for R and maybe QGIS.

1. [Are you hunting for ghosts?] Use Monte-Carlo test to be sure that your rocks are not randomly distributed across the study area (regardless of human activity). To perform it see the corresponding section of tutorial or this post. In this case o-win object (it is the boundaries of the study area, see tutorial) will be the whole area that you study.
2. Use Monte-Carlo test to determine whether rocks are randomly distributed across human activity sites (in this case o-win object will be polygon shp-file of human activity).
3. Get the empirical graph of the dependency of intensity of rocks distribution on the distance to the human activity sites. If humans determine rocks location then intensity of spatial distribution of rocks will depend on the distance to the areas of human activities. Use rhohat function for it. If distribution of rocks indeed depends on human activity sites then the graph in your case should look similar to this: Here o-win object again will be the whole study area. Use a distfun function for a polygons of human activity as a covariate in rhohat.

Here is some example code for the third step from one of my projects:

library(maptools)
library(rgdal)
library(spatstat)
spatstat.options(gpclib=TRUE)
gpclibPermit()

# load point shp-file for analysis
S <- readShapePoints("/dumps.shp", 
                      proj4string= CRS("your_+proj_sting_here"))
SP <- as(S, "SpatialPoints")
P <- as(SP, "ppp")

# load boundary layer and make it o-win object
Z <- readShapePoly("/boundaries.shp", 
            proj4string= CRS("your_+proj_sting_here"))
Z1 <- as(Z, "SpatialPolygons")
W <- as.owin(Z1)
P <- P[W]

# in my case covariate polygons were quite small so I loaded them as lines
# to avoid spatstat issue with polygons. You should represent your
# human activity polygons as points as will be described below
c <- readShapeLines("/ccovariate.shp", 
        proj4string= CRS("your_+proj_sting_here"))

cr <- as.psp(c)
cr <- cr[W]

# create a distance function 
сrdist <- distfun(cr)

# create and plot your graph
plot(rhohat(P, сrdist, covname="quarry"), 
         xlab= "Расстояние до карьера, м",
         legendpos = "topright", ##see help(plot.fv)
         main = NULL)

That's pretty much it. Now some important details.

• If you dont't have a point shp-files for rocks (in this case you should skip first two steps cause it will be pointless) you can recreate it using your polygon layer. Use QGIS for this. Go Vector -> Research Tools -> Random points. Here choose your polygon layer of rocks as an input layer and set up one of the option for the Individual Polygons (density or the number of points per polygon).
• If you have issues with human activity polygons as a covariate for the rhohat function (spatstat sometimes not working well with the polygons). You can replace that polygons with the point layer just as it was suggested in the previous paragraph, but using Regular points instead of Random points.

P.S. You may take some other approach (get inspiration from the tutorial for spatatat) but it is essential to use point pattern analysis in this case, not some polygons of locations of rocks.