How to get feature values from multiple polygons of shapefile while using extract function on raster

Question

I have a shapefile which contains arround 800 circular polygons with diameter of 60 meters each. Raster is 20 meter resolution so I would have multiple cell values extracted for each polygon. I am using an extract function on multiple raster layers in a for loop, and am putting the results in the dataframe. I know i can get the cellnumbers within the extract function, but I would also like to obtain the feature values (in my case a column in attribute table containing specific polygon code) directly from the shapefile for each pixel value extracted. Thus I would have in one row a pixel cell number, its extracted value and the polygon code.

my code:

dirs=list.dirs(full.names=TRUE)
dirs=dirs[grepl("R20",dirs)]
dirs=as.list(dirs)

DataToWrite=list()
sumirana=data.frame()
for (j in 1:length(dirs)) {
  setwd(paste("E:/",as.character(dirs[j]),sep="/"))
  files=list.files(full.names=FALSE,pattern="*20m*.tif$")
  for (i in 1:length(files)) {
      curRaster=raster(paste(getwd(),files[i],sep="/"))
      rasterOut=(na.omit(extract(curRaster,azo_plohe, cellnumbers=TRUE, 
      weights=TRUE, df=TRUE)))
      if((length(rasterOut) > 1)) {

         DataToWrite=c(DataToWrite,rasterOut[3])

         } else {

         lista=list(rep("x",15))
         tablica_folder=as.data.frame(t(unlist(lista)))      
         colnames(tablica_folder)=c(1:15)
         }
   }
   if((length(rasterOut) > 1)) {
      mylist=as.data.frame(DataToWrite)
      naziv_fajla=as.vector(rep((substr((names(DataToWrite)
      [1]),50,55)),nrow(mylist)))  
      tablica_folder=as.data.frame(c(rasterOut[2],naziv_fajla[1], 
      rasterOut[4],mylist),col.names=c(1:15))
      sumirana=rbind(sumirana,tablica_folder)
      DataToWrite=list()

      } else {
      sumirana=rbind(sumirana,tablica_folder)
      }
}

Answer 1

You have to unlist() the output of extract, and maintain that list-level grouping to know which object the value is from - this is a pain, and can be tricky for non-matches, so I put the workflow into package 'tabularaster' https://CRAN.R-project.org/package=tabularaster

You might find this is enough, where r is a raster and poly is a polygondataframe:

library(tabularaster)
cells <- cellnumbers(r, poly)
cells$value <- extract(r, cells$cell_)

So value is the pixel value, and cell_ is the pixel index, and object_ is the row-number of poly.

Then use object_ to get a value from the query layer:

## replace "polyID" with whatever column you want
cells$polyID <- poly$polyID[cells$object_]

For a case like this, making the cell_ explicit is overkill, but it works well for general extraction (i.e. time series coming in the door).

Answer 2

This would be easier to answer with a reproducible set of data. I've created some here:

library(rgdal)
library(raster)
library(sp)
library(rgeos)
#Create fake raster data
r1 <- r2 <- r3 <- raster(res=20)
values(r1) <- runif(ncell(r1))
values(r2) <- runif(ncell(r2))
values(r3) <- runif(ncell(r3))
s <- stack(r1, r2, r3)
#name the rasters
names(s) <- c("raster1","vegetation", "raster3")
plot(s)
#set random sample for 10 points to be converted to polygons
set.seed(1)
rn<-runif(10, min=-150, max=150)
set.seed(999)
rn2<-runif(10, min=-150, max=150)
df<-data.frame(x=rn, y=rn2, id=seq(10))
cor<-df[,c(1,2)]
#convert to points
p<-SpatialPointsDataFrame(data=df, coords=cor)
#buffer points to create polygons
pb<-gBuffer(p, width=5, byid=TRUE)

plot(s[[1]])
plot(pb,add=T)

Once you have your rasters identified from your directory structure
the code below will take values from your polygons and put them into a data.frame with the name & value of the raster and the cell #. Since your rasters may be in different extents, stacking may not work but see this question "Handing multiple extent problem to create raster stack in r"

  #extract values from all rasters for each polygon
  #add new field with "code values"
  pb$code<-sample(10:20,10)
  ext_poly<-extract(s,pb,cellnumbers=TRUE, df=TRUE)
  ext_poly
#join polygon data to extracted values
  m<-merge(ext_poly,pb@data)

Here's the output with raster cell, raster cell values, and the joined polygon data (X,Y and the new field created above with "codes")

ID cell    raster1 vegetation   raster3          x         y code
1   1  114 0.47134842 0.11526912 0.7113075  -70.34740 -33.27859   16
2   1   96 0.71419069 0.64174118 0.8939626  -70.34740 -33.27859   16
3   2   62 0.56123275 0.97922882 0.6705206  -38.36283  24.91822   15
4   2   61 0.32030622 0.16661933 0.8966194  -38.36283  24.91822   15
5   5    5 0.48085049 0.84698577 0.5859292  -89.49542  86.02403   19
6   5   NA         NA         NA        NA  -89.49542  86.02403   19
7   7   52 0.64265299 0.59516900 0.8185208  133.40258  31.93410   17
8   7   70 0.14084638 0.08905376 0.1137092  133.40258  31.93410   17
9   9  120 0.33380777 0.32627085 0.6227753   38.73421 -32.76827   11
10  9  119 0.57976869 0.11692465 0.3070686   38.73421 -32.76827   11
11  9  101 0.76347539 0.94087637 0.4202294   38.73421 -32.76827   11
12  9  102 0.31022532 0.17236697 0.5379526   38.73421 -32.76827   11
13 10   39 0.02508786 0.63041189 0.5683433 -131.46412  35.84172   10

Also in your code, I would clean up where you set your workspace (this should be outside of the loop) and don't use getwd in your loop either, set a variable for the path you need.