#!/usr/bin/env python

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


if 1:
	sectionals_path = '/var/www/aviationtoolbox/raw_data/FAA_sectionals/current/'
	sectional_filenames = glob.glob('%s* North.tif' % (sectionals_path))
	sectional_filenames.sort()

	horizontal_skip_min = 950
	horizontal_skip_max = 2000

	
	for base_sectional_filename in sectional_filenames:
		sectional_region = base_sectional_filename[len(sectionals_path):-13]
		sectional_version = base_sectional_filename[-12:-10]

	

		for side in ['North', 'South']:
			sectional_filename = '%s%s %s %s.tif' % (sectionals_path, sectional_region, sectional_version, side)
	#		#sectional_filename = 'Chicago_excerpt North.tif'
	#	for sectional_filename in sys.argv[1:]:

			sectional_image = gdal.Open(sectional_filename, gdalconst.GA_ReadOnly)
			sectional_band = sectional_image.GetRasterBand(1)
			sectional_color_table = 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] = 2

			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)

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

		

			# Find left edge (line of longitude).
			geomatrix = src_ds.GetGeoTransform()
			srs = osr.SpatialReference()
			srs.ImportFromWkt(src_ds.GetProjection())

			srsLatLong = srs.CloneGeogCS()
			ct = osr.CoordinateTransformation(srs, srsLatLong)

			left_longitudes = []
			for iY in range(top_line, src_ds.RasterYSize, 100):
				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]
						
				inmap_horizontal = 0
				check_area_size = 30
				threshold_high = check_area_size * 0.8
				threshold_low = check_area_size * 0.2
				horizontal_beyond_white = 0
				threshold_white_band = check_area_size * 0.6

				for iXo in range(horizontal_skip_max - horizontal_skip_min + check_area_size):
					check_area = src_data_indexed[iXo:iXo + check_area_size]
					check_area_sum = sum(check_area)

					# print iY, iXo, check_area_sum, horizontal_beyond_white

					if not horizontal_beyond_white:
						if check_area_sum <= threshold_white_band:
							horizontal_beyond_white = 1
						continue	

					if not inmap_horizontal and check_area_sum > threshold_high:
						inmap_horizontal = iXo
					elif inmap_horizontal and check_area_sum < threshold_low:
						inmap_horizontal = 0

				horizontal_start = inmap_horizontal + horizontal_skip_min + check_area_size - int(threshold_high) 

				X = geomatrix[0] + geomatrix[1] * horizontal_start + geomatrix[2] * iY
				Y = geomatrix[3] + geomatrix[4] * horizontal_start + geomatrix[5] * iY
				(long, lat, height) = ct.TransformPoint(X, Y)
				# print iY, horizontal_start, long, lat
				left_longitudes.append(long)

				dst_data = src_data

			left_longitudes.sort()
			left_longitudes_outliers_count = int(len(left_longitudes) * .2)
			left_longitudes_middle = left_longitudes[left_longitudes_outliers_count:-left_longitudes_outliers_count]
			left_edge_lon = round(sum(left_longitudes_middle) / len(left_longitudes_middle) * 10) / 10.0


			if side == 'South':
				# Find bottom edge (line of latitude) if there is one.
		
				x = src_ds.RasterXSize / 2
				y = src_ds.RasterYSize / 2
				start_lat = lat
		
				vertical_space = 0
				for iY in range(int(src_ds.RasterYSize * 0.8), src_ds.RasterYSize):
					for iX in range(horizontal_skip_min, src_ds.RasterXSize):
						X = geomatrix[0] + geomatrix[1] * iX + geomatrix[2] * iY
						Y = geomatrix[3] + geomatrix[4] * iX + geomatrix[5] * iY
						(long, lat, height) = ct.TransformPoint(X, Y)
						if long >= left_edge_lon:
							break
		
					src_data = src_band.ReadAsArray(
						iX, iY,
						check_area_size, 1
					)
					src_data_indexed = Numeric.take(colorish_index, src_data)[0]
					check_area_score = sum(src_data_indexed) / float(check_area_size)
					# print iY, iX + horizontal_skip_min, check_area_score
					if check_area_score < 0.8:
						if vertical_space == 0:
							vertical_space_start = (iX, iY)
						vertical_space += 1
						if vertical_space > 60:
							break
					elif check_area_score > 0.9:
						vertical_space = 0
						vertical_space_start = None
		
					continue

				if vertical_space_start:
					(iX, iY) = vertical_space_start
					X = geomatrix[0] + geomatrix[1] * iX + geomatrix[2] * iY
					Y = geomatrix[3] + geomatrix[4] * iX + geomatrix[5] * iY
					(long, lat, height) = ct.TransformPoint(X, Y)
					bottom_edge_lat = round(lat * 10) / 10.0
				else:
					bottom_edge_lat = None
			else:
				bottom_edge_lat = None

			if side == 'North':
				print sectional_region, left_edge_lon,
				left_edge_lon_north = left_edge_lon

			else:
				if not left_edge_lon == left_edge_lon_north:
					print '(%f)' % (left_edge_lon),
				if bottom_edge_lat is None:
					print '-',
				else:
					print bottom_edge_lat,

				print

		#sys.exit(0)