#!/usr/bin/env python

import gdal
import osr
import glob
import string

if 0:
	lon = -86.9368989
	lat = 40.4123169
	area = 'dvd9/area05'

	ds = gdal.Open('/var/www/aviationtoolbox/raw_data/USGS/NED/current/%s/neddata' % (area))

	size = (ds.RasterXSize, ds.RasterYSize)
	projection = ds.GetProjection()
	geotransform = ds.GetGeoTransform()

	print 'size=(%d, %d)' % size

	srs = osr.SpatialReference()
	srs.ImportFromWkt(projection)

	srs_lonlat = srs.CloneGeogCS()

	ct_to_lonlat = osr.CoordinateTransformation(srs, srs_lonlat)
	ct_from_lonlat = osr.CoordinateTransformation(srs_lonlat, srs)

	(X, Y, Z) = ct_from_lonlat.TransformPoint(lon, lat)

	print X, Y, Z

	(a, b, c, d, e, f) = geotransform

	x = int((c * (Y - d) - f * (X - a)) / (c * e - b * f))
	y = int((b * (Y - d) - e * (X - a)) / (b * f - c * e))

	print x, y

	elevation_band = ds.GetRasterBand(1)

	elevations = elevation_band.ReadAsArray(x, y, 1, 1)

	# Print as feet (not meters).
	print elevations * 3.2808399

if 1:
	#(lon, lat) = (-86.9368989, 40.4123169) # LAF (606)
	#(lon, lat) = (-104.6731775, 39.8616564) # DEN (5431)
	#(lon, lat) = (-107.9084800, 37.9537586) # TEX (9078)
	area_paths = glob.glob('/var/www/aviationtoolbox/raw_data/USGS/NED/current/dvd*/area*')

airports_file = open('../raw_data/FAA_ATA-100/APT.del', 'r')

for line in airports_file.readlines():
	fields = string.split(line, '|')
	if not fields[0] == 'APT':
		continue
	if not fields[2] == 'AIRPORT':
		continue

	landing_facility = fields[2]
	identifier = fields[3]
	county = fields[9]
	state = fields[10]

	if state in ['AK', 'HI']:
		continue

	lat = float(fields[24][:-1]) / 3600.0
	lon = float(fields[26][:-1]) / -3600.0
	rp_det = fields[27]  # No reference point is surveyed.
	elev = float(fields[28])
	elev_det = fields[29]

	if not elev_det == 'S':
		continue

	#if not identifier == 'LAF':
	#	continue

	found = 0
	for area_path in area_paths:
		#print area_path
		ds = gdal.Open(area_path + '/neddata/')
		size = (ds.RasterXSize, ds.RasterYSize)
		projection = ds.GetProjection()
		geotransform = ds.GetGeoTransform()

		srs = osr.SpatialReference()
		srs.ImportFromWkt(projection)


		srs_lonlat = srs.CloneGeogCS()

		ct_to_lonlat = osr.CoordinateTransformation(srs, srs_lonlat)
		#(X, Y, Z) = ct_from_lonlat.TransformPoint(lon, lat)

		corners_xy = ( (0, 0), size )
		corners_XY = []
		for (x, y) in corners_xy:
			X = geotransform[0] + geotransform[1] * x + geotransform[2] * y
			Y = geotransform[3] + geotransform[4] * x + geotransform[5] * y

			corners_XY.append(ct_to_lonlat.TransformPoint(X, Y)[:2])
 
		#print corners_XY

		if (
			corners_XY[0][0] < lon and corners_XY[1][0] > lon
			and
			corners_XY[0][1] > lat and corners_XY[1][1] < lat
		):
			found = 1
			break

	if not found:
		print '%s %f %f missing coverage' % (identifier, lon, lat)
		continue

	#print area_path
	
	ct_from_lonlat = osr.CoordinateTransformation(srs_lonlat, srs)

	(X, Y, Z) = ct_from_lonlat.TransformPoint(lon, lat)

	#print X, Y, Z

	(a, b, c, d, e, f) = geotransform

	x = int((c * (Y - d) - f * (X - a)) / (c * e - b * f))
	y = int((b * (Y - d) - e * (X - a)) / (b * f - c * e))

	elevation_band = ds.GetRasterBand(1)

	elevation = elevation_band.ReadAsArray(x, y, 1, 1)[0][0] * 3.2808399

	if elevation < -1000:
		print '%s %f %f missing information' % (identifier, lon, lat)
		continue

	print identifier, lon, lat, elev, elevation, elev - elevation #, projection[8:13]