Get Point WKT from a point shapefile using ogr
from osgeo import ogr
import os
shapefile = "../points.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
for feature in layer:
geom = feature.GetGeometryRef()
print geom.Centroid().ExportToWkt()
POINT (624501.568181818 5628082.59695052) POINT (656596.878308401 5630027.76726122) POINT (694388.75863061 5628499.41915995) POINT (718425.506041427 5639336.79660529) POINT (725233.602128884 5664901.89211738)
Use rasterstats (pip install rasterstats) to get the value at a point location.
from rasterstats import point_query
point = "POINT (725233.602128884 5664901.89211738)"
point_query(point, "../NDVI_example_1.tif")
[0.6670739601440762]
Loop over a point shapefile and extract values on a single band raster (but a multi band one could be used)
#from osgeo import ogr
#import os
shapefile = "D:/rasterStats/points.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
for feature in layer:
geom = feature.GetGeometryRef()
point = geom.Centroid().ExportToWkt()
print point_query(point, "..NDVI_example_1.tif")
[-0.40849066211442864] [0.652933816283167] [-0.4690519447620269] [0.3553161198672801] [0.6670739601440762]
Get mean and max values from a series of polygons from a shapefile with multiple features
from rasterstats import zonal_stats
stats = zonal_stats("../polys.shp", "../NDVI_example_1.tif")
print stats[1].keys()
#['count', 'min', 'max', 'mean']
print [f['mean'] for f in stats]
print [f['max'] for f in stats]
['count', 'max', 'mean', 'min'] [0.17415108794477352, 0.6178150208015752, 0.2820635024507548, 0.4951216243519921] [0.3866455554962158, 0.683485209941864, 0.45673757791519165, 0.6485721468925476]
Walk through all your .tif images in a directory and extract values
for root, dirs, files in os.walk("../"):
for file in files:
if file.endswith(".tif"):
raster = (os.path.join(root, file))
stats = zonal_stats("../polys.shp", raster)
print "filename: ", file
print [f['mean'] for f in stats]
filename: NDVI_example_1.tif [0.17415108794477352, 0.6178150208015752, 0.2820635024507548, 0.4951216243519921] filename: NDVI_example_2.tif [0.7579913160863959, 0.6720901117092226, 0.7985192637380338, 0.6334189456158745]