How to itterate over the bands raster and extract pixel value with gdal python

Question

I have this code from here:

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

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

I have a .vrt that has 7 bands. I want to iterate over each band and to extract the pixel values for each band and to save it to a nested list.

pixel_values=[['fid_1',value1,value2,value3,value4,value6,value7],['fid_2', value1,value2,value3,value4,value6,value7].. ..]

After that I want to export it to a csv.

Example:

pixel_values =[]
for band in range(1,8):
    pixel_value.append(intval[0])

pixel_value.to_csv('data.csv')

My question is that: How can I extract the pixel value based on a point .shp file from a raster with multiple bands?

Answer 1

I think the following should do it.

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 = []  # A list to hold all the raster bands
for band_num in range(1, 8):
    rb.append(src_ds.GetRasterBand(band_num))

ds = ogr.Open(shp_filename)
lyr = ds.GetLayer()

pixel_values = []  # List to store pixel values along with feature ID

for feat in lyr:
    fid = feat.GetFID()  # Get the feature ID
    geom = feat.GetGeometryRef()
    mx, my = geom.GetX(), geom.GetY()  # Coordinates 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

    # Read pixel values for all bands and store in a list along with feature ID
    pixel_values_for_feature = [fid]
    for band_num in range(1, 8):
        intval = rb[band_num - 1].ReadAsArray(px, py, 1, 1)
        pixel_values_for_feature.append(intval[0][0])

    pixel_values.append(pixel_values_for_feature)

# Print the resulting array
for row in pixel_values:
    print(row)

There are three major modifications. Firstly, changing rb into an array of bands. Secondly, retrieving the feature ID and as last change reading the value from each band and put it into an array.

Note: I haven't executed the code, but I think it should it.