source: inundation/pyvolution/interpolate_sww.py @ 1884

Last change on this file since 1884 was 1814, checked in by duncan, 19 years ago

added more error checking to files loaded

File size: 6.8 KB
Line 
1""" Interpolation of a sww file.
2Used to interpolate height information from sww files.
3
4When using as a stand alone application,
5Input;
6 - The sww file with stage and bed elevation information
7 - An xya file specifying the points (x,y) where the height values
8   (w.r.t. time) are needed.
9
10Ouput;
11An xya file with x,y and height values w.r.t. time.
12
13The first row of the output xya file has the time value for each height
14Each following row has the format x,y,[height,]
15
16NOTE: stage = bed elevation + height
17
18   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
19   Geoscience Australia, 2004.   
20"""
21
22##FIXME (DSG-DSG)  no sww file? give a better error message.
23
24from Numeric import transpose
25from least_squares import Interpolation
26
27DEFAULT_QUANTITY = "depth"
28
29def  interpolate_sww2xya(sww_file, quantity_name,point_file_in,
30                         point_file_out, display_errors = True):
31    """
32    This function catches exceptions.
33    """
34    try:
35        interp = Interpolate_sww(sww_file, quantity_name)
36        interp.interpolate_xya(point_file_in)
37        interp.write_depth_xya(point_file_out)
38    except IOError,e: #need to convert key error to ioerror
39        if display_errors:
40            print "Could not load bad file. ", e
41        import sys; sys.exit()
42
43        # FIXME (DSG-DSG): how are bad quantities caught?
44        #try:
45        #    interp = Interpolate_sww(sww_file, quantity_name)
46        #except KeyError:
47        #    print "Error: Unknown quantity"
48        #    sys.exit(1)
49
50        interp.interpolate_xya(point_file_in)
51        interp.write_depth_xya(point_file_out)
52
53
54
55class Interpolate_sww:
56    def __init__(self, file_name, quantity_name):
57
58        #It's bad to have the quantity_name passed in here.
59        # But it works with how the program is currently used.
60        # Refactor when necessary. - DSG
61       
62        x, y, volumes, time, quantity = self.read_sww(file_name, quantity_name)
63        vertex_coordinates = transpose([x,y])
64       
65        if False:
66            print "****************************"
67            print "x "
68            print x
69            print "****************************"
70            print "Y "
71            print y
72            print "****************************"
73            print "V "
74            print volumes
75            print "****************************"
76            print "Time "
77            print time
78            print "****************************"
79            print "quantity "
80            print quantity
81            print "****************************"
82            print "vertex_coordinates"
83            print vertex_coordinates
84            print "****************************"
85       
86        self.interp = Interpolation(vertex_coordinates, volumes, alpha=0)
87        self.time = time
88
89        self.quantity_name = quantity_name
90        self.quantity = quantity
91       
92    def interpolate_xya(self, file_name):
93        """
94        Given a point file, interpolate the height w.r.t. time at the points
95        specified in the point file.
96
97        Input;
98        file_name - the xya file
99        """
100       
101        from load_mesh.loadASCII import import_points_file
102       
103        point_dict = import_points_file(file_name)
104        self.point_coordinates = point_dict['pointlist']
105        self.interp.build_interpolation_matrix_A(self.point_coordinates)
106        self.interpolated_quantity_raw = self.interp.interpolate(transpose(self.quantity))
107        self.interpolated_quantity = {}
108        for i,time_slice in enumerate(self.time):
109            self.interpolated_quantity[str(time_slice)] = self.interpolated_quantity_raw[:,i]
110       
111           
112    def read_sww(self, file_name, quantity_name):
113        """
114        Read in an sww file.
115
116        Input;
117        file_name - the sww file
118
119        Output;
120        x - Vector of x values
121        y - Vector of y values
122        z - Vector of bed elevation
123        volumes - Array.  Each row has 3 values, representing
124                  the vertices that define the volume
125        time - Vector of the times where there is stage information
126        stage - array with respect to time and vertices (x,y)
127        """
128       
129        #FIXME Have this reader as part of data_manager?
130
131        from Scientific.IO.NetCDF import NetCDFFile     
132        import tempfile
133        import sys
134        import os
135           
136        #Check contents
137        #Get NetCDF
138   
139        # see if the file is there.  Throw a QUIET IO error if it isn't
140        fd = open(file_name,'r')
141        fd.close()
142
143        #throws prints to screen if file not present
144       
145        junk = tempfile.mktemp(".txt")
146        fd = open(junk,'w')
147        stdout = sys.stdout
148        sys.stdout = fd
149        fid = NetCDFFile(file_name, 'r') 
150        sys.stdout = stdout
151        fd.close()
152        #clean up
153        os.remove(junk)         
154     
155        # Get the variables
156        x = fid.variables['x'][:]
157        y = fid.variables['y'][:]
158        volumes = fid.variables['volumes'][:] 
159        time = fid.variables['time'][:]
160        try:
161            if quantity_name == 'depth':
162                z = fid.variables['elevation'][:]
163                stage = fid.variables['stage'][:,:]
164                quantity = stage - z  # 2D, using broadcasting
165                #print quantity
166            else:
167                quantity = fid.variables[quantity_name][:,:]   # 2D
168        except (KeyError, IndexError),e:
169            fid.close()
170            raise KeyError
171           
172        fid.close()
173        return x, y, volumes, time, quantity
174
175    def write_depth_xya(self,point_file, delimiter = ','):
176        """
177        pre condition:
178        point attributes have been determined
179        (interpolate_xya has been called)
180        The time list is defined
181        """
182        from load_mesh.loadASCII import export_points_file
183
184        xya_dict = {}
185        xya_dict['pointlist'] = self.point_coordinates
186        xya_dict['attributelist'] = self.interpolated_quantity
187        export_points_file(point_file, xya_dict)
188       
189       
190#-------------------------------------------------------------
191if __name__ == "__main__":
192    """
193    Load in an sww file and an xya file and return an xya file
194    """
195    import os, sys
196    usage = "usage: %s pyvolution_results.sww points.xya depth.xya [depth|stage|(other quantities)]" %         os.path.basename(sys.argv[0])
197    if len(sys.argv) < 4:
198        print usage
199    else:
200        sww_file = sys.argv[1]
201        point_file_in = sys.argv[2]
202        point_file_out = sys.argv[3]
203        if len(sys.argv) == 5:
204            quantity_name = sys.argv[4]
205        else:
206            quantity_name = DEFAULT_QUANTITY   
207        #print "quantity",quantity
208
209
210
211
212
213
214
215
216
Note: See TracBrowser for help on using the repository browser.