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 math import floor
from osgeo import gdal,ogr
import struct
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
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.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #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. It uses the built-in gdal.ApplyGeoTransform() function in order to deal with axes rotations.
from math import floor
from osgeo import gdal,ogr
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
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.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #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
You could try only doing the Within() test against points that are within the polygon features envelopes, i.e.
def PointsInPolygons(pointsLayer, polysLayer):
nbInside = 0
polyFeat=polysLayer.GetNextFeature()
while polyFeat:
polyGeo = polyFeat.GetGeometryRef()
pointsLayer.SetSpatialFilter(polyGeo) #<----Only test points within feature envelope
pointFeat = pointsLayer.GetNextFeature()
while pointFeat:
pointGeo = pointFeat.GetGeometryRef()
if pointGeo.Intersects(polyGeo):
nbInside+=1
pointFeat = pointsLayer.GetNextFeature()
polyFeat=polysLayer.GetNextFeature()
return nbInside
Best Answer
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.