5

I want to use GDAL/OGR to extract pixel values from a raster file. The extraction should be based on polygons, stored in a shapefile. My thought was to iterate over the features in the shapefile, extract the pixel values of those pixels that have there center point within the polygon and calculate a percentile from those values. When this is done I want to write the calculated percentile back to the polygon and go on with the next polygon.

Iterating over the polygons is easy, but what I do not get is how to extract the pixel values from the pixels, covered by the polygon.

I'm using Python.

Could somebody help me at this point?

Improve this question
6
  • If I iterate over the features and while one is 'selected', can I call the GetExtent function of the shapefile layer and will the returned extent only represent the extent of the current feature?
    – Thomas Becker
    Commented Aug 13, 2013 at 11:37
  • 1
    It sounds as though you are wanting to do zonal statistics. Have a look at gis.stackexchange.com/questions/43748/…"
    – MappaGnosis
    Commented Aug 13, 2013 at 11:56
  • @MappaGnosis That version of zonal stats appears extremely limited; I see no evidence from the documentation that it computes percentiles. Do you know of a zonal stats solution for QGIS that does compute specified percentiles by zone?
    – whuber
    Commented Aug 13, 2013 at 15:13
  • 1
    alright, the percentile is no problem for me since I had to define a function for it anyway. My issue goes with the Python coding using GDAL/OGR... I'm used to do the stuff using ArcGIS and their Python classes. However, it takes ages and in the end the script is crashing. So, I'm on my way to re-write the script without arcpy. Can anybody point me to some code snippet where I can learn on how to do the stuff?
    – Thomas Becker
    Commented Aug 13, 2013 at 16:04
  • 1
    @WHuber Offhand I don't know of a more extensive Zonal Stats Plugin for QGIS. However, my link was meant to point to the programmed Python + GDAL solution in my reply in that thread. My example there is just looking at getting the minimum value, but being a programmed solution, it shouldn't be too hard to extend it a bit to calculate percentiles of the zone data. Instead of going for the 'easy statistics' options there is also a scipy.ndimage.filters.percentile_filter function available.
    – MappaGnosis
    Commented Aug 14, 2013 at 15:07

1 Answer 1

Reset to default
3

Your question defines a complex process which is difficult to help with all of it. Here, I will suggest something that might help with the part getting the raster data using the vector.

Refer to this: http://machinalis8.rssing.com/chan-51920729/latest.php#item3

In that blog post I rasterize vector data to subset pixels from an image.

The idea there is that you convert the vector data into a raster, like a mask. Then, using numpy, you use that mask to filter or choose the pixels in the raster.

You may need to check the part about "the center of the pixel" fallin within the poligon: it will depend on GDAL's implementation of RasterizeLayer.

Next, you can do the rest of the calculations.

Improve this answer
3
  • 3
    +1 Nice use of random forests in your example. It would be helpful to include an example on this site in the event that the link goes away.
    – Aaron ♦
    Commented Mar 15, 2016 at 23:15
  • That link doesn't work anymore sadly
    – Cam
    Commented Mar 28, 2019 at 12:45
  • 1
    @cam link fixed to an archive version
    – Akhorus
    Commented Apr 1, 2019 at 21:21

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.