Ignore:
Timestamp:
Dec 12, 2008, 10:35:06 AM (15 years ago)
Author:
rwilson
Message:

Fixed bug in remove_lone_verts(), added test, did pep8 changes.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r6070 r6072  
    1 #12345678901234567890123456789012345678901234567890123456789012345678901234567890
    21"""This module contains various auxiliary function used by pyvolution.
    32
    … …  
    17521751        # verts=take(verts,X)  to Remove the loners from verts
    17531752        # but I think it would use more memory
    1754         new_i = lone_start      # point at first loner - first 'shuffle down' target
     1753        new_i = lone_start      # point at first loner - 'shuffle down' target
    17551754        for i in range(lone_start, N):
    17561755            if loners[i] >= N:  # [i] is a loner, leave alone
    … …  
    17691768    return verts, triangles
    17701769
    1771  
     1770
     1771##
     1772# @brief Compute centroid values from vertex values
     1773# @param x Values at vertices of triangular mesh.
     1774# @param triangles Nx3 integer array pointing to vertex information.
     1775# @return [N] array of centroid values.
    17721776def get_centroid_values(x, triangles):
    17731777    """Compute centroid values from vertex values
    … …  
    17781782    indices into x
    17791783    """
    1780 
    17811784       
    17821785    xc = zeros(triangles.shape[0], Float) # Space for centroid info
    … …  
    17901793        xc[k] = (x[i0] + x[i1] + x[i2])/3
    17911794
    1792 
    17931795    return xc
     1796
    17941797
    17951798# @note TEMP
    … …  
    18021805                                extra_plot_name='test' ):
    18031806
    1804     msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs ',
    1805     msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import csv2timeseries_graphs"'
    1806 
     1807    msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs '
     1808    msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import ' \
     1809           'csv2timeseries_graphs"'
    18071810    raise Exception, msg
     1811
    18081812    return csv2timeseries_graphs(directories_dic,
    18091813                                 output_dir,
    … …  
    18141818                                 assess_all_csv_files
    18151819                                 )
    1816    
     1820
     1821
     1822##
     1823# @brief Plot time series from CSV files.
     1824# @param directories_dic
     1825# @param output_dir
     1826# @param base_name
     1827# @param plot_numbers
     1828# @param quantities
     1829# @param extra_plot_name
     1830# @param assess_all_csv_files
     1831# @param create_latex
     1832# @param verbose
     1833# @note Assumes that 'elevation' is in the CSV file(s).
    18171834def csv2timeseries_graphs(directories_dic={},
    1818                             output_dir='',
    1819                             base_name=None,
    1820                             plot_numbers='',
    1821                             quantities=['stage'],
    1822                             extra_plot_name='',
    1823                             assess_all_csv_files=True,                           
    1824                             create_latex=False,
    1825                             verbose=False):
     1835                          output_dir='',
     1836                          base_name=None,
     1837                          plot_numbers='',
     1838                          quantities=['stage'],
     1839                          extra_plot_name='',
     1840                          assess_all_csv_files=True,
     1841                          create_latex=False,
     1842                          verbose=False):
    18261843                               
    18271844    """
    … …  
    18881905                                  verbose=True)   
    18891906       
    1890             This will produce one plot for each quantity (therefore 3) in the current directory,
    1891             each plot will have 2 lines on them. The first plot named 'new' will have the time
    1892             offseted by 20secs and the stage height adjusted by -0.1m
     1907            This will produce one plot for each quantity (therefore 3) in the
     1908            current directory, each plot will have 2 lines on them. The first
     1909            plot named 'new' will have the time offseted by 20secs and the stage
     1910            height adjusted by -0.1m
    18931911       
    18941912    Inputs:
    … …  
    18971915                         the time series) and the (value to add to
    18981916                         stage if needed). For example
    1899                         {dir1:['Anuga_ons',5000, 0],
    1900                          dir2:['b_emoth',5000,1.5],
    1901                          dir3:['b_ons',5000,1.5]}
     1917                         {dir1:['Anuga_ons',5000, 0],
     1918                          dir2:['b_emoth',5000,1.5],
     1919                          dir3:['b_ons',5000,1.5]}
    19021920                         Having multiple directories defined will plot them on
    1903                          one plot, therefore there will be 3 lines on each of these
    1904                          plot. If you only want one line per plot call csv2timeseries_graph
    1905                          separately for each directory, eg only have one directory in the
    1906                          'directories_dic' in each call.
     1921                         one plot, therefore there will be 3 lines on each of
     1922                         these plot. If you only want one line per plot call
     1923                         csv2timeseries_graph separately for each directory,
     1924                         eg only have one directory in the 'directories_dic' in
     1925                         each call.
    19071926                         
    1908         output_dir: directory for the plot outputs. Only important to define
    1909                     when you have more than one directory in your directories_dic,
    1910                     if you have not defined it and you have multiple directories in
    1911                     'directories_dic' there will be plots in each directory however only
    1912                     one directory will contain the complete plot/graphs.
    1913        
    1914         base_name: Is used a couple of times. 1) to find the csv files to be plotted
    1915                    if there is no 'plot_numbers' then csv files with 'base_name' are
    1916                    plotted and 2) in the title of the plots, the lenght of base_name is
    1917                    removed from the front of the filename to be used in the title.
     1927        output_dir: directory for the plot outputs. Only important to define when
     1928                    you have more than one directory in your directories_dic, if
     1929                    you have not defined it and you have multiple directories in
     1930                    'directories_dic' there will be plots in each directory,
     1931                    however only one directory will contain the complete
     1932                    plot/graphs.
     1933       
     1934        base_name: Is used a couple of times.
     1935                   1) to find the csv files to be plotted if there is no
     1936                      'plot_numbers' then csv files with 'base_name' are plotted
     1937                   2) in the title of the plots, the length of base_name is
     1938                      removed from the front of the filename to be used in the
     1939                      title.
    19181940                   This could be changed if needed.
    19191941                   Note is ignored if assess_all_csv_files=True
    … …  
    19211943        plot_numbers: a String list of numbers to plot. For example
    19221944                      [0-4,10,15-17] will read and attempt to plot
    1923                        the follow 0,1,2,3,4,10,15,16,17
    1924                        NOTE: if no plot numbers this will create
    1925                        one plot per quantity, per gauge
    1926         quantities: Will get available quantities from the header in the csv file.
    1927                     quantities must be one of these.
     1945                      the follow 0,1,2,3,4,10,15,16,17
     1946                      NOTE: if no plot numbers this will create one plot per
     1947                            quantity, per gauge
     1948
     1949        quantities: Will get available quantities from the header in the csv
     1950                    file.  Quantities must be one of these.
    19281951                    NOTE: ALL QUANTITY NAMES MUST BE lower case!
    19291952                   
    … …  
    19421965    OUTPUTS: saves the plots to
    19431966              <output_dir><base_name><plot_number><extra_plot_name>.png
    1944              
    1945      
    19461967    """
     1968
    19471969    try:
    19481970        import pylab
    … …  
    19721994    quantities_label['bearing'] = 'degrees (o)'
    19731995    quantities_label['elevation'] = 'elevation (m)'
    1974 
    19751996   
    19761997    if extra_plot_name != '':
    1977         extra_plot_name='_'+extra_plot_name
     1998        extra_plot_name = '_' + extra_plot_name
    19781999
    19792000    new_plot_numbers=[]
    … …  
    19832004        if '-' in num_string:
    19842005            start = int(num_string[:num_string.rfind('-')])
    1985             end = int(num_string[num_string.rfind('-')+1:])+1
     2006            end = int(num_string[num_string.rfind('-') + 1:]) + 1
    19862007            for x in range(start, end):
    19872008                new_plot_numbers.append(str(x))
    … …  
    19952016    if verbose: print 'Determining files to access for axes ranges \n'
    19962017   
    1997     #print directories_dic.keys(), base_name
    1998  
    19992018    for i,directory in enumerate(directories_dic.keys()):
    20002019        all_csv_filenames.append(get_all_files_with_extension(directory,
    2001                                   base_name,'.csv'))
     2020                                                              base_name, '.csv'))
    20022021
    20032022        filenames=[]
    20042023        if plot_numbers == '':
    20052024            list_filenames.append(get_all_files_with_extension(directory,
    2006                                   base_name,'.csv'))
     2025                                                               base_name,'.csv'))
    20072026        else:
    20082027            for number in new_plot_numbers:
    2009 #                print 'number!!!', base_name, number
    2010                 filenames.append(base_name+number)
    2011 #                print filenames
     2028                filenames.append(base_name + number)
    20122029            list_filenames.append(filenames)
    2013 
    2014 
    2015 
    2016     #print "list_filenames", list_filenames
    20172030
    20182031    #use all the files to get the values for the plot axis
    … …  
    20202033    min_start_time = 100000
    20212034   
    2022    
    20232035    if verbose: print 'Determining uniform axes \n'
     2036
    20242037    #this entire loop is to determine the min and max range for the
    20252038    #axes of the plots
    … …  
    20332046    max_quantity_value={}
    20342047
    2035    
    20362048    for i, directory in enumerate(directories_dic.keys()):
    2037         filename_quantity_value={}
    2038         if assess_all_csv_files==False:
     2049        filename_quantity_value = {}
     2050        if assess_all_csv_files == False:
    20392051            which_csv_to_assess = list_filenames[i]
    20402052        else:
    20412053            #gets list of filenames for directory "i"
    20422054            which_csv_to_assess = all_csv_filenames[i]
    2043 #        print'IN DIR', list_filenames[i]
    2044        
    2045 
    2046 
    20472055       
    20482056        for j, filename in enumerate(which_csv_to_assess):
    2049             quantity_value={}
    2050 
    2051             dir_filename=join(directory,filename)
    2052             attribute_dic, title_index_dic = csv2dict(dir_filename+
    2053                                                        '.csv')
     2057            quantity_value = {}
     2058
     2059            dir_filename = join(directory,filename)
     2060            attribute_dic, title_index_dic = csv2dict(dir_filename + '.csv')
    20542061            directory_start_time = directories_dic[directory][1]
    20552062            directory_add_tide = directories_dic[directory][2]
    20562063
    20572064            if verbose: print 'reading: %s.csv' %dir_filename
    2058 #            print 'keys',attribute_dic.keys()
     2065
    20592066            #add time to get values
    20602067            for k, quantity in enumerate(quantities):
    2061                 quantity_value[quantity] = [float(x) for x in attribute_dic[quantity]]
     2068                quantity_value[quantity] = [float(x) for
     2069                                                x in attribute_dic[quantity]]
    20622070
    20632071                #add tide to stage if provided
    20642072                if quantity == 'stage':
    2065                      quantity_value[quantity]=array(quantity_value[quantity])+directory_add_tide
     2073                     quantity_value[quantity] = array(quantity_value[quantity]) \
     2074                                                      + directory_add_tide
    20662075
    20672076                #condition to find max and mins for all the plots
    … …  
    20692078                # then compare to the other values to determine abs max and min
    20702079                if i==0 and j==0:
    2071 
    20722080                    min_quantity_value[quantity], \
    2073                     max_quantity_value[quantity] = get_min_max_values(quantity_value[quantity])
     2081                        max_quantity_value[quantity] = \
     2082                            get_min_max_values(quantity_value[quantity])
    20742083
    20752084                    if quantity != 'time':
    2076                         min_quantity_value[quantity] = min_quantity_value[quantity] *1.1
    2077                         max_quantity_value[quantity] = max_quantity_value[quantity] *1.1
    2078 
    2079                    
    2080 #                    print '1 min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename
     2085                        min_quantity_value[quantity] = \
     2086                            min_quantity_value[quantity] *1.1
     2087                        max_quantity_value[quantity] = \
     2088                            max_quantity_value[quantity] *1.1
    20812089                else:
    2082 #                    print 'min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename
    20832090                    min, max = get_min_max_values(quantity_value[quantity])
    2084 #                    print "MIN",min, max
    20852091               
    2086                     #min and max are multipled by "1+increase_axis" to get axes that are slighty bigger
    2087                     # than the max and mins so the plots look good.
     2092                    # min and max are multipled by "1+increase_axis" to get axes
     2093                    # that are slighty bigger than the max and mins
     2094                    # so the plots look good.
    20882095
    20892096                    increase_axis = (max-min)*0.05
    2090 #                    print quantity, "MIN MAX", max, min
    2091                     if min<=min_quantity_value[quantity]:
     2097                    if min <= min_quantity_value[quantity]:
    20922098                        if quantity == 'time':
    2093                             min_quantity_value[quantity]=min
     2099                            min_quantity_value[quantity] = min
    20942100                        else:
    20952101                            if round(min,2) == 0.00:
    2096                                 min_quantity_value[quantity]=-increase_axis
    2097 #                                min_quantity_value[quantity]=-2.
    2098                                 #min_quantity_value[quantity]= -max_quantity_value[quantity]*increase_axis
     2102                                min_quantity_value[quantity] = -increase_axis
     2103#                                min_quantity_value[quantity] = -2.
     2104                                #min_quantity_value[quantity] = \
     2105                                #    -max_quantity_value[quantity]*increase_axis
    20992106                            else:
    2100 #                                min_quantity_value[quantity]=min*(1+increase_axis)
     2107#                                min_quantity_value[quantity] = \
     2108#                                    min*(1+increase_axis)
    21012109                                min_quantity_value[quantity]=min-increase_axis
    2102 #                        print quantity, min_quantity_value[quantity]
    21032110                   
    2104                     if max>max_quantity_value[quantity]:
     2111                    if max > max_quantity_value[quantity]:
    21052112                        if quantity == 'time':
    2106                             max_quantity_value[quantity]=max
     2113                            max_quantity_value[quantity] = max
    21072114                        else:
    2108                             max_quantity_value[quantity]=max+increase_axis
     2115                            max_quantity_value[quantity] = max + increase_axis
    21092116#                            max_quantity_value[quantity]=max*(1+increase_axis)
    2110 #                        print quantity, max_quantity_value[quantity],increase_axis
    2111                
    2112 #                print 'min,maj',quantity, min_quantity_value[quantity],max_quantity_value[quantity]
    2113 
    2114            
    2115            
    2116            
     2117
    21172118            #set the time... ???
    21182119            if min_start_time > directory_start_time:
    … …  
    21202121            if max_start_time < directory_start_time:
    21212122                max_start_time = directory_start_time
    2122             #print 'start_time' , max_start_time, min_start_time
    21232123           
    21242124            filename_quantity_value[filename]=quantity_value
    … …  
    21302130   
    21312131    for i, quantity in enumerate(quantities):
    2132         quantities_axis[quantity] = (float(min_start_time)/float(seconds_in_minutes),
    2133                                         (float(max_quantity_value['time'])+float(max_start_time))\
    2134                                               /float(seconds_in_minutes),
    2135                                          min_quantity_value[quantity],
    2136                                          max_quantity_value[quantity])
     2132        quantities_axis[quantity] = (float(min_start_time) \
     2133                                         / float(seconds_in_minutes),
     2134                                     (float(max_quantity_value['time']) \
     2135                                          + float(max_start_time)) \
     2136                                              / float(seconds_in_minutes),
     2137                                     min_quantity_value[quantity],
     2138                                     max_quantity_value[quantity])
     2139
    21372140        if verbose and (quantity != 'time' and quantity != 'elevation'):
    2138             print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' %(quantity,
    2139                                quantities_axis[quantity][0],
    2140                                quantities_axis[quantity][1],
    2141                                quantities_label['time'],
    2142                                quantities_axis[quantity][2],
    2143                                quantities_axis[quantity][3],
    2144                                quantities_label[quantity])
    2145    
    2146         #print  quantities_axis[quantity]
     2141            print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' \
     2142                  % (quantity,
     2143                     quantities_axis[quantity][0],
     2144                     quantities_axis[quantity][1],
     2145                     quantities_label['time'],
     2146                     quantities_axis[quantity][2],
     2147                     quantities_axis[quantity][3],
     2148                     quantities_label[quantity])
    21472149
    21482150    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
    … …  
    21512153   
    21522154    i_max = len(directories_dic.keys())
    2153     legend_list_dic={}
    2154     legend_list =[]
     2155    legend_list_dic = {}
     2156    legend_list = []
    21552157    for i, directory in enumerate(directories_dic.keys()):
    2156        
    2157         if verbose: print'Plotting in %s %s' %(directory, new_plot_numbers)
    2158 #        print 'LIST',list_filenames
    2159         #FIXME THIS SORT IS VERY IMPORTANT, without it the assigned plot numbers may not work correctly
    2160         #there must be a better way
     2158        if verbose: print 'Plotting in %s %s' % (directory, new_plot_numbers)
     2159
     2160        # FIXME THIS SORT IS VERY IMPORTANT
     2161        # Without it the assigned plot numbers may not work correctly
     2162        # there must be a better way
    21612163        list_filenames[i].sort()
    21622164        for j, filename in enumerate(list_filenames[i]):
    2163 #            print'IN plot', list_filenames[i]
    2164            
    2165             if verbose: print'Starting %s' %filename 
     2165            if verbose: print 'Starting %s' % filename 
     2166
    21662167            directory_name = directories_dic[directory][0]
    21672168            directory_start_time = directories_dic[directory][1]
    21682169            directory_add_tide = directories_dic[directory][2]
    21692170           
    2170             #create an if about the start time and tide hieght
    2171             #if they don't exist
    2172             #print 'i %s,j %s, number %s, file %s' %(i,j,number,file)
    2173             attribute_dic, title_index_dic = csv2dict(directory+sep+filename+'.csv')
     2171            # create an if about the start time and tide height if don't exist
     2172            attribute_dic, title_index_dic = csv2dict(directory + sep
     2173                                                      + filename + '.csv')
    21742174            #get data from dict in to list
    21752175            #do maths to list by changing to array
    2176             t=(array(directory_quantity_value[directory][filename]['time'])+directory_start_time)/seconds_in_minutes
     2176            t = (array(directory_quantity_value[directory][filename]['time'])
     2177                     + directory_start_time) / seconds_in_minutes
    21772178
    21782179            #finds the maximum elevation, used only as a test
    … …  
    21872188                print "Note! Elevation changes in %s" %dir_filename
    21882189
    2189             # creates a dictionary with keys that is the filename and attributes are a list of
    2190             # lists containing 'directory_name' and 'elevation'. This is used to make the contents
    2191             # for the legends in the graphs, this is the name of the model and the elevation.
    2192             # All in this great one liner from DG. If the key 'filename' doesn't exist it creates the
    2193             # entry if the entry exist it appends to the key.
    2194            
    2195             legend_list_dic.setdefault(filename,[]).append([directory_name,round(max_ele,3)])
    2196 
    2197             # creates a LIST for the legend on the last iteration of the directories
    2198             # which is when "legend_list_dic" has been fully populated. Creates a list of strings
    2199             # which is used in the legend
     2190            # creates a dictionary with keys that is the filename and attributes
     2191            # are a list of lists containing 'directory_name' and 'elevation'.
     2192            # This is used to make the contents for the legends in the graphs,
     2193            # this is the name of the model and the elevation.  All in this
     2194            # great one liner from DG. If the key 'filename' doesn't exist it
     2195            # creates the entry if the entry exist it appends to the key.
     2196
     2197            legend_list_dic.setdefault(filename,[]) \
     2198                .append([directory_name, round(max_ele, 3)])
     2199
     2200            # creates a LIST for the legend on the last iteration of the
     2201            # directories which is when "legend_list_dic" has been fully
     2202            # populated. Creates a list of strings which is used in the legend
    22002203            # only runs on the last iteration for all the gauges(csv) files
    22012204            # empties the list before creating it
    2202             if i==i_max-1:
    2203                 legend_list=[]
    2204    
    2205                 #print 'DIC',legend_list_dic
     2205
     2206            if i == i_max - 1:
     2207                legend_list = []
     2208   
    22062209                for name_and_elevation in legend_list_dic[filename]:
    2207                     legend_list.append('%s (elevation = %sm)'%(name_and_elevation[0],name_and_elevation[1]))
     2210                    legend_list.append('%s (elevation = %sm)'\
     2211                                       % (name_and_elevation[0],
     2212                                          name_and_elevation[1]))
    22082213           
    2209             #print 'filename',filename, quantities
    22102214            #skip time and elevation so it is not plotted!
    22112215            for k, quantity in enumerate(quantities):
    22122216                if quantity != 'time' and quantity != 'elevation':
    2213                    
    2214                     num=int(k*100+j)
     2217                    num = int(k*100+j)
    22152218                    pylab.figure(num)
    22162219                    pylab.ylabel(quantities_label[quantity])
    2217                     pylab.plot(t, directory_quantity_value[directory][filename][quantity], c = cstr[i], linewidth=1)
     2220                    pylab.plot(t,
     2221                               directory_quantity_value[directory]\
     2222                                                       [filename][quantity],
     2223                               c = cstr[i], linewidth=1)
    22182224                    pylab.xlabel(quantities_label['time'])
    22192225                    pylab.axis(quantities_axis[quantity])
    22202226                    pylab.legend(legend_list,loc='upper right')
    22212227                   
    2222                     pylab.title('%s at %s gauge' %(quantity,filename[len(base_name):]))
     2228                    pylab.title('%s at %s gauge'
     2229                                % (quantity, filename[len(base_name):]))
     2230
    22232231                    if output_dir == '':
    2224                         figname = '%s%s%s_%s%s.png' %(directory,sep,
    2225                                             filename,
    2226                                             quantity,
    2227                                             extra_plot_name)
     2232                        figname = '%s%s%s_%s%s.png' \
     2233                                  % (directory, sep, filename, quantity,
     2234                                     extra_plot_name)
    22282235                    else:
    2229                         figname = '%s%s%s_%s%s.png' %(output_dir,sep,
    2230                                             filename,
    2231                                             quantity,
    2232                                             extra_plot_name)
     2236                        figname = '%s%s%s_%s%s.png' \
     2237                                  % (output_dir, sep, filename, quantity,
     2238                                     extra_plot_name)
     2239
    22332240                    if verbose: print 'saving figure here %s' %figname
     2241
    22342242                    pylab.savefig(figname)
    22352243           
    22362244    if verbose: print 'Closing all plots'
     2245
    22372246    pylab.close('all')
    22382247    del pylab
     2248
    22392249    if verbose: print 'Finished closing plots'
    22402250
     2251##
     2252# @brief Return min and max of an iterable.
     2253# @param list The iterable to return min & max of.
     2254# @return (min, max) of 'list'.
    22412255def get_min_max_values(list=None):
    22422256    """
    22432257    Returns the min and max of the list it was provided.
    22442258    """
     2259
    22452260    if list == None: print 'List must be provided'
    22462261       
    … …  
    22482263
    22492264
    2250 
     2265##
     2266# @brief Get runup around a point in a CSV file.
     2267# @param gauge_filename gauge file name.
     2268# @param sww_filename SWW file name.
     2269# @param runup_filename Name of file to report into.
     2270# @param size ??
     2271# @param verbose ??
    22512272def get_runup_data_for_locations_from_file(gauge_filename,
    22522273                                           sww_filename,
    … …  
    22542275                                           size=10,
    22552276                                           verbose=False):
    2256 
    22572277    """this will read a csv file with the header x,y. Then look in a square
    22582278    'size'x2 around this position for the 'max_inundaiton_height' in the
    2259     'sww_filename' and report the findings in the 'runup_filename
     2279    'sww_filename' and report the findings in the 'runup_filename'.
    22602280   
    22612281    WARNING: NO TESTS!
    … …  
    22662286                                                 csv2dict
    22672287                                                 
    2268     file = open(runup_filename,"w")
     2288    file = open(runup_filename, "w")
    22692289    file.write("easting,northing,runup \n ")
    22702290    file.close()
    … …  
    22792299    runup_locations=[]
    22802300    for i, x in enumerate(northing):
    2281 #        print 'easting,northing',i,easting[i],northing[i]
    22822301        poly = [[int(easting[i]+size),int(northing[i]+size)],
    22832302                [int(easting[i]+size),int(northing[i]-size)],
    … …  
    22862305       
    22872306        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
    2288                                               polygon=poly,
    2289                                               verbose=False)
     2307                                                  polygon=poly,
     2308                                                  verbose=False)
     2309
    22902310        #if no runup will return 0 instead of NONE
    22912311        if run_up==None: run_up=0
    … …  
    22962316       
    22972317        #writes to file
    2298         file = open(runup_filename,"a")
    2299         temp = '%s,%s,%s \n' %(x_y[0], x_y[1], run_up)
     2318        file = open(runup_filename, "a")
     2319        temp = '%s,%s,%s \n' % (x_y[0], x_y[1], run_up)
    23002320        file.write(temp)
    23012321        file.close()
    2302        
    2303 
     2322
     2323
     2324##
     2325# @brief ??
     2326# @param sww_file ??
     2327# @param gauge_file ??
     2328# @param out_name ??
     2329# @param quantities ??
     2330# @param verbose ??
     2331# @param use_cache ??
    23042332def sww2csv_gauges(sww_file,
    23052333                   gauge_file,
    23062334                   out_name='gauge_',
    2307                    quantities = ['stage', 'depth', 'elevation',
    2308                                  'xmomentum', 'ymomentum'],
     2335                   quantities=['stage', 'depth', 'elevation',
     2336                               'xmomentum', 'ymomentum'],
    23092337                   verbose=False,
    23102338                   use_cache = True):
    … …  
    23552383
    23562384    This is really returning speed, not velocity.
    2357      
    2358    
    23592385    """
    23602386   
    … …  
    23642390    import string
    23652391    from anuga.shallow_water.data_manager import get_all_swwfiles
     2392
    23662393#    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
    23672394    #print "points",points
    23682395
    2369     assert type(gauge_file) == type(''),\
    2370            'Gauge filename must be a string'
    2371            
    2372     assert type(out_name) == type(''),\
    2373            'Output filename prefix must be a string'
     2396    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
     2397    assert type(out_name) == type(''), 'Output filename prefix must be a string'
    23742398   
    23752399    try:
    2376     #    fid = open(gauge_file)
    23772400        point_reader = reader(file(gauge_file))
    2378 
    23792401    except Exception, e:
    2380         msg = 'File "%s" could not be opened: Error="%s"'\
    2381                   %(gauge_file, e)
     2402        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
    23822403        raise msg
     2404
    23832405    if verbose: print '\n Gauges obtained from: %s \n' %gauge_file
    2384    
    23852406   
    23862407    point_reader = reader(file(gauge_file))
    … …  
    23902411    #read point info from file
    23912412    for i,row in enumerate(point_reader):
    2392 #        print 'i',i,'row',row
    23932413        #read header and determine the column numbers to read correcty.
    23942414        if i==0:
    23952415            for j,value in enumerate(row):
    2396 #                print 'j',j,value, row
    23972416                if value.strip()=='easting':easting=j
    23982417                if value.strip()=='northing':northing=j
    … …  
    24002419                if value.strip()=='elevation':elevation=j
    24012420        else:
    2402 #            print i,'easting',easting,'northing',northing, row[easting]
    24032421            points.append([float(row[easting]),float(row[northing])])
    24042422            point_name.append(row[name])
    … …  
    24102428
    24112429    dir_name, base = os.path.split(sww_file)   
    2412     #print 'dirname',dir_name, base
     2430
    24132431    #need to get current directory so when path and file
    24142432    #are "joined" below the directory is correct
    … …  
    24192437        if verbose: print 'File %s exists' %(sww_file)
    24202438    else:
    2421         msg = 'File "%s" could not be opened: no read permission'\
    2422                %(sww_file)
     2439        msg = 'File "%s" could not be opened: no read permission' % sww_file
    24232440        raise msg
    24242441
    … …  
    24362453   
    24372454    for sww_file in sww_files:
    2438    
    2439 #        print 'sww_file',sww_file
    24402455        sww_file = join(dir_name, sww_file+'.sww')
    2441 #        print 'sww_file',sww_file, core_quantities
    2442        
    24432456        callable_sww = file_function(sww_file,
    2444                                  quantities=core_quantities,
    2445                                  interpolation_points=points_array,
    2446                                  verbose=verbose,
    2447                                  use_cache=use_cache)
    2448        
     2457                                     quantities=core_quantities,
     2458                                     interpolation_points=points_array,
     2459                                     verbose=verbose,
     2460                                     use_cache=use_cache)
    24492461    gauge_file = out_name
    24502462
    24512463    heading = [quantity for quantity in quantities]
    24522464    heading.insert(0,'time')
    2453 
    2454 #    print 'start time', callable_sww.starttime, heading, quantities
    24552465
    24562466    #create a list of csv writers for all the points and write header
    24572467    points_writer = []
    24582468    for i,point in enumerate(points):
    2459         #print 'gauge file:',dir_name+sep+'gauge_'+point_name[i]+'.csv'
    2460         points_writer.append(writer(file(dir_name+sep+gauge_file+point_name[i]+'.csv', "wb")))
     2469        points_writer.append(writer(file(dir_name + sep + gauge_file
     2470                                         + point_name[i] + '.csv', "wb")))
    24612471        points_writer[i].writerow(heading)
    24622472   
    … …  
    24642474
    24652475    for time in callable_sww.get_time():
    2466 
    24672476        for point_i, point in enumerate(points_array):
    24682477            #add domain starttime to relative time.
    2469             points_list = [time+callable_sww.starttime]
    2470 #            print'time',time,'point_i',point_i,point, points_array
     2478            points_list = [time + callable_sww.starttime]
    24712479            point_quantities = callable_sww(time,point_i)
    2472 #            print "quantities", point_quantities
    24732480           
    24742481            for quantity in quantities:
    2475                 if quantity==NAN:
    2476                     print 'quantity does not exist in' %callable_sww.get_name
     2482                if quantity == NAN:
     2483                    print 'quantity does not exist in' % callable_sww.get_name
    24772484                else:
    24782485                    if quantity == 'stage':
    … …  
    24892496                       
    24902497                    if quantity == 'depth':
    2491                         points_list.append(point_quantities[0] - point_quantities[1])
    2492                        
    2493                        
     2498                        points_list.append(point_quantities[0]
     2499                                           - point_quantities[1])
     2500
    24942501                    if quantity == 'momentum':
    2495                         momentum = sqrt(point_quantities[2]**2 +\
    2496                                         point_quantities[3]**2)
     2502                        momentum = sqrt(point_quantities[2]**2
     2503                                        + point_quantities[3]**2)
    24972504                        points_list.append(momentum)
    24982505                       
    … …  
    25032510                        else:
    25042511                            if point_quantities[2] < 1.0e6:
    2505                                 momentum = sqrt(point_quantities[2]**2 +\
    2506                                                 point_quantities[3]**2)
     2512                                momentum = sqrt(point_quantities[2]**2
     2513                                                + point_quantities[3]**2)
    25072514    #                            vel = momentum/depth             
    2508                                 vel = momentum/(point_quantities[0] - point_quantities[1])
     2515                                vel = momentum / (point_quantities[0]
     2516                                                  - point_quantities[1])
    25092517    #                            vel = momentum/(depth + 1.e-6/depth)
    25102518                            else:
    … …  
    25172525                        points_list.append(calc_bearing(point_quantities[2],
    25182526                                                        point_quantities[3]))
    2519                        
    2520                 #print 'list',points_list
    2521                
     2527
    25222528            points_writer[point_i].writerow(points_list)
    25232529       
    25242530
    2525 def greens_law(d1,d2,h1,verbose=False):
    2526     """
     2531##
     2532# @brief Get a wave height at a certain depth given wave height at another depth.
     2533# @param d1 The first depth.
     2534# @param d2 The second depth.
     2535# @param h1 Wave ampitude at d1
     2536# @param verbose True if this function is to be verbose.
     2537# @return The wave height at d2.
     2538def greens_law(d1, d2, h1, verbose=False):
     2539    """Green's Law
    25272540
    25282541    Green's Law allows an approximation of wave amplitude at
    … …  
    25472560
    25482561    if d1 <= 0.0:
    2549         msg = 'the first depth, d1 (%f), must be strictly positive' %(d1)
     2562        msg = 'the first depth, d1 (%f), must be strictly positive' % (d1)
    25502563        raise Exception(msg)
    25512564
    25522565    if d2 <= 0.0:
    2553         msg = 'the second depth, d2 (%f), must be strictly positive' %(d2)
     2566        msg = 'the second depth, d2 (%f), must be strictly positive' % (d2)
    25542567        raise Exception(msg)
    25552568   
    25562569    if h1 <= 0.0:
    2557         msg = 'the wave amplitude, h1 (%f), must be strictly positive' %(h1)
     2570        msg = 'the wave amplitude, h1 (%f), must be strictly positive' % (h1)
    25582571        raise Exception(msg)
    25592572   
    … …  
    25652578       
    25662579
     2580##
     2581# @brief Get the square-root of a value.
     2582# @param s The value to get the square-root of.
     2583# @return The square-root of 's'.
    25672584def square_root(s):
    25682585    return sqrt(s)
     2586
     2587
Note: See TracChangeset for help on using the changeset viewer.