4

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)
  • You may find you answer from this link. mucing.github.io/NDVI-for-Large-Raster-Python – Lexi_GIS-RS May 3 '16 at 16:31
  • It looks potentially helpful but links can go away sometimes. Perhaps you could provide some commentary about what the link provides? – jbchurchill May 3 '16 at 17:02
  • 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 – jbchurchill May 3 '16 at 17:03
2

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 ) 
  • The OP is writing to Float32 not Byte tiffs, so there's no need to multiply. – user2856 Apr 15 '14 at 23:48
  • 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. – Doug Apr 16 '14 at 2:46
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.

  • +1 This has caught more than a few people out, myself included... Should probably be added to the GDAL python gotchas wiki page. – user2856 Apr 16 '14 at 7:43
  • 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. – Doug Apr 16 '14 at 13:41
  • 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. – sgillies Apr 16 '14 at 16:09
1

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 )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.