I have a raster tiff file containing vegetation index & also a shapefile of fields. I want to extract the field wise values from that raster file using GDAL. How can I do this?

I can read the shapefile in python using gdal, but don't know how to extract its coordinates and how to extract that polygon from the raster tiff file.

import optparse
from osgeo import gdal, ogr, osr
import numpy as np
g_nd = gdal.Open('ndvi.tif')
nd_array = g_nd.ReadAsArray()
nd_array = np.array(nd_array, dtype=float)
myShapefile = r"/home/cr.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(Shapefile, 0)
dataSource = ogr.Open(Shapefile)

2 Answers 2


If you want GDAL you can use this (I found this code in another post and it works nice):

rs = r'path\to\your\raster'

def get_point_values(rs):

    src_shp = r'path\to\your\shape.shp' 

    for feat in lyr:
        geom = feat.GetGeometryRef()
        mx,my=geom.GetX(), geom.GetY()
        print mx, my
        px = int((mx - gt[0]) / gt[1]) 
        py = int((my - gt[3]) / gt[5]) 

        print intval[0] 

I prefer use rasterio with fiona to get the coordinates and use the sample tool:

with rasterio.open(raster) as src:
    for val in src.sample([(x, y)]): #You can loop through shape coordinates 
        print float(val)

Both solutions are just for points!


There is a python module called rasterstats that does compute zonal statistics. For each Polygon, your fields, there will be the underlying rastercells accumulated by SUM, MEAN, MIN or MAX.

You find more infos and examples here:


Otherwise you find an example without this module here:


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.