[GIS] Extract raster pixel values from shapefile in QGIS

Question

I have been encountering difficulties trying to extract the raster pixel values from a polygon shapefile that has been assigned the same CRS.

The problem under inspection is how to get a list or a table with coordinates and pixel values from a polygon shapefile over a singleband raster for further investigation of values.

What is the easiest and most efficient way to perform this otherwise seemingly very simple task?

Answer 1

I would prefer using R to do this, is really easy with extract() function. In QGIS is a long process, but result:

Suppose you have both layers:

Clip (or crop) raster layer to have a small input to work:

Create an empty raster with same extent/resolution:

Convert vector to raster using a unique ID by geometry:

Merge both raster layers:

Convert each layer to a csv file or use gdal2xyz -band 1 -band 2 -csv /path/to/file.csv in cmd/terminal windows:

You'll obtain two .csv files (or one using command line):

Filter table by polygon ID:

---

R approach:

library(raster)

raster <- stack('~/Downloads/S029W072/AVERAGE/S029W072_AVE_DSM.tif')
poly <- shapefile('~/Desktop/eliminar/Poly.shp')

val <- extract(raster,poly)

The result is a list with n slots, each slot represents a polygon feature used to extract.

summary(val)
     Length Class  Mode   
[1,] 22667  -none- numeric
[2,]  9190  -none- numeric
[3,]  8212  -none- numeric

head(val[[1]])
     S029W072_AVE_DSM
[1,]              593
[2,]              588
[3,]              598
[4,]              607
[5,]              577
[6,]              586

---

Saving output in a csv file:

# I will use a field from my vector to create an idintifier (use unique values)
listnames <- poly$id

# create a empty list to save data frames to export
valList <- list()

# create as many data frames as features used to extract
for(i in 1:length(val)){
  valList[[i]] <- data.frame(ID=listnames[[i]],Value = val[[i]][,1])
}

# join all data frames and save it to an csv file
write.csv(do.call(rbind.data.frame,valList),"test.csv")

The output file: