[R-sig-Geo] Union/Overlay single spatialPolygonsDataFrame containing overlapping polygonsto create polygons with all unique combinations of attributes
Robert J. Hijmans
r.hijmans at gmail.com
Wed Aug 12 19:06:35 CEST 2015
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]]
>> >>
>> >> _______________________________________________
>> >> 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
>
>
More information about the R-sig-Geo
mailing list