Union/Overlay single spatialPolygonsDataFrame containing overlapping polygonsto create polygons with all unique combinations of attributes

Nevil,

raster::union does not _require_ two SpatialPolygons*, it also works
with one. See ?raster::union,  this is the (valid) case where 'y' is
missing.

library(raster)
u <- union(p4)
data.frame(u)

  ID.1 ID.2 ID.3 ID.4 count
1    1    0    0    0     1
2    0    1    1    0     2
3    0    0    0    1     1

u

class       : SpatialPolygonsDataFrame
features    : 3
extent      : 1, 5, 1, 5  (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables   : 5
names       : ID.1, ID.2, ID.3, ID.4, count
min values  :    0,    0,    0,    0,     1
max values  :    1,    1,    1,    1,     2

plot(u, col=rainbow(3))

Robert

On Tue, Aug 11, 2015 at 8:14 PM, nevil amos <nevil.amos at gmail.com> wrote:

Thanks for the suggestions.

unfortunately raster::union does not get the desired result,

First it requires two, and only two spatialPolygons objects as inputs. -

I

think the function I need must either have a single input containing all

the

polygons, or one input for each polygon to work in more complex cases

with

many polygons with more than tow overlapping polygons

Second in my example it does indeed produce four polygons, however these
comprise the two original square polygons and two polygons from the
overlapping areas:

#what union::(SPDF2,SPDF2) acheives

par(mfrow=c(2,4))
myUnion<-raster::union(mySPDF2,mySPDF2)
myUnion$Area<-gArea(myUnion, byid=TRUE)
for (i in 1:4){
  plot(myUnion[i,],col=cols[i],axes=T,xlim=c(0,5),ylim=c(0,5))
}
print(myUnion at data)
print("Areas myUnion:")
print(gArea(myUnion, byid=TRUE))


#what is required:

p4<-SpatialPolygons(list(Polygons(list(Polygon(cbind(c(1,4,4,3,3,1,1),c(1,1,3,3,4,4,1)),hole

= F)), "1_ "),

Polygons(list(Polygon(cbind(c(3,4,4,3,3),c(3,3,4,4,3)),hole = F)),

"1_2"),

Polygons(list(Polygon(cbind(c(3,4,4,3,3),c(3,3,4,4,3)),hole = F)),

"2_1"),

Polygons(list(Polygon(cbind(c(4,4,5,5,3,3,4),c(4,3,3,5,5,4,4)),hole =

F)),

"2_")))
for (i in 1:4){
  plot(p4[i,],col=cols[i],axes=T,xlim=c(0,5),ylim=c(0,5))
}
print("Areas p4:")
print(gArea(p4, byid=TRUE))


On Wed, Aug 12, 2015 at 11:02 AM, Robert J. Hijmans <r.hijmans at gmail.com

wrote:

I would think raster::union(polygons) should get you there.
Robert

On Tue, Aug 11, 2015 at 7:58 AM, Michael Sumner <mdsumner at gmail.com>
wrote:

On Mon, 10 Aug 2015 at 17:52 nevil amos <nevil.amos at gmail.com> wrote:

I am attempting to combine overlapping polygons in a
single spatialPolygonsDataFrame so that the output is a set of

polygons

such that where there is no overlap one polygon, of the
non-intersecting
area will be formed, and retain the attributes of the original

polygon.

 where polygons intersect  new polygons, of the intersecting area

only

will
be formed, one with the attributes of each of the original polygons
intersecting in that area.

The example includes a single overlap, the real data may contain many
tens
or hundreds of intersecting polygons, from among tens of thousands of
polygons.

the step that I cannot work out is "myUnion"  it is aiming to achieve
the
equivalent of a self- union in ESRI geoprocessing.

I realize that I can achieve a similar result with raster but would
prefer
to stay with polygons for this process.

Is raster::cover closer to what you want?

library(raster)

plot(mySPDF2)
 plot(cover(mySPDF2[1,], mySPDF2[,2]), add = TRUE, col = rgb(c(0, 1),
c(0,
0), c(1, 0), c(0.5)))
 plot(mySPDF2)
 plot(cover(mySPDF2[2,], mySPDF2[1,]), add = TRUE, col = rgb(c(0, 0),
c(0,
1), c(1, 0), c(0.5)))

You'll still need to take the non-overlapping piece from the

complement

overlay. I don't know if there's an easier way.
See ?"raster-package" for a complete overview.

HTH, Mike.


many thanks in advance

Nevil Amos

Trivial example:
library(rgeos)
library(maptools)
library (tidyr)
library(plyr)
rm(list = ls())
par(mfrow=c(2,2))
cols=c("red", "black", "blue", "yellow")
cols<-adjustcolor(cols,alpha.f=.5)

p2<-

SpatialPolygons(list(Polygons(list(Polygon(cbind(c(1,4,4,1,1),c(1,1,4,4,1)),hole

= F)), 1 ),

 Polygons(list(Polygon(cbind(c(3,5,5,3,3),c(3,3,5,5,3)),hole = F)),
2)))

ID<-c(1,2)
YEAR<-c(1999,2002)
Type<-c(1,2)
mydf2<-data.frame(ID,YEAR,Type)
mySPDF2<-SpatialPolygonsDataFrame(p2,mydf2)
mySPDF2$Area<-gArea(mySPDF2, byid=TRUE)
print("mySPDF2 at data")
print(mySPDF2 at data)
plot(mySPDF2,axes=T,col=cols)
mtext("mySPDF2",cex=2)
mtext(paste("colors =","red", "yellow", "blue", "green,"), outer =
TRUE,
cex = 1.5)
#this is the step that does not give the desired result. It has the
same
two polygons ( areas 9 and 4) as the input, no intersection and
clipping of
the polygons occurs   I am trying to acheive the equivalent of a

"self

union in ESRI ARCGIS.  The result I would like be four polygons: two

"L

shaped plygons formed from the non- intersecting parts of the input
polygons (Areas 8 and 3 respecively), and two square polygons (
coincident
both of Area = )  I manually create my desired output in the next
section
to illustrate what I would like to achieve
myUnion<-unionSpatialPolygons(mySPDF2,IDs=mySPDF2$ID)
print("myUnion Areas:")
print (gArea(myUnion, byid=TRUE))
plot(myUnion,col=cols,axes=T)
mtext("myUnion",cex=2)
#Create example of desired output of "union" operation

p4<-SpatialPolygons(list(Polygons(list(Polygon(cbind(c(1,4,4,3,3,1,1),c(1,1,3,3,4,4,1)),hole

= F)), 1),

 Polygons(list(Polygon(cbind(c(4,4,5,5,3,3,4),c(4,3,3,5,5,4,4)),hole

=

F)),
2),

 Polygons(list(Polygon(cbind(c(3,4,4,3,3),c(3,3,4,4,3)),hole = F)),

3),

 Polygons(list(Polygon(cbind(c(3,4,4,3,3),c(3,3,4,4,3)),hole = F)),
4)))
ID<-c(1,2,3,4)
YEAR<-c(1999,2002,1999,2002)
Type<-c(1,2,1,2)
mydf4<- data.frame(ID,YEAR,Type)
mySPDF4<-SpatialPolygonsDataFrame(p4,mydf4)


centroids<-coordinates(mySPDF4)
mySPDF4$X_Centroid<-centroids[,1]
mySPDF4$Y_Centroid<-centroids[,2]
mySPDF4$MergeID<-paste(mySPDF4$X_Centroid,mySPDF4$Y_Centroid)
mySPDF4$MergeID<-unclass(as.factor(mySPDF4$MergeID))
mySPDF4$Area<-gArea(mySPDF4, byid=TRUE)
print("mySPDF4 at data")
print(mySPDF4 at data)
plot(mySPDF4, col=cols,axes=T)
mtext("mySPDF4",cex=2)
# the ultimate aim is to get a Flat ( ie non-overlapping polygons)
with
attributes that can then be joined back wide form table of to the

YEAR

and
Type of mySPDF4

#Flatten on "MergeID"
myUnion2<-unionSpatialPolygons(mySPDF4,IDs=mySPDF4$MergeID)
myDF<-as.data.frame(mySPDF4)
myDF<-arrange(myDF,MergeID,desc(YEAR))

#add sequnce number for subsequent conversion to wide form
for( row in 1:nrow(myDF)){
  thisVal = myDF$MergeID[row]
  if (row ==1) Seq = 1 else {if (thisVal==lastVal) Seq =Seq+1 else

Seq

=1}
  myDF$Sequence[row]<-Seq
  lastVal=thisVal
}
#make YEAR and Type into wide form
yearDF<-spread(myDF[,c("MergeID","Sequence",'YEAR')],Sequence,YEAR)
names(yearDF)<-c("ID","YEAR1","YEAR2")
typeDF<-spread(myDF[,c("MergeID","Sequence",'Type')],Sequence,Type)
names(typeDF)<-c("ID","Type1","Type2")
flatDF<-join(yearDF,typeDF)
rownames(flatDF)<-flatDF[,"ID"]
#append data frame to myUnion2 spatialpolygons
myUnion2<-SpatialPolygonsDataFrame(myUnion2,flatDF)
plot(myUnion2,col=cols,axes=T)
mtext("myUnion2",cex=2)
myUnion2$Area<-gArea(myUnion2, byid=TRUE)
print("myUnion2 at data")
print(myUnion2 at data)

        [[alternative HTML version deleted]]

        [[alternative HTML version deleted]]