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 can calculate a buffer directly with ogr2ogr using GDAL >= 1.10 with SQLite and SpatiaLite support. If the raster spatial resolution is e.g. 1m:
ogr2ogr -dialect SQLite -sql "SELECT ST_Buffer(geometry,2.5) FROM lines" buffer_lines.shp lines.shp
Finally, you can adapt and apply the first example of gdal_rasterize documentation:
gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 0 -burn 0 -l buffer_lines buffer_lines.shp work.tif
Best Answer
GDAL version <3.2
To set the pixel size in
gdal_grid, you have to play with the following options:because:
New in GDAL version 3.2.