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

attribute-tableextractrrastershapefile

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)
      }
}

Best Answer

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).

Related Question