Load Raster Values into Shapefile

Question

I have a shapefile containing many points.

Some of these points fall within an area for which I have raster data.

I would like to load the values of the raster which underlie each point into the shapefile such that they appear as a column in its DBF.

Is there a way to do this with, say, GDAL?

It is possible in Python with the GDAL/OGR module – gene Nov 08 '13 at 17:53 — gene, Nov 08 '13 at 17:53

score 2 · Answer 1 · edited Apr 13 '17 at 12:34

Starting from:

...and gluing all the pieces together:

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,1)
lyr=ds.GetLayer()
# Add a new field
new_field = ogr.FieldDefn('VALUE', ogr.OFTInteger)
lyr.CreateField(new_field)

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.
    #If raster is rotated, see http://code.google.com/p/metageta/source/browse/trunk/metageta/geometry.py#493
    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
    feat.SetField('VALUE', intval[0]) 
    # trigger the update
    lyr.SetFeature(feat)

# Close the Shapefile
ds = None

Load Raster Values into Shapefile

1 Answers1