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