[R-sig-Geo] Create SpatialPolygonsDataFrame indicating the polygons in a single source SpatialPolygonsDataFrame that overlap with the new polygons?

[Roger Bivand]

On Fri, 4 Jul 2014, nevil amos wrote:

> I figured out how to achieve this  in raster, however for the real data the
> raster files and array exceed my availible 16gb memory.
>
> I would still be interested in any suggestions on how to achieve directly
> with the SpatialPolygons

Yes, use the rgeos package. In fact your problem isn't very well formed, 
because the outcome is 9 intersections, 3 of which are equal:

gII <- gIntersection(SDF, SDF, byid=c(TRUE, TRUE))
gEquals(gII, gII, byid=TRUE)

Retrieve the IDs from the row.names() of GII, but you'll need to copy back 
the off-diagonal equalities.

You can set up gIntersection() for large n by first using gIntersects(..., 
returnDense=FALSE) to return only an n-list of predicate TRUE indices, or 
query a binary STRtree to generate a similar list. See for example 
suggestions in:

http://spatial.nhh.no/misc/ogrs14/ogrs140612.zip

in ogrs_140612.R (there line/polygon intersections)

This should put you on a feasible track, I hope.

Roger



>
> For the record here is the raster/ array approach I figured out:
> rm(list=ls())
> library(rgdal)
> library(raster)
> library(sp)
> library (plyr)
> p1<-Polygons(list(Polygon(cbind(x=c(0,2,2,0,0),y=c(0,0,2,2,0)))),1)
> p2<-Polygons(list(Polygon(cbind(x=c(1,3,3,1,1),y=c(1,1,3,3,1)))),2)
> p3<-Polygons(list(Polygon(cbind(x=c(0,2,2,0,0),y=c(1,1,2,2,1)))),3)
>
> mypolys<-list(p1,p2,p3)
>
>
> SDF<-SpatialPolygons(mypolys)
> plot(SDF)
> ext<-(extent(SDF))
>
> xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
> n <- 1
> r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)
> for(i in 1:length(SDF)){
>  rr<-rasterize(SDF[i,],r)
>  if(exists("S")) S<-stack(S,rr) else S<-stack(rr)
> }
> A<-as.array(S)
>
> C<-apply(A,3,as.numeric)
> U<-unique(C)
> C<-data.frame(C)
> U<-data.frame(U)
> ID<-as.numeric(rownames(U))
> U<-cbind(ID,U)
> Vals<-matrix(join(C,U)$ID,dim(S)[1],dim(S)[2])
> Vals<-setValues(r,Vals)
> Vals<-ratify(Vals)
> levels(Vals)[[1]]<-U
> names(levels(Vals)[[1]])<-c("ID","p1","p2","p3")
> plot(Vals)
> plot(SDF,add=T)
> print(Vals)
>
>
>
> On Thu, Jul 3, 2014 at 4:17 PM, nevil amos <nevil.amos at gmail.com> wrote:
>
>> I am trying to create a polygon dataset that creates a new polygon for
>> each unique combination of location and overlap of polygons in an input
>> polygon dataset, and is attributed with the input polygons underlying each
>> of the derived polygons, and th output does not contain any overlapping
>> polygons.
>>
>>
>> #trivial example
>> #single layer with three partially overlapping polygons
>> library(sp)
>> p1<-Polygons(list(Polygon(cbind(x=c(0,2,2,0,0),y=c(0,0,2,2,0)))),1)
>> p2<-Polygons(list(Polygon(cbind(x=c(1,3,3,1,1),y=c(1,1,3,3,1)))),2)
>>
>> p3<-Polygons(list(Polygon(cbind(x=c(1.5,3,3,1.5,1.5),y=c(1.5,1.5,3,3,1.5)))),3)
>> mypolys<-list(p1,p2,p3)
>> mydata<-data.frame(Name=c("p1","p2","p3"),row.names=c(1,2,3))
>>
>> SDF<-SpatialPolygonsDataFrame(SpatialPolygons(mypolys),mydata)
>> par(mfrow=c(3,4))
>> plot(SDF)
>>
>>
>> #the desired result is a SpatialPolygonsDataFrame with a single separate
>> polygon for each of the six areas deifined by the interection of the
>> boundaries of the three input polygons
>> text(.5,.5,1)
>> text(1.25,1.25,2)
>> text(1.75,1.75,3)
>> text(2.5,2.5,4)
>> text(1.25,2.5,5)
>> text(2.5,1.25,6)
>>
>> # such that the data frame lists the input polygons (p1,p2,p3) that
>> overlap at the location of each of the output polygons (1:6) and the output
>> has no overlap between the polygons
>> # the results would be the six unique polygons shown in the plot and so
>> the attributes would look like:
>> data.frame(polys=c("p1","p1 p2","p1 p2 p3","p2 p3","p2","p2"))
>> # or alternatively
>>
>> data.frame(cbind(c("p1","p1","p1","p2","p2","p2"),c(NA,"p2","p2","p3",NA,NA),c(NA,NA,"p3",NA,NA,NA)))
>> # or alternatively
>> data.frame(cbind(p1=c(1,1,1,0,0,0),p2=c(0,1,1,1,1,1),p3=c(0,0,1,1,0,0)))
>>
>>
>> #This last plot shows the "apparent polygons" I want in separate colors (
>> because you cannot see the overlapping portions)
>> Outpoly<-union(SDF,SDF)
>> plot(OutPoly,col=1:9)
>>
>> plot(SpatialPolygons(list(Polygons(list(Polygon(cbind(x=c(2,2,3,3,2),y=c(1,1.5,1.5,1,1)))),1))),col=12,add=T)
>> text(.5,.5,1)
>> text(1.25,1.25,2)
>> text(1.75,1.75,3)
>> text(2.5,2.5,4)
>> text(1.25,2.5,5)
>> text(2.5,1.25,6)
>> #The nearest I could get was raster::union however this requires two ( and
>> only 2) SpatialDataFrames as input , what I need is union with a single
>> input union(SDF,SDF) results in 9 polygons not 6, these are all simple
>> rectangles of overlapping pairs (union(union(SDF,SDF)SDF)) results in 27
>> polygons - all again just the simple rectangular overlaps and original
>> polygons
>> union(SDF,SDF)@data
>>
>>
>>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no