Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am working on a computational model of the abundance of wild pollinators across a landscape. The model itself is complete, and I am now struggling with a post-processing step.

I have my GDAL pollinator supply raster that looks something like this (lighter colors mean higher pollinator visitation to a pixel):

Greyscale raster representing pollinator supply on a landscape

And I have an OGR shapefile of points representing sample locations on the landscape:

enter image description here

I'm trying to run some analysis on the pixels under these points, but to do so, I need to be able to extract the value of a pixel under a point.

Is it possible to extract the value of a pixel under a point using only OGR and GDAL through Python? I would prefer to avoid reading the entire raster into memory through ReadAsArray(), as my output rasters are very, very large (too large to fit into memory).

I noticed this post, which is similar, but requires a command-line call.

share|improve this question
1  
What about ReadAsArray() and only reading in the point? So only read the single cell that you are interested in? You would need to convert from the point coords to pixel space and extract the necessary cell. –  Jay Laura Jan 11 '13 at 22:14
1  
Look at the code for gdalsrsinfo, it shows you how to use GDALInvertGeoTransform() and switch between geographic space and pixel space: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp –  kyle Jan 11 '13 at 23:18
    
If you don't mind using PostGIS, see this. It is extremely fast and is just 1 SQL line. –  mlt Jul 16 '13 at 18:57
    
I'll keep that in mind if I come across this problem and have access to a PostGIS database! I didn't for this particular problem, so the GDAL solution below did the trick. Thanks, though! –  James Aug 5 '13 at 19:03
    
@kyle I don't know if things have changed but it looks like it is GDALInvGeoTransform not invert and this is an example. –  mlt Sep 19 '13 at 22:29

1 Answer 1

up vote 17 down vote accepted

You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module.

See also the QGIS "Point Sampling Tool" plugin for a GUI way of doing this.

An example:

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    #If raster is rotated, see http://code.google.com/p/metageta/source/browse/trunk/metageta/geometry.py#493
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
    intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

    print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
share|improve this answer
    
Awesome. Works like a charm! Thanks! –  James Jan 15 '13 at 1:12
    
This is amazing. Thanks! –  bananafish Jan 30 '13 at 23:30

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.