Thursday, 18 August 2016

python - Getting pixel value of GDAL raster under OGR point without NumPy?


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.



Answer



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.


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.
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


Alternatively, since the reason you gave for not using numpy was to avoid reading the entire array in using ReadAsArray(), below is an example that uses numpy and does not read the entire raster in.


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.

px = int((mx - gt[0]) / gt[1]) #x pixel
py = int((my - gt[3]) / gt[5]) #y pixel

intval=rb.ReadAsArray(px,py,1,1)
print intval[0] #intval is a numpy array, length=1 as we only asked for 1 pixel value

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...