[R-sig-Geo] Union/Overlay single spatialPolygonsDataFrame containing overlapping polygonsto create polygons with all unique combinations of attributes

[nevil amos]

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]]
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> >         [[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
>

	[[alternative HTML version deleted]]