Union/Overlay single spatialPolygonsDataFrame containing overlapping polygonsto create polygons with all unique combinations of attributes
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