source: trunk/anuga_core/source/anuga/shallow_water/most2nc.py @ 7844

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

Replaced 'print' statements with log.critical() calls.

File size: 2.6 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
11from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
12import anuga.utilities.log as log
13
14
15##
16# @brief Convert a MOST file to NetCDF format.
17# @param input_file The input file to convert.
18# @param output_file The name of the oputput NetCDF file to create,
19# @param inverted_bathymetry ??
20# @param verbose True if the function is to be verbose.
21def most2nc(input_file, output_file, inverted_bathymetry=False, verbose=True):
22    # variable names
23    long_name = 'LON'
24    lat_name = 'LAT'
25    elev_name = 'ELEVATION'
26
27    # set up bathymetry
28    if inverted_bathymetry:
29        up = -1.
30    else:
31        up = +1.
32
33    # read data from the MOST file
34    in_file = open(input_file,'r')
35
36    if verbose: log.critical('reading header')
37
38    nx_ny_str = in_file.readline()
39    nx_str,ny_str = nx_ny_str.split()
40    nx = int(nx_str)
41    ny = int(ny_str)
42    h1_list=[]
43    for i in range(nx):
44        h1_list.append(float(in_file.readline()))
45
46    h2_list=[]
47    for j in range(ny):
48        h2_list.append(float(in_file.readline()))
49
50    h2_list.reverse()
51
52    if verbose: log.critical('reading depths')
53
54    in_depth_list = in_file.readlines()
55    in_file.close()
56
57    out_depth_list = [[]]
58
59    if verbose: log.critical('processing depths')
60
61    k=1
62    for in_line in in_depth_list:
63        for string in in_line.split():
64            #j = k/nx
65            out_depth_list[(k-1)/nx].append(float(string)*up)
66            if k==nx*ny:
67                break
68            if k-(k/nx)*nx ==0:
69                out_depth_list.append([])
70            k+=1
71
72    in_file.close()
73    out_depth_list.reverse()
74    depth_list = out_depth_list
75
76    # write the NetCDF file
77    if verbose: log.critical('writing results')
78
79    out_file = NetCDFFile(output_file, netcdf_mode_w)
80
81    out_file.createDimension(long_name,nx)
82
83    out_file.createVariable(long_name,'d',(long_name,))
84    out_file.variables[long_name].point_spacing='uneven'
85    out_file.variables[long_name].units='degrees_east'
86    out_file.variables[long_name].assignValue(h1_list)
87
88    out_file.createDimension(lat_name,ny)
89    out_file.createVariable(lat_name,'d',(lat_name,))
90    out_file.variables[lat_name].point_spacing='uneven'
91    out_file.variables[lat_name].units='degrees_north'
92    out_file.variables[lat_name].assignValue(h2_list)
93
94    out_file.createVariable(elev_name,'d',(lat_name,long_name))
95    out_file.variables[elev_name].point_spacing='uneven'
96    out_file.variables[elev_name].units='meters'
97    out_file.variables[elev_name].assignValue(depth_list)
98
99    out_file.close()
Note: See TracBrowser for help on using the repository browser.