# Extraction of raster with shapefile python

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 = np.array(nd_array, dtype=float)
osgeo.ogr.UseExceptions()
myShapefile = r"/home/cr.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(Shapefile, 0)
dataSource = ogr.Open(Shapefile)


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.

https://github.com/perrygeo/python-rasterstats

Otherwise you find an example without this module here:

https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

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'
src_ds=gdal.Open(rs)
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(src_shp)
lyr=ds.GetLayer()
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])


with rasterio.open(raster) as src: