Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Simple script to calculate NDVI using pure GDAL
from osgeo import gdal
import numpy
g = gdal.Open ("baikal_subset.tif")
if g is None:
raise IOError, "Couldn't open baikal_subset.tif"
b3 = g.GetRasterBand(3).ReadAsArray().astype(np.float32)
b4 = g.GetRasterBand(4).ReadAsArray().astype(np.float32)
ndvi = (b4 - b3)/(b4 + b3)
drv = gdal.GetDriverByName ( "GTiff" )
dst_ds = drv.Create ( "output_file.tif", g.RasterXSize, g.RasterYSize, 1,
gdal.GDT_Float32, options=["COMPRESS=LZW"] )
dst_ds.GetRasterBand(1).WriteArray ( dst_ds.astype (np.float32) )
dst_ds = None
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.