source: anuga_core/source/anuga/shallow_water/most2nc.py @ 6058

Last change on this file since 6058 was 6058, checked in by rwilson, 15 years ago

Removed old 'run from commandline' code, added comments, pep8'd it.

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