from anuga.geospatial_data.geospatial_data import *
import project
from os.path import exists
from os import sep

def get_zone(filename):
    from anuga.coordinate_transforms.redfearn import redfearn
    fr=open(filename+'.xyz','r')
    line =fr.readline()
    parts = line.rstrip().split()
    zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
    fr.close
    print easting,northing
    return zone

def clip_csv(filename,poly):
    fr=open(filename+'.csv','r')
    fw=open(filename+'_clipped'+'.csv','w')
    fw.write('x,y,elevation\n')
    #extremes = [East,West,North,South] #May only work for positive values
    extremes=[-1.0e10,1.0e10,-1.0e10,1.0e10]
    for points in poly:
        if points[0]>extremes[0]:
            extremes[0]=points[0]
        elif points[0]<extremes[1]:
            extremes[1]=points[0]
        if points[1]>extremes[2]:
            extremes[2]=points[1]
        elif points[1]<extremes[3]:
            extremes[3]=points[1]
    print "polygon_extent= ", extremes  

    start = True
    for line in fr:
        if start:
            #Read in Header
            start=False
        else:
            parts = line.rstrip().split(',')
            if float(parts[0])<extremes[0] and float(parts[0])>extremes[1] and float(parts[1])<extremes[2] and float(parts[1])>extremes[3]:
                fw.write(','.join([parts[0],parts[1],parts[2]])+'\n')
    fr.close
    fw.close
    
def xyz2csv(filename):
    """Convert .xyz file with an arbitrary number of space delimeters to .csv
    format.

    File format
    xyz:
        (No Header)
        longitude, latitude, elevation
    
    csv:
        (Header x,y,elevation)
        easting, northing, elevation
    """
    from anuga.coordinate_transforms.redfearn import redfearn
    
    #Write header
    fr=open(filename+'.xyz','r')
    fw=open(filename+'.csv','w')
    fw.write('x,y,elevation\n')
    for line in fr:
        parts = line.rstrip().split()
        zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
        fw.write(','.join([str(easting),str(northing),parts[2]])+'\n')
        fr.close
        fw.close

filenames=['patong_bay_1s_grid','andaman_3s_grid']

dir='data'
print 'Converting .xya files to .csv format'
for filename in filenames:
    if not exists(dir+sep+filename+'.csv'):
        xyz2csv(dir+sep+filename)
    else:
        print dir+sep+filename+'.csv already exists ensure the .xya changed has not changed'

print 'Clipping .csv files'
filename='andaman_3s_grid'
if not exists(dir+sep+filename+'_clipped'+'.csv'):
    clip_csv(dir+sep+filename,project.bounding_polygon)
else:
    print dir+sep+filename+'_clipped'+'.csv already exists ensure the .csv file has not changed'

print 'Creating geospatial data'
if exists(project.good_combined_dir_name +'.pts'):
    print project.good_combined_dir_name +'.pts'+' already exists ensure the input .csv files have not changed'
else:
    G0 = Geospatial_data(file_name=dir+sep+filenames[0]+'.csv')
    print 'read bay'
    G1=Geospatial_data(file_name=dir+sep+filenames[1]+'_clipped.csv')
    print 'read andaman'
    G1.clip(project.bounding_polygon)
    print 'clipped andaman'
    G1.clip_outside(G0)

    G = G0 + G1

    print 'Creating .pts file'
    G.export_points_file(project.good_combined_dir_name +'.pts')