I am having some difficulties in R trying to calculate the area of overlap between polygons. As an example, for Switzerland I have data on the administrative boundaries of the cantons and the ethnic groups in the country.
## Libraries
library(sp)
library(spdep)
library(maptools)
library(rgeos)
## Data
print(load(url("http://gadm.org/data/rda/CHE_adm1.RData")))
greg<-readShapeSpatial("raw_data/GREG.shp",
proj4string=CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
switzerland<-greg[greg$COW==225,]
## Identical projections
proj<-CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
switzerland<-spTransform(switzerland,proj)
Plotting the data and we get something that looks like this:
## Plot data
par(mar=c(1,1,1,1))
plot(gadm,col="grey80")
plot(switzerland,add=TRUE,lwd=2,border="red")
We can see that the boundaries for the ethnic groups don't follow the national borders entirely, but good enough. What I am trying to do is for each of the cantons calculate the number of ethnic groups, taking into account the area of the ethnic group within the canton. So for Graubünden in the East I want to know the area occupied by German Swiss, Italian Swiss, Rhaetoromaniens, etc..
Reading some of the similar questions here on stackoverflow I thought that gIntersection
is the command to use, but this gives me the following error:
int<-gIntersection(gadm,switzerland,byid=TRUE) # Updated
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_intersection") :
TopologyException: no outgoing dirEdge found at 7.3306802801147688 47.439399101791921
In addition: Warning message:
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_intersection") :
spgeom1 and spgeom2 have different proj4 strings
Not entirely sure what causes the second warning mission though since
identicalCRS(gadm,switzerland)
[1] TRUE
Does anyone have a suggestion on how I could calculate the overlap between the cantons and the ethnic groups?
Update: Possible solution
This might be a possible solution, although the warning on different proj4 strings persists. Also note that there is some measurement error (in Aargua for instance) due to the fact that the ethnic groups shape file does not follow the national borders correctly.
## Possible solution
int<-gIntersection(gadm,switzerland,byid=T) # Multiple polygons per province
n<-names(int)
n<-data.frame(t(data.frame(strsplit(n," ",fixed=TRUE))))
colnames(n)[1:2]<-c("id0","ethnic.group")
n$area<-sapply(int@polygons, function(x) x@area)
a<-data.frame(id0=row.names(gadm),total.area=gadm$Shape_Area)
df<-merge(n,a,all.x=TRUE)
df$share.area<-df$area/df$total.area*100