2

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.

Forest area polygons (in blue) overlayed with ordinance survey grid cells (in pink). Both layers are polygon shapefiles.

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

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"])
  • 2
    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

Not the answer you're looking for? Browse other questions tagged or ask your own question.