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

Last change on this file since 5597 was 3564, checked in by ole, 18 years ago

Moved more supporting files into shallow_water.
Moved old least_squares to obsolete_code

File size: 2.5 KB
RevLine 
[3563]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 MOSTs output.
3
4 $Author: Peter Row
5 
6"""
7
8
[1123]9def most2nc(input_file=None,output_file=None,inverted_bathymetry = False,\
10            verbose = True):
[1122]11    #input_file = 'small.txt'
12    #output_file = 'small_e.nc'
[1120]13
[1122]14    long_name = 'LON'
15    lat_name = 'LAT'
16    if inverted_bathymetry:
[1124]17        up = -1.
[1123]18    else:
[1124]19        up = +1.
20
[1122]21    from Scientific.IO.NetCDF import NetCDFFile
22    import sys
[1124]23
[1122]24    try:
[1124]25        if input_file is None:
26            input_file = sys.argv[1]
27        if output_file is None:
28            output_file = sys.argv[2]
[1122]29    except:
[1124]30        raise 'usage is: most2nc input output'
31
[1122]32    in_file = open(input_file,'r')
[1123]33    if verbose: print 'reading header'
[1122]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()))
[1120]41
[1122]42    h2_list=[]
43    for j in range(ny):
44        h2_list.append(float(in_file.readline()))
[1120]45
[1122]46    h2_list.reverse()
[1120]47
[1123]48    if verbose: print 'reading depths'
[1124]49
[1122]50    in_depth_list = in_file.readlines()
51    in_file.close()
[1120]52
[1122]53    out_depth_list = [[]]
[1120]54
[1123]55    if verbose: print 'processing depths'
[1124]56
[1122]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
[1120]68
[1122]69    in_file.close()
70    out_depth_list.reverse()
71    depth_list = out_depth_list
[1120]72
[1123]73    if verbose: print 'writing results'
[1124]74
[1122]75    out_file = NetCDFFile(output_file,'w')
[1120]76
[1122]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)
[1120]82
[1122]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)
[1120]88
[1122]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)
[1120]93
[1122]94    out_file.close()
Note: See TracBrowser for help on using the repository browser.