1

I have two polygon shapefiles and I want to see to what extent the observed amount of overlap is due to chance. I'd like to randomly move the polygons in one shapefile to new locations, measure the total area of intersection, repeat this a thousand times or so to get a distribution of areas of intersection, and then measure how likely the observed area is compared to the distribution of random areas.

I'm stuck on moving the polygons: how can I shift polygons to a random new location in a bounding box? I'm using R, here's an attempt at a reproducible example:

require(sp)
require(splancs)
require(rgeos)
require(maptools)

# simple example polygons
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)

Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
plot(SpP, col = 1:3, pbg="white")

# make bounding box to contain random relocations of polygons
bbx <- bboxx(bbox(SpP))
bbx <- rbind(bbx, bbx[1,])
bbx <- SpatialPolygons(list(Polygons(list(Polygon(bbx)), "bbx")))

# make random points in bounding box to use as centroids to sample from
rnd <- spsample(bbx, 100, type = "random")

# get centroids of polygons
cnt <- gCentroid(SpP, byid = TRUE)

Here's where I'm stuck...

# replace centroids of polygons with random points in bounding box
# in the hope that this will shift the polygons
out <- SpP # to store output
for(i in 1:length(SpP)) {
    # create a random offset to shift the polygon by
    offset <- unname(cnt[i,]@coords) - unname(rnd@coords[sample(1:length(rnd),1),]) 
    # shift the polygon by this random offset
    # this doesn't seem to be doing what I want...
    out[i,] <-  elide(SpP[i,], shift = offset)
}
plot(out) # not the desired result...

# check to see that something is happening...
sapply(1:length(SpP), function(k) SpP@polygons[[k]]@labpt)

The desired output is a shapefile with the polygons in new random locations in the bounding box. My general approach is based on this. What am I doing wrong and how can I fix it?

  • You have to shift not centroids of the polygons, but every vertex instead.BTW I'm not a great statistician, but to my mind your general approach is flawed. But I can not suggest anything due to lack of information. – SS_Rebelious Jun 22 '13 at 20:20
  • Thanks for the detail about vertices, I'll give that a shot (feel free to supply code in more detailed answer!). Your comment has prompted me to get some feeback on the statistical method, so I've put more detail about my specific use-case in a question at CV – Ben Jun 23 '13 at 3:31
0

Seems that I have to move every coordinate stored in the polygon, vertices, centroids and bounding box coords. Corrections and suggestions for improvement welcome!

This is what I've found to work, where rocks comes from a polygon shapefile and rnd are the random points similar to in my question:

rocks_rnd <- vector("list", length = length(rocks))
# randomly relocate every rock in the shapefile
for(i in 1:length(rocks)) { 
# this is where we move each rock, one by one
# get a rock from our observed rocks to shift
ri <- rocks[i,]
# we have to shift all of the coords relating to this rock...
# get coords of vertices (outline of rock)
cds <- slot(slot(slot(ri, "polygons")[[1]], "Polygons")[[1]], "coords")
# get coords of labpt, labpt, bbox also (centrepoints and bounding box)
l1 <- slot(slot(slot(ri, "polygons")[[1]], "Polygons")[[1]], "labpt")
l2 <- slot(slot(ri,  "polygons")[[1]], "labpt")
b1 <- slot(ri,  "bbox")
# get a random point in the study area to shift to
rn <- unname(rnd@coords[sample(1:length(rnd),1),])
# shift all the vertices coords
ofst <- cds[1,] - rn
newcds <- t(apply(cds, 1, function(x) x - ofst))
# shift lapt1 (centrepoint)
ofst <- l1 - rn
newl1 <- l1 - ofst
# shift lapt2 (centrepoint)
ofst <- l2 - rn
newl2 <- l2 - ofst
# shift bbox (bounding box)
ofst <- t(apply(b1, 2,  function(x) x-rn))
newb1 <- b1 - ofst
# put these shifted points back into the polygon object
slot(slot(slot(ri, "polygons")[[1]], "Polygons")[[1]], "coords") <- newcds
slot(slot(slot(ri, "polygons")[[1]], "Polygons")[[1]], "labpt") <- newl1
slot(slot(ri,  "polygons")[[1]], "labpt") <- newl2
slot(ri, "bbox") <- newb1
# assign this rock to the list of rocks in the shapefile
rocks_rnd[[i]] <- ri 
}
# make SpatialPolygonsDataFrame out of the list
# this is the object with the shifted polygons
r_r <- do.call( rbind, rocks_list)
0

It seems now you can use the shift() function in the {raster} package to move SpatialPolygons...

http://www.inside-r.org/packages/cran/raster/docs/shift

Actually it didn't work for me. So I used the elide() function in the {maptools} package...

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.