Finding the area of polygons within other polygons in R

Question

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 03 '15 at 17:04

score 1 · Answer 1 · answered Dec 05 '15 at 00:38

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 06 '15 at 01:27
I know that I'm just kind of philosophically opposed in general. I'll check it out ;) — mdsumner, Dec 06 '15 at 06:09

Finding the area of polygons within other polygons in R

1 Answers1