l8tc.py

import rasterio as rio
import numpy as np
import os
import sys

# Tasseled Cap with Baig et al. 2014 coefficients.
# Use: l8tc.py bundle/*B{2..7}.TIF tass.tif

coefficients = [
  [0.3029, 0.2786, 0.4733 ,0.5599, 0.508, 0.1872],
  [-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608],
  [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559],
  [-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773],
  [-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085],
  [ 0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252]
]

with rio.open(sys.argv[1]) as for_meta:
  meta = for_meta.meta

meta.update({
  'count': 6,
  'transform': None
})

paths = sys.argv[1:-1]

stack = np.array(
  [rio.open(path).read(1).astype(np.float32) for path in paths]
)

with rio.open(sys.argv[-1], 'w', **meta) as dst:
  for component in range(1, 7):
    image = np.zeros_like(stack[0])
    for bidx in range(1, 7):
      image += stack[bidx-1] * coefficients[component-1][bidx-1]
      
    image = np.clip((image/2) + 2**15, 0, (2**16)-1).astype(np.uint16)
    
    dst.write(image, component)