I would like to get the areas related to the polygons resulting from the intersection between two SpatialPolygonsDataFrames: 'A' and 'B'. SpatialPolygonsDataFrame 'B' has 3 different 'classes', which represent polygons with different shape and size. In reality, 'A' is PRIO-GRID data (regular polygons) and 'B' is GREG ethnic database (with 3 layers that indicate different ethnic groups). In ArcGIS, the command is 'Tabulate Intersection" with option 'classes' where we can put the different class layers'. Is there any alternative in R?
-
1Could you explain more what you want to do. Maybe providing a reproducible example of your data? And the desired output? For example: is this questions related to what you are trying to do: gis.stackexchange.com/questions/64537/…? Because it might happen folks working with R do not know how the 'Tabulate Intersection' tool from ArcGIS works.– Andre SilvaApr 6, 2016 at 12:49
Add a comment
|
1 Answer
Here is how you can do that, I think. I am following this example:
library(raster)
library(rgeos)
library(rgdal)
# example data
p <- shapefile(system.file("external/lux.shp", package="raster"))[, 1]
p$Color <- rep(c('blue', 'green', 'red'), 4)
p <- p[,2]
z <- raster(p, nrow=2, ncol=2, vals=1:4)
names(z) <- 'Zone'
z <- as(z, 'SpatialPolygonsDataFrame')
# inspect
p
z
plot(p)
plot(z, add=TRUE, border='blue')
# intersect
i <- intersect(p, z)
# compute area
i$area <- abs(area(i))/1000000
# get the attribute table
d <- data.frame(i)
# aggregate and sum the areas
a <- aggregate(d[, 'area', drop=FALSE], d[, c('Color', 'Zone')], sum)
# get the total area by zone
aa <- aggregate(d[, 'area', drop=FALSE], d[, 'Zone', drop=FALSE], sum)
colnames(aa)[2] <- 'zonearea'
# merge that to the data
m <- merge(a, aa)
# compute percentage
m$percentage <- 100 * m$area / m$zonearea
# drop intermediate variable
m$zonearea <- NULL
m
## Zone Color area percentage
##1 1 blue 329.80536 39.52691
##2 1 green 391.08141 46.87080
##3 1 red 113.49501 13.60229
##4 2 blue 58.67825 33.56140
##5 2 green 79.69006 45.57924
##6 2 red 36.47019 20.85936
##7 3 blue 224.19645 28.88015
##8 3 green 79.20874 10.20338
##9 3 red 472.89422 60.91647
##10 4 blue 156.11708 20.03322
##11 4 green 379.33556 48.67702
##12 4 red 243.83826 31.28976