Changeset 6070
- Timestamp:
- Dec 11, 2008, 4:32:08 PM (16 years ago)
- Location:
- anuga_core/source/anuga/abstract_2d_finite_volumes
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
r5897 r6070 1263 1263 1264 1264 check_list(['stage','xmomentum']) 1265 1266 def test_remove_lone_verts_d(self): 1265 1266 ###### 1267 # Test the remove_lone_verts() function 1268 ###### 1269 1270 def test_remove_lone_verts_a(self): 1267 1271 verts = [[0,0],[1,0],[0,1]] 1268 1272 tris = [[0,1,2]] 1269 1273 new_verts, new_tris = remove_lone_verts(verts, tris) 1270 assert new_verts == verts 1271 assert new_tris == tris 1272 1273 1274 def test_remove_lone_verts_e(self): 1274 self.failUnless(new_verts.tolist() == verts) 1275 self.failUnless(new_tris.tolist() == tris) 1276 1277 def test_remove_lone_verts_b(self): 1275 1278 verts = [[0,0],[1,0],[0,1],[99,99]] 1276 1279 tris = [[0,1,2]] 1277 1280 new_verts, new_tris = remove_lone_verts(verts, tris) 1278 assert new_verts == verts[0:3]1279 assert new_tris == tris1280 1281 def test_remove_lone_verts_ a(self):1281 self.failUnless(new_verts.tolist() == verts[0:3]) 1282 self.failUnless(new_tris.tolist() == tris) 1283 1284 def test_remove_lone_verts_c(self): 1282 1285 verts = [[99,99],[0,0],[1,0],[99,99],[0,1],[99,99]] 1283 1286 tris = [[1,2,4]] 1284 1287 new_verts, new_tris = remove_lone_verts(verts, tris) 1285 #print "new_verts", new_verts 1286 assert new_verts == [[0,0],[1,0],[0,1]] 1287 assert new_tris == [[0,1,2]] 1288 self.failUnless(new_verts.tolist() == [[0,0],[1,0],[0,1]]) 1289 self.failUnless(new_tris.tolist() == [[0,1,2]]) 1288 1290 1289 def test_remove_lone_verts_ c(self):1291 def test_remove_lone_verts_d(self): 1290 1292 verts = [[0,0],[1,0],[99,99],[0,1]] 1291 1293 tris = [[0,1,3]] 1292 1294 new_verts, new_tris = remove_lone_verts(verts, tris) 1293 #print "new_verts", new_verts 1294 assert new_verts == [[0,0],[1,0],[0,1]] 1295 assert new_tris == [[0,1,2]] 1296 1297 def test_remove_lone_verts_b(self): 1295 self.failUnless(new_verts.tolist() == [[0,0],[1,0],[0,1]]) 1296 self.failUnless(new_tris.tolist() == [[0,1,2]]) 1297 1298 def test_remove_lone_verts_e(self): 1298 1299 verts = [[0,0],[1,0],[0,1],[99,99],[99,99],[99,99]] 1299 1300 tris = [[0,1,2]] 1300 1301 new_verts, new_tris = remove_lone_verts(verts, tris) 1301 assert new_verts == verts[0:3]1302 assert new_tris == tris1302 self.failUnless(new_verts.tolist() == verts[0:3]) 1303 self.failUnless(new_tris.tolist() == tris) 1303 1304 1304 1305 def test_remove_lone_verts_e(self): 1306 verts = [[0,0],[1,0],[0,1],[99,99]] 1307 tris = [[0,1,2]] 1305 def test_remove_lone_verts_f(self): 1306 verts = [[0,0],[1,0],[99,99],[0,1],[99,99],[1,1],[99,99]] 1307 tris = [[0,1,3],[0,1,5]] 1308 1308 new_verts, new_tris = remove_lone_verts(verts, tris) 1309 assert new_verts == verts[0:3] 1310 assert new_tris == tris 1309 self.failUnless(new_verts.tolist() == [[0,0],[1,0],[0,1],[1,1]]) 1310 self.failUnless(new_tris.tolist() == [[0,1,2],[0,1,3]]) 1311 1312 ###### 1313 # 1314 ###### 1311 1315 1312 1316 def test_get_min_max_values(self): … … 1829 1833 #------------------------------------------------------------- 1830 1834 if __name__ == "__main__": 1831 suite = unittest.makeSuite(Test_Util,'test')1832 # suite = unittest.makeSuite(Test_Util,'test_sww2csv')1835 # suite = unittest.makeSuite(Test_Util,'test') 1836 suite = unittest.makeSuite(Test_Util,'test_remove_lone_verts') 1833 1837 # runner = unittest.TextTestRunner(verbosity=2) 1834 1838 runner = unittest.TextTestRunner(verbosity=1) -
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r6035 r6070 1 #12345678901234567890123456789012345678901234567890123456789012345678901234567890 1 2 """This module contains various auxiliary function used by pyvolution. 2 3 … … 15 16 16 17 from anuga.utilities.numerical_tools import ensure_numeric 18 from Scientific.IO.NetCDF import NetCDFFile 17 19 from Numeric import arange, choose, zeros, Float, array, allclose, take, compress 18 20 … … 20 22 from math import sqrt, atan, degrees 21 23 22 23 24 # FIXME (Ole): Temporary short cuts - remove and update scripts where they are used 24 # FIXME (Ole): Temporary short cuts - 25 # FIXME (Ole): remove and update scripts where they are used 25 26 from anuga.utilities.system_tools import get_revision_number 26 27 from anuga.utilities.system_tools import store_version_info 27 28 28 29 30 ## 31 # @brief Read time history of data from NetCDF file, return callable object. 32 # @param filename Name of .sww or .tms file. 33 # @param domain Associated domain object. 34 # @param quantities Name of quantity to be interpolated or a list of names. 35 # @param interpolation_points List of absolute UTM coordinates for points 36 # (N x 2) or geospatial object or 37 # points file name at which values are sought. 38 # @param time_thinning 39 # @param verbose True if this function is to be verbose. 40 # @param use_cache True means that caching of intermediate result is attempted. 41 # @param boundary_polygon 42 # @return A callable object. 29 43 def file_function(filename, 30 44 domain=None, … … 69 83 interpolation_points - list of absolute UTM coordinates for points (N x 2) 70 84 or geospatial object or points file name at which values are sought 71 85 86 time_thinning - 87 88 verbose - 89 72 90 use_cache: True means that caching of intermediate result of 73 91 Interpolation_function is attempted 74 92 75 76 See Interpolation function in anuga.fit_interpolate.interpolation for further documentation 93 boundary_polygon - 94 95 96 See Interpolation function in anuga.fit_interpolate.interpolation for 97 further documentation 77 98 """ 78 79 99 80 100 # FIXME (OLE): Should check origin of domain against that of file … … 90 110 if verbose: 91 111 msg = 'Quantities specified in file_function are None,' 92 msg += ' so I will use stage, xmomentum, and ymomentum in that order .'112 msg += ' so I will use stage, xmomentum, and ymomentum in that order' 93 113 print msg 94 95 114 quantities = ['stage', 'xmomentum', 'ymomentum'] 96 97 115 98 116 # Use domain's startime if available … … 102 120 domain_starttime = None 103 121 104 105 122 # Build arguments and keyword arguments for use with caching or apply. 106 123 args = (filename,) 107 108 124 109 125 # FIXME (Ole): Caching this function will not work well … … 118 134 'boundary_polygon': boundary_polygon} 119 135 120 121 136 # Call underlying engine with or without caching 122 137 if use_cache is True: … … 133 148 compression=False, 134 149 verbose=verbose) 135 136 150 else: 137 151 f, starttime = apply(_file_function, 138 152 args, kwargs) 139 140 153 141 154 #FIXME (Ole): Pass cache arguments, such as compression, in some sort of … … 155 168 156 169 if verbose: print msg 170 157 171 domain.set_starttime(starttime) #Modifying model time 172 158 173 if verbose: print 'Domain starttime is now set to %f'\ 159 174 %domain.starttime 160 161 175 return f 162 176 163 177 164 178 ## 179 # @brief ?? 180 # @param filename Name of .sww or .tms file. 181 # @param domain Associated domain object. 182 # @param quantities Name of quantity to be interpolated or a list of names. 183 # @param interpolation_points List of absolute UTM coordinates for points 184 # (N x 2) or geospatial object or 185 # points file name at which values are sought. 186 # @param time_thinning 187 # @param verbose True if this function is to be verbose. 188 # @param use_cache True means that caching of intermediate result is attempted. 189 # @param boundary_polygon 165 190 def _file_function(filename, 166 191 quantities=None, … … 174 199 See file_function for documentatiton 175 200 """ 176 177 201 178 202 assert type(filename) == type(''),\ … … 182 206 fid = open(filename) 183 207 except Exception, e: 184 msg = 'File "%s" could not be opened: Error="%s"'\ 185 %(filename, e) 208 msg = 'File "%s" could not be opened: Error="%s"' % (filename, e) 186 209 raise msg 187 210 211 # read first line of file, guess file type 188 212 line = fid.readline() 189 213 fid.close() 190 191 214 192 215 if line[:3] == 'CDF': … … 199 222 boundary_polygon=boundary_polygon) 200 223 else: 201 # FIXME (Ole): Could add csv file here to address Ted Rigby's suggestion about reading hydrographs. 224 # FIXME (Ole): Could add csv file here to address Ted Rigby's 225 # suggestion about reading hydrographs. 202 226 # This may also deal with the gist of ticket:289 203 227 raise 'Must be a NetCDF File' 204 228 205 229 206 230 ## 231 # @brief ?? 232 # @param filename Name of .sww or .tms file. 233 # @param quantity_names Name of quantity to be interpolated or a list of names. 234 # @param interpolation_points List of absolute UTM coordinates for points 235 # (N x 2) or geospatial object or 236 # points file name at which values are sought. 237 # @param domain_starttime Start time from domain object. 238 # @param time_thinning ?? 239 # @param verbose True if this function is to be verbose. 240 # @param boundary_polygon ?? 241 # @return A callable object. 207 242 def get_netcdf_file_function(filename, 208 243 quantity_names=None, … … 222 257 223 258 See Interpolation function for further documetation 224 225 259 """ 226 227 260 228 261 # FIXME: Check that model origin is the same as file's origin 229 262 # (both in UTM coordinates) … … 232 265 # Take this code from e.g. dem2pts in data_manager.py 233 266 # FIXME: Use geo_reference to read and write xllcorner... 234 235 267 236 268 import time, calendar, types 237 269 from anuga.config import time_format 238 from Scientific.IO.NetCDF import NetCDFFile239 270 from Numeric import array, zeros, Float, alltrue, concatenate, reshape 240 271 241 272 # Open NetCDF file 242 273 if verbose: print 'Reading', filename 274 243 275 fid = NetCDFFile(filename, 'r') 244 276 … … 249 281 msg = 'No quantities are specified in file_function' 250 282 raise Exception, msg 251 252 283 253 284 if interpolation_points is not None: 254 285 interpolation_points = ensure_absolute(interpolation_points) 255 msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1]286 msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1] 256 287 assert interpolation_points.shape[1] == 2, msg 257 258 288 259 289 # Now assert that requested quantitites (and the independent ones) … … 266 296 if len(missing) > 0: 267 297 msg = 'Quantities %s could not be found in file %s'\ 268 % (str(missing), filename)298 % (str(missing), filename) 269 299 fid.close() 270 300 raise Exception, msg … … 346 376 else: 347 377 gauge_neighbour_id.append(-1) 348 if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 and boundary_id[0]==0: 378 if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \ 379 and boundary_id[0]==0: 349 380 gauge_neighbour_id.append(0) 350 381 else: 351 382 gauge_neighbour_id.append(-1) 352 383 gauge_neighbour_id=ensure_numeric(gauge_neighbour_id) 353 if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1: 384 385 if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \ 386 != len(temp)-1: 354 387 msg='incorrect number of segments' 355 388 raise msg … … 365 398 interpolation_points[:,0] -= xllcorner 366 399 interpolation_points[:,1] -= yllcorner 367 368 400 else: 369 401 gauge_neighbour_id=None 370 402 371 403 if domain_starttime is not None: 372 373 404 # If domain_startime is *later* than starttime, 374 405 # move time back - relative to domain's time … … 414 445 # Return Interpolation_function instance as well as 415 446 # starttime for use to possible modify that of domain 416 return Interpolation_function(time, 417 quantities, 418 quantity_names, 419 vertex_coordinates, 420 triangles, 421 interpolation_points, 422 time_thinning=time_thinning, 423 verbose=verbose,gauge_neighbour_id=gauge_neighbour_id), starttime 447 return (Interpolation_function(time, 448 quantities, 449 quantity_names, 450 vertex_coordinates, 451 triangles, 452 interpolation_points, 453 time_thinning=time_thinning, 454 verbose=verbose, 455 gauge_neighbour_id=gauge_neighbour_id), 456 starttime) 424 457 425 458 # NOTE (Ole): Caching Interpolation function is too slow as … … 427 460 428 461 429 430 462 ## 463 # @brief Replace multiple substrings in a string. 464 # @param text The string to operate on. 465 # @param dictionary A dict containing replacements, key->value. 466 # @return The new string. 431 467 def multiple_replace(text, dictionary): 432 468 """Multiple replace of words in text … … 436 472 437 473 Python Cookbook 3.14 page 88 and page 90 474 http://code.activestate.com/recipes/81330/ 438 475 """ 439 476 … … 450 487 451 488 452 453 454 def apply_expression_to_dictionary(expression, dictionary):#dictionary): 489 ## 490 # @brief Apply arbitrary expressions to the values of a dict. 491 # @param expression A string expression to apply. 492 # @param dictionary The dictionary to apply the expression to. 493 def apply_expression_to_dictionary(expression, dictionary): 455 494 """Apply arbitrary expression to values of dictionary 456 495 457 496 Given an expression in terms of the keys, replace key by the 458 497 corresponding values and evaluate. 459 460 498 461 499 expression: Arbitrary, e.g. arithmetric, expression relating keys … … 467 505 expression e.g. by overloading 468 506 469 due to a limitation with Numeric, this can not evaluate 0/0507 Due to a limitation with Numeric, this can not evaluate 0/0 470 508 In general, the user can fix by adding 1e-30 to the numerator. 471 509 SciPy core can handle this situation. … … 481 519 D = {} 482 520 for key in dictionary: 483 D[key] = 'dictionary["%s"]' % key521 D[key] = 'dictionary["%s"]' % key 484 522 485 523 #Perform substitution of variables … … 490 528 return eval(expression) 491 529 except NameError, e: 492 msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)530 msg = 'Expression "%s" could not be evaluated: %s' % (expression, e) 493 531 raise NameError, msg 494 532 except ValueError, e: 495 msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)533 msg = 'Expression "%s" could not be evaluated: %s' % (expression, e) 496 534 raise ValueError, msg 497 498 535 536 537 ## 538 # @brief Format a float into a string. 539 # @param value Float value to format. 540 # @param format The format to use (%.2f is default). 541 # @return The formatted float as a string. 499 542 def get_textual_float(value, format = '%.2f'): 500 543 """Get textual representation of floating point numbers … … 521 564 raise 'Illegal input to get_textual_float:', value 522 565 else: 523 return format % float(value)524 525 526 527 # ###################################528 #### OBSOLETE STUFF529 530 566 return format % float(value) 567 568 569 ################################################################################# 570 # OBSOLETE STUFF 571 ################################################################################# 572 573 # @note TEMP 531 574 def angle(v1, v2): 532 575 """Temporary Interface to new location""" … … 539 582 540 583 return NT.angle(v1, v2) 541 584 585 586 # @note TEMP 542 587 def anglediff(v0, v1): 543 588 """Temporary Interface to new location""" … … 552 597 553 598 599 # @note TEMP 554 600 def mean(x): 555 601 """Temporary Interface to new location""" … … 563 609 return NT.mean(x) 564 610 611 612 # @note TEMP 565 613 def point_on_line(*args, **kwargs): 566 614 """Temporary Interface to new location""" … … 571 619 572 620 return utilities.polygon.point_on_line(*args, **kwargs) 573 621 622 623 # @note TEMP 574 624 def inside_polygon(*args, **kwargs): 575 625 """Temporary Interface to new location""" … … 579 629 580 630 return utilities.polygon.inside_polygon(*args, **kwargs) 581 631 632 633 # @note TEMP 582 634 def outside_polygon(*args, **kwargs): 583 635 """Temporary Interface to new location""" … … 589 641 590 642 643 # @note TEMP 591 644 def separate_points_by_polygon(*args, **kwargs): 592 645 """Temporary Interface to new location""" 593 646 594 647 print 'separate_points_by_polygon has moved from util.py. ', 595 print 'Please use "from anuga.utilities.polygon import separate_points_by_polygon"' 648 print 'Please use "from anuga.utilities.polygon import ' \ 649 'separate_points_by_polygon"' 596 650 597 651 return utilities.polygon.separate_points_by_polygon(*args, **kwargs) 598 652 599 653 600 654 # @note TEMP 601 655 def read_polygon(*args, **kwargs): 602 656 """Temporary Interface to new location""" … … 608 662 609 663 664 # @note TEMP 610 665 def populate_polygon(*args, **kwargs): 611 666 """Temporary Interface to new location""" … … 616 671 return utilities.polygon.populate_polygon(*args, **kwargs) 617 672 618 ##################### end of obsolete stuff ? ############ 619 673 674 ################################################################################# 675 # End of obsolete stuff ? 676 ################################################################################# 677 678 # @note TEMP 620 679 def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', 621 680 verbose=False): 622 681 """Temporary Interface to new location""" 623 from anuga.shallow_water.data_manager import start_screen_catcher as dm_start_screen_catcher 682 from anuga.shallow_water.data_manager import start_screen_catcher \ 683 as dm_start_screen_catcher 624 684 625 685 print 'start_screen_catcher has moved from util.py. ', 626 print 'Please use "from anuga.shallow_water.data_manager import start_screen_catcher"' 627 628 return dm_start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', 629 verbose=False) 630 631 686 print 'Please use "from anuga.shallow_water.data_manager import ' \ 687 'start_screen_catcher"' 688 689 return dm_start_screen_catcher(dir_name, myid='', numprocs='', 690 extra_info='', verbose=False) 691 692 693 ## 694 # @brief Read a .sww file and plot the time series. 695 # @param swwfiles Dictionary of .sww files. 696 # @param gauge_filename Name of gauge file. 697 # @param production_dirs ?? 698 # @param report If True, write figures to report directory. 699 # @param reportname Name of generated report (if any). 700 # @param plot_quantity List containing quantities to plot. 701 # @param generate_fig If True, generate figures as well as CSV files. 702 # @param surface If True, then generate solution surface with 3d plot. 703 # @param time_min Beginning of user defined time range for plotting purposes. 704 # @param time_max End of user defined time range for plotting purposes. 705 # @param time_thinning ?? 706 # @param time_unit ?? 707 # @param title_on If True, export standard graphics with title. 708 # @param use_cache If True, use caching. 709 # @param verbose If True, this function is verbose. 632 710 def sww2timeseries(swwfiles, 633 711 gauge_filename, 634 712 production_dirs, 635 report = None, 636 reportname = None, 637 plot_quantity = None, 638 generate_fig = False, 639 surface = None, 640 time_min = None, 641 time_max = None, 642 time_thinning = 1, 643 time_unit = None, 644 title_on = None, 645 use_cache = False, 646 verbose = False): 647 713 report=None, 714 reportname=None, 715 plot_quantity=None, 716 generate_fig=False, 717 surface=None, 718 time_min=None, 719 time_max=None, 720 time_thinning=1, 721 time_unit=None, 722 title_on=None, 723 use_cache=False, 724 verbose=False): 648 725 """ Read sww file and plot the time series for the 649 726 prescribed quantities at defined gauge locations and … … 655 732 generating latex output. It will be part of 656 733 the directory name of file_loc (typically the timestamp). 657 Helps to differentiate latex files for different simulations658 for a particular scenario.734 Helps to differentiate latex files for different 735 simulations for a particular scenario. 659 736 - assume that all conserved quantities have been stored 660 737 - assume each sww file has been simulated with same timestep … … 668 745 production_dirs - A list of list, example {20061101_121212: '1 in 10000', 669 746 'boundaries': 'urs boundary'} 670 this will use the second part as the label and the first part671 as the ?747 this will use the second part as the label and the 748 first part as the ? 672 749 #FIXME: Is it a list or a dictionary 673 750 # This is probably obsolete by now … … 696 773 - default will be ['stage', 'speed', 'bearing'] 697 774 698 generate_fig - if True, generate figures as well as csv file s of quantities775 generate_fig - if True, generate figures as well as csv file 699 776 - if False, csv files created only 700 777 … … 713 790 714 791 715 716 792 Output: 717 793 … … 743 819 """ 744 820 745 msg = 'NOTE: A new function is available to create csv files from sww files called'746 msg += ' sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'747 msg += ' PLUS another new function to create graphs from csv files called'821 msg = 'NOTE: A new function is available to create csv files from sww ' 822 msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util' 823 msg += ' PLUS another new function to create graphs from csv files called ' 748 824 msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util' 749 825 print msg … … 764 840 use_cache, 765 841 verbose) 766 767 842 return k 768 843 844 845 ## 846 # @brief Read a .sww file and plot the time series. 847 # @param swwfiles Dictionary of .sww files. 848 # @param gauge_filename Name of gauge file. 849 # @param production_dirs ?? 850 # @param report If True, write figures to report directory. 851 # @param reportname Name of generated report (if any). 852 # @param plot_quantity List containing quantities to plot. 853 # @param generate_fig If True, generate figures as well as CSV files. 854 # @param surface If True, then generate solution surface with 3d plot. 855 # @param time_min Beginning of user defined time range for plotting purposes. 856 # @param time_max End of user defined time range for plotting purposes. 857 # @param time_thinning ?? 858 # @param time_unit ?? 859 # @param title_on If True, export standard graphics with title. 860 # @param use_cache If True, use caching. 861 # @param verbose If True, this function is verbose. 769 862 def _sww2timeseries(swwfiles, 770 863 gauge_filename, … … 783 876 verbose = False): 784 877 785 assert type(gauge_filename) == type(''),\ 786 'Gauge filename must be a string' 878 assert type(gauge_filename) == type(''), 'Gauge filename must be a string' 787 879 788 880 try: 789 881 fid = open(gauge_filename) 790 882 except Exception, e: 791 msg = 'File "%s" could not be opened: Error="%s"'\ 792 %(gauge_filename, e) 883 msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e) 793 884 raise msg 794 885 … … 799 890 plot_quantity = ['depth', 'speed'] 800 891 else: 801 assert type(plot_quantity) == list,\ 802 'plot_quantity must be a list' 892 assert type(plot_quantity) == list, 'plot_quantity must be a list' 803 893 check_list(plot_quantity) 804 894 … … 813 903 814 904 if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename 905 815 906 gauges, locations, elev = get_gauges_from_file(gauge_filename) 816 907 … … 825 916 826 917 for swwfile in swwfiles.keys(): 827 828 918 try: 829 919 fid = open(swwfile) 830 920 except Exception, e: 831 msg = 'File "%s" could not be opened: Error="%s"'\ 832 %(swwfile, e) 921 msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e) 833 922 raise msg 834 923 … … 841 930 #print 'label', label 842 931 leg_label.append(label) 843 844 845 932 846 933 f = file_function(swwfile, … … 872 959 file_loc.append(swwfile[:index+1]) 873 960 label_id.append(swwfiles[swwfile]) 874 875 961 876 962 f_list.append(f) … … 894 980 raise Exception, msg 895 981 896 if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures' 982 if verbose and len(gauge_index) > 0: 983 print 'Inputs OK - going to generate figures' 897 984 898 985 if len(gauge_index) <> 0: 899 texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface, 900 leg_label, f_list, gauges, locations, elev, gauge_index, 901 production_dirs, time_min, time_max, time_unit, 902 title_on, label_id, generate_fig, verbose) 986 texfile, elev_output = \ 987 generate_figures(plot_quantity, file_loc, report, reportname, 988 surface, leg_label, f_list, gauges, locations, 989 elev, gauge_index, production_dirs, time_min, 990 time_max, time_unit, title_on, label_id, 991 generate_fig, verbose) 903 992 else: 904 993 texfile = '' … … 906 995 907 996 return texfile, elev_output 908 997 998 999 ## 1000 # @brief Read gauge info from a file. 1001 # @param filename The name of the file to read. 1002 # @return A (gauges, gaugelocation, elev) tuple. 909 1003 def get_gauges_from_file(filename): 910 1004 """ Read in gauge information from file 911 1005 """ 1006 912 1007 from os import sep, getcwd, access, F_OK, mkdir 1008 1009 # Get data from the gauge file 913 1010 fid = open(filename) 914 1011 lines = fid.readlines() … … 923 1020 line11 = line1.split(',') 924 1021 925 if isinstance(line11[0], str) is True:1022 if isinstance(line11[0], str) is True: 926 1023 # We have found text in the first line 927 1024 east_index = None … … 929 1026 name_index = None 930 1027 elev_index = None 1028 931 1029 for i in range(len(line11)): 932 if line11[i].strip( '\n').strip('\r').strip(' ').lower() == 'easting':east_index = i933 if line11[i].strip( '\n').strip('\r').strip(' ').lower() == 'northing':north_index = i934 if line11[i].strip( '\n').strip('\r').strip(' ').lower() == 'name':name_index = i935 if line11[i].strip( '\n').strip('\r').strip(' ').lower() == 'elevation': elev_index = i1030 if line11[i].strip().lower() == 'easting': east_index = i 1031 if line11[i].strip().lower() == 'northing': north_index = i 1032 if line11[i].strip().lower() == 'name': name_index = i 1033 if line11[i].strip().lower() == 'elevation': elev_index = i 936 1034 937 1035 if east_index < len(line11) and north_index < len(line11): 938 1036 pass 939 1037 else: 940 msg = 'WARNING: %s does not contain correct header information' %(filename) 1038 msg = 'WARNING: %s does not contain correct header information' \ 1039 % filename 941 1040 msg += 'The header must be: easting, northing, name, elevation' 942 1041 raise Exception, msg … … 952 1051 # No header, assume that this is a simple easting, northing file 953 1052 954 msg = 'There was no header in file %s and the number of columns is %d' %(filename, len(line11)) 955 msg += '- I was assuming two columns corresponding to Easting and Northing' 1053 msg = 'There was no header in file %s and the number of columns is %d' \ 1054 % (filename, len(line11)) 1055 msg += '- was assuming two columns corresponding to Easting and Northing' 956 1056 assert len(line11) == 2, msg 957 1057 … … 977 1077 978 1078 979 1079 ## 1080 # @brief Check that input quantities in quantity list are legal. 1081 # @param quantity Quantity list to check. 1082 # @note Raises an exception of list is not legal. 980 1083 def check_list(quantity): 981 1084 """ Check that input quantities in quantity list are possible 982 1085 """ 1086 import sys 1087 from sets import Set as set 1088 983 1089 all_quantity = ['stage', 'depth', 'momentum', 'xmomentum', 984 1090 'ymomentum', 'speed', 'bearing', 'elevation'] 985 1091 986 987 988 import sys 989 #if not sys.version.startswith('2.4'): 990 # Backwards compatibility 991 # from sets import Set as set 992 993 from sets import Set as set 994 1092 # convert all quanitiy names to lowercase 995 1093 for i,j in enumerate(quantity): 996 1094 quantity[i] = quantity[i].lower() 1095 1096 # check that all names in 'quantity' appear in 'all_quantity' 997 1097 p = list(set(quantity).difference(set(all_quantity))) 998 if len(p) <>0:1098 if len(p) != 0: 999 1099 msg = 'Quantities %s do not exist - please try again' %p 1000 1100 raise Exception, msg 1001 1002 return 1003 1101 1102 1103 ## 1104 # @brief Calculate velocity bearing from North. 1105 # @param uh ?? 1106 # @param vh ?? 1107 # @return The calculated bearing. 1004 1108 def calc_bearing(uh, vh): 1005 1109 """ Calculate velocity bearing from North … … 1014 1118 1015 1119 angle = degrees(atan(vh/(uh+1.e-15))) 1120 1016 1121 if (0 < angle < 90.0): 1017 1122 if vh > 0: … … 1029 1134 return bearing 1030 1135 1136 1137 ## 1138 # @brief Generate figures from quantities and gauges for each sww file. 1139 # @param plot_quantity ?? 1140 # @param file_loc ?? 1141 # @param report ?? 1142 # @param reportname ?? 1143 # @param surface ?? 1144 # @param leg_label ?? 1145 # @param f_list ?? 1146 # @param gauges ?? 1147 # @param locations ?? 1148 # @param elev ?? 1149 # @param gauge_index ?? 1150 # @param production_dirs ?? 1151 # @param time_min ?? 1152 # @param time_max ?? 1153 # @param time_unit ?? 1154 # @param title_on ?? 1155 # @param label_id ?? 1156 # @param generate_fig ?? 1157 # @param verbose?? 1158 # @return (texfile2, elev_output) 1031 1159 def generate_figures(plot_quantity, file_loc, report, reportname, surface, 1032 1160 leg_label, f_list, gauges, locations, elev, gauge_index, … … 1036 1164 each sww file 1037 1165 """ 1038 # from math import sqrt, atan, degrees1039 1166 from Numeric import ones, allclose, zeros, Float, ravel 1040 1167 from os import sep, altsep, getcwd, mkdir, access, F_OK, environ … … 1055 1182 label_id1 = label_id[0].replace(sep,'') 1056 1183 label_id2 = label_id1.replace('_','') 1057 texfile = texdir +reportname+'%s' %(label_id2)1058 texfile2 = reportname +'%s' %(label_id2)1184 texfile = texdir + reportname + '%s' % label_id2 1185 texfile2 = reportname + '%s' % label_id2 1059 1186 texfilename = texfile + '.tex' 1187 fid = open(texfilename, 'w') 1188 1060 1189 if verbose: print '\n Latex output printed to %s \n' %texfilename 1061 fid = open(texfilename, 'w')1062 1190 else: 1063 1191 texfile = texdir+reportname 1064 1192 texfile2 = reportname 1065 1193 texfilename = texfile + '.tex' 1194 fid = open(texfilename, 'w') 1195 1066 1196 if verbose: print '\n Latex output printed to %s \n' %texfilename 1067 fid = open(texfilename, 'w')1068 1197 else: 1069 1198 texfile = '' … … 1078 1207 n0 = int(n0) 1079 1208 m = len(locations) 1080 model_time = zeros((n0, m,p), Float)1081 stages = zeros((n0, m,p), Float)1082 elevations = zeros((n0, m,p), Float)1083 momenta = zeros((n0, m,p), Float)1084 xmom = zeros((n0, m,p), Float)1085 ymom = zeros((n0, m,p), Float)1086 speed = zeros((n0, m,p), Float)1087 bearings = zeros((n0, m,p), Float)1088 due_east = 90.0*ones((n0, 1), Float)1089 due_west = 270.0*ones((n0, 1), Float)1090 depths = zeros((n0, m,p), Float)1091 eastings = zeros((n0, m,p), Float)1209 model_time = zeros((n0, m, p), Float) 1210 stages = zeros((n0, m, p), Float) 1211 elevations = zeros((n0, m, p), Float) 1212 momenta = zeros((n0, m, p), Float) 1213 xmom = zeros((n0, m, p), Float) 1214 ymom = zeros((n0, m, p), Float) 1215 speed = zeros((n0, m, p), Float) 1216 bearings = zeros((n0, m, p), Float) 1217 due_east = 90.0*ones((n0, 1), Float) 1218 due_west = 270.0*ones((n0, 1), Float) 1219 depths = zeros((n0, m, p), Float) 1220 eastings = zeros((n0, m, p), Float) 1092 1221 min_stages = [] 1093 1222 max_stages = [] … … 1101 1230 min_speeds = [] 1102 1231 max_depths = [] 1103 model_time_plot3d = zeros((n0, m), Float)1104 stages_plot3d = zeros((n0, m), Float)1105 eastings_plot3d = zeros((n0, m),Float)1232 model_time_plot3d = zeros((n0, m), Float) 1233 stages_plot3d = zeros((n0, m), Float) 1234 eastings_plot3d = zeros((n0, m),Float) 1106 1235 if time_unit is 'mins': scale = 60.0 1107 1236 if time_unit is 'hours': scale = 3600.0 1237 1108 1238 ##### loop over each swwfile ##### 1109 1239 for j, f in enumerate(f_list): 1240 if verbose: print 'swwfile %d of %d' % (j, len(f_list)) 1241 1110 1242 starttime = f.starttime 1111 if verbose: print 'swwfile %d of %d' %(j, len(f_list)) 1112 comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv' 1243 comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv' 1113 1244 fid_compare = open(comparefile, 'w') 1114 file0 = file_loc[j] +'gauges_t0.csv'1245 file0 = file_loc[j] + 'gauges_t0.csv' 1115 1246 fid_0 = open(file0, 'w') 1247 1116 1248 ##### loop over each gauge ##### 1117 1249 for k in gauge_index: 1250 if verbose: print 'Gauge %d of %d' % (k, len(gauges)) 1251 1118 1252 g = gauges[k] 1119 if verbose: print 'Gauge %d of %d' %(k, len(gauges))1120 1253 min_stage = 10 1121 1254 max_stage = 0 … … 1126 1259 max_depth = 0 1127 1260 gaugeloc = str(locations[k]) 1128 thisfile = file_loc[j] +sep+'gauges_time_series'+'_'\1129 + gaugeloc+'.csv'1261 thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \ 1262 + gaugeloc + '.csv' 1130 1263 fid_out = open(thisfile, 'w') 1131 1264 s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' 1132 1265 fid_out.write(s) 1266 1133 1267 #### generate quantities ####### 1134 1268 for i, t in enumerate(f.get_time()): … … 1156 1290 thisgauge = gauges[k] 1157 1291 eastings[i,k,j] = thisgauge[0] 1158 s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' %(t, w, m, vel, z, uh, vh, bearing) 1292 s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \ 1293 % (t, w, m, vel, z, uh, vh, bearing) 1159 1294 fid_out.write(s) 1160 1295 if t == 0: 1161 s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)1296 s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w) 1162 1297 fid_0.write(s) 1163 1298 if t/60.0 <= 13920: tindex = i … … 1175 1310 1176 1311 1177 s = '%.2f, %.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, thisgauge[0], leg_label[j]) 1312 s = '%.2f, %.2f, %.2f, %.2f, %s\n' \ 1313 % (max_stage, min_stage, z, thisgauge[0], leg_label[j]) 1178 1314 fid_compare.write(s) 1179 1315 max_stages.append(max_stage) … … 1194 1330 eastings_plot3d[:,] = eastings[:,:,j] 1195 1331 1196 if surface is 1332 if surface is True: 1197 1333 print 'Printing surface figure' 1198 1334 for i in range(2): … … 1200 1336 ax = p3.Axes3D(fig) 1201 1337 if len(gauges) > 80: 1202 ax.plot_surface(model_time[:,:,j],eastings[:,:,j],stages[:,:,j]) 1338 ax.plot_surface(model_time[:,:,j], 1339 eastings[:,:,j], 1340 stages[:,:,j]) 1203 1341 else: 1204 #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j]) 1205 ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j])) 1342 ax.plot3D(ravel(eastings[:,:,j]), 1343 ravel(model_time[:,:,j]), 1344 ravel(stages[:,:,j])) 1206 1345 ax.set_xlabel('time') 1207 1346 ax.set_ylabel('x') … … 1209 1348 fig.add_axes(ax) 1210 1349 p1.show() 1211 surfacefig = 'solution_surface%s' % leg_label[j]1350 surfacefig = 'solution_surface%s' % leg_label[j] 1212 1351 p1.savefig(surfacefig) 1213 1352 p1.close() … … 1218 1357 if surface is True: 1219 1358 figure(11) 1220 plot(eastings[tindex,:,j], stages[tindex,:,j])1359 plot(eastings[tindex,:,j], stages[tindex,:,j]) 1221 1360 xlabel('x') 1222 1361 ylabel('stage') … … 1226 1365 elev_output = [] 1227 1366 if generate_fig is True: 1228 depth_axis = axis([starttime/scale, time_max/scale, -0.1, max(max_depths)*1.1]) 1229 stage_axis = axis([starttime/scale, time_max/scale, min(min_stages), max(max_stages)*1.1]) 1230 vel_axis = axis([starttime/scale, time_max/scale, min(min_speeds), max(max_speeds)*1.1]) 1231 mom_axis = axis([starttime/scale, time_max/scale, min(min_momentums), max(max_momentums)*1.1]) 1232 xmom_axis = axis([starttime/scale, time_max/scale, min(min_xmomentums), max(max_xmomentums)*1.1]) 1233 ymom_axis = axis([starttime/scale, time_max/scale, min(min_ymomentums), max(max_ymomentums)*1.1]) 1367 depth_axis = axis([starttime/scale, time_max/scale, -0.1, 1368 max(max_depths)*1.1]) 1369 stage_axis = axis([starttime/scale, time_max/scale, 1370 min(min_stages), max(max_stages)*1.1]) 1371 vel_axis = axis([starttime/scale, time_max/scale, 1372 min(min_speeds), max(max_speeds)*1.1]) 1373 mom_axis = axis([starttime/scale, time_max/scale, 1374 min(min_momentums), max(max_momentums)*1.1]) 1375 xmom_axis = axis([starttime/scale, time_max/scale, 1376 min(min_xmomentums), max(max_xmomentums)*1.1]) 1377 ymom_axis = axis([starttime/scale, time_max/scale, 1378 min(min_ymomentums), max(max_ymomentums)*1.1]) 1234 1379 cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] 1235 1380 nn = len(plot_quantity) … … 1244 1389 count1 = 0 1245 1390 if report == True and len(label_id) > 1: 1246 s = '\\begin{figure}[ht] \n \\centering \n \\begin{tabular}{cc} \n' 1391 s = '\\begin{figure}[ht] \n' \ 1392 '\\centering \n' \ 1393 '\\begin{tabular}{cc} \n' 1247 1394 fid.write(s) 1248 1395 if len(label_id) > 1: graphname_report = [] 1396 1249 1397 #### generate figures for each gauge #### 1250 1398 for j, f in enumerate(f_list): … … 1256 1404 word_quantity = '' 1257 1405 if report == True and len(label_id) == 1: 1258 s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n' 1406 s = '\\begin{figure}[hbt] \n' \ 1407 '\\centering \n' \ 1408 '\\begin{tabular}{cc} \n' 1259 1409 fid.write(s) 1260 1410 … … 1264 1414 figure(count, frameon = False) 1265 1415 if which_quantity == 'depth': 1266 plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) 1416 plot(model_time[0:n[j]-1,k,j], 1417 depths[0:n[j]-1,k,j], '-', c = cstr[j]) 1267 1418 units = 'm' 1268 1419 axis(depth_axis) 1269 1420 if which_quantity == 'stage': 1270 1421 if elevations[0,k,j] <= 0: 1271 plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j]) 1422 plot(model_time[0:n[j]-1,k,j], 1423 stages[0:n[j]-1,k,j], '-', c = cstr[j]) 1272 1424 axis(stage_axis) 1273 1425 else: 1274 plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) 1426 plot(model_time[0:n[j]-1,k,j], 1427 depths[0:n[j]-1,k,j], '-', c = cstr[j]) 1275 1428 #axis(depth_axis) 1276 1429 units = 'm' 1277 1430 if which_quantity == 'momentum': 1278 plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j]) 1431 plot(model_time[0:n[j]-1,k,j], 1432 momenta[0:n[j]-1,k,j], '-', c = cstr[j]) 1279 1433 axis(mom_axis) 1280 1434 units = 'm^2 / sec' 1281 1435 if which_quantity == 'xmomentum': 1282 plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j]) 1436 plot(model_time[0:n[j]-1,k,j], 1437 xmom[0:n[j]-1,k,j], '-', c = cstr[j]) 1283 1438 axis(xmom_axis) 1284 1439 units = 'm^2 / sec' 1285 1440 if which_quantity == 'ymomentum': 1286 plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j]) 1441 plot(model_time[0:n[j]-1,k,j], 1442 ymom[0:n[j]-1,k,j], '-', c = cstr[j]) 1287 1443 axis(ymom_axis) 1288 1444 units = 'm^2 / sec' 1289 1445 if which_quantity == 'speed': 1290 plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j]) 1446 plot(model_time[0:n[j]-1,k,j], 1447 speed[0:n[j]-1,k,j], '-', c = cstr[j]) 1291 1448 axis(vel_axis) 1292 1449 units = 'm / sec' 1293 1450 if which_quantity == 'bearing': 1294 plot(model_time[0:n[j]-1,k,j], bearings[0:n[j]-1,k,j], '-',1451 plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-', 1295 1452 model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 1296 1453 model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.') … … 1301 1458 if time_unit is 'mins': xlabel('time (mins)') 1302 1459 if time_unit is 'hours': xlabel('time (hours)') 1303 #if which_quantity == 'stage' and elevations[0:n[j]-1,k,j] > 0: 1460 #if which_quantity == 'stage' \ 1461 # and elevations[0:n[j]-1,k,j] > 0: 1304 1462 # ylabel('%s (%s)' %('depth', units)) 1305 1463 #else: … … 1321 1479 mkdir (figdir) 1322 1480 latex_file_loc = figdir.replace(sep,altsep) 1323 graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory 1324 graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file 1481 # storing files in production directory 1482 graphname_latex = '%sgauge%s%s' \ 1483 % (latex_file_loc, gaugeloc2, 1484 which_quantity) 1485 # giving location in latex output file 1486 graphname_report_input = '%sgauge%s%s' % \ 1487 ('..' + altsep + 1488 'report_figures' + altsep, 1489 gaugeloc2, which_quantity) 1325 1490 graphname_report.append(graphname_report_input) 1326 1491 1327 savefig(graphname_latex) # save figures in production directory for report generation 1492 # save figures in production directory for report 1493 savefig(graphname_latex) 1328 1494 1329 1495 if report == True: 1330 1331 figdir = getcwd()+sep+'report_figures'+sep 1332 if access(figdir,F_OK) == 0 : 1333 mkdir (figdir) 1496 figdir = getcwd() + sep + 'report_figures' + sep 1497 if access(figdir,F_OK) == 0: 1498 mkdir(figdir) 1334 1499 latex_file_loc = figdir.replace(sep,altsep) 1335 1500 1336 1501 if len(label_id) == 1: 1337 graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory 1338 graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file 1339 s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png') 1502 # storing files in production directory 1503 graphname_latex = '%sgauge%s%s%s' % \ 1504 (latex_file_loc, gaugeloc2, 1505 which_quantity, label_id2) 1506 # giving location in latex output file 1507 graphname_report = '%sgauge%s%s%s' % \ 1508 ('..' + altsep + 1509 'report_figures' + altsep, 1510 gaugeloc2, which_quantity, 1511 label_id2) 1512 s = '\includegraphics' \ 1513 '[width=0.49\linewidth, height=50mm]{%s%s}' % \ 1514 (graphname_report, '.png') 1340 1515 fid.write(s) 1341 1516 if where1 % 2 == 0: … … 1348 1523 1349 1524 if title_on == True: 1350 title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2)) 1351 #title('Gauge %s (MOST elevation %.2f, ANUGA elevation %.2f)' %(gaugeloc2, elevations[10,k,0], elevations[10,k,1] )) 1525 title('%s scenario: %s at %s gauge' % \ 1526 (label_id, which_quantity, gaugeloc2)) 1527 #title('Gauge %s (MOST elevation %.2f, ' \ 1528 # 'ANUGA elevation %.2f)' % \ 1529 # (gaugeloc2, elevations[10,k,0], 1530 # elevations[10,k,1])) 1352 1531 1353 1532 savefig(graphname) # save figures with sww file … … 1356 1535 for i in range(nn-1): 1357 1536 if nn > 2: 1358 if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: 1537 if plot_quantity[i] == 'stage' \ 1538 and elevations[0,k,j] > 0: 1359 1539 word_quantity += 'depth' + ', ' 1360 1540 else: 1361 1541 word_quantity += plot_quantity[i] + ', ' 1362 1542 else: 1363 if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: 1543 if plot_quantity[i] == 'stage' \ 1544 and elevations[0,k,j] > 0: 1364 1545 word_quantity += 'depth' + ', ' 1365 1546 else: … … 1370 1551 else: 1371 1552 word_quantity += ' and ' + plot_quantity[nn-1] 1372 caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' ')) 1553 caption = 'Time series for %s at %s location ' \ 1554 '(elevation %.2fm)' % \ 1555 (word_quantity, locations[k], elev[k]) 1373 1556 if elev[k] == 0.0: 1374 caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j]) 1557 caption = 'Time series for %s at %s location ' \ 1558 '(elevation %.2fm)' % \ 1559 (word_quantity, locations[k], 1560 elevations[0,k,j]) 1375 1561 east = gauges[0] 1376 1562 north = gauges[1] 1377 elev_output.append([locations[k],east,north,elevations[0,k,j]]) 1378 label = '%sgauge%s' %(label_id2, gaugeloc2) 1379 s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label) 1563 elev_output.append([locations[k], east, north, 1564 elevations[0,k,j]]) 1565 label = '%sgauge%s' % (label_id2, gaugeloc2) 1566 s = '\end{tabular} \n' \ 1567 '\\caption{%s} \n' \ 1568 '\label{fig:%s} \n' \ 1569 '\end{figure} \n \n' % (caption, label) 1380 1570 fid.write(s) 1381 1571 cc += 1 … … 1400 1590 for which_quantity in plot_quantity: 1401 1591 where1 += 1 1402 s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report[index], '.png') 1592 s = '\includegraphics' \ 1593 '[width=0.49\linewidth, height=50mm]{%s%s}' % \ 1594 (graphname_report[index], '.png') 1403 1595 index += 1 1404 1596 fid.write(s) … … 1411 1603 word_quantity += ' and ' + plot_quantity[nn-1] 1412 1604 label = 'gauge%s' %(gaugeloc2) 1413 caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) 1605 caption = 'Time series for %s at %s location ' \ 1606 '(elevation %.2fm)' % \ 1607 (word_quantity, locations[k], elev[k]) 1414 1608 if elev[k] == 0.0: 1415 caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j]) 1609 caption = 'Time series for %s at %s location ' \ 1610 '(elevation %.2fm)' % \ 1611 (word_quantity, locations[k], 1612 elevations[0,k,j]) 1416 1613 thisgauge = gauges[k] 1417 1614 east = thisgauge[0] 1418 1615 north = thisgauge[1] 1419 elev_output.append([locations[k],east,north,elevations[0,k,j]]) 1616 elev_output.append([locations[k], east, north, 1617 elevations[0,k,j]]) 1420 1618 1421 s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label) 1619 s = '\end{tabular} \n' \ 1620 '\\caption{%s} \n' \ 1621 '\label{fig:%s} \n' \ 1622 '\end{figure} \n \n' % (caption, label) 1422 1623 fid.write(s) 1423 1624 if float((k+1)/div - pp) == 0.: 1424 1625 fid.write('\\clearpage \n') 1425 1626 pp += 1 1426 1427 1627 #### finished generating figures ### 1428 1628 … … 1430 1630 1431 1631 return texfile2, elev_output 1632 1432 1633 1433 1634 # FIXME (DSG): Add unit test, make general, not just 2 files, 1434 1635 # but any number of files. 1636 ## 1637 # @brief ?? 1638 # @param dir_name ?? 1639 # @param filename1 ?? 1640 # @param filename2 ?? 1641 # @return ?? 1642 # @note TEMP 1435 1643 def copy_code_files(dir_name, filename1, filename2): 1436 1644 """Temporary Interface to new location""" … … 1445 1653 1446 1654 1655 ## 1656 # @brief Create a nested sub-directory path. 1657 # @param root_directory The base diretory path. 1658 # @param directories An iterable of sub-directory names. 1659 # @return The final joined directory path. 1660 # @note If each sub-directory doesn't exist, it will be created. 1447 1661 def add_directories(root_directory, directories): 1448 1662 """ 1449 Add the first directory in directories to root_directory. 1450 Then add the second 1451 directory to the first directory and so on. 1663 Add the first sub-directory in 'directories' to root_directory. 1664 Then add the second sub-directory to the accumulating path and so on. 1452 1665 1453 1666 Return the path of the final directory. 1454 1667 1455 This is handy for specifying and creating a directory 1456 where data will go. 1668 This is handy for specifying and creating a directory where data will go. 1457 1669 """ 1458 1670 dir = root_directory … … 1460 1672 dir = os.path.join(dir, new_dir) 1461 1673 if not access(dir,F_OK): 1462 mkdir 1674 mkdir(dir) 1463 1675 return dir 1464 1676 1465 def get_data_from_file(filename,separator_value = ','): 1677 1678 ## 1679 # @brief 1680 # @param filename 1681 # @param separator_value 1682 # @return 1683 # @note TEMP 1684 def get_data_from_file(filename, separator_value=','): 1466 1685 """Temporary Interface to new location""" 1467 1686 from anuga.shallow_water.data_manager import \ … … 1473 1692 return dm_get_data_from_file(filename,separator_value = ',') 1474 1693 1694 1695 ## 1696 # @brief 1697 # @param verbose 1698 # @param kwargs 1699 # @return 1700 # @note TEMP 1475 1701 def store_parameters(verbose=False,**kwargs): 1476 1702 """Temporary Interface to new location""" … … 1484 1710 return dm_store_parameters(verbose=False,**kwargs) 1485 1711 1712 1713 ## 1714 # @brief Remove vertices that are not associated with any triangle. 1715 # @param verts An iterable (or array) of points. 1716 # @param triangles An iterable of 3 element tuples. 1717 # @param number_of_full_nodes ?? 1718 # @return (verts, triangles) where 'verts' has been updated. 1486 1719 def remove_lone_verts(verts, triangles, number_of_full_nodes=None): 1487 """ 1488 Removes vertices that are not associated with any triangles. 1489 1490 verts is a list/array of points 1491 triangles is a list of 3 element tuples. 1492 Each tuple represents a triangle. 1493 1720 """Removes vertices that are not associated with any triangles. 1721 1722 verts is a list/array of points. 1723 triangles is a list of 3 element tuples. Each tuple represents a triangle. 1494 1724 number_of_full_nodes relate to parallelism when a mesh has an 1495 1725 extra layer of ghost points. 1496 1497 1726 """ 1727 1498 1728 verts = ensure_numeric(verts) 1499 1729 triangles = ensure_numeric(triangles) … … 1502 1732 1503 1733 # initialise the array to easily find the index of the first loner 1504 loners=arange(2*N, N, -1) # if N=3 [6,5,4] 1734 # ie, if N=3 -> [6,5,4] 1735 loners=arange(2*N, N, -1) 1505 1736 for t in triangles: 1506 1737 for vert in t: 1507 1738 loners[vert]= vert # all non-loners will have loners[i]=i 1508 #print loners1509 1739 1510 1740 lone_start = 2*N - max(loners) # The index of the first loner … … 1522 1752 # verts=take(verts,X) to Remove the loners from verts 1523 1753 # but I think it would use more memory 1524 new_i = 01754 new_i = lone_start # point at first loner - first 'shuffle down' target 1525 1755 for i in range(lone_start, N): 1526 if loners[i] >= N: 1527 # Loner! 1756 if loners[i] >= N: # [i] is a loner, leave alone 1528 1757 pass 1529 else: 1758 else: # a non-loner, move down 1530 1759 loners[i] = new_i 1531 1760 verts[new_i] = verts[i] … … 1564 1793 return xc 1565 1794 1795 # @note TEMP 1566 1796 def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]}, 1567 1797 output_dir='',
Note: See TracChangeset
for help on using the changeset viewer.