[GIS] Unable to create NDVI image using Python gdal package

Question

I have Python 3 and all GDAL bindings installed on my Windows machine, so that when I run this:

import numpy as np
from sys import argv
from osgeo import gdal, gdalconst

I get no error messages. What I want to achieve now is to create NDVI image from two Landsat 8 images. To implement this, I'm using this nice looking script – https://gist.github.com/celoyd/7443646 . I really like how simple it is, however, when I run it like so:

C:\LC81750202017091LGN00>python ndvi.py LC81750202017091LGN00_B4.TIF LC81750202017091LGN00_B5.TIF output-ndvi.TIF
ERROR 1: Deleting output-ndvi.TIF failed:
Permission denied
ERROR 1: TIFFResetField:output-ndvi.TIF: Could not find tag 273.
Traceback (most recent call last):
  File "ndvi.py", line 31, in <module>
    ndvi = (n - r)/(n + r)
MemoryError

I tried to debug it and see that the first error message

ERROR 1: Deleting output-ndvi.TIF failed:

comes from these lines of code (19-23):

output = geotiff.Create(
   argv[3], 
   red.RasterXSize, red.RasterYSize, 
   1, 
   gdal.GDT_UInt16)

How can I fix it?

Answer 1

I complemented the code with a few lines, to set spatial reference and same geotransform parameters of original rasters, and it works perfectly in a Linux system (python 2.7.x and GDAL 1.11.x). In a Windows 8.x system (python 2.7.x and GDAL 1.11.x), it too worked but it printed these console errors (first one is similar to your error but I got a resulting raster):

with next command:

python ndvi_github.py b3.ND.tif b4.ND.tif output-ndvi.tif

Apparently, resulting raster hasn't any problem.

Complete used code was:

#!/usr/bin/env/python

# ndvi.py red.tif nir.tif output-ndvi.tif
# Calculate NDVI (see Wikipedia). Assumes atmospheric correction.
# (Although I use it without all the time for quick experiments.)

import numpy as np
from sys import argv
from osgeo import gdal, gdalconst, osr

#type for internal calculations
t = np.float32

red, nir = map(gdal.Open, argv[1:3])

print argv[1:3]

transform = red.GetGeoTransform()

geotiff = gdal.GetDriverByName('GTiff')
output = geotiff.CreateCopy(argv[3], red, 0)

output = geotiff.Create(argv[3], 
                       red.RasterXSize, 
                       red.RasterYSize, 
                       1,
                       gdal.GDT_Float32)

# Ugly syntax, but fast:
r = red.GetRasterBand(1).ReadAsArray(0,0,red.RasterXSize,red.RasterYSize)
n = nir.GetRasterBand(1).ReadAsArray(0,0,nir.RasterXSize,nir.RasterYSize)

# Convert the 16-bit Landsat 8 values to floats for the division operation:
r = r.astype(t)
n = n.astype(t)

# Tell numpy not to complain about division by 0:
np.seterr(invalid='ignore')

# Here's the meat of this whole thing, the actual NDVI formula:
ndvi = (n - r)/(n + r)

# The ndvi value is in the range -1..1, but we want it to be displayable, so:
# Make the value positive and scale it back up to the 16-bit range:
ndvi = (ndvi + 1) * (2**15 - 1)

# And do the type conversion back:
ndvi = ndvi.astype(np.uint16)

output.GetRasterBand(1).WriteArray(ndvi)
output.SetGeoTransform(transform)

wkt = red.GetProjection()

#setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
output.SetProjection( srs.ExportToWkt() )

red = None
nir = None
output = None

and next image was resulting raster, obtained in Windows 8 system, loaded in QGIS.