I have one shape file containing lots of SpatialPolygons which are partly overlapping. These polygons belong to the application of a fungicide on a field and each polygon has an associated application rate as attribute.
What I want to obtain is to correct AsApplied map taking into account the overlapping areas meaning that the rate should be summed up and merged if two (or more) polygons are overlapping.
The following example code creates a SpatialPolygonsDataFrame simplifying the problem:
library(raster)
library(sp)
p<-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_")))
pid <- sapply(slot(p, "polygons"), function(x) slot(x, "ID"))
p.df <- data.frame( ID=1:length(p), row.names = pid)
p <- SpatialPolygonsDataFrame(p, p.df)
p$Rate <- c(100, 100, 100, 100)
crs(p) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
plot(p)
You can see two squares from four polygons which are partly overlapping. Each polygon has an associated rate of 100. What I would like to have is three polygons. The two non overlapping ones should have the rate 100 and the two overlapping ones should be joined to one polygon having a value of 200.
I already tried the union or intersect functions of the raster package but was only able to get the information which polygons overlap but not the summing and merging. In addition I am seeking explicitly for a solution in R.
Any help solving this problem is highly appreciated.
Update: The solution provided by RobertH provided below works for my simple example. Thank you very much already!
However when switching to my real usecase I am getting to following kind of errors and warnings:
Error in if (is.numeric(i) && i < 0) { :
missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
Too few points in geometry component at or near point 8.3634020800000002 50.056772690000003
...
An example shape file is uploaded here: (obsolete)
Any ideas how to deal with this problem?
Update #2 Using the current development version 2.5-10 indeed fixes the warnings in RGEOSUnaryPredFunc. However, if polygons are only overlapping very very little I am still getting the error:
Error in if (is.numeric(i) && i < 0) { :
missing value where TRUE/FALSE needed
An example shape file for which this is happening is uploaded here: http://www.share-online.biz/dl/O4ZIVH8OBW. More precise the field looks like the following:
The two polygons marked in red cause the error and if one of the two is removed the union works fine.
Thank you very much already for your great help!
rgeos
? I would be happy to look at other cases, but not via that annoying fileshare. You can email me if you want. – Robert Hijmansrgeos
but thesp
package I had to update. So thank you very much for your help! :) – Aludra