ndvi.py

#!/usr/bin/env python
# Very simple rasterio NDVI (or general normalized difference ratio) demo.
# python ndvi.py $landsat_scene/*B{4,5}.TIF ndvi.tif

import rasterio as rio
from sys import argv
import numpy as np

math_type = np.float32 # type for internal calculations

red_file, nir_file = (rio.open(img) for img in argv[1:3])
red, nir = (f.read(1).astype(math_type) for f in (red_file, nir_file))

dst_type = red_file.dtypes[0]
largest = np.iinfo(dst_type).max # for scaling

with rio.open(argv[3], 'w', **red_file.meta) as dst:
  np.seterr(divide='ignore')
  # NDVI = (nir - red)/(nir + red), which is in -1..1,
  # which we output as 0..(largest representable value)
  scale_ndvi = (nir * largest)/(red + nir)
  dst.write_band(1, scale_ndvi.astype(dst_type))