#!/usr/bin/env python

import gdal
import gdalconst
import sys
import Numeric
import os.path
import osr
import colorsys
import Image
import glob
import string
import os
import pickle

class image_info:
	def __init__(self):
		# color indices used
		self.histogram = [0] * 256
		
		self.y_min = self.y_max = self.x_min = self.x_max = None


def transform_pixel(ct, geomatrix, pixel):
	(pixel_x, pixel_y) = pixel
	X = geomatrix[0] + geomatrix[1] * pixel_x + geomatrix[2] * pixel_y
	Y = geomatrix[3] + geomatrix[4] * pixel_x + geomatrix[5] * pixel_y
	(lon, lat, height) = ct.TransformPoint(X, Y)
	return(lon, lat)

def transform_latlon(ct, geomatrix, pixel):
	(a, b, c, d, e, f) = geomatrix
	(lon, lat) = pixel
	(X, Y, Z) = ct.TransformPoint(lon, lat)
	x = int((c * (Y - d) - f * (X - a)) / (c * e - b * f))
	y = int((b * (Y - d) - e * (X - a)) / (b * f - c * e))
	return(x, y)

if 1:
	masked_sectionals_path = './masked/'
	reference_filename = glob.glob('../raw_data/FAA_sectionals/current/Chicago*South.tif')[0]
	reference_image = gdal.Open(reference_filename, gdalconst.GA_ReadOnly)
	reference_projection = reference_image.GetProjection()
	reference_srs = osr.SpatialReference()
	reference_srs.ImportFromWkt(reference_projection)

	#for sectional_filename in glob.glob('%s*.tif' % (masked_sectionals_path)):
	for sectional_filename in glob.glob('%sHalifax North.tif' % (masked_sectionals_path)):
		print sectional_filename
		sectional_region = sectional_filename[len(masked_sectionals_path):-10]
		sectional_side = sectional_filename[-9:-4]
		print sectional_region, sectional_side

		output_filename = './warped/%s %s.tif' % (sectional_region, sectional_side)
		info_filename = './masked/%s %s.info' % (sectional_region, sectional_side)

		info = pickle.load(open(info_filename, 'r'))

		if 0 and os.path.exists(output_filename):
			print 'already have %s, skipping' % (output_filename)
			continue

		sectional_image = gdal.Open(sectional_filename, gdalconst.GA_ReadOnly)
		sectional_band = sectional_image.GetRasterBand(1)
		sectional_color_table = sectional_band.GetRasterColorTable()
		sectional_geomatrix = sectional_image.GetGeoTransform()
		sectional_srs = osr.SpatialReference()
		sectional_projection = sectional_image.GetProjection()
		sectional_srs.ImportFromWkt(sectional_projection)
		srs_latlon = sectional_srs.CloneGeogCS()
		sectional_ct_to_latlon = osr.CoordinateTransformation(sectional_srs, srs_latlon)
		sectional_ct_to_pixel = osr.CoordinateTransformation(srs_latlon, sectional_srs)

		sectional_to_reference_ct = osr.CoordinateTransformation(sectional_srs, reference_srs)

		print info.x_min, info.x_max
		print info.y_min, info.y_max

		output_x_pixel_size = 64.0
		output_y_pixel_size = 64.0

		print sectional_geomatrix
		
		corners = [
			(info.x_min, info.y_min),
			(info.x_min, info.y_max),
			(info.x_max, info.y_min),
			(info.x_max, info.y_max),
		]

		print corners

		targ_xmin = None
		for corner in corners:
			(x, y) = corner
			X = sectional_geomatrix[0] + sectional_geomatrix[1] * x + sectional_geomatrix[2] * y
			Y = sectional_geomatrix[3] + sectional_geomatrix[4] * x + sectional_geomatrix[5] * y
			(x, y, z) = sectional_to_reference_ct.TransformPoint(X, Y)
			print corner, x,y
			if targ_xmin is None:
				targ_xmin = targ_xmax = x
				targ_ymin = targ_ymax = y
			else:
				targ_xmin = min(targ_xmin, x)
				targ_xmax = max(targ_xmax, x)
				targ_ymin = min(targ_ymin, y)
				targ_ymax = max(targ_ymax, y)

		output_xsize = (targ_xmax - targ_xmin) / output_x_pixel_size + 1
		output_ysize = (targ_ymax - targ_ymin) / output_y_pixel_size + 1
		print output_xsize, output_ysize
		print info.x_max - info.x_min, info.y_max - info.y_min

		output_gt = (targ_xmin, output_x_pixel_size, 0.0, targ_ymax, 0.0, - output_y_pixel_size)

		output_image = gdal.GetDriverByName('GTiff').Create(
			output_filename,
			output_xsize, output_ysize,
			3, gdal.GDT_Byte, # RGB
			[ 'TILED=YES' ]
		)

		output_image.SetProjection( reference_srs.ExportToWkt() )
		output_image.SetGeoTransform( output_gt )
		# output_image.GetRasterBand(1).SetRasterColorTable(sectional_image.GetRasterBand(1).GetRasterColorTable())
		output_image.GetRasterBand(1)
		output_image.GetRasterBand(2)
		output_image.GetRasterBand(3)