pca_multi.py

# hackety hack

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

with rio.open(argv[1]) as src:
  
  count = src.meta['count']
  dtype = src.meta['dtype']
  bidxs = [x+1 for x in range(count)]
  
  # imagine this but faster and clearer:
  pixels = (src.read(b) for b in bidxs)
  pixels = np.dstack([c.ravel() for c in pixels])[0]
  
  pca = decomposition.PCA(n_components=count, whiten=True)
  pca.fit(pixels)
  
  print pca.components_
  
  out = pca.transform(pixels)
  
  print count, range(count)
  print out.shape
  print pixels[0].shape
  
  xyz = np.array([
    out[:,c].reshape(src.meta['height'], src.meta['width'])
    for c in range(count)
  ])
  
  # scale up to fill out the bit depth:
  xyz = np.array([
    ((c - np.amin(c))/(np.amax(c) - np.amin(c)))*np.iinfo(dtype).max
    for c in xyz
  ])
  
  xyz = xyz.astype(dtype)
  
  with rio.open(argv[2],
    'w',
    driver='GTiff',
    width=src.meta['width'],
    height=src.meta['height'],
    count=count,
    dtype=dtype,
    meta=src.meta,
    photometric='rgb') as dst:
    
    for bidx in bidxs:
      dst.write_band(bidx, xyz[bidx-1])