1

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?

  • Can you expand this a little bit please? What have you tried? What statistic from the raster are you trying to calculate for your polygons (average/max etc)? What error messages are you having? – Ed Rollason Jul 25 '17 at 13:54
  • Do you need to extract each pixel value or statistics? – aldo_tapia Jul 25 '17 at 13:54
  • @aldo_tapia I need to extract each pixel value, not statistics. – ODstuck Jul 25 '17 at 14:00
  • @EdRollason I cannot really explain here all the tries I have made and the walls I have hit, really.. Most of them have been huge transformations from raster to vector and reverse, difiiculties with CRS recognition inside QGIS for certain tools etc. I just need something solid to try. – ODstuck Jul 25 '17 at 14:03
3

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:

enter image description here

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

enter image description here

Create an empty raster with same extent/resolution:

enter image description here

Convert vector to raster using a unique ID by geometry:

enter image description here

Merge both raster layers:

enter image description here

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

enter image description here

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

enter image description here

Filter table by polygon ID:

enter image description here


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: enter image description here

  • Thank you very much on the detailed answer sir, however I seem to have different tools. I am currently working under QGIS 2.18.11. I will try to work my way around and maybe try extract() by R as you proposed which may be a lot simpler indeed – ODstuck Jul 25 '17 at 14:58
  • 1
    I edit my answer adding the code for doing the same in R – aldo_tapia Jul 25 '17 at 15:09
  • I installed R, added the libraries needed, and tried the following script (the same as yours but with the addition of initial input and output in order to be able to save it) and I believe it breaks for some reason. ##poly=vector ##val=output table library(raster) raster <- stack('_____/clipExtentNewLayer.tif') poly <- shapefile('_____/Polygons.shp') val <- extract(raster,poly) However, I tried firstly ending the code with >val in order to see the results and it loaded the hundreds of values (I believe correctly) until QGIS just collapsed. Do you maybe have any insight on this? – ODstuck Jul 25 '17 at 17:56
  • Oh, and no output saved at location set – ODstuck Jul 25 '17 at 17:59
  • @ODstuck ok, I put you the code to export values. You can check them in excel or another software – aldo_tapia Jul 25 '17 at 19:22
0

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.

  • Thank you for the answer . Of course it will not work under different CRS. However I do not need the Sampling tool, I need all the values under the polygon. – ODstuck Jul 25 '17 at 14:09
0

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.

Raster statistics for polygons window

  • Thank you, but this again is statistics output and I need all the pixel values... – ODstuck Jul 25 '17 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.