Changeset 7767
 Timestamp:
 Jun 2, 2010, 3:00:14 PM (13 years ago)
 Location:
 trunk/anuga_core/source/anuga/file
 Files:

 1 added
 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/file/sww.py
r7762 r7767 984 984 log.critical('') 985 985 986 987 988 989 ## 990 # @brief Get the extents of a NetCDF data file. 991 # @param file_name The path to the SWW file. 992 # @return A list of x, y, z and stage limits (min, max). 993 def extent_sww(file_name): 994 """Read in an sww file. 995 996 Input: 997 file_name  the sww file 998 999 Output: 1000 A list: [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)] 1001 """ 1002 1003 from Scientific.IO.NetCDF import NetCDFFile 1004 1005 #Get NetCDF 1006 fid = NetCDFFile(file_name, netcdf_mode_r) 1007 1008 # Get the variables 1009 x = fid.variables['x'][:] 1010 y = fid.variables['y'][:] 1011 stage = fid.variables['stage'][:] 1012 1013 fid.close() 1014 1015 return [min(x), max(x), min(y), max(y), num.min(stage), num.max(stage)] 1016 1017 1018 ## 1019 # @brief 1020 # @param filename 1021 # @param boundary 1022 # @param t 1023 # @param fail_if_NaN 1024 # @param NaN_filler 1025 # @param verbose 1026 # @param very_verbose 1027 # @return 1028 def load_sww_as_domain(filename, boundary=None, t=None, 1029 fail_if_NaN=True, NaN_filler=0, 1030 verbose=False, very_verbose=False): 1031 """ 1032 Usage: domain = load_sww_as_domain('file.sww', 1033 t=time (default = last time in file)) 1034 1035 Boundary is not recommended if domain.smooth is not selected, as it 1036 uses unique coordinates, but not unique boundaries. This means that 1037 the boundary file will not be compatable with the coordinates, and will 1038 give a different final boundary, or crash. 1039 """ 1040 1041 from Scientific.IO.NetCDF import NetCDFFile 1042 from shallow_water import Domain 1043 1044 # initialise NaN. 1045 NaN = 9.969209968386869e+036 1046 1047 if verbose: log.critical('Reading from %s' % filename) 1048 1049 fid = NetCDFFile(filename, netcdf_mode_r) # Open existing file for read 1050 time = fid.variables['time'] # Timesteps 1051 if t is None: 1052 t = time[1] 1053 time_interp = get_time_interp(time,t) 1054 1055 # Get the variables as numeric arrays 1056 x = fid.variables['x'][:] # xcoordinates of vertices 1057 y = fid.variables['y'][:] # ycoordinates of vertices 1058 elevation = fid.variables['elevation'] # Elevation 1059 stage = fid.variables['stage'] # Water level 1060 xmomentum = fid.variables['xmomentum'] # Momentum in the xdirection 1061 ymomentum = fid.variables['ymomentum'] # Momentum in the ydirection 1062 1063 starttime = fid.starttime[0] 1064 volumes = fid.variables['volumes'][:] # Connectivity 1065 coordinates = num.transpose(num.asarray([x.tolist(), y.tolist()])) 1066 # FIXME (Ole): Something like this might be better: 1067 # concatenate((x, y), axis=1) 1068 # or concatenate((x[:,num.newaxis], x[:,num.newaxis]), axis=1) 1069 1070 conserved_quantities = [] 1071 interpolated_quantities = {} 1072 other_quantities = [] 1073 1074 # get geo_reference 1075 try: # sww files don't have to have a geo_ref 1076 geo_reference = Geo_reference(NetCDFObject=fid) 1077 except: # AttributeError, e: 1078 geo_reference = None 1079 1080 if verbose: log.critical(' getting quantities') 1081 1082 for quantity in fid.variables.keys(): 1083 dimensions = fid.variables[quantity].dimensions 1084 if 'number_of_timesteps' in dimensions: 1085 conserved_quantities.append(quantity) 1086 interpolated_quantities[quantity] = \ 1087 interpolated_quantity(fid.variables[quantity][:], time_interp) 1088 else: 1089 other_quantities.append(quantity) 1090 1091 other_quantities.remove('x') 1092 other_quantities.remove('y') 1093 #other_quantities.remove('z') 1094 other_quantities.remove('volumes') 1095 try: 1096 other_quantities.remove('stage_range') 1097 other_quantities.remove('xmomentum_range') 1098 other_quantities.remove('ymomentum_range') 1099 other_quantities.remove('elevation_range') 1100 except: 1101 pass 1102 1103 conserved_quantities.remove('time') 1104 1105 if verbose: log.critical(' building domain') 1106 1107 # From domain.Domain: 1108 # domain = Domain(coordinates, volumes,\ 1109 # conserved_quantities = conserved_quantities,\ 1110 # other_quantities = other_quantities,zone=zone,\ 1111 # xllcorner=xllcorner, yllcorner=yllcorner) 1112 1113 # From shallow_water.Domain: 1114 coordinates = coordinates.tolist() 1115 volumes = volumes.tolist() 1116 # FIXME:should this be in mesh? (peter row) 1117 if fid.smoothing == 'Yes': 1118 unique = False 1119 else: 1120 unique = True 1121 if unique: 1122 coordinates, volumes, boundary = weed(coordinates, volumes,boundary) 1123 1124 1125 1126 try: 1127 domain = Domain(coordinates, volumes, boundary) 1128 except AssertionError, e: 1129 fid.close() 1130 msg = 'Domain could not be created: %s. ' \ 1131 'Perhaps use "fail_if_NaN=False and NaN_filler = ..."' % e 1132 raise DataDomainError, msg 1133 1134 if not boundary is None: 1135 domain.boundary = boundary 1136 1137 domain.geo_reference = geo_reference 1138 1139 domain.starttime = float(starttime) + float(t) 1140 domain.time = 0.0 1141 1142 for quantity in other_quantities: 1143 try: 1144 NaN = fid.variables[quantity].missing_value 1145 except: 1146 pass # quantity has no missing_value number 1147 X = fid.variables[quantity][:] 1148 if very_verbose: 1149 log.critical(' %s' % str(quantity)) 1150 log.critical(' NaN = %s' % str(NaN)) 1151 log.critical(' max(X)') 1152 log.critical(' %s' % str(max(X))) 1153 log.critical(' max(X)==NaN') 1154 log.critical(' %s' % str(max(X)==NaN)) 1155 log.critical('') 1156 if max(X) == NaN or min(X) == NaN: 1157 if fail_if_NaN: 1158 msg = 'quantity "%s" contains no_data entry' % quantity 1159 raise DataMissingValuesError, msg 1160 else: 1161 data = (X != NaN) 1162 X = (X*data) + (data==0)*NaN_filler 1163 if unique: 1164 X = num.resize(X, (len(X)/3, 3)) 1165 domain.set_quantity(quantity, X) 1166 # 1167 for quantity in conserved_quantities: 1168 try: 1169 NaN = fid.variables[quantity].missing_value 1170 except: 1171 pass # quantity has no missing_value number 1172 X = interpolated_quantities[quantity] 1173 if very_verbose: 1174 log.critical(' %s' % str(quantity)) 1175 log.critical(' NaN = %s' % str(NaN)) 1176 log.critical(' max(X)') 1177 log.critical(' %s' % str(max(X))) 1178 log.critical(' max(X)==NaN') 1179 log.critical(' %s' % str(max(X)==NaN)) 1180 log.critical('') 1181 if max(X) == NaN or min(X) == NaN: 1182 if fail_if_NaN: 1183 msg = 'quantity "%s" contains no_data entry' % quantity 1184 raise DataMissingValuesError, msg 1185 else: 1186 data = (X != NaN) 1187 X = (X*data) + (data==0)*NaN_filler 1188 if unique: 1189 X = num.resize(X, (X.shape[0]/3, 3)) 1190 domain.set_quantity(quantity, X) 1191 1192 fid.close() 1193 1194 return domain 1195 1196 1197 ## 1198 # @brief 1199 # @parm time 1200 # @param t 1201 # @return An (index, ration) tuple. 1202 def get_time_interp(time, t=None): 1203 #Finds the ratio and index for time interpolation. 1204 #It is borrowed from previous abstract_2d_finite_volumes code. 1205 if t is None: 1206 t=time[1] 1207 index = 1 1208 ratio = 0. 1209 else: 1210 T = time 1211 tau = t 1212 index=0 1213 msg = 'Time interval derived from file %s [%s:%s]' \ 1214 % ('FIXMEfilename', T[0], T[1]) 1215 msg += ' does not match model time: %s' % tau 1216 if tau < time[0]: raise DataTimeError, msg 1217 if tau > time[1]: raise DataTimeError, msg 1218 while tau > time[index]: index += 1 1219 while tau < time[index]: index = 1 1220 if tau == time[index]: 1221 #Protect against case where tau == time[1] (last time) 1222 #  also works in general when tau == time[i] 1223 ratio = 0 1224 else: 1225 #t is now between index and index+1 1226 ratio = (tau  time[index])/(time[index+1]  time[index]) 1227 1228 return (index, ratio) 1229 1230 1231 1232 ## 1233 # @brief Interpolate a quantity wrt time. 1234 # @param saved_quantity The quantity to interpolate. 1235 # @param time_interp (index, ratio) 1236 # @return A vector of interpolated values. 1237 def interpolated_quantity(saved_quantity, time_interp): 1238 '''Given an index and ratio, interpolate quantity with respect to time.''' 1239 1240 index, ratio = time_interp 1241 1242 Q = saved_quantity 1243 1244 if ratio > 0: 1245 q = (1ratio)*Q[index] + ratio*Q[index+1] 1246 else: 1247 q = Q[index] 1248 1249 #Return vector of interpolated values 1250 return q 1251 1252 1253 1254 1255 ## 1256 # @brief 1257 # @param coordinates 1258 # @param volumes 1259 # @param boundary 1260 def weed(coordinates, volumes, boundary=None): 1261 if isinstance(coordinates, num.ndarray): 1262 coordinates = coordinates.tolist() 1263 if isinstance(volumes, num.ndarray): 1264 volumes = volumes.tolist() 1265 1266 unique = False 1267 point_dict = {} 1268 same_point = {} 1269 for i in range(len(coordinates)): 1270 point = tuple(coordinates[i]) 1271 if point_dict.has_key(point): 1272 unique = True 1273 same_point[i] = point 1274 #to change all point i references to point j 1275 else: 1276 point_dict[point] = i 1277 same_point[i] = point 1278 1279 coordinates = [] 1280 i = 0 1281 for point in point_dict.keys(): 1282 point = tuple(point) 1283 coordinates.append(list(point)) 1284 point_dict[point] = i 1285 i += 1 1286 1287 for volume in volumes: 1288 for i in range(len(volume)): 1289 index = volume[i] 1290 if index > 1: 1291 volume[i] = point_dict[same_point[index]] 1292 1293 new_boundary = {} 1294 if not boundary is None: 1295 for segment in boundary.keys(): 1296 point0 = point_dict[same_point[segment[0]]] 1297 point1 = point_dict[same_point[segment[1]]] 1298 label = boundary[segment] 1299 #FIXME should the bounday attributes be concaterated 1300 #('exterior, pond') or replaced ('pond')(peter row) 1301 1302 if new_boundary.has_key((point0, point1)): 1303 new_boundary[(point0,point1)] = new_boundary[(point0, point1)] 1304 1305 elif new_boundary.has_key((point1, point0)): 1306 new_boundary[(point1,point0)] = new_boundary[(point1, point0)] 1307 else: new_boundary[(point0, point1)] = label 1308 1309 boundary = new_boundary 1310 1311 return coordinates, volumes, boundary
Note: See TracChangeset
for help on using the changeset viewer.