1 | from anuga.geospatial_data.geospatial_data import * |
---|
2 | import project |
---|
3 | from os.path import exists |
---|
4 | from os import sep |
---|
5 | |
---|
6 | def 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 | |
---|
16 | def 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 | |
---|
45 | def 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 | |
---|
71 | filenames=['patong_bay_1s_grid','andaman_3s_grid'] |
---|
72 | |
---|
73 | dir='data' |
---|
74 | print 'Converting .xya files to .csv format' |
---|
75 | for 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 | |
---|
81 | print 'Clipping .csv files' |
---|
82 | filename='andaman_3s_grid' |
---|
83 | if not exists(dir+sep+filename+'_clipped'+'.csv'): |
---|
84 | clip_csv(dir+sep+filename,project.bounding_polygon) |
---|
85 | |
---|
86 | print 'Creating geospatial data' |
---|
87 | if exists(project.good_combined_dir_name +'.pts'): |
---|
88 | print project.good_combined_dir_name +'.pts'+' already exists ensure the input .csv files have not changed' |
---|
89 | else: |
---|
90 | G0 = Geospatial_data(file_name=dir+sep+filenames[0]+'.csv') |
---|
91 | G1=Geospatial_data(file_name=dir+sep+filenames[1]+'_clipped.csv') |
---|
92 | G1.clip(project.bounding_polygon) |
---|
93 | G1.clip_outside(G0) |
---|
94 | |
---|
95 | G = G0 + G1 |
---|
96 | |
---|
97 | print 'Creating .pts file' |
---|
98 | G.export_points_file(project.good_combined_dir_name +'.pts') |
---|