Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

PS: I'm using Python.

Could somebody help me at this point?

Cheers Thomas

share|improve this question
    
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 Aug 13 '13 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 Aug 13 '13 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 Aug 13 '13 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 Aug 13 '13 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 Aug 14 '13 at 15:07

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://www.machinalis.com/blog/python-for-geospatial-data-processing/

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.

Hope it helps.

share|improve this answer

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

    
+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 Mar 15 at 23:15

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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