source: trunk/anuga_work/development/boxingday08/combine_good_data.py @ 7924

from anuga.geospatial_data.geospatial_data import *
from anuga.utilities.numerical_tools import ensure_numeric
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 .xyz 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 .xyz 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'
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'
bounding_polygon=ensure_numeric(project.bounding_polygon)
G1 = G1.clip(bounding_polygon)
print 'clipped andaman'
bay_grid_poly=ensure_numeric(project.bay_grid_poly)
G1 = G1.clip_outside(bay_grid_poly)

# Plot the two point sets
data_3s=G1.get_data_points()
data_1s=G0.get_data_points()

from pylab import *
hold(True)
plot(data_3s[:,0],data_3s[:,1],'r.')
#plot(data_1s[:,0],data_1s[:,1],'b.')
#plot(bounding_polygon[:,0],bounding_polygon[:,1],'k-')
show()


G = G0 + G1

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