[GIS] Extraction of raster with shapefile python

gdalgeotiff-tiffpythonrastershapefile

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)

Best Answer

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!