- Timestamp:
- Oct 31, 2008, 4:54:40 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5880 r5882 5533 5533 return base_name, files 5534 5534 5535 def write_mux(self,lat_long_points, time_step_count, time_step, 5536 depth=None, ha=None, ua=None, va=None 5537 ): 5535 def write_mux(self, lat_long_points, time_step_count, time_step, 5536 depth=None, ha=None, ua=None, va=None): 5538 5537 """ 5539 5538 This will write 3 non-gridded mux files, for testing. … … 5542 5541 Depth and ua will be the Northing value. 5543 5542 """ 5543 5544 5544 #print "lat_long_points", lat_long_points 5545 5545 #print "time_step_count",time_step_count 5546 5546 #print "time_step", 5547 5547 5548 5548 5549 points_num = len(lat_long_points) … … 5583 5584 os.close(file_handle) 5584 5585 os.remove(base_name) 5585 5586 5586 5587 files = [] 5587 for i, q in enumerate(quantities):5588 for i, q in enumerate(quantities): 5588 5589 quantities_init[i] = ensure_numeric(quantities_init[i]) 5589 5590 #print "HA_init", HA_init … … 6063 6064 6064 6065 files = [] 6065 for i, q in enumerate(quantities):6066 for i, q in enumerate(quantities): 6066 6067 q_time = zeros((time_step_count, points_num), Float) 6067 6068 quantities_init[i] = ensure_numeric(quantities_init[i]) 6068 6069 for time in range(time_step_count): 6069 q_time[time,:]=quantities_init[i][:,time] 6070 #print i, q, time, quantities_init[i][:,time] 6071 q_time[time,:] = quantities_init[i][:,time] 6072 #print i, q, time, q_time[time, :] 6070 6073 6071 6074 #Write C files 6072 6075 columns = 3 # long, lat , depth 6073 6076 file = base_name + mux_names[i] 6074 #print "base_name file",file 6077 6078 #print "base_name file", file 6075 6079 f = open(file, 'wb') 6076 6080 files.append(file) … … 6106 6110 # Find when first station starts recording 6107 6111 min_tstep = min(first_tstep) 6108 # Find when all stations have stop ed recording6112 # Find when all stations have stopped recording 6109 6113 max_tstep = max(last_tstep) 6110 6114 … … 6114 6118 for point_i in range(points_num): 6115 6119 if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]: 6116 f.write(pack('f',q_time[time,point_i])) 6120 #print 'writing', time, point_i, q_time[time, point_i] 6121 f.write(pack('f', q_time[time, point_i])) 6122 6123 f.close() 6117 6124 6118 6125 return base_name, files … … 6304 6311 assert allclose(yvelocity,va) 6305 6312 6306 6307 def test_read_mux_platform_problem(self): 6308 """test_read_mux_platform_problem 6313 6314 6315 def test_read_mux_platform_problem1(self): 6316 """test_read_mux_platform_problem1 6309 6317 6310 6318 This is to test a situation where read_mux returned 6311 6319 wrong values Win32 6320 6321 This test passes but test_read_mux_platform_problem2 does not 6312 6322 """ 6313 6323 … … 6320 6330 time_step_count = 10 6321 6331 time_step = 0.2 6322 6323 6332 times_ref = arange(0, time_step_count*time_step, time_step) 6324 #print 'time vector', times_ref 6325 6333 6326 6334 lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)] 6327 6335 n = len(lat_long_points) 6328 6336 6329 6330 6337 # Create different timeseries starting and ending at different times 6331 first_tstep=ones(n, Int)6338 first_tstep=ones(n, Int) 6332 6339 first_tstep[0]+=2 # Point 0 starts at 2 6333 6340 first_tstep[1]+=4 # Point 1 starts at 4 … … 6338 6345 last_tstep[1]-=2 # Point 1 ends 2 steps early 6339 6346 last_tstep[4]-=3 # Point 4 ends 3 steps early 6340 6341 6342 6347 6343 6348 # Create varying elevation data (positive values for seafloor) … … 6346 6351 gauge_depth[i] += i**2 6347 6352 6348 #print 'gauge_depth', gauge_depth6349 6350 6353 # Create data to be written to first mux file 6351 6354 ha0=2*ones((n,time_step_count),Float) … … 6361 6364 for i in range(n): 6362 6365 # For each point 6363 6364 6366 for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count): 6365 6367 # For timesteps before and after recording range 6366 6368 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0 6367 6368 6369 6369 6370 6370 # Write first mux file to be combined by urs2sts … … 6377 6377 va=va0) 6378 6378 6379 6380 6379 # Create ordering file 6380 permutation = ensure_numeric([4,0,2]) 6381 6382 _, ordering_filename = tempfile.mkstemp('') 6383 order_fid = open(ordering_filename, 'w') 6384 order_fid.write('index, longitude, latitude\n') 6385 for index in permutation: 6386 order_fid.write('%d, %f, %f\n' %(index, 6387 lat_long_points[index][1], 6388 lat_long_points[index][0])) 6389 order_fid.close() 6390 6391 6392 6393 # ------------------------------------- 6394 # Now read files back and check values 6395 weights = ensure_numeric([1.0]) 6396 6397 # For each quantity read the associated list of source mux2 file with 6398 # extention associated with that quantity 6399 file_params=-1*ones(3,Float) #[nsta,dt,nt] 6400 OFFSET = 5 6401 6402 for j, file in enumerate(filesI): 6403 data = read_mux2(1, [file], weights, file_params, permutation, verbose) 6404 6405 number_of_selected_stations = data.shape[0] 6406 6407 # Index where data ends and parameters begin 6408 parameters_index = data.shape[1]-OFFSET 6409 6410 for i in range(number_of_selected_stations): 6411 if j == 0: assert allclose(data[i][:parameters_index], ha0[permutation[i], :]) 6412 if j == 1: assert allclose(data[i][:parameters_index], ua0[permutation[i], :]) 6413 if j == 2: assert allclose(data[i][:parameters_index], va0[permutation[i], :]) 6414 6415 6416 6417 6418 6419 def test_read_mux_platform_problem2(self): 6420 """test_read_mux_platform_problem2 6421 6422 This is to test a situation where read_mux returned 6423 wrong values Win32 6424 6425 This test does not pass but test_read_mux_platform_problem1 does 6426 """ 6427 6428 from Numeric import sin, cos 6429 from urs_ext import read_mux2 6430 6431 verbose = False 6432 6433 tide = 1.5 6434 time_step_count = 10 6435 time_step = 0.2 6436 6437 times_ref = arange(0, time_step_count*time_step, time_step) 6438 6439 lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)] 6440 n = len(lat_long_points) 6441 6442 # Create different timeseries starting and ending at different times 6443 first_tstep=ones(n,Int) 6444 first_tstep[0]+=2 # Point 0 starts at 2 6445 first_tstep[1]+=4 # Point 1 starts at 4 6446 first_tstep[2]+=3 # Point 2 starts at 3 6447 6448 last_tstep=(time_step_count)*ones(n,Int) 6449 last_tstep[0]-=1 # Point 0 ends 1 step early 6450 last_tstep[1]-=2 # Point 1 ends 2 steps early 6451 last_tstep[4]-=3 # Point 4 ends 3 steps early 6452 6453 # Create varying elevation data (positive values for seafloor) 6454 gauge_depth=20*ones(n,Float) 6455 for i in range(n): 6456 gauge_depth[i] += i**2 6457 6381 6458 # Create data to be written to second mux file 6382 6459 ha1=ones((n,time_step_count),Float) … … 6398 6475 va1[3]=2*sin(times_ref-0.71) 6399 6476 6400 6401 6477 # Ensure data used to write mux file to be zero when gauges are 6402 6478 # not recording 6403 6479 for i in range(n): 6404 6480 # For each point 6405 6406 6481 for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count): 6407 6482 # For timesteps before and after recording range … … 6409 6484 6410 6485 6411 6486 #print 'Second station to be written to MUX' 6487 #print 'ha', ha1[2,:] 6488 #print 'ua', ua1[2,:] 6489 #print 'va', va1[2,:] 6490 6412 6491 # Write second mux file to be combined by urs2sts 6413 6492 base_nameII, filesII = self.write_mux2(lat_long_points, … … 6440 6519 # For each quantity read the associated list of source mux2 file with 6441 6520 # extention associated with that quantity 6442 file_params=-1*ones(3,Float) # [nsta,dt,nt]6521 file_params=-1*ones(3,Float) # [nsta,dt,nt] 6443 6522 OFFSET = 5 6444 6523 6445 # FILE I 6446 for j, file in enumerate(filesI): 6524 for j, file in enumerate(filesII): 6447 6525 data = read_mux2(1, [file], weights, file_params, permutation, verbose) 6448 6526 6449 nsta=int(file_params[0]) 6450 dt=file_params[1] 6451 6527 #print 6528 #print 'j:', j 6529 #print data[2][:] 6530 #print file_params 6531 #print 6532 6452 6533 number_of_selected_stations = data.shape[0] 6453 6534 … … 6455 6536 parameters_index = data.shape[1]-OFFSET 6456 6537 6457 times=dt*arange(parameters_index)6458 latitudes=zeros(number_of_selected_stations, Float)6459 longitudes=zeros(number_of_selected_stations, Float)6460 elevation=zeros(number_of_selected_stations, Float)6461 6538 quantity=zeros((number_of_selected_stations, parameters_index), Float) 6462 6539 6463 starttime=1e166464 6540 6465 6541 for i in range(number_of_selected_stations): 6466 6542 quantity[i][:]=data[i][:parameters_index] 6467 6543 6468 #print i, j, parameters_index6469 #print quantity[i][:]6470 6471 if j == 0:6472 # HA6473 if i == 0:6474 assert allclose(quantity[i][:],6475 [2., 2., 2., 2., 2., 2., 2., 0., 0., 0.])6476 if i == 1:6477 assert allclose(quantity[i][:],6478 [0., 0., 2., 3., 4., 5., 6., 7., 8., 0.])6479 if i == 2:6480 assert allclose(quantity[i][:],6481 [0., 0., 0., 23., 24., 25., 26., 27., 28., 29.])6482 6483 if j == 1:6484 # UA6485 if i == 0:6486 assert allclose(quantity[i][:],6487 [5., 5., 5., 5., 5., 5., 5., 0., 0., 0.])6488 if i == 1:6489 assert allclose(quantity[i][:],6490 [0., 0., 5., 5., 5., 5., 5., 5., 5., 0.])6491 if i == 2:6492 assert allclose(quantity[i][:],6493 [0., 0., 0., 5., 5., 5., 5., 5., 5., 5.])6494 if j == 2:6495 # VA6496 if i == 0:6497 assert allclose(quantity[i][:],6498 [-10., -10., -10., -10., -10., -10., -10., 0., 0., 0.])6499 if i == 1:6500 assert allclose(quantity[i][:],6501 [0., 0., -10., -10., -10., -10., -10., -10., -10., 0.])6502 if i == 2:6503 assert allclose(quantity[i][:],6504 [0., 0., 0., -10., -10., -10., -10., -10., -10., -10.])6505 6506 6507 # FILE II6508 for j, file in enumerate(filesII):6509 data = read_mux2(1, [file], weights, file_params, permutation, verbose)6510 6511 nsta=int(file_params[0])6512 dt=file_params[1]6513 6514 number_of_selected_stations = data.shape[0]6515 6516 # Index where data ends and parameters begin6517 parameters_index = data.shape[1]-OFFSET6518 6519 times=dt*arange(parameters_index)6520 latitudes=zeros(number_of_selected_stations, Float)6521 longitudes=zeros(number_of_selected_stations, Float)6522 elevation=zeros(number_of_selected_stations, Float)6523 quantity=zeros((number_of_selected_stations, parameters_index), Float)6524 6525 starttime=1e166526 6527 for i in range(number_of_selected_stations):6528 quantity[i][:]=data[i][:parameters_index]6529 6530 6544 #print i, parameters_index 6531 6545 #print quantity[i][:] 6532 6533 6534 if j == 0: 6535 # HA 6536 if i == 0: 6537 assert allclose(quantity[i][:], 6538 [-0.64421767, -0.29552022, 0.09983341, 0.47942555, 6539 0.78332692, 0.9635582, 0.99166483, 0., 0., 0.]) 6540 if i == 1: 6541 assert allclose(quantity[i][:], 6542 [0., 0., 0.38941833, 0.56464249, 0.71735609, 0.84147096, 6543 0.93203908, 0.98544973, 0.99957359, 0.]) 6544 if i == 2: 6545 assert allclose(quantity[i][:], 6546 [0., 0., 0., 3.377316, -0.29187071, -3.78401256, 6547 -4.98082304, -3.15633321, 0.58274603, 3.9683392]) 6548 6549 6550 if j == 1: 6551 # UA 6552 if i == 0: 6553 assert allclose(quantity[i][:], 6554 [2., 2., 2., 2., 2., 2., 2., 0., 0., 0.]) 6555 if i == 1: 6556 assert allclose(quantity[i][:], 6557 [0., 0., 2.76318288, 2.47600675, 2.09012008, 1.62090695, 6558 1.08707321, 0.5099014, -0.08759857, 0.]) 6559 if i == 2: 6560 assert allclose(quantity[i][:], 6561 [0., 0., 0., 33., 34., 35., 36., 37., 38., 39.]) 6562 6546 6547 6548 if j == 0: assert allclose(data[i][:parameters_index], ha1[permutation[i], :]) 6549 if j == 1: assert allclose(data[i][:parameters_index], ua1[permutation[i], :]) 6563 6550 if j == 2: 6564 # VA 6565 if i == 0: 6566 assert allclose(quantity[i][:], 6567 [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) 6568 if i == 1: 6569 assert allclose(quantity[i][:], 6570 [0., 0., 1.78313661, 1.92754185, 1.99510205, 1.98312378, 6571 1.89208472, 1.72561419, 1.49034882, 0.]) 6572 if i == 2: 6573 assert allclose(quantity[i][:], 6574 [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) 6575 6576 6577 6578 6551 # FIXME (Ole): This is where the output is wrong on Win32 6552 6553 #print 6554 #print j, i 6555 #print data[i][:parameters_index] 6556 #print va1[permutation[i], :] 6557 6558 if j == 2 and i == 1: 6559 pass 6560 # Skip assert for this combination for now as the second error is more obvious 6561 else: 6562 assert allclose(data[i][:parameters_index], va1[permutation[i], :]) 6563 6564 6579 6565 def test_urs2sts0(self): 6580 6566 """ … … 8031 8017 #print ymomentum 8032 8018 8033 from sys import platform 8034 if platform == 'win32': 8035 # FIXME (Ole) - one array element differs on windoze. Why? 8036 pass 8037 else: 8038 # It works fine on Linux 32 and 64 bit platforms. 8039 assert allclose(transpose(va_ref*depth_ref), ymomentum) 8019 assert allclose(transpose(va_ref*depth_ref), ymomentum) 8040 8020 8041 8021 # check the elevation values. … … 10933 10913 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources') 10934 10914 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources') 10935 #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_prob ')10915 #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem') 10936 10916 10937 10917
Note: See TracChangeset
for help on using the changeset viewer.
