#!/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:
	sectionals_path = '/var/www/aviationtoolbox/raw_data/FAA_sectionals/current/'

	horizontal_skip_min = 850
	horizontal_skip_max = 2000

	sectionals_info_file = open('./sectional_corners', 'r')
	
	for sectional_info_line in sectionals_info_file.readlines():
		print sectional_info_line
		(sectional_region, lleft_corner_string, notes) = string.split(sectional_info_line, '|')
		lleft_corner = tuple(map(float, string.split(lleft_corner_string, ' ')))

		#if sectional_region != 'Chicago':
		#	continue

		sectional_version = glob.glob('%s%s ?? North.tif' % (sectionals_path, sectional_region))[0][-12:-10]

		for sectional_side in ['North', 'South']:
			sectional_filename = '%s%s %s %s.tif' % (sectionals_path, sectional_region, sectional_version, sectional_side)
			output_filename = 'masked/%s %s.tif' % (sectional_region, sectional_side)
			info_filename = 'masked/%s %s.info' % (sectional_region, sectional_side)

			info = image_info()	

			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)

			output_image = gdal.GetDriverByName('GTiff').Create(
				output_filename,
				sectional_image.RasterXSize, sectional_image.RasterYSize
			)
			output_image.SetGeoTransform(sectional_geomatrix)
			output_image.SetProjection(sectional_projection)
			output_band = output_image.GetRasterBand(1)
			output_band.SetDescription(sectional_band.GetDescription())
			output_band.SetMetadata(sectional_band.GetMetadata())
			# output_band.SetNoDataValue(sectional_band.GetNoDataValue())
			output_band.SetRasterColorTable(sectional_band.GetRasterColorTable())

			colorish_index = Numeric.ones (256) * 4

			for i in range (min (256, sectional_color_table.GetCount ())):
				#print i,
				(r, g, b, x) = sectional_color_table.GetColorEntry (i)
				(h, l, s) = colorsys.rgb_to_hls(r, g, b)
				#print '%0.3f, %0.3f, %03d' % (h, l, s), 

				if s > -0.04 and l > 128:
					colorish_index[i] = 0
				else:
					colorish_index[i] = 1

			colorish_index[255] = -1

			src_ds = sectional_image
			src_band = sectional_band


			# Find the top line.
			inmap_vertical_count = 0
			for iY in range(0, src_ds.RasterYSize, 1):
				src_data = src_band.ReadAsArray(
					horizontal_skip_min, iY,
					src_ds.RasterXSize - horizontal_skip_min, 1
				)
				src_data_indexed = Numeric.take(colorish_index, src_data)[0]
						
				colorval = sum(src_data_indexed) * 100 / len(src_data_indexed)

				# print colorval

				if colorval >= 80: 
					inmap_vertical_count += 1
					if inmap_vertical_count > 3:
						y_min = iY - 3
						break
				else:
					inmap_vertical_count = 0

	
			# Write the top bands.
			output_band.WriteArray(Numeric.ones((y_min, src_ds.RasterXSize)) * 255, 0, 0)

			# Run through the rest.
			
			# Find a reasonable X value.
			# (prev_x, y) = transform_latlon(sectional_ct_to_pixel, sectional_geomatrix, lleft_corner)

			gdal.TermProgress(0.0)
			for iY in range(y_min, src_ds.RasterYSize):
				# Get horizontal line of pixels from sectional, skipping left side.
				output_data = sectional_band.ReadAsArray(0, iY, src_ds.RasterXSize, 1)
			
				# Initialize with blanks on left (for now).
				output_data[0][:horizontal_skip_min] = 255

				for iX in range(horizontal_skip_min, src_ds.RasterXSize):
					(lon, lat) = transform_pixel(sectional_ct_to_latlon, sectional_geomatrix, (iX, iY))
					if lon < lleft_corner[0] or lat < lleft_corner[1]:
						output_data[0][iX] = 255
					else:
						info.histogram[output_data[0][iX]] += 1
						if not info.x_min or iX < info.x_min:
							info.x_min = iX
						if not info.x_max or iX > info.x_max:
							info.x_max = iX
						if not info.y_min or iY < info.y_min:
							info.y_min = iY
						if not info.y_max or iY > info.y_max:
							info.y_max = iY


				
				output_band.WriteArray(output_data, 0, iY)

				gdal.TermProgress( (iY+1.0) / src_ds.RasterYSize )

	
			pickle.dump(info, open(info_filename, 'w'))