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

[Bernd Vogelgesang]

Hi,

I don't know if that is of any help, but right at the moment I am building  
a model in QGIS which should achieve sth similar:
I have overlapping polygons with minute values per square meter, which I  
self-intersect, calculate the new area and then dissolve the intersecting  
parts while summing up the minute values.
Though this seems to me like a quite common task in GIS, its really hard  
to implement somehow.


The self-intersection in my model is done through v.clean break from  
GRASS. This "cuts" all polygons with the borders of all other overlapping  
polygons.
The dissolve is done with a gdal function which looks like this:
ogr2ogr output input -dialect sqlite -sql "SELECT  
ST_Union(geometry),dissolve_field, SUM(sum_field) AS sum_diss GROUP BY  
dissolve_field"

I haven't tested any of your code, so this is just a hint in case you  
can't find any native possibilities. The grass function might be quite  
easy to implement, but I have no idea about the dissolve from gdal.

Good luck
Bernd

Am 10.08.2015, 09:51 Uhr, schrieb nevil amos <nevil.amos at gmail.com>:

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


-- 
Bernd Vogelgesang
Siedlerstraße 2
91083 Baiersdorf/Igelsdorf
Tel: 09133-825374