Changeset 6070


Ignore:
Timestamp:
Dec 11, 2008, 4:32:08 PM (16 years ago)
Author:
rwilson
Message:

Fixed problem in remove_lone_verts() and added an extra test case.

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  
    12631263
    12641264        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):
    12671271        verts = [[0,0],[1,0],[0,1]]
    12681272        tris = [[0,1,2]]
    12691273        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):
    12751278        verts = [[0,0],[1,0],[0,1],[99,99]]
    12761279        tris = [[0,1,2]]
    12771280        new_verts, new_tris = remove_lone_verts(verts, tris)
    1278         assert new_verts == verts[0:3]
    1279         assert new_tris == tris
    1280        
    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):
    12821285        verts = [[99,99],[0,0],[1,0],[99,99],[0,1],[99,99]]
    12831286        tris = [[1,2,4]]
    12841287        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]])
    12881290     
    1289     def test_remove_lone_verts_c(self):
     1291    def test_remove_lone_verts_d(self):
    12901292        verts = [[0,0],[1,0],[99,99],[0,1]]
    12911293        tris = [[0,1,3]]
    12921294        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):
    12981299        verts = [[0,0],[1,0],[0,1],[99,99],[99,99],[99,99]]
    12991300        tris = [[0,1,2]]
    13001301        new_verts, new_tris = remove_lone_verts(verts, tris)
    1301         assert new_verts == verts[0:3]
    1302         assert new_tris == tris
     1302        self.failUnless(new_verts.tolist() == verts[0:3])
     1303        self.failUnless(new_tris.tolist() == tris)
    13031304     
    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]]
    13081308        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######
    13111315       
    13121316    def test_get_min_max_values(self):
     
    18291833#-------------------------------------------------------------
    18301834if __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')
    18331837#    runner = unittest.TextTestRunner(verbosity=2)
    18341838    runner = unittest.TextTestRunner(verbosity=1)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r6035 r6070  
     1#12345678901234567890123456789012345678901234567890123456789012345678901234567890
    12"""This module contains various auxiliary function used by pyvolution.
    23
     
    1516
    1617from anuga.utilities.numerical_tools import ensure_numeric
     18from Scientific.IO.NetCDF import NetCDFFile
    1719from Numeric import arange, choose, zeros, Float, array, allclose, take, compress
    1820   
     
    2022from math import sqrt, atan, degrees
    2123
    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
    2526from anuga.utilities.system_tools import get_revision_number
    2627from anuga.utilities.system_tools import store_version_info
    2728
    2829
     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.
    2943def file_function(filename,
    3044                  domain=None,
     
    6983    interpolation_points - list of absolute UTM coordinates for points (N x 2)
    7084    or geospatial object or points file name at which values are sought
    71    
     85
     86    time_thinning -
     87
     88    verbose -
     89
    7290    use_cache: True means that caching of intermediate result of
    7391               Interpolation_function is attempted
    7492
    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
    7798    """
    78 
    7999
    80100    # FIXME (OLE): Should check origin of domain against that of file
     
    90110        if verbose:
    91111            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'
    93113            print msg
    94 
    95114        quantities = ['stage', 'xmomentum', 'ymomentum']
    96            
    97115
    98116    # Use domain's startime if available
     
    102120        domain_starttime = None
    103121
    104 
    105122    # Build arguments and keyword arguments for use with caching or apply.
    106123    args = (filename,)
    107 
    108124
    109125    # FIXME (Ole): Caching this function will not work well
     
    118134              'boundary_polygon': boundary_polygon}
    119135
    120 
    121136    # Call underlying engine with or without caching
    122137    if use_cache is True:
     
    133148                             compression=False,                 
    134149                             verbose=verbose)
    135 
    136150    else:
    137151        f, starttime = apply(_file_function,
    138152                             args, kwargs)
    139 
    140153
    141154    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
     
    155168           
    156169            if verbose: print msg
     170
    157171            domain.set_starttime(starttime) #Modifying model time
     172
    158173            if verbose: print 'Domain starttime is now set to %f'\
    159174               %domain.starttime
    160 
    161175    return f
    162176
    163177
    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
    165190def _file_function(filename,
    166191                   quantities=None,
     
    174199    See file_function for documentatiton
    175200    """
    176    
    177201
    178202    assert type(filename) == type(''),\
     
    182206        fid = open(filename)
    183207    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)
    186209        raise msg
    187210
     211    # read first line of file, guess file type
    188212    line = fid.readline()
    189213    fid.close()
    190 
    191214
    192215    if line[:3] == 'CDF':
     
    199222                                        boundary_polygon=boundary_polygon)
    200223    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.
    202226        # This may also deal with the gist of ticket:289
    203227        raise 'Must be a NetCDF File'
    204228
    205229
    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.
    207242def get_netcdf_file_function(filename,
    208243                             quantity_names=None,
     
    222257
    223258    See Interpolation function for further documetation
    224    
    225259    """
    226    
    227    
     260
    228261    # FIXME: Check that model origin is the same as file's origin
    229262    # (both in UTM coordinates)
     
    232265    # Take this code from e.g. dem2pts in data_manager.py
    233266    # FIXME: Use geo_reference to read and write xllcorner...
    234        
    235267
    236268    import time, calendar, types
    237269    from anuga.config import time_format
    238     from Scientific.IO.NetCDF import NetCDFFile
    239270    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
    240271
    241272    # Open NetCDF file
    242273    if verbose: print 'Reading', filename
     274
    243275    fid = NetCDFFile(filename, 'r')
    244276
     
    249281        msg = 'No quantities are specified in file_function'
    250282        raise Exception, msg
    251 
    252283 
    253284    if interpolation_points is not None:
    254285        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]
    256287        assert interpolation_points.shape[1] == 2, msg
    257 
    258288
    259289    # Now assert that requested quantitites (and the independent ones)
     
    266296    if len(missing) > 0:
    267297        msg = 'Quantities %s could not be found in file %s'\
    268               %(str(missing), filename)
     298              % (str(missing), filename)
    269299        fid.close()
    270300        raise Exception, msg
     
    346376                else:
    347377                    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:
    349380                gauge_neighbour_id.append(0)
    350381            else:
    351382                gauge_neighbour_id.append(-1)
    352383            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:
    354387                msg='incorrect number of segments'
    355388                raise msg
     
    365398            interpolation_points[:,0] -= xllcorner
    366399            interpolation_points[:,1] -= yllcorner       
    367    
    368400    else:
    369401        gauge_neighbour_id=None
    370402       
    371403    if domain_starttime is not None:
    372 
    373404        # If domain_startime is *later* than starttime,
    374405        # move time back - relative to domain's time
     
    414445    # Return Interpolation_function instance as well as
    415446    # 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)
    424457
    425458    # NOTE (Ole): Caching Interpolation function is too slow as
     
    427460
    428461
    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.
    431467def multiple_replace(text, dictionary):
    432468    """Multiple replace of words in text
     
    436472
    437473    Python Cookbook 3.14 page 88 and page 90
     474    http://code.activestate.com/recipes/81330/
    438475    """
    439476
     
    450487
    451488
    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.
     493def apply_expression_to_dictionary(expression, dictionary):
    455494    """Apply arbitrary expression to values of dictionary
    456495
    457496    Given an expression in terms of the keys, replace key by the
    458497    corresponding values and evaluate.
    459    
    460498
    461499    expression: Arbitrary, e.g. arithmetric, expression relating keys
     
    467505                expression e.g. by overloading
    468506
    469     due to a limitation with Numeric, this can not evaluate 0/0
     507    Due to a limitation with Numeric, this can not evaluate 0/0
    470508    In general, the user can fix by adding 1e-30 to the numerator.
    471509    SciPy core can handle this situation.
     
    481519    D = {}
    482520    for key in dictionary:
    483         D[key] = 'dictionary["%s"]' %key
     521        D[key] = 'dictionary["%s"]' % key
    484522
    485523    #Perform substitution of variables   
     
    490528        return eval(expression)
    491529    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)
    493531        raise NameError, msg
    494532    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)
    496534        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.
    499542def get_textual_float(value, format = '%.2f'):
    500543    """Get textual representation of floating point numbers
     
    521564                raise 'Illegal input to get_textual_float:', value
    522565        else:
    523             return format %float(value)
    524 
    525 
    526 
    527 ####################################
    528 ####OBSOLETE STUFF
    529 
    530 
     566            return format % float(value)
     567
     568
     569#################################################################################
     570# OBSOLETE STUFF
     571#################################################################################
     572
     573# @note TEMP
    531574def angle(v1, v2):
    532575    """Temporary Interface to new location"""
     
    539582
    540583    return NT.angle(v1, v2)
    541    
     584
     585   
     586# @note TEMP
    542587def anglediff(v0, v1):
    543588    """Temporary Interface to new location"""
     
    552597
    553598   
     599# @note TEMP
    554600def mean(x):
    555601    """Temporary Interface to new location"""
     
    563609    return NT.mean(x)   
    564610
     611
     612# @note TEMP
    565613def point_on_line(*args, **kwargs):
    566614    """Temporary Interface to new location"""
     
    571619
    572620    return utilities.polygon.point_on_line(*args, **kwargs)     
    573    
     621
     622   
     623# @note TEMP
    574624def inside_polygon(*args, **kwargs):
    575625    """Temporary Interface to new location"""
     
    579629
    580630    return utilities.polygon.inside_polygon(*args, **kwargs)   
    581    
     631
     632   
     633# @note TEMP
    582634def outside_polygon(*args, **kwargs):
    583635    """Temporary Interface to new location"""
     
    589641
    590642
     643# @note TEMP
    591644def separate_points_by_polygon(*args, **kwargs):
    592645    """Temporary Interface to new location"""
    593646
    594647    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"'
    596650
    597651    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)   
    598652
    599653
    600 
     654# @note TEMP
    601655def read_polygon(*args, **kwargs):
    602656    """Temporary Interface to new location"""
     
    608662
    609663
     664# @note TEMP
    610665def populate_polygon(*args, **kwargs):
    611666    """Temporary Interface to new location"""
     
    616671    return utilities.polygon.populate_polygon(*args, **kwargs)   
    617672
    618 ##################### end of obsolete stuff ? ############
    619 
     673
     674#################################################################################
     675# End of obsolete stuff ?
     676#################################################################################
     677
     678# @note TEMP
    620679def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
    621680                         verbose=False):
    622681    """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
    624684
    625685    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.
    632710def sww2timeseries(swwfiles,
    633711                   gauge_filename,
    634712                   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):
    648725    """ Read sww file and plot the time series for the
    649726    prescribed quantities at defined gauge locations and
     
    655732                      generating latex output. It will be part of
    656733                      the directory name of file_loc (typically the timestamp).
    657                       Helps to differentiate latex files for different simulations
    658                       for a particular scenario. 
     734                      Helps to differentiate latex files for different
     735                      simulations for a particular scenario. 
    659736                    - assume that all conserved quantities have been stored
    660737                    - assume each sww file has been simulated with same timestep
     
    668745    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
    669746                                                'boundaries': 'urs boundary'}
    670                       this will use the second part as the label and the first part
    671                       as the ?
     747                      this will use the second part as the label and the
     748                      first part as the ?
    672749                      #FIXME: Is it a list or a dictionary
    673750                      # This is probably obsolete by now
     
    696773                    - default will be ['stage', 'speed', 'bearing']
    697774
    698     generate_fig     - if True, generate figures as well as csv files of quantities
     775    generate_fig     - if True, generate figures as well as csv file
    699776                     - if False, csv files created only
    700777                     
     
    713790
    714791
    715                      
    716792    Output:
    717793   
     
    743819    """
    744820
    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 '
    748824    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
    749825    print msg
     
    764840                        use_cache,
    765841                        verbose)
    766 
    767842    return k
    768843
     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.
    769862def _sww2timeseries(swwfiles,
    770863                    gauge_filename,
     
    783876                    verbose = False):   
    784877       
    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'
    787879   
    788880    try:
    789881        fid = open(gauge_filename)
    790882    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)
    793884        raise msg
    794885
     
    799890        plot_quantity = ['depth', 'speed']
    800891    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'
    803893        check_list(plot_quantity)
    804894
     
    813903   
    814904    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
     905
    815906    gauges, locations, elev = get_gauges_from_file(gauge_filename)
    816907
     
    825916
    826917    for swwfile in swwfiles.keys():
    827 
    828918        try:
    829919            fid = open(swwfile)
    830920        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)
    833922            raise msg
    834923
     
    841930        #print 'label', label
    842931        leg_label.append(label)
    843 
    844        
    845932
    846933        f = file_function(swwfile,
     
    872959        file_loc.append(swwfile[:index+1])
    873960        label_id.append(swwfiles[swwfile])
    874 
    875961       
    876962        f_list.append(f)
     
    894980            raise Exception, msg
    895981
    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'
    897984
    898985    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)
    903992    else:
    904993        texfile = ''
     
    906995
    907996    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.
    9091003def get_gauges_from_file(filename):
    9101004    """ Read in gauge information from file
    9111005    """
     1006
    9121007    from os import sep, getcwd, access, F_OK, mkdir
     1008
     1009    # Get data from the gauge file
    9131010    fid = open(filename)
    9141011    lines = fid.readlines()
     
    9231020    line11 = line1.split(',')
    9241021
    925     if isinstance(line11[0],str) is True:
     1022    if isinstance(line11[0], str) is True:
    9261023        # We have found text in the first line
    9271024        east_index = None
     
    9291026        name_index = None
    9301027        elev_index = None
     1028
    9311029        for i in range(len(line11)):
    932             if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'easting': east_index = i
    933             if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'northing': north_index = i
    934             if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'name': name_index = i
    935             if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'elevation': elev_index = i
     1030            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
    9361034
    9371035        if east_index < len(line11) and north_index < len(line11):
    9381036            pass
    9391037        else:
    940             msg = 'WARNING: %s does not contain correct header information' %(filename)
     1038            msg = 'WARNING: %s does not contain correct header information' \
     1039                  % filename
    9411040            msg += 'The header must be: easting, northing, name, elevation'
    9421041            raise Exception, msg
     
    9521051        # No header, assume that this is a simple easting, northing file
    9531052
    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'
    9561056        assert len(line11) == 2, msg
    9571057
     
    9771077
    9781078
    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.
    9801083def check_list(quantity):
    9811084    """ Check that input quantities in quantity list are possible
    9821085    """
     1086    import sys
     1087    from sets import Set as set
     1088
    9831089    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
    9841090                    'ymomentum', 'speed', 'bearing', 'elevation']
    9851091
    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
    9951093    for i,j in enumerate(quantity):
    9961094        quantity[i] = quantity[i].lower()
     1095
     1096    # check that all names in 'quantity' appear in 'all_quantity'
    9971097    p = list(set(quantity).difference(set(all_quantity)))
    998     if len(p) <> 0:
     1098    if len(p) != 0:
    9991099        msg = 'Quantities %s do not exist - please try again' %p
    10001100        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.
    10041108def calc_bearing(uh, vh):
    10051109    """ Calculate velocity bearing from North
     
    10141118       
    10151119    angle = degrees(atan(vh/(uh+1.e-15)))
     1120
    10161121    if (0 < angle < 90.0):
    10171122        if vh > 0:
     
    10291134    return bearing
    10301135
     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)
    10311159def generate_figures(plot_quantity, file_loc, report, reportname, surface,
    10321160                     leg_label, f_list, gauges, locations, elev, gauge_index,
     
    10361164    each sww file
    10371165    """
    1038 #    from math import sqrt, atan, degrees
    10391166    from Numeric import ones, allclose, zeros, Float, ravel
    10401167    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
     
    10551182            label_id1 = label_id[0].replace(sep,'')
    10561183            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
    10591186            texfilename = texfile + '.tex'
     1187            fid = open(texfilename, 'w')
     1188
    10601189            if verbose: print '\n Latex output printed to %s \n' %texfilename
    1061             fid = open(texfilename, 'w')
    10621190        else:
    10631191            texfile = texdir+reportname
    10641192            texfile2 = reportname
    10651193            texfilename = texfile + '.tex'
     1194            fid = open(texfilename, 'w')
     1195
    10661196            if verbose: print '\n Latex output printed to %s \n' %texfilename
    1067             fid = open(texfilename, 'w')
    10681197    else:
    10691198        texfile = ''
     
    10781207    n0 = int(n0)
    10791208    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)
    10921221    min_stages = []
    10931222    max_stages = []
     
    11011230    min_speeds = []   
    11021231    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)
    11061235    if time_unit is 'mins': scale = 60.0
    11071236    if time_unit is 'hours': scale = 3600.0
     1237
    11081238    ##### loop over each swwfile #####
    11091239    for j, f in enumerate(f_list):
     1240        if verbose: print 'swwfile %d of %d' % (j, len(f_list))
     1241
    11101242        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'
    11131244        fid_compare = open(comparefile, 'w')
    1114         file0 = file_loc[j]+'gauges_t0.csv'
     1245        file0 = file_loc[j] + 'gauges_t0.csv'
    11151246        fid_0 = open(file0, 'w')
     1247
    11161248        ##### loop over each gauge #####
    11171249        for k in gauge_index:
     1250            if verbose: print 'Gauge %d of %d' % (k, len(gauges))
     1251
    11181252            g = gauges[k]
    1119             if verbose: print 'Gauge %d of %d' %(k, len(gauges))
    11201253            min_stage = 10
    11211254            max_stage = 0
     
    11261259            max_depth = 0           
    11271260            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'
    11301263            fid_out = open(thisfile, 'w')
    11311264            s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
    11321265            fid_out.write(s)
     1266
    11331267            #### generate quantities #######
    11341268            for i, t in enumerate(f.get_time()):
     
    11561290                    thisgauge = gauges[k]
    11571291                    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)
    11591294                    fid_out.write(s)
    11601295                    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)
    11621297                        fid_0.write(s)
    11631298                    if t/60.0 <= 13920: tindex = i
     
    11751310                   
    11761311                   
    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])
    11781314            fid_compare.write(s)
    11791315            max_stages.append(max_stage)
     
    11941330        eastings_plot3d[:,] = eastings[:,:,j]
    11951331           
    1196         if surface is  True:
     1332        if surface is True:
    11971333            print 'Printing surface figure'
    11981334            for i in range(2):
     
    12001336                ax = p3.Axes3D(fig)
    12011337                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])
    12031341                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]))
    12061345                ax.set_xlabel('time')
    12071346                ax.set_ylabel('x')
     
    12091348                fig.add_axes(ax)
    12101349                p1.show()
    1211                 surfacefig = 'solution_surface%s' %leg_label[j]
     1350                surfacefig = 'solution_surface%s' % leg_label[j]
    12121351                p1.savefig(surfacefig)
    12131352                p1.close()
     
    12181357    if surface is True:
    12191358        figure(11)
    1220         plot(eastings[tindex,:,j],stages[tindex,:,j])
     1359        plot(eastings[tindex,:,j], stages[tindex,:,j])
    12211360        xlabel('x')
    12221361        ylabel('stage')
     
    12261365    elev_output = []
    12271366    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])
    12341379        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
    12351380        nn = len(plot_quantity)
     
    12441389            count1 = 0
    12451390            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'
    12471394                fid.write(s)
    12481395            if len(label_id) > 1: graphname_report = []
     1396
    12491397            #### generate figures for each gauge ####
    12501398            for j, f in enumerate(f_list):
     
    12561404                word_quantity = ''
    12571405                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'
    12591409                    fid.write(s)
    12601410                   
     
    12641414                    figure(count, frameon = False)
    12651415                    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])
    12671418                        units = 'm'
    12681419                        axis(depth_axis)
    12691420                    if which_quantity == 'stage':
    12701421                        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])
    12721424                            axis(stage_axis)
    12731425                        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])
    12751428                            #axis(depth_axis)                 
    12761429                        units = 'm'
    12771430                    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])
    12791433                        axis(mom_axis)
    12801434                        units = 'm^2 / sec'
    12811435                    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])
    12831438                        axis(xmom_axis)
    12841439                        units = 'm^2 / sec'
    12851440                    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])
    12871443                        axis(ymom_axis)
    12881444                        units = 'm^2 / sec'
    12891445                    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])
    12911448                        axis(vel_axis)
    12921449                        units = 'm / sec'
    12931450                    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],'-',
    12951452                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.',
    12961453                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
     
    13011458                    if time_unit is 'mins': xlabel('time (mins)')
    13021459                    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:
    13041462                    #    ylabel('%s (%s)' %('depth', units))
    13051463                    #else:
     
    13211479                            mkdir (figdir)
    13221480                        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)
    13251490                        graphname_report.append(graphname_report_input)
    13261491                       
    1327                         savefig(graphname_latex) # save figures in production directory for report generation
     1492                        # save figures in production directory for report
     1493                        savefig(graphname_latex)
    13281494
    13291495                    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)
    13341499                        latex_file_loc = figdir.replace(sep,altsep)   
    13351500
    13361501                        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')
    13401515                            fid.write(s)
    13411516                            if where1 % 2 == 0:
     
    13481523                   
    13491524                    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]))
    13521531
    13531532                    savefig(graphname) # save figures with sww file
     
    13561535                    for i in range(nn-1):
    13571536                        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:
    13591539                                word_quantity += 'depth' + ', '
    13601540                            else:
    13611541                                word_quantity += plot_quantity[i] + ', '
    13621542                        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:
    13641545                                word_quantity += 'depth' + ', '
    13651546                            else:
     
    13701551                    else:
    13711552                        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])
    13731556                    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])
    13751561                        east = gauges[0]
    13761562                        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)
    13801570                    fid.write(s)
    13811571                    cc += 1
     
    14001590                    for which_quantity in plot_quantity:
    14011591                        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')
    14031595                        index += 1
    14041596                        fid.write(s)
     
    14111603                word_quantity += ' and ' + plot_quantity[nn-1]           
    14121604                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])
    14141608                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])
    14161613                        thisgauge = gauges[k]
    14171614                        east = thisgauge[0]
    14181615                        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]])
    14201618                       
    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)
    14221623                fid.write(s)
    14231624                if float((k+1)/div - pp) == 0.:
    14241625                    fid.write('\\clearpage \n')
    14251626                    pp += 1
    1426                
    14271627                #### finished generating figures ###
    14281628
     
    14301630       
    14311631    return texfile2, elev_output
     1632
    14321633
    14331634# FIXME (DSG): Add unit test, make general, not just 2 files,
    14341635# but any number of files.
     1636##
     1637# @brief ??
     1638# @param dir_name ??
     1639# @param filename1 ??
     1640# @param filename2 ??
     1641# @return ??
     1642# @note TEMP
    14351643def copy_code_files(dir_name, filename1, filename2):
    14361644    """Temporary Interface to new location"""
     
    14451653
    14461654
     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.
    14471661def add_directories(root_directory, directories):
    14481662    """
    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.
    14521665
    14531666    Return the path of the final directory.
    14541667
    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.
    14571669    """
    14581670    dir = root_directory
     
    14601672        dir = os.path.join(dir, new_dir)
    14611673        if not access(dir,F_OK):
    1462             mkdir (dir)
     1674            mkdir(dir)
    14631675    return dir
    14641676
    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
     1684def get_data_from_file(filename, separator_value=','):
    14661685    """Temporary Interface to new location"""
    14671686    from anuga.shallow_water.data_manager import \
     
    14731692    return dm_get_data_from_file(filename,separator_value = ',')
    14741693
     1694
     1695##
     1696# @brief
     1697# @param verbose
     1698# @param kwargs
     1699# @return
     1700# @note TEMP
    14751701def store_parameters(verbose=False,**kwargs):
    14761702    """Temporary Interface to new location"""
     
    14841710    return dm_store_parameters(verbose=False,**kwargs)
    14851711
     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.
    14861719def 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.
    14941724    number_of_full_nodes relate to parallelism when a mesh has an
    14951725        extra layer of ghost points.
    1496        
    14971726    """
     1727
    14981728    verts = ensure_numeric(verts)
    14991729    triangles = ensure_numeric(triangles)
     
    15021732   
    15031733    # 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)
    15051736    for t in triangles:
    15061737        for vert in t:
    15071738            loners[vert]= vert # all non-loners will have loners[i]=i
    1508     #print loners
    15091739
    15101740    lone_start = 2*N - max(loners) # The index of the first loner
     
    15221752        # verts=take(verts,X)  to Remove the loners from verts
    15231753        # but I think it would use more memory
    1524         new_i = 0
     1754        new_i = lone_start      # point at first loner - first 'shuffle down' target
    15251755        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
    15281757                pass
    1529             else:
     1758            else:               # a non-loner, move down
    15301759                loners[i] = new_i
    15311760                verts[new_i] = verts[i]
     
    15641793    return xc
    15651794
     1795# @note TEMP
    15661796def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]},
    15671797                                output_dir='',
Note: See TracChangeset for help on using the changeset viewer.