1

I am new to R, and I am trying to use it to calculate the areas of polygons that overlap with each gridcell, in order to assign the attributes of the largest overlapping polygon area to the corresponding cell and build a dataframe, in which each gridcell represents a row, and attributes represent columns.

Gridcell is overlayed on top of polygons representing indigenous lands.

As shown in the figure above, I have a polygon shapefile representing indegenous lands and an overlayed gridcell (polygons).

Is there anyway I could do this? I have seen this question/answer: How to measure the area of polygons within a fishnet but for some reason the intersect function in QGIS (which is the GIS software I am using) does not generate a new layer.

Just in case, here you could download my spatial data: https://www.dropbox.com/sh/hdbpcq80uzd1ben/AADdG5Lg-BBbvFwA3QRLycSRa?dl=0

1

Here is a possible approach I came up with a while ago in response to a very similar question asked here on GIS StackExchange as Convert ESRI Shapefile to Raster using R - biggest polygon area to raster cell value? Although the latter was about rasterizing a shapefile based on the subpolygon covering the largest area inside each grid cell, only little changes are required to meet your demands. Please note that the sample data you provided is stored in the subfolder 'data' inside my current working directory.

## load packages
library(rgdal)
library(raster)
library(rgeos)

## import data
spy_poly <- readOGR(dsn = "data", layer = "CCNN_PNCB_Toy")
spy_grid <- readOGR(dsn = "data", layer = "grid_toy")

## per grid cell, identify ID covering largest area
lst <- lapply(1:length(spy_grid), function(i) {

  # create polygons subset based on current grid cell
  spy_poly_crp <- crop(spy_poly, spy_grid[i, ])

  # case 1: no polygon intersects with current cell
  if (is.null(spy_poly_crp)) {
    out <- data.frame(matrix(ncol = ncol(spy_poly), nrow = 1))
    names(out) <- names(spy_poly)
    return(out)
  # case 2: one polygon intersects with current cell  
  } else if (nrow(spy_poly_crp@data) == 1)  {
    return(spy_poly_crp@data) 
  # case 3: multiple polygons intersect with current cell
  # -> choose sub-polygon with largest area
  } else {
    areas <- gArea(spy_poly_crp, byid = TRUE)
    index <- which.max(areas)
    return(spy_poly_crp@data[index, ])
  }
})

## to 'data.frame'
do.call("rbind", lst)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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