Changeset 6553 for branches/numpy/anuga/abstract_2d_finite_volumes/util.py
- Timestamp:
- Mar 19, 2009, 1:43:34 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/util.py
r6533 r6553 925 925 verbose = False): 926 926 927 # FIXME(Ole): Shouldn't print statements here be governed by verbose? 927 928 assert type(gauge_filename) == type(''), 'Gauge filename must be a string' 928 929 … … 971 972 raise msg 972 973 973 print 'swwfile', swwfile 974 if verbose: 975 print 'swwfile', swwfile 974 976 975 977 # Extract parent dir name and use as label … … 1309 1311 thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \ 1310 1312 + gaugeloc + '.csv' 1311 fid_out = open(thisfile, 'w') 1312 s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' 1313 fid_out.write(s) 1313 if j == 0: 1314 fid_out = open(thisfile, 'w') 1315 s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' 1316 fid_out.write(s) 1314 1317 1315 1318 #### generate quantities ####### … … 2368 2371 file.close() 2369 2372 2370 2371 2373 ## 2372 2374 # @brief ?? 2373 # @param sww_file??2375 # @param ?? 2374 2376 # @param gauge_file ?? 2375 2377 # @param out_name ?? … … 2383 2385 'xmomentum', 'ymomentum'], 2384 2386 verbose=False, 2385 use_cache =True):2387 use_cache=True): 2386 2388 """ 2387 2389 2388 2390 Inputs: 2389 2390 2391 NOTE: if using csv2timeseries_graphs after creating csv file, 2391 2392 it is essential to export quantities 'depth' and 'elevation'. … … 2412 2413 myfile_2_point1.csv if <out_name> ='myfile_2_' 2413 2414 2414 2415 2415 They will all have a header 2416 2416 … … 2436 2436 import string 2437 2437 from anuga.shallow_water.data_manager import get_all_swwfiles 2438 2439 # quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']2440 #print "points",points2441 2438 2442 2439 assert type(gauge_file) == type(''), 'Gauge filename must be a string' … … 2455 2452 point_name = [] 2456 2453 2457 # read point info from file2454 # read point info from file 2458 2455 for i,row in enumerate(point_reader): 2459 # read header and determine the column numbers to read correcty.2456 # read header and determine the column numbers to read correctly. 2460 2457 if i==0: 2461 2458 for j,value in enumerate(row): … … 2489 2486 base_name=base, 2490 2487 verbose=verbose) 2488 #print 'sww files just after get_all_swwfiles()', sww_files 2489 # fudge to get SWW files in 'correct' order, oldest on the left 2490 sww_files.sort() 2491 2492 if verbose: 2493 print 'sww files', sww_files 2491 2494 2492 2495 #to make all the quantities lower case for file_function … … 2497 2500 2498 2501 core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] 2499 2502 2503 gauge_file = out_name 2504 2505 heading = [quantity for quantity in quantities] 2506 heading.insert(0,'time') 2507 heading.insert(1,'hours') 2508 2509 #create a list of csv writers for all the points and write header 2510 points_writer = [] 2511 for point_i,point in enumerate(points): 2512 points_writer.append(writer(file(dir_name + sep + gauge_file 2513 + point_name[point_i] + '.csv', "wb"))) 2514 points_writer[point_i].writerow(heading) 2515 2516 if verbose: print 'Writing csv files' 2517 2518 quake_offset_time = None 2519 2500 2520 for sww_file in sww_files: 2501 2521 sww_file = join(dir_name, sww_file+'.sww') … … 2505 2525 verbose=verbose, 2506 2526 use_cache=use_cache) 2507 gauge_file = out_name 2508 2509 heading = [quantity for quantity in quantities] 2510 heading.insert(0,'time') 2511 2512 #create a list of csv writers for all the points and write header 2513 points_writer = [] 2514 for i,point in enumerate(points): 2515 points_writer.append(writer(file(dir_name + sep + gauge_file 2516 + point_name[i] + '.csv', "wb"))) 2517 points_writer[i].writerow(heading) 2518 2519 if verbose: print 'Writing csv files' 2520 2521 for time in callable_sww.get_time(): 2522 for point_i, point in enumerate(points_array): 2523 #add domain starttime to relative time. 2524 points_list = [time + callable_sww.starttime] 2525 point_quantities = callable_sww(time,point_i) 2526 2527 for quantity in quantities: 2528 if quantity == NAN: 2529 print 'quantity does not exist in' % callable_sww.get_name 2530 else: 2531 if quantity == 'stage': 2532 points_list.append(point_quantities[0]) 2533 2534 if quantity == 'elevation': 2535 points_list.append(point_quantities[1]) 2536 2537 if quantity == 'xmomentum': 2538 points_list.append(point_quantities[2]) 2539 2540 if quantity == 'ymomentum': 2541 points_list.append(point_quantities[3]) 2542 2543 if quantity == 'depth': 2544 points_list.append(point_quantities[0] 2545 - point_quantities[1]) 2546 2547 if quantity == 'momentum': 2548 momentum = sqrt(point_quantities[2]**2 2549 + point_quantities[3]**2) 2550 points_list.append(momentum) 2551 2552 if quantity == 'speed': 2553 #if depth is less than 0.001 then speed = 0.0 2554 if point_quantities[0] - point_quantities[1] < 0.001: 2555 vel = 0.0 2556 else: 2557 if point_quantities[2] < 1.0e6: 2558 momentum = sqrt(point_quantities[2]**2 2559 + point_quantities[3]**2) 2560 # vel = momentum/depth 2561 vel = momentum / (point_quantities[0] 2562 - point_quantities[1]) 2563 # vel = momentum/(depth + 1.e-6/depth) 2527 2528 if quake_offset_time is None: 2529 quake_offset_time = callable_sww.starttime 2530 2531 for time in callable_sww.get_time(): 2532 for point_i, point in enumerate(points_array): 2533 # print 'gauge_file = ', str(point_name[point_i]) 2534 #print 'point_i = ', str(point_i), ' point is = ', str(point) 2535 #add domain starttime to relative time. 2536 quake_time = time + quake_offset_time 2537 points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug 2538 #print 'point list = ', str(points_list) 2539 point_quantities = callable_sww(time,point_i) 2540 #print 'point quantities = ', str(point_quantities) 2541 2542 for quantity in quantities: 2543 if quantity == NAN: 2544 print 'quantity does not exist in' % callable_sww.get_name 2545 else: 2546 if quantity == 'stage': 2547 points_list.append(point_quantities[0]) 2548 2549 if quantity == 'elevation': 2550 points_list.append(point_quantities[1]) 2551 2552 if quantity == 'xmomentum': 2553 points_list.append(point_quantities[2]) 2554 2555 if quantity == 'ymomentum': 2556 points_list.append(point_quantities[3]) 2557 2558 if quantity == 'depth': 2559 points_list.append(point_quantities[0] 2560 - point_quantities[1]) 2561 2562 if quantity == 'momentum': 2563 momentum = sqrt(point_quantities[2]**2 2564 + point_quantities[3]**2) 2565 points_list.append(momentum) 2566 2567 if quantity == 'speed': 2568 #if depth is less than 0.001 then speed = 0.0 2569 if point_quantities[0] - point_quantities[1] < 0.001: 2570 vel = 0.0 2564 2571 else: 2565 momentum = 0 2566 vel = 0 2572 if point_quantities[2] < 1.0e6: 2573 momentum = sqrt(point_quantities[2]**2 2574 + point_quantities[3]**2) 2575 vel = momentum / (point_quantities[0] 2576 - point_quantities[1]) 2577 else: 2578 momentum = 0 2579 vel = 0 2580 2581 points_list.append(vel) 2567 2582 2568 points_list.append(vel) 2569 2570 if quantity == 'bearing': 2571 points_list.append(calc_bearing(point_quantities[2], 2572 point_quantities[3])) 2573 2574 points_writer[point_i].writerow(points_list) 2575 2583 if quantity == 'bearing': 2584 points_list.append(calc_bearing(point_quantities[2], 2585 point_quantities[3])) 2586 2587 points_writer[point_i].writerow(points_list) 2576 2588 2577 2589 ##
Note: See TracChangeset
for help on using the changeset viewer.
