-1

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')

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

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Duststar
  • 11
  • 3
  • 1
    So what's your question? GIS SE is a problem-solving site. If your problem statement doesn't contain a problem, just a list of tasks which may or may not be addressed by the provided code (which uses a non-standard indent scheme), then it's difficult to anwser. – Vince Jul 16 '23 at 18:16
  • 1
    Please do not keep re-posting closed questions. Fix your original one. – user2856 Jul 17 '23 at 01:49
  • Please merge your two accounts and do not use them to ask the same question twice. – PolyGeo Sep 21 '23 at 11:02

1 Answers1

0

If you do not definitely want to write your own code, I would recommend to use the GDAL utility "gdallocationinfo" https://gdal.org/programs/gdallocationinfo.html.

This utility is intended to provide a variety of information about a pixel. Currently it reports:

The location of the pixel in pixel/line space.

The result of a LocationInfo metadata query against the datasource. This is implement for VRT files which will report the file(s) used to satisfy requests for that pixel, and by the MBTiles driver

The raster pixel value of that pixel for all or a subset of the bands.

The unscaled pixel value if a Scale and/or Offset apply to the band.

The pixel selected is requested by x/y coordinate on the command line, or read from stdin. More than one coordinate pair can be supplied when reading coordinates from stdin. By default pixel/line coordinates are expected. However with use of the -geoloc, -wgs84, or -l_srs switches it is possible to specify the location in other coordinate systems.

The default report is in a human readable text format. It is possible to instead request xml output with the -xml switch.

For scripting purposes, the -valonly and -lifonly switches are provided to restrict output to the actual pixel values, or the LocationInfo files identified for the pixel.

user30184
  • 65,331
  • 4
  • 65
  • 118
  • I'm new to programming and still I don't know to write my own code. Thank you for your answer but a code answer would help me a lot. – Duststar Jul 16 '23 at 19:01
  • 1
    Gdallocationinfo is C++, but maybe the source code can still be useful https://github.com/OSGeo/gdal/blob/master/apps/gdallocationinfo.cpp. And actually there is also a Python version https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py. – user30184 Jul 16 '23 at 19:37