extract() while keeping polygon ID in output

Question

im trying to extract the mean from a satellite imagery. The point is, that when Im using this code (simplified for this example):

    > semana21 <- raster("Semana21.tif")
    > poly <- st_read("shapefile.shp")

Reading layer `shapefile' from data source `shapefile.shp' using driver `ESRI Shapefile'
Simple feature collection with 43 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -6.18 ymin: 38.89 xmax: -2.72 ymax: 42.90
epsg (SRID):    NA
proj4string:    NA

    > extract(semana18, st_zm(poly), fun=sum, na.rm=TRUE, df=TRUE)

Extract seems to auto-generate ID in extract() output dataframe. So, extract() generate a new row for each polygon, but extract() row position does not correspond with the shapefile row position. So I dont have way to know which output corresponds to each element in the shapefile. So, Is there any way to keep "name" in the generated extract() dataframe? Or another way to know which result correspond to each polygon?

I want a dataframe output like: (name, extract_calculation) That way, I can then know which extract operation belong to each polygon in file

Answer 1

It looks like your problem is that the polygon data is multipart geometry. This means that you have multiple features (polygons) associated with single rows (attributes). Even if raster::extract works, this makes very little sense from a results standpoint. For your data to match, you need to explode your geometry into single part.

Here is an example, you will need the spData package (on CRAN) for the example multipart polygon data.

library(sf)
library(raster)

Read the multipart polygon (MULTIPOLYGON) data and verify geometry.

poly <- st_read(system.file("shapes/sids.shp", package = "spData")[1]) 
  unique(as.character(sf::st_geometry_type(poly)))

Create an associated raster for extracting values.

r <- raster(extent(poly), nrow=200, ncol=200)
  r[] <- runif(ncell(r))

Now, lets extract the raster values and output the mean. You will see that the number of resulting mean values is 100.

pm <- as.numeric(extract(r, poly, fun=mean))
  length(pm)

The number of features in the poly data matches the number of extracted mean values. However, if we explode the multipart data into singlepart you will see that there are, in fact, 108 polygons. Just ignore the warring issued from st_cast.

dim(poly) 
poly.single <- sf::st_cast(poly, "POLYGON")  
dim(poly.single)  

If we calculate the raster mean vlaues on the singlepart polygons we now get a one-to-one match to the features.

pm <- as.numeric(extract(r, poly.single, fun=mean))
  length(pm)

Since the results from extract are ordered, we can now assign the raster mean values to the polygon data.

poly.single$rmean <- as.numeric(extract(r, poly.single, fun=mean)) 
  head(poly.single)

Answer 2

You could use exactextractr::exact_extract for this. It will work with multi-part polygon inputs.

library(raster)
library(sf)
library(exactextractr)

semana21 <- raster("Semana21.tif")
poly <- st_read("shapefile.shp")

poly$sum_semana21 <- exact_extract(semana21, poly, 'sum')

Note that these results will differ from raster::extract because the summation will include the values of all pixels wholly or partially covered by the polygon, whereas raster::extract will include only pixels whose center is covered by the polygon.