source: anuga_core/source/obsolete_code/interpolate_sww.py @ 7585

Last change on this file since 7585 was 4163, checked in by duncan, 18 years ago

moving obsolete code

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