Extraction of raster with shapefile python

Question

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

Answer 1

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

        intval=rb.ReadAsArray(px,py,1,1)
        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!

Answer 2

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:

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