[R-sig-Geo] Union/Overlay single spatialPolygonsDataFrame containing overlapping polygonsto create polygons with all unique combinations of attributes
nevil amos
nevil.amos at gmail.com
Thu Aug 13 05:23:05 CEST 2015
Thanks Robert,
I was getting an error when I ran it, and no mention of the case in the
help.
I then realised that my raster package was not up to date.
I have now updated and it works fine.
On Thu, Aug 13, 2015 at 3:06 AM, Robert J. Hijmans <r.hijmans at gmail.com>
wrote:
> 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
> >
> >
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list