 Timestamp:
 Jul 4, 2008, 11:48:37 AM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/test_data_manager.py
r5442 r5462 6319 6319 va=va) 6320 6320 6321 urs2sts(base_name,mean_stage=tide,verbose=False) 6321 urs2sts(base_name, 6322 basename_out=base_name, 6323 mean_stage=tide,verbose=False) 6322 6324 6323 6325 # now I want to check the sts file ... … … 6441 6443 6442 6444 # Call urs2sts with multiple mux files 6443 urs2sts([base_nameI, base_nameII], weights=[1.0,1.0], 6445 urs2sts([base_nameI, base_nameII], 6446 basename_out=base_nameI, 6447 weights=[1.0,1.0], 6444 6448 mean_stage=tide, 6445 6449 verbose=False) … … 6534 6538 self.delete_mux(filesII) 6535 6539 os.remove(sts_file) 6540 6541 6542 def test_urs2sts_ordering(self): 6543 """Test multiple sources with ordering file 6544 """ 6545 6546 tide=0 6547 time_step_count = 3 6548 time_step = 2 6549 lat_long_points =[(21.5,114.5),(21,114.5),(21.5,115), (21.,115.)] 6550 n=len(lat_long_points) 6551 first_tstep=ones(n,Int) 6552 first_tstep[0]+=1 6553 first_tstep[2]+=1 6554 last_tstep=(time_step_count)*ones(n,Int) 6555 last_tstep[0]=1 6556 6557 gauge_depth=20*ones(n,Float) 6558 ha=2*ones((n,time_step_count),Float) 6559 ha[0]=arange(0,time_step_count) 6560 ha[1]=arange(time_step_count,2*time_step_count) 6561 ha[2]=arange(2*time_step_count,3*time_step_count) 6562 ha[3]=arange(3*time_step_count,4*time_step_count) 6563 ua=5*ones((n,time_step_count),Float) 6564 va=10*ones((n,time_step_count),Float) 6565 6566 # Create two identical mux files to be combined by urs2sts 6567 base_nameI, filesI = self.write_mux2(lat_long_points, 6568 time_step_count, time_step, 6569 first_tstep, last_tstep, 6570 depth=gauge_depth, 6571 ha=ha, 6572 ua=ua, 6573 va=va) 6574 6575 base_nameII, filesII = self.write_mux2(lat_long_points, 6576 time_step_count, time_step, 6577 first_tstep, last_tstep, 6578 depth=gauge_depth, 6579 ha=ha, 6580 ua=ua, 6581 va=va) 6582 6583 6584 # Create ordering file 6585 permutation = [3,0,2] 6586 6587 # Do it wrongly at first and check that exception is being raised 6588 _, ordering_filename = tempfile.mkstemp('') 6589 order_fid = open(ordering_filename, 'w') 6590 order_fid.write('index, longitude, latitude\n') 6591 for index in permutation: 6592 order_fid.write('%d, %f, %f\n' %(index, 6593 lat_long_points[index][0], 6594 lat_long_points[index][1])) 6595 order_fid.close() 6596 6597 try: 6598 urs2sts([base_nameI, base_nameII], 6599 basename_out=base_nameI, 6600 ordering_filename=ordering_filename, 6601 weights=[1.0,1.0], 6602 mean_stage=tide, 6603 verbose=False) 6604 os.remove(ordering_filename) 6605 except: 6606 pass 6607 else: 6608 msg = 'Should have caught wrong lat longs' 6609 raise Exception, msg 6610 6611 6612 6613 6614 # Then do it right 6615 _, ordering_filename = tempfile.mkstemp('') 6616 order_fid = open(ordering_filename, 'w') 6617 order_fid.write('index, longitude, latitude\n') 6618 for index in permutation: 6619 order_fid.write('%d, %f, %f\n' %(index, 6620 lat_long_points[index][1], 6621 lat_long_points[index][0])) 6622 order_fid.close() 6623 6624 6625 6626 6627 # Call urs2sts with multiple mux files 6628 urs2sts([base_nameI, base_nameII], 6629 basename_out=base_nameI, 6630 ordering_filename=ordering_filename, 6631 weights=[1.0,1.0], 6632 mean_stage=tide, 6633 verbose=False) 6634 6635 # now I want to check the sts file ... 6636 sts_file = base_nameI + '.sts' 6637 6638 #Let's interrogate the sts file 6639 # Note, the sts info is not gridded. It is point data. 6640 fid = NetCDFFile(sts_file) 6641 6642 # Make x and y absolute 6643 x = fid.variables['x'][:] 6644 y = fid.variables['y'][:] 6645 6646 geo_reference = Geo_reference(NetCDFObject=fid) 6647 points = geo_reference.get_absolute(map(None, x, y)) 6648 points = ensure_numeric(points) 6649 6650 x = points[:,0] 6651 y = points[:,1] 6652 6653 for i, index in enumerate(permutation): 6654 # Check that STS points are stored in the correct order 6655 6656 # Work out the UTM coordinates sts point i 6657 zone, e, n = redfearn(lat_long_points[index][0], 6658 lat_long_points[index][1]) 6659 6660 assert allclose([x[i],y[i]], [e,n]) 6661 6662 6663 # Check the time vector 6664 times = fid.variables['time'][:] 6665 6666 times_actual = [] 6667 for i in range(time_step_count): 6668 times_actual.append(time_step * i) 6669 6670 assert allclose(ensure_numeric(times), 6671 ensure_numeric(times_actual)) 6672 6673 6674 # Check sts values 6675 stage = fid.variables['stage'][:] 6676 xmomentum = fid.variables['xmomentum'][:] 6677 ymomentum = fid.variables['ymomentum'][:] 6678 elevation = fid.variables['elevation'][:] 6679 6680 # Set original data used to write mux file to be zero when gauges are 6681 # not recdoring 6682 6683 ha[0][0]=0.0 6684 ha[0][time_step_count1]=0.0 6685 ha[2][0]=0.0 6686 ua[0][0]=0.0 6687 ua[0][time_step_count1]=0.0 6688 ua[2][0]=0.0 6689 va[0][0]=0.0 6690 va[0][time_step_count1]=0.0 6691 va[2][0]=0.0; 6692 6693 6694 # The stage stored in the .sts file should be the sum of the stage 6695 # in the two mux2 files because both have weights = 1. In this case 6696 # the mux2 files are the same so stage == 2.0 * ha 6697 #print 2.0*transpose(ha)  stage 6698 6699 ha_permutation = take(ha, permutation) 6700 ua_permutation = take(ua, permutation) 6701 va_permutation = take(va, permutation) 6702 gauge_depth_permutation = take(gauge_depth, permutation) 6703 6704 6705 assert allclose(2.0*transpose(ha_permutation), stage) # Meters 6706 6707 #Check the momentums  ua 6708 #momentum = velocity*(stageelevation) 6709 # elevation =  depth 6710 #momentum = velocity_ua *(stage+depth) 6711 6712 depth=zeros((len(lat_long_points),time_step_count),Float) 6713 for i in range(len(lat_long_points)): 6714 depth[i]=gauge_depth[i]+tide+2.0*ha[i] 6715 #2.0*ha necessary because using two files with weights=1 are used 6716 6717 depth_permutation = take(depth, permutation) 6718 6719 6720 # The xmomentum stored in the .sts file should be the sum of the ua 6721 # in the two mux2 files multiplied by the depth. 6722 assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum) 6723 6724 #Check the momentums  va 6725 #momentum = velocity*(stageelevation) 6726 # elevation =  depth 6727 #momentum = velocity_va *(stage+depth) 6728 6729 # The ymomentum stored in the .sts file should be the sum of the va 6730 # in the two mux2 files multiplied by the depth. 6731 assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum) 6732 6733 # check the elevation values. 6734 # ve since urs measures depth, sww meshers height, 6735 assert allclose(gauge_depth_permutation, elevation) #Meters 6736 6737 fid.close() 6738 self.delete_mux(filesI) 6739 self.delete_mux(filesII) 6740 os.remove(sts_file) 6741 6536 6742 6537 6743 def test_file_boundary_stsI(self): … … 6563 6769 6564 6770 sts_file=base_name 6565 urs2sts(base_name,sts_file,mean_stage=tide,verbose=False) 6771 urs2sts(base_name, 6772 sts_file, 6773 mean_stage=tide,verbose=False) 6566 6774 self.delete_mux(files) 6567 6775
Note: See TracChangeset
for help on using the changeset viewer.