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

[Michael Sumner]

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]]