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?

• xunilk: I tried out the code, slightly modified by me (see my answer), and I got similar errors in my Windows 8 system but, with a correct resulting raster. In my case, I used python 2.7.x and GDAL 1.11.x (no python 3.x).

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):

enter image description here

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.

enter image description here