# Finding the area of polygons within other polygons in R [duplicate]

Is there a function in R which will loop through a layer of zonal polygons (i.e. a polygon grid) and return the area of a second polygon layer that is contained in each grid cell? I have a layer of polygons representing areas of forest (in blue), and a second layer representing ordinance survey grid cells (in pink - each grid cell has a unique "name"). Both layers are polygon shapefiles. I would like to calculate how much forest is present in each grid cell but I cannot find a tool that will do this in R. I'm looking for something equivalent to tabulate area in ArcGIS. As you can see, the forest polygons are complex and in many cases there are several of them per grid cell, which is making it trickier. • rgeos::gIntersection with byid=TRUE or raster::intersect both will give you what your after. To return area you can do something like: sapply(slot(mypolys, "polygons"), function(i) slot(i, "area")) – Jeffrey Evans Dec 3 '15 at 17:04
• or use rgeos::gArea – mdsumner Dec 5 '15 at 0:22

Here's a way.

PLEASE NOTE: this calculates area in degrees, which makes no sense - but otherwise doesn't impact the illustration. You should calculate area in a coordinate system that is sensible for your data. (Use rgdal to project to a local Lambert Azimuthal Equal Area projection or similar if your data are in longitude/latitude, or otherwise unsuitable coordinate system).

`````` library(maptools)
data(wrld_simpl)
library(rgeos)

## make a grid of polygons on the space
library(raster)
g <- as(raster(extent(-180, 180, -90, 90), nrow = 9, ncol = 12, crs = projection(wrld_simpl)), "SpatialPolygonsDataFrame")
g\$pieces_area <- NA_real_
print(nrow(g))
for (i in seq(nrow(g))) {
asub <- gOverlaps(g[i, ], wrld_simpl, byid = TRUE)[,1]
if (any(asub)) {
x <- gIntersection(g[i, ], wrld_simpl[asub, ])
if (!is.null(x)) {
g\$pieces_area[i] <- gArea(x)
print(i)
}
}
}

## simple plot of the result
#spplot(g["pieces_area"])
``````
• To compute the area of lon/lat polygons you can use `'geosphere::areaPolygon` – Robert Hijmans Dec 6 '15 at 1:27
• I know that I'm just kind of philosophically opposed in general. I'll check it out ;) – mdsumner Dec 6 '15 at 6:09