source: anuga_work/development/boxingday08/combine_good_data.py @ 5393

Last change on this file since 5393 was 5393, checked in by jakeman, 16 years ago

updating scripts for boxingday validation

File size: 3.2 KB
Line 
1from anuga.geospatial_data.geospatial_data import *
2import project
3from os.path import exists
4from os import sep
5
6def get_zone(filename):
7    from anuga.coordinate_transforms.redfearn import redfearn
8    fr=open(filename+'.xyz','r')
9    line =fr.readline()
10    parts = line.rstrip().split()
11    zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
12    fr.close
13    print easting,northing
14    return zone
15
16def clip_csv(filename,poly):
17    fr=open(filename+'.csv','r')
18    fw=open(filename+'_clipped'+'.csv','w')
19    fw.write('x,y,elevation\n')
20    #extremes = [East,West,North,South] #May only work for positive values
21    extremes=[-1.0e10,1.0e10,-1.0e10,1.0e10]
22    for points in poly:
23        if points[0]>extremes[0]:
24            extremes[0]=points[0]
25        elif points[0]<extremes[1]:
26            extremes[1]=points[0]
27        if points[1]>extremes[2]:
28            extremes[2]=points[1]
29        elif points[1]<extremes[3]:
30            extremes[3]=points[1]
31    print "polygon_extent= ", extremes 
32
33    start = True
34    for line in fr:
35        if start:
36            #Read in Header
37            start=False
38        else:
39            parts = line.rstrip().split(',')
40            if float(parts[0])<extremes[0] and float(parts[0])>extremes[1] and float(parts[1])<extremes[2] and float(parts[1])>extremes[3]:
41                fw.write(','.join([parts[0],parts[1],parts[2]])+'\n')
42    fr.close
43    fw.close
44   
45def xyz2csv(filename):
46    """Convert .xyz file with an arbitrary number of space delimeters to .csv
47    format.
48
49    File format
50    xyz:
51        (No Header)
52        longitude, latitude, elevation
53   
54    csv:
55        (Header x,y,elevation)
56        easting, northing, elevation
57    """
58    from anuga.coordinate_transforms.redfearn import redfearn
59   
60    #Write header
61    fr=open(filename+'.xyz','r')
62    fw=open(filename+'.csv','w')
63    fw.write('x,y,elevation\n')
64    for line in fr:
65        parts = line.rstrip().split()
66        zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
67        fw.write(','.join([str(easting),str(northing),parts[2]])+'\n')
68        fr.close
69        fw.close
70
71filenames=['patong_bay_1s_grid','andaman_3s_grid']
72
73dir='data'
74print 'Converting .xya files to .csv format'
75for filename in filenames:
76    if not exists(dir+sep+filename+'.csv'):
77        xyz2csv(dir+sep+filename)
78    else:
79        print dir+sep+filename+'.csv already exists ensure the .xya changed has not changed'
80
81print 'Clipping .csv files'
82filename='andaman_3s_grid'
83if not exists(dir+sep+filename+'_clipped'+'.csv'):
84    clip_csv(dir+sep+filename,project.bounding_polygon)
85else:
86    print dir+sep+filename+'_clipped'+'.csv already exists ensure the .csv file has not changed'
87
88print 'Creating geospatial data'
89if exists(project.good_combined_dir_name +'.pts'):
90    print project.good_combined_dir_name +'.pts'+' already exists ensure the input .csv files have not changed'
91else:
92    G0 = Geospatial_data(file_name=dir+sep+filenames[0]+'.csv')
93    print 'read bay'
94    G1=Geospatial_data(file_name=dir+sep+filenames[1]+'_clipped.csv')
95    print 'read andaman'
96    G1.clip(project.bounding_polygon)
97    print 'clipped andaman'
98    G1.clip_outside(G0)
99
100    G = G0 + G1
101
102    print 'Creating .pts file'
103    G.export_points_file(project.good_combined_dir_name +'.pts')
Note: See TracBrowser for help on using the repository browser.