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:

Answer 2

Using polygon shapefile you cannot extract each pixel value within that polygon. You can only get the statistics using Zonal Statistics tool from Rater -> Zonal Statistics.

If you want to extract each pixel value, you need to use point shapefile not polygon and then you can use Point Sampling tool to extract the pixel value under each point. The CRS for both Raster and Point data should be same. It will not work using different CRS.

As a workaround, you can rasterize the polygon with reasonable pixel size, and then convert the rasterized polygon into points and use Point Sampling tool mentioned above.

Answer 3

You can use one of the SAGA functions in QGIS or in SAGA itself. In QGIS, go to your Processing Toolbox and look for "Raster statistics for polygons". Then, you set your raster files (yes, you can do more than one at a time, ideal for time series) and your shapefile (polygons). Finally, you choose how to aggregate the extract data, let's say, the mean pixel values within a polygon. The processing is wickedly fast, I tried other methods like python-rasterstats, R-extract and even SQL... SAGA stills being the fastest. Just make sure that you are processing a large number of rasters at the same time, otherwise. your computer can run out of memory, so, don't be greedy.