4

I want to be able to write pixel values of each pixel falling under a polygon onto a polygon attribute table or a CSV. I have about 1200 polygons and a georefernced GeoTiff file with NDVI values. While I understand there are many tools available to calculate zonal stats and I have tried zonal statistics, LecoS and rasterstats plugin (all QGIS), they all provide only summary statistics (mean, median, min, max, etc.) but not individual pixel values. So pretty much are of no use here. I also took a look at rasterstats source code but could not find a quick way to modify the zonal output parameters to suit my needs. So I'm wondering if there is a GDAL, QGIS, python based solution out there?

TL;DR: I would like to write pixel values of all the pixels falling under a poylgon onto a CSV or a shapefile attribute table (1200 polygons).

Improve this question
5
  • Are you prepared to script this? The Esri ASCII format has what you want with a few header rows, some simple file/text manipulation will do it. Do you already have some pyQGIS code? I have done something similar using GDAL_Translate in a subprocess then file to file - how does that sound to you?
    – Michael Stimson
    Commented Sep 22, 2015 at 1:20
  • @MichaelMiles-Stimson I do have a basic understanding of pyQGIS but haven't played too much with it. In case it is needed I would be open to giving it a shot. Right now I am having difficulty in getting the correct pixel values using the export to xyz format option.
    – SinghD
    Commented Sep 24, 2015 at 18:55
  • How are you exporting to XYZ? Raster translate would be my pick (it's a GUI for gdal_translate) with the output format XYZ. Out of interest what about the cells outside the polygons but inside the extent? Rasters must be rectangular so there will be situations where cells are in the extent of a polygon but not inside. When I did this previously I used the rectangular extent of each polygon to generate the rasters as outside cells were considered valid.
    – Michael Stimson
    Commented Sep 24, 2015 at 22:52
  • I am using raster > conversion > translate. I used a small raster clip this time and it worked fine. But now the problem is with spatial join. Join attributes by selection doesn't show xyz layer in the option. Strangely, MMQGIS plugin is also crashing on me.
    – SinghD
    Commented Sep 24, 2015 at 23:05
  • You need to make your XYZ data into a points feature type (step 1) have a look at gis.stackexchange.com/questions/73804/…
    – Michael Stimson
    Commented Sep 24, 2015 at 23:14

1 Answer 1

Reset to default
3

If you want to get all pixel values as value per row there is one way to do it easily though might be pretty high demanding on memory / data processing time if you have huge grid dataset.

The way is to convert raster to points and then select points by polygons.

  1. Save raster as xyz format and load as delimited text values as points (also can be converted from raster to points through processing toolbox - Grass / r.to.vect)
  2. use Vector / Data management Tools / Join attributes by location to assign polygon ID value into points table
  3. delete points with empty (Null) polygon ID column and save points layer as CSV or select all in points table and paste into your favorite spreadsheet editor
Improve this answer
5
  • 1
    That was my first thought too, but I used Raster->Point and then iterative clip but it took for ever, it would probably help to divide the input geometries into batches and run one batch per CPU thread... In the end I got the extents of each polygon, used GDAL_Translate -projwin to format AAIGRID then opened each grid, skipped the header lines and formatted the data with replace(" ",",") to CSV, this was much faster.
    – Michael Stimson
    Commented Sep 22, 2015 at 3:03
  • 1
    @MichaelMiles-Stimson Using python / gdal can definitely help here a lot talking about processing time, but if the grid is not huge and the number of polygons 1200, it should be ok to do it even without programming skills. With python / gdal I would loop through polygons, for every polygon clip the source grid and export clipped non-null values into text file. All depends on skills and how big data we are talking about.
    – Miro
    Commented Sep 22, 2015 at 3:17
  • 100% Miro! It depends on whether the effort of scripting is going to be less than geoprocessing, whether the OP has the understanding and skills and, of course, whether it's acceptable to leave a process going overnight/few days or if the result is needed sooner. That is always the cost/benefit of scripting something that can be geoprocessed. Your method works, which is why you've got my upvote.
    – Michael Stimson
    Commented Sep 22, 2015 at 3:26
  • @Miro When I export raster to xyz, it shows the pixel values as negative numbers which is supposed to be 0 to 1 for NDVI. Any thoughts ?
    – SinghD
    Commented Sep 24, 2015 at 18:48
  • @plant_GIS By negative you mean proper raster * -1? Is it possible that your raster has more bands? If yes then while exporting to xyz you need to edit gdal_translate and use -b as explained here: gdal.org/gdal_translate.html
    – Miro
    Commented Sep 25, 2015 at 0:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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