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

This question already has an answer here:

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. ## marked as duplicate by Andre Silva, PolyGeo♦Dec 13 '16 at 0:06

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

• 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

## 1 Answer

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