I need a script that can calculate NDVI from two separate input Geotiff files, then output the results as a Geotiff. I wrote following script, most of which is borrowed from various on-line sources (I've never had any training with python and only minimal experience). I know the calculation part works, because if I print out the ndvi array, I get reasonable values. If I calculate statistics on the array, I also get reasonable results. The script does create an output file called ndvi.tif, but all the values in the file are 0 ( and not 0.0).
I need to convert this into a batch script that can do this calculate on a large number of input files, but right now I can't figure out why it won't output calculated values to a file. Any advice?
+++++++++++++++++++++++++++++++++++++++++++++ from osgeo import gdal import numpy as np from numpy import * g = gdal.Open("fhmod1.tif") red = g.ReadAsArray() g = gdal.Open("fhmod2.tif") nir = g.ReadAsArray() red = array(red, dtype = float) nir = array(nir, dtype = float) check = np.logical_and ( red > 1, nir > 1 ) ndvi = np.where ( check, (nir - red ) / ( nir + red ), -999 ) geo = g.GetGeoTransform() proj = g.GetProjection() shape = red.shape driver = gdal.GetDriverByName("GTiff") dst_ds = driver.Create( "ndvi.tif", shape, shape, 1, gdal.GDT_Float32) dst_ds.SetGeoTransform( geo ) dst_ds.SetProjection( proj ) dst_ds.GetRasterBand(1).WriteArray(ndvi)