using Python to calculate NDVI, output as a Geotiff

Question

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)

• Lexi_GIS-RS: You may find you answer from this link. mucing.github.io/NDVI-for-Large-Raster-Python
• jbchurchill: It looks potentially helpful but links can go away sometimes. Perhaps you could provide some commentary about what the link provides?
• jbchurchill: While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review

Answer 1

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

ndvi = np.where ( check, (nir - red ) / ( nir + red ) * 100, -999 ) 

• user2856: The OP is writing to Float32 not Byte tiffs, so there's no need to multiply.
• Doug: thanks for the suggestion, but the problem isn't in the ND calculation. If I print out the NDVI array, it has values as expected. The problem seems to be in the WriteArray command, since it's the TIFF file that has zeroes in it.

Answer 2

It doesn't look like your code properly saves/closes the dataset. To do this, add this to the end:

dst_ds = None # save, close

Also, although it looks like you want to use -999 for NODATA, this needs to be set to the resulting band.

If you want to learn more about raster processing with Python, check out rasterio.

• user2856: +1 This has caught more than a few people out, myself included... Should probably be added to the GDAL python gotchas wiki page.
• Doug: Mike T -- that did it! Thanks! Also, thanks for the pointers to rasterio and the gotcha wiki. Both look like they will be a big help.
• sgillies: Doug, like Mike said, rasterio's datasets are designed to be just like familiar Python file objects. There's an open() method and 'r' or 'w' modes and to close them and save, you call a close() method. No gotchas.

Answer 3

Try:

 red = np.array(red, dtype = float)
 nir = np.array(nir, dtype = float)

Besides I think you mean:

check = np.logical_or ( red > 0, nir > 0 )