1 | """This has to do with creating elevation data files for use with ferret2sww. |
---|
2 | It reads a bathymetry ascii file and creates a NetCDF (nc) file similar to MOSTs output. |
---|
3 | |
---|
4 | $Author: Peter Row |
---|
5 | |
---|
6 | """ |
---|
7 | |
---|
8 | |
---|
9 | def most2nc(input_file=None,output_file=None,inverted_bathymetry = False,\ |
---|
10 | verbose = True): |
---|
11 | #input_file = 'small.txt' |
---|
12 | #output_file = 'small_e.nc' |
---|
13 | |
---|
14 | long_name = 'LON' |
---|
15 | lat_name = 'LAT' |
---|
16 | if inverted_bathymetry: |
---|
17 | up = -1. |
---|
18 | else: |
---|
19 | up = +1. |
---|
20 | |
---|
21 | from Scientific.IO.NetCDF import NetCDFFile |
---|
22 | import sys |
---|
23 | |
---|
24 | try: |
---|
25 | if input_file is None: |
---|
26 | input_file = sys.argv[1] |
---|
27 | if output_file is None: |
---|
28 | output_file = sys.argv[2] |
---|
29 | except: |
---|
30 | raise 'usage is: most2nc input output' |
---|
31 | |
---|
32 | in_file = open(input_file,'r') |
---|
33 | if verbose: print 'reading header' |
---|
34 | nx_ny_str = in_file.readline() |
---|
35 | nx_str,ny_str = nx_ny_str.split() |
---|
36 | nx = int(nx_str) |
---|
37 | ny = int(ny_str) |
---|
38 | h1_list=[] |
---|
39 | for i in range(nx): |
---|
40 | h1_list.append(float(in_file.readline())) |
---|
41 | |
---|
42 | h2_list=[] |
---|
43 | for j in range(ny): |
---|
44 | h2_list.append(float(in_file.readline())) |
---|
45 | |
---|
46 | h2_list.reverse() |
---|
47 | |
---|
48 | if verbose: print 'reading depths' |
---|
49 | |
---|
50 | in_depth_list = in_file.readlines() |
---|
51 | in_file.close() |
---|
52 | |
---|
53 | out_depth_list = [[]] |
---|
54 | |
---|
55 | if verbose: print 'processing depths' |
---|
56 | |
---|
57 | k=1 |
---|
58 | for in_line in in_depth_list: |
---|
59 | for string in in_line.split(): |
---|
60 | #j = k/nx |
---|
61 | out_depth_list[(k-1)/nx].append(float(string)*up) |
---|
62 | #print k,len(out_depth_list),(k-1)/nx,out_depth_list[(k-1)/nx][-1],len(out_depth_list[(k-1)/nx]) |
---|
63 | if k==nx*ny: |
---|
64 | break |
---|
65 | if k-(k/nx)*nx ==0: |
---|
66 | out_depth_list.append([]) |
---|
67 | k+=1 |
---|
68 | |
---|
69 | in_file.close() |
---|
70 | out_depth_list.reverse() |
---|
71 | depth_list = out_depth_list |
---|
72 | |
---|
73 | if verbose: print 'writing results' |
---|
74 | |
---|
75 | out_file = NetCDFFile(output_file,'w') |
---|
76 | |
---|
77 | out_file.createDimension(long_name,nx) |
---|
78 | out_file.createVariable(long_name,'d',(long_name,)) |
---|
79 | out_file.variables[long_name].point_spacing='uneven' |
---|
80 | out_file.variables[long_name].units='degrees_east' |
---|
81 | out_file.variables[long_name].assignValue(h1_list) |
---|
82 | |
---|
83 | out_file.createDimension(lat_name,ny) |
---|
84 | out_file.createVariable(lat_name,'d',(lat_name,)) |
---|
85 | out_file.variables[lat_name].point_spacing='uneven' |
---|
86 | out_file.variables[lat_name].units='degrees_north' |
---|
87 | out_file.variables[lat_name].assignValue(h2_list) |
---|
88 | |
---|
89 | out_file.createVariable('ELEVATION','d',(lat_name,long_name)) |
---|
90 | out_file.variables['ELEVATION'].point_spacing='uneven' |
---|
91 | out_file.variables['ELEVATION'].units='meters' |
---|
92 | out_file.variables['ELEVATION'].assignValue(depth_list) |
---|
93 | |
---|
94 | out_file.close() |
---|