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[1], shape[0], 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform( geo )
dst_ds.SetProjection( proj )
dst_ds.GetRasterBand(1).WriteArray(ndvi)
Best Answer
Without having gone through your code just a suggestion:
NDVI is an index and results are between 0 and 1.
you will probably work with 8-bit Tiff which stores values between 0 and 255
so if you multiply your results by 100 it should work