do not use this

simple-l8-pan-draft.py

# PRE-ALPHA!

import asyncio
import rasterio as rio
import numpy as np
from rasterio.warp import reproject, RESAMPLING

'''
./pan.py print/LC81100362014076LGN00/*B{4,3,2,8}.TIF out.tif

This is broken. If you're looking for a simple Landsat 8 pansharpener,
here's an imperfect but working one: https://gist.github.com/celoyd/2e7beed82951d22b9b90

- fix horizontal null bands appearing in larger images -- rounding/off-by-one?
- fix everything to do with windowing
- explain that the weird windows are because
  Landsat 8 bands come as Nx1
- parallelize
- ^ to that end, have the pansharpen function do the writing?
- fix overflow at bright values

'''

yblocksize=2**14
xblocksize=256

matht = np.float32

def pansharpen(rgb, color_meta, p, pan_meta):
	np.seterr(invalid='ignore', divide='ignore')
	for i, c in enumerate(rgb):
		x = np.empty(p.shape, matht)
		reproject(
			c, x,
			src_transform=color_meta['transform'],
			src_crs=color_meta['crs'],
			dst_transform=pan_meta['transform'],
			dst_crs=pan_meta['crs'],
			resampling=RESAMPLING.bilinear)
		rgb[i] = x
	ratio = (p / ((rgb[0] + rgb[1] + rgb[2])/3)) # this overflows or something
	return (x*ratio for x in rgb)

def half_window(window):
	# same spatial window for perfectly a aligned image at half linear res
	((fr, tr), (fc, tc)) = window
	return ((int(fr/2), int(tr/2)), (int(fc/2), int(tc/2)))

def make_windows(w, h):
 for x in range(int(w/yblocksize+1)): # well, this got away from me fast
 	for y in range(int(h/xblocksize+1)):
 		yield (
			(y*xblocksize, min(h, (y*xblocksize)+xblocksize)),
			(x*yblocksize, min(w, (x*yblocksize)+yblocksize)))

def main(red, green, blue, pan, out):
	r, g, b, p = (rio.open(f) for f in [red, green, blue, pan])
	rgb = (r, g, b)
	color_meta, pan_meta = b.meta, p.meta
	w, h = pan_meta['width'], pan_meta['height']
	print(w, h)
	pan_meta.update(compress='lzw', photometric='rgb', count=3)
	with rio.open(out, 'w', **pan_meta) as dest:
		n = 1
		ws = list(make_windows(w, h)) # Why did I make this a generator if I just ...
		for w in ws:
			print("%s/%s" % (n, len(ws))) # oh, fine.
			n += 1
			print(w, half_window(w))
			pan_data = p.read_band(1, window=w).astype(matht)
			rgb_data = list(x.read_band(1, window=half_window(w)).astype(matht) for x in rgb)
			for channel, data in enumerate(pansharpen(rgb_data, color_meta, pan_data, pan_meta)):
				dest.write_band(channel+1, data.astype(np.uint16), window=w)
	print('that was cool and fun')

if __name__ == '__main__':
	import sys
	main(*sys.argv[1:])