Changeset 7690

Timestamp:

Apr 21, 2010, 1:48:09 PM (14 years ago)

Author:

hudson

Message:

Check for complex polygons raises an exception.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/file_function.py (modified) (2 diffs)
• abstract_2d_finite_volumes/gauge.py (modified) (3 diffs)
• abstract_2d_finite_volumes/test_gauge.py (modified) (21 diffs)
• abstract_2d_finite_volumes/util.py (modified) (1 diff)
• fit_interpolate/general_fit_interpolate.py (modified) (1 diff)
• utilities/mainland_only.csv (modified) (7 diffs)
• utilities/numerical_tools.py (modified) (1 diff)
• utilities/polygon.py (modified) (2 diffs)
• utilities/test_polygon.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/file_function.py

@@ -17 +17 @@
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import ensure_numeric
+
+import anuga.utilities.log as log
+
 
 
@@ -192 +195 @@
                    verbose=False,
                    boundary_polygon=None,
-                                   output_centroids=False):
+                   output_centroids=False):
     """Internal function
     

anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

@@ -12 +12 @@
 
 from anuga.geospatial_data.geospatial_data import ensure_absolute
-from util import check_list, generate_figures
+from util import check_list
 from file_function import file_function
 
@@ -548 +548 @@
     if len(gauge_index) <> 0:
         texfile, elev_output = \
-            generate_figures(plot_quantity, file_loc, report, reportname,
+            _generate_figures(plot_quantity, file_loc, report, reportname,
                              surface, leg_label, f_list, gauges, locations,
                              elev, gauge_index, production_dirs, time_min,
@@ -638 +638 @@
 
     return gauges, gaugelocation, elev
+
+    
+##
+# @brief Generate figures from quantities and gauges for each sww file.
+#        This is a legacy function, used only by sww2timeseries
+# @param plot_quantity  ??
+# @param file_loc ??
+# @param report ??
+# @param reportname ??
+# @param surface ??
+# @param leg_label ??
+# @param f_list ??
+# @param gauges ??
+# @param locations ??
+# @param elev ??
+# @param gauge_index ??
+# @param production_dirs ??
+# @param time_min ??
+# @param time_max ??
+# @param time_unit ??
+# @param title_on ??
+# @param label_id ??
+# @param generate_fig ??
+# @param verbose??
+# @return (texfile2, elev_output)
+def _generate_figures(plot_quantity, file_loc, report, reportname, surface,
+                     leg_label, f_list, gauges, locations, elev, gauge_index,
+                     production_dirs, time_min, time_max, time_unit,
+                     title_on, label_id, generate_fig, verbose):
+    """ Generate figures based on required quantities and gauges for
+    each sww file
+    """
+    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
+
+    if generate_fig is True:
+        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
+             xlabel, ylabel, title, close, subplot
+    
+        if surface is True:
+            import pylab as p1
+            import mpl3d.mplot3d as p3
+        
+    if report == True:    
+        texdir = getcwd()+sep+'report'+sep
+        if access(texdir,F_OK) == 0:
+            mkdir (texdir)
+        if len(label_id) == 1:
+            label_id1 = label_id[0].replace(sep,'')
+            label_id2 = label_id1.replace('_','')
+            texfile = texdir + reportname + '%s' % label_id2
+            texfile2 = reportname + '%s' % label_id2
+            texfilename = texfile + '.tex'
+            fid = open(texfilename, 'w')
+
+            if verbose: log.critical('Latex output printed to %s' % texfilename)
+        else:
+            texfile = texdir+reportname 
+            texfile2 = reportname
+            texfilename = texfile + '.tex' 
+            fid = open(texfilename, 'w')
+
+            if verbose: log.critical('Latex output printed to %s' % texfilename)
+    else:
+        texfile = ''
+        texfile2 = ''
+
+    p = len(f_list)
+    n = []
+    n0 = 0
+    for i in range(len(f_list)):
+        n.append(len(f_list[i].get_time()))
+        if n[i] > n0: n0 = n[i]  
+    n0 = int(n0)
+    m = len(locations)
+    model_time = num.zeros((n0, m, p), num.float) 
+    stages = num.zeros((n0, m, p), num.float)
+    elevations = num.zeros((n0, m, p), num.float) 
+    momenta = num.zeros((n0, m, p), num.float)
+    xmom = num.zeros((n0, m, p), num.float)
+    ymom = num.zeros((n0, m, p), num.float)
+    speed = num.zeros((n0, m, p), num.float)
+    bearings = num.zeros((n0, m, p), num.float)
+    due_east = 90.0*num.ones((n0, 1), num.float)
+    due_west = 270.0*num.ones((n0, 1), num.float)
+    depths = num.zeros((n0, m, p), num.float)
+    eastings = num.zeros((n0, m, p), num.float)
+    min_stages = []
+    max_stages = []
+    min_momentums = []    
+    max_momentums = []
+    max_xmomentums = []
+    max_ymomentums = []
+    min_xmomentums = []
+    min_ymomentums = []
+    max_speeds = []
+    min_speeds = []    
+    max_depths = []
+    model_time_plot3d = num.zeros((n0, m), num.float)
+    stages_plot3d = num.zeros((n0, m), num.float)
+    eastings_plot3d = num.zeros((n0, m),num.float)
+    if time_unit is 'mins': scale = 60.0
+    if time_unit is 'hours': scale = 3600.0
+
+    ##### loop over each swwfile #####
+    for j, f in enumerate(f_list):
+        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
+
+        starttime = f.starttime
+        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
+        fid_compare = open(comparefile, 'w')
+        file0 = file_loc[j] + 'gauges_t0.csv'
+        fid_0 = open(file0, 'w')
+
+        ##### loop over each gauge #####
+        for k in gauge_index:
+            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
+
+            g = gauges[k]
+            min_stage = 10
+            max_stage = 0
+            max_momentum = max_xmomentum = max_ymomentum = 0
+            min_momentum = min_xmomentum = min_ymomentum = 100
+            max_speed = 0
+            min_speed = 0            
+            max_depth = 0            
+            gaugeloc = str(locations[k])
+            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
+                       + gaugeloc + '.csv'
+            if j == 0:
+                fid_out = open(thisfile, 'w')
+                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
+                fid_out.write(s)            
+
+            #### generate quantities #######
+            for i, t in enumerate(f.get_time()):
+                if time_min <= t <= time_max:
+                    w = f(t, point_id = k)[0]
+                    z = f(t, point_id = k)[1]
+                    uh = f(t, point_id = k)[2]
+                    vh = f(t, point_id = k)[3]
+                    depth = w-z      
+                    m = sqrt(uh*uh + vh*vh)
+                    if depth < 0.001:
+                        vel = 0.0
+                    else:
+                        vel = m / (depth + 1.e-6/depth) 
+                    bearing = calc_bearing(uh, vh)                    
+                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
+                    stages[i,k,j] = w
+                    elevations[i,k,j] = z 
+                    xmom[i,k,j] = uh 
+                    ymom[i,k,j] = vh 
+                    momenta[i,k,j] = m 
+                    speed[i,k,j] = vel 
+                    bearings[i,k,j] = bearing 
+                    depths[i,k,j] = depth
+                    thisgauge = gauges[k]
+                    eastings[i,k,j] = thisgauge[0]
+                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
+                            % (t, w, m, vel, z, uh, vh, bearing)
+                    fid_out.write(s)
+                    if t == 0:
+                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
+                        fid_0.write(s)
+                    if t/60.0 <= 13920: tindex = i
+                    if w > max_stage: max_stage = w
+                    if w < min_stage: min_stage = w
+                    if m > max_momentum: max_momentum = m
+                    if m < min_momentum: min_momentum = m                    
+                    if uh > max_xmomentum: max_xmomentum = uh
+                    if vh > max_ymomentum: max_ymomentum = vh
+                    if uh < min_xmomentum: min_xmomentum = uh
+                    if vh < min_ymomentum: min_ymomentum = vh
+                    if vel > max_speed: max_speed = vel
+                    if vel < min_speed: min_speed = vel                    
+                    if z > 0 and depth > max_depth: max_depth = depth
+                    
+                    
+            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
+                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
+            fid_compare.write(s)
+            max_stages.append(max_stage)
+            min_stages.append(min_stage)
+            max_momentums.append(max_momentum)
+            max_xmomentums.append(max_xmomentum)
+            max_ymomentums.append(max_ymomentum)
+            min_xmomentums.append(min_xmomentum)
+            min_ymomentums.append(min_ymomentum)
+            min_momentums.append(min_momentum)            
+            max_depths.append(max_depth)
+            max_speeds.append(max_speed)
+            min_speeds.append(min_speed)            
+            #### finished generating quantities for each swwfile #####
+        
+        model_time_plot3d[:,:] = model_time[:,:,j]
+        stages_plot3d[:,:] = stages[:,:,j]
+        eastings_plot3d[:,] = eastings[:,:,j]
+            
+        if surface is True:
+            log.critical('Printing surface figure')
+            for i in range(2):
+                fig = p1.figure(10)
+                ax = p3.Axes3D(fig)
+                if len(gauges) > 80:
+                    ax.plot_surface(model_time[:,:,j],
+                                    eastings[:,:,j],
+                                    stages[:,:,j])
+                else:
+                    ax.plot3D(num.ravel(eastings[:,:,j]),
+                              num.ravel(model_time[:,:,j]),
+                              num.ravel(stages[:,:,j]))
+                ax.set_xlabel('time')
+                ax.set_ylabel('x')
+                ax.set_zlabel('stage')
+                fig.add_axes(ax)
+                p1.show()
+                surfacefig = 'solution_surface%s' % leg_label[j]
+                p1.savefig(surfacefig)
+                p1.close()
+            
+    #### finished generating quantities for all swwfiles #####
+
+    # x profile for given time
+    if surface is True:
+        figure(11)
+        plot(eastings[tindex,:,j], stages[tindex,:,j])
+        xlabel('x')
+        ylabel('stage')
+        profilefig = 'solution_xprofile' 
+        savefig('profilefig')
+
+    elev_output = []
+    if generate_fig is True:
+        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
+                           max(max_depths)*1.1])
+        stage_axis = axis([starttime/scale, time_max/scale,
+                           min(min_stages), max(max_stages)*1.1])
+        vel_axis = axis([starttime/scale, time_max/scale,
+                         min(min_speeds), max(max_speeds)*1.1])
+        mom_axis = axis([starttime/scale, time_max/scale,
+                         min(min_momentums), max(max_momentums)*1.1])
+        xmom_axis = axis([starttime/scale, time_max/scale,
+                          min(min_xmomentums), max(max_xmomentums)*1.1])
+        ymom_axis = axis([starttime/scale, time_max/scale,
+                          min(min_ymomentums), max(max_ymomentums)*1.1])
+        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
+        nn = len(plot_quantity)
+        no_cols = 2
+        
+        if len(label_id) > 1: graphname_report = []
+        pp = 1
+        div = 11.
+        cc = 0
+        for k in gauge_index:
+            g = gauges[k]
+            count1 = 0
+            if report == True and len(label_id) > 1:
+                s = '\\begin{figure}[ht] \n' \
+                    '\\centering \n' \
+                    '\\begin{tabular}{cc} \n'
+                fid.write(s)
+            if len(label_id) > 1: graphname_report = []
+
+            #### generate figures for each gauge ####
+            for j, f in enumerate(f_list):
+                ion()
+                hold(True)
+                count = 0
+                where1 = 0
+                where2 = 0
+                word_quantity = ''
+                if report == True and len(label_id) == 1:
+                    s = '\\begin{figure}[hbt] \n' \
+                        '\\centering \n' \
+                        '\\begin{tabular}{cc} \n'
+                    fid.write(s)
+                    
+                for which_quantity in plot_quantity:
+                    count += 1
+                    where1 += 1
+                    figure(count, frameon = False)
+                    if which_quantity == 'depth':
+                        plot(model_time[0:n[j]-1,k,j],
+                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
+                        units = 'm'
+                        axis(depth_axis)
+                    if which_quantity == 'stage':
+                        if elevations[0,k,j] <= 0:
+                            plot(model_time[0:n[j]-1,k,j],
+                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
+                            axis(stage_axis)
+                        else:
+                            plot(model_time[0:n[j]-1,k,j],
+                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
+                            #axis(depth_axis)                 
+                        units = 'm'
+                    if which_quantity == 'momentum':
+                        plot(model_time[0:n[j]-1,k,j],
+                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
+                        axis(mom_axis)
+                        units = 'm^2 / sec'
+                    if which_quantity == 'xmomentum':
+                        plot(model_time[0:n[j]-1,k,j],
+                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
+                        axis(xmom_axis)
+                        units = 'm^2 / sec'
+                    if which_quantity == 'ymomentum':
+                        plot(model_time[0:n[j]-1,k,j],
+                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
+                        axis(ymom_axis)
+                        units = 'm^2 / sec'
+                    if which_quantity == 'speed':
+                        plot(model_time[0:n[j]-1,k,j],
+                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
+                        axis(vel_axis)
+                        units = 'm / sec'
+                    if which_quantity == 'bearing':
+                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
+                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
+                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
+                        units = 'degrees from North'
+                        #ax = axis([time_min, time_max, 0.0, 360.0])
+                        legend(('Bearing','West','East'))
+
+                    if time_unit is 'mins': xlabel('time (mins)')
+                    if time_unit is 'hours': xlabel('time (hours)')
+                    #if which_quantity == 'stage' \
+                    #   and elevations[0:n[j]-1,k,j] > 0:
+                    #    ylabel('%s (%s)' %('depth', units))
+                    #else:
+                    #    ylabel('%s (%s)' %(which_quantity, units))
+                        #ylabel('%s (%s)' %('wave height', units))
+                    ylabel('%s (%s)' %(which_quantity, units))
+                    if len(label_id) > 1: legend((leg_label),loc='upper right')
+
+                    #gaugeloc1 = gaugeloc.replace(' ','')
+                    #gaugeloc2 = gaugeloc1.replace('_','')
+                    gaugeloc2 = str(locations[k]).replace(' ','')
+                    graphname = '%sgauge%s_%s' %(file_loc[j],
+                                                 gaugeloc2,
+                                                 which_quantity)
+
+                    if report == True and len(label_id) > 1:
+                        figdir = getcwd()+sep+'report_figures'+sep
+                        if access(figdir,F_OK) == 0 :
+                            mkdir (figdir)
+                        latex_file_loc = figdir.replace(sep,altsep) 
+                        # storing files in production directory    
+                        graphname_latex = '%sgauge%s%s' \
+                                          % (latex_file_loc, gaugeloc2,
+                                             which_quantity)
+                        # giving location in latex output file
+                        graphname_report_input = '%sgauge%s%s' % \
+                                                 ('..' + altsep + 
+                                                      'report_figures' + altsep,
+                                                  gaugeloc2, which_quantity)
+                        graphname_report.append(graphname_report_input)
+                        
+                        # save figures in production directory for report
+                        savefig(graphname_latex)
+
+                    if report == True:
+                        figdir = getcwd() + sep + 'report_figures' + sep
+                        if access(figdir,F_OK) == 0:
+                            mkdir(figdir)
+                        latex_file_loc = figdir.replace(sep,altsep)    
+
+                        if len(label_id) == 1: 
+                            # storing files in production directory  
+                            graphname_latex = '%sgauge%s%s%s' % \
+                                              (latex_file_loc, gaugeloc2,
+                                               which_quantity, label_id2)
+                            # giving location in latex output file
+                            graphname_report = '%sgauge%s%s%s' % \
+                                               ('..' + altsep +
+                                                    'report_figures' + altsep,
+                                                gaugeloc2, which_quantity,
+                                                label_id2)
+                            s = '\includegraphics' \
+                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
+                                (graphname_report, '.png')
+                            fid.write(s)
+                            if where1 % 2 == 0:
+                                s = '\\\\ \n'
+                                where1 = 0
+                            else:
+                                s = '& \n'
+                            fid.write(s)
+                            savefig(graphname_latex)
+                    
+                    if title_on == True:
+                        title('%s scenario: %s at %s gauge' % \
+                              (label_id, which_quantity, gaugeloc2))
+                        #title('Gauge %s (MOST elevation %.2f, ' \
+                        #      'ANUGA elevation %.2f)' % \
+                        #      (gaugeloc2, elevations[10,k,0],
+                        #       elevations[10,k,1]))
+
+                    savefig(graphname) # save figures with sww file
+
+                if report == True and len(label_id) == 1:
+                    for i in range(nn-1):
+                        if nn > 2:
+                            if plot_quantity[i] == 'stage' \
+                               and elevations[0,k,j] > 0:
+                                word_quantity += 'depth' + ', '
+                            else:
+                                word_quantity += plot_quantity[i] + ', '
+                        else:
+                            if plot_quantity[i] == 'stage' \
+                               and elevations[0,k,j] > 0:
+                                word_quantity += 'depth' + ', '
+                            else:
+                                word_quantity += plot_quantity[i]
+                        
+                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
+                        word_quantity += ' and ' + 'depth'
+                    else:
+                        word_quantity += ' and ' + plot_quantity[nn-1]
+                    caption = 'Time series for %s at %s location ' \
+                              '(elevation %.2fm)' % \
+                              (word_quantity, locations[k], elev[k])
+                    if elev[k] == 0.0:
+                        caption = 'Time series for %s at %s location ' \
+                                  '(elevation %.2fm)' % \
+                                  (word_quantity, locations[k],
+                                   elevations[0,k,j])
+                        east = gauges[0]
+                        north = gauges[1]
+                        elev_output.append([locations[k], east, north,
+                                            elevations[0,k,j]])
+                    label = '%sgauge%s' % (label_id2, gaugeloc2)
+                    s = '\end{tabular} \n' \
+                        '\\caption{%s} \n' \
+                        '\label{fig:%s} \n' \
+                        '\end{figure} \n \n' % (caption, label)
+                    fid.write(s)
+                    cc += 1
+                    if cc % 6 == 0: fid.write('\\clearpage \n')
+                    savefig(graphname_latex)               
+                    
+            if report == True and len(label_id) > 1:
+                for i in range(nn-1):
+                    if nn > 2:
+                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
+                            word_quantity += 'depth' + ','
+                        else:
+                            word_quantity += plot_quantity[i] + ', '
+                    else:
+                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
+                            word_quantity += 'depth'
+                        else:
+                            word_quantity += plot_quantity[i]
+                    where1 = 0
+                    count1 += 1
+                    index = j*len(plot_quantity)
+                    for which_quantity in plot_quantity:
+                        where1 += 1
+                        s = '\includegraphics' \
+                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
+                            (graphname_report[index], '.png')
+                        index += 1
+                        fid.write(s)
+                        if where1 % 2 == 0:
+                            s = '\\\\ \n'
+                            where1 = 0
+                        else:
+                            s = '& \n'
+                        fid.write(s)
+                word_quantity += ' and ' + plot_quantity[nn-1]            
+                label = 'gauge%s' %(gaugeloc2) 
+                caption = 'Time series for %s at %s location ' \
+                          '(elevation %.2fm)' % \
+                          (word_quantity, locations[k], elev[k])
+                if elev[k] == 0.0:
+                        caption = 'Time series for %s at %s location ' \
+                                  '(elevation %.2fm)' % \
+                                  (word_quantity, locations[k],
+                                   elevations[0,k,j])
+                        thisgauge = gauges[k]
+                        east = thisgauge[0]
+                        north = thisgauge[1]
+                        elev_output.append([locations[k], east, north,
+                                            elevations[0,k,j]])
+                        
+                s = '\end{tabular} \n' \
+                    '\\caption{%s} \n' \
+                    '\label{fig:%s} \n' \
+                    '\end{figure} \n \n' % (caption, label)
+                fid.write(s)
+                if float((k+1)/div - pp) == 0.:
+                    fid.write('\\clearpage \n')
+                    pp += 1
+                #### finished generating figures ###
+
+            close('all')
+        
+    return texfile2, elev_output

anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

@@ -29 +29 @@
     return x+y
 
-class Test_Util(unittest.TestCase):
+class Test_Gauge(unittest.TestCase):
     def setUp(self):
-        pass
-
-    def tearDown(self):
-        pass
-
-    def test_sww2csv(self):
 
         def elevation_function(x, y):
             return -x
         
-        """Most of this test was copied from test_interpolate
-        test_interpole_sww2csv
-        
-        This is testing the sww2csv_gauges function, by creating a sww file and
-        then exporting the gauges and checking the results.
-        """
-        
-        # Create mesh
+        """ Setup for all tests. """
+        
         mesh_file = tempfile.mktemp(".tsh")    
         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
@@ -60 +48 @@
         domain = Domain(mesh_file)
         os.remove(mesh_file)
-        
-        domain.default_order=2
-        
+ 
+        domain.default_order = 2
+
         # This test was made before tight_slope_limiters were introduced
         # Since were are testing interpolation values this is OK
-        domain.tight_slope_limiters = 0 
-        
-
+        domain.tight_slope_limiters = 0
+                
         # Set some field values
         domain.set_quantity('elevation', elevation_function)
@@ -83 +70 @@
         domain.set_quantity('stage', 1.0)
 
-
         domain.set_name('datatest' + str(time.time()))
         domain.smooth = True
         domain.reduction = mean
-
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        domain.set_quantity('stage', 10.0) # This is automatically limited
+        
+        self.domain = domain
+        
+        
+    def tearDown(self):
+        """Called at end of each test."""
+        if self.sww:
+            os.remove(self.sww.filename)
+
+    def _create_sww(self):
+        self.sww = SWW_file(self.domain)
+        self.sww.store_connectivity()
+        self.sww.store_timestep()
+        self.domain.set_quantity('stage', 10.0) # This is automatically limited
         # so it will not be less than the elevation
-        domain.time = 2.
-        sww.store_timestep()
-
+        self.domain.time = 2.
+        self.sww.store_timestep()
+        
+        
+    def test_sww2csv(self):
+
+        """Most of this test was copied from test_interpolate
+        test_interpole_sww2csv
+        
+        This is testing the sww2csv_gauges function, by creating a sww file and
+        then exporting the gauges and checking the results.
+        """
+        
+        domain = self.domain
+        self._create_sww()
+        
         # test the function
         points = [[5.0,1.],[0.5,2.]]
@@ -109 +116 @@
 
         
-        sww2csv_gauges(sww.filename, 
+        sww2csv_gauges(self.sww.filename, 
                        points_file,
                        verbose=False,
@@ -146 +153 @@
         point1_handle.close()
         point2_handle.close()
-        #print "sww.filename",sww.filename 
-        os.remove(sww.filename)
         os.remove(points_file)
         os.remove(point1_filename)
@@ -159 +164 @@
         import time
         import string
-
-        def elevation_function(x, y):
-            return -x
         
         """Most of this test was copied from test_interpolate
@@ -173 +175 @@
         """
         
-        # Create mesh
-        mesh_file = tempfile.mktemp(".tsh")    
-        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
-        m = Mesh()
-        m.add_vertices(points)
-        m.auto_segment()
-        m.generate_mesh(verbose=False)
-        m.export_mesh_file(mesh_file)
-        
-        # Create shallow water domain
-        domain = Domain(mesh_file)
-        os.remove(mesh_file)
-        
-        domain.default_order=2
-
-        # Set some field values
-        domain.set_quantity('elevation', elevation_function)
-        domain.set_quantity('friction', 0.03)
-        domain.set_quantity('xmomentum', 3.0)
-        domain.set_quantity('ymomentum', 4.0)
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        # This call mangles the stage values.
-        domain.distribute_to_vertices_and_edges()
-        domain.set_quantity('stage', 1.0)
-
-
-        domain.set_name('datatest' + str(time.time()))
-        domain.smooth = True
-        domain.reduction = mean
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        domain.set_quantity('stage', 10.0) # This is automatically limited
-        # so it will not be less than the elevation
-        domain.time = 2.
-        sww.store_timestep()
-
+        domain = self.domain
+        self._create_sww()
+        
         # test the function
         points = [[5.0,1.],[0.5,2.]]
@@ -227 +189 @@
         file_id.close()
 
-        sww2csv_gauges(sww.filename, 
+        sww2csv_gauges(self.sww.filename, 
                             points_file,
                             quantities=['stage', 'elevation'],
@@ -263 +225 @@
         # clean up
         point1_handle.close()
-        point2_handle.close()
-        #print "sww.filename",sww.filename 
-        os.remove(sww.filename)
+        point2_handle.close() 
         os.remove(points_file)
         os.remove(point1_filename)
-        os.remove(point2_filename)
-
+        os.remove(point2_filename)        
+        
 
     def test_sww2csv_gauges2(self):
-
-        def elevation_function(x, y):
-            return -x
         
         """Most of this test was copied from test_interpolate
@@ -286 +243 @@
         """
         
-        # Create mesh
-        mesh_file = tempfile.mktemp(".tsh")    
-        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
-        m = Mesh()
-        m.add_vertices(points)
-        m.auto_segment()
-        m.generate_mesh(verbose=False)
-        m.export_mesh_file(mesh_file)
-        
-        # Create shallow water domain
-        domain = Domain(mesh_file)
-        os.remove(mesh_file)
-        
-        domain.default_order=2
-
-        # This test was made before tight_slope_limiters were introduced
-        # Since were are testing interpolation values this is OK
-        domain.tight_slope_limiters = 0         
-
-        # Set some field values
-        domain.set_quantity('elevation', elevation_function)
-        domain.set_quantity('friction', 0.03)
-        domain.set_quantity('xmomentum', 3.0)
-        domain.set_quantity('ymomentum', 4.0)
+        domain = self.domain
         domain.set_starttime(5)
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        # This call mangles the stage values.
-        domain.distribute_to_vertices_and_edges()
-        domain.set_quantity('stage', 1.0)
-        
-
-
-        domain.set_name('datatest' + str(time.time()))
-        domain.smooth = True
-        domain.reduction = mean
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        domain.set_quantity('stage', 10.0) # This is automatically limited
-        # so it will not be less than the elevation
-        domain.time = 2.
-        sww.store_timestep()
-
+        self._create_sww()
+        
         # test the function
         points = [[5.0,1.],[0.5,2.]]
@@ -346 +258 @@
         file_id.close()
 
-        
-        sww2csv_gauges(sww.filename, 
+        sww2csv_gauges(self.sww.filename, 
                             points_file,
                             verbose=False,
@@ -364 +275 @@
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
                          float(row[4]), float(row[5]), float(row[6])])
-            #print 'assert line',line[i],'point1',point1_answers_array[i]
+            #print 'assert line',line[i],'answer',point1_answers_array[i]
             assert num.allclose(line[i], point1_answers_array[i])
 
@@ -384 +295 @@
         point1_handle.close()
         point2_handle.close()
-        #print "sww.filename",sww.filename 
-        os.remove(sww.filename)
         os.remove(points_file)
         os.remove(point1_filename)
         os.remove(point2_filename)
-        
-
+
+
+       
+    # def test_sww2csv_gauge_point_off_mesh(self):
+        # from anuga.pmesh.mesh import Mesh
+        # from anuga.shallow_water import Domain, Transmissive_boundary
+        # from csv import reader,writer
+        # import time
+        # import string
+        
+        # """Most of this test was copied from test_interpolate
+        # test_interpole_sww2csv
+        
+        # This is testing the sww2csv_gauges function with one gauge off the mesh, by creating a sww file and
+        # then exporting the gauges and checking the results.
+        
+        # This tests the correct values for when a gauge is off the mesh, which is important for parallel.
+        # """
+        
+        # domain = self.domain
+        # self._create_sww()     
+     
+        # # test the function
+        # points = [[50.0,1.],[50.5,-20.25]]
+
+# #        points_file = tempfile.mktemp(".csv")
+        # points_file = 'test_point.csv'
+        # file_id = open(points_file,"w")
+        # file_id.write("name,easting,northing \n\
+# point1, 50.0, 1.0\n\
+# point2, 50.5, 20.25\n")
+        # file_id.close()
+
+        # sww2csv_gauges(sww.filename, 
+                            # points_file,
+                            # quantities=['stage', 'elevation'],
+                            # use_cache=False,
+                            # verbose=False)
+
+        # point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
+        # point1_filename = 'gauge_point1.csv'
+        # point1_handle = file(point1_filename)
+        # point1_reader = reader(point1_handle)
+        # point1_reader.next()
+
+        # line=[]
+        # for i,row in enumerate(point1_reader):
+# #            print 'i',i,'row',row
+            # # note the 'hole' (element 1) below - skip the new 'hours' field
+            # line.append([float(row[0]),float(row[2]),float(row[3])])
+            # #print 'line',line[i],'point1',point1_answers_array[i]
+            # assert num.allclose(line[i], point1_answers_array[i])
+
+        # point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
+        # point2_filename = 'gauge_point2.csv' 
+        # point2_handle = file(point2_filename)
+        # point2_reader = reader(point2_handle)
+        # point2_reader.next()
+                        
+        # line=[]
+        # for i,row in enumerate(point2_reader):
+# #            print 'i',i,'row',row
+            # # note the 'hole' (element 1) below - skip the new 'hours' field
+            # line.append([float(row[0]),float(row[2]),float(row[3])])
+# #            print 'line',line[i],'point1',point1_answers_array[i]
+            # assert num.allclose(line[i], point2_answers_array[i])
+                         
+        # # clean up
+        # point1_handle.close()
+        # point2_handle.close()
+        # os.remove(points_file)
+# #        os.remove(point1_filename)
+# #        os.remove(point2_filename)
+        
+        
     def test_sww2csv_centroid(self):
-
-        def elevation_function(x, y):
-            return -x
         
         """Check sww2csv timeseries at centroid.
@@ -402 +381 @@
         """
         
-        # Create rectangular mesh
-        mesh_file = tempfile.mktemp(".tsh")    
-        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
-        m = Mesh()
-        m.add_vertices(points)
-        m.auto_segment()
-        m.generate_mesh(verbose=False)
-        m.export_mesh_file(mesh_file)
-        
-        # Create shallow water domain
-        domain = Domain(mesh_file)
-        
-        domain.default_order=2
-        
-        # This test was made before tight_slope_limiters were introduced
-        # Since were are testing interpolation values this is OK
-        domain.tight_slope_limiters = 0 
-        
-
-        # Set some field values
-        domain.set_quantity('elevation', elevation_function)
-        domain.set_quantity('friction', 0.03)
-        domain.set_quantity('xmomentum', 3.0)
-        domain.set_quantity('ymomentum', 4.0)
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        # This call mangles the stage values.
-        domain.distribute_to_vertices_and_edges()
-        domain.set_quantity('stage', 1.0)
-
-
-        domain.set_name('datatest' + str(time.time()))
-        domain.smooth = True
-        domain.reduction = mean
-
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        domain.set_quantity('stage', 10.0) # This is automatically limited
-        # so it will not be less than the elevation
-        domain.time = 2.
-        sww.store_timestep()
-
+        domain = self.domain
+        sww = self._create_sww()
+        
         # create a csv file containing our gauge points
         points_file = tempfile.mktemp(".csv")
@@ -466 +400 @@
         file_id.close()
 
-        sww2csv_gauges(sww.filename, 
+        sww2csv_gauges(self.sww.filename, 
                        points_file,
                        verbose=False,
@@ -501 +435 @@
         point1_handle.close()
         point2_handle.close()
-        os.remove(mesh_file)        
-        os.remove(sww.filename)
         os.remove(points_file)
         os.remove(point1_filename)
@@ -509 +441 @@
 
     def test_sww2csv_output_centroid_attribute(self):
-
-        def elevation_function(x, y):
-            return -x
         
         """Check sww2csv timeseries at centroid, then output the centroid coordinates.
@@ -519 +448 @@
         """
         
-        # Create rectangular mesh
-        mesh_file = tempfile.mktemp(".tsh")    
-        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
-        m = Mesh()
-        m.add_vertices(points)
-        m.auto_segment()
-        m.generate_mesh(verbose=False)
-        m.export_mesh_file(mesh_file)
-        
-        # Create shallow water domain
-        domain = Domain(mesh_file)
-        
-        domain.default_order=2
-        
-        # This test was made before tight_slope_limiters were introduced
-        # Since were are testing interpolation values this is OK
-        domain.tight_slope_limiters = 0 
-        
-
-        # Set some field values
-        domain.set_quantity('elevation', elevation_function)
-        domain.set_quantity('friction', 0.03)
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        # This call mangles the stage values.
-        domain.distribute_to_vertices_and_edges()
-        domain.set_quantity('stage', 1.0)
-
-
-        domain.set_name('datatest' + str(time.time()))
-        domain.smooth = True
-        domain.reduction = mean
-
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        domain.set_quantity('stage', 10.0) # This is automatically limited
-        # so it will not be less than the elevation
-        domain.time = 2.
-        sww.store_timestep()
-
+        domain = self.domain        
+        self._create_sww()
+        
         # create a csv file containing our gauge points
         points_file = tempfile.mktemp(".csv")
@@ -575 +461 @@
         file_id.close()
 
-        sww2csv_gauges(sww.filename, 
+        sww2csv_gauges(self.sww.filename, 
                        points_file,
                        quantities=['stage', 'xcentroid', 'ycentroid'],
@@ -595 +481 @@
 
         # clean up
-        point1_handle.close()
-        os.remove(mesh_file)        
-        os.remove(sww.filename)
+        point1_handle.close()        
         os.remove(points_file)
         os.remove(point1_filename)
@@ -606 +490 @@
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Util, 'test')
+    suite = unittest.makeSuite(Test_Gauge, 'test')
 #    runner = unittest.TextTestRunner(verbosity=2)
     runner = unittest.TextTestRunner(verbosity=1)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -317 +317 @@
 
     return bearing
-
-
-##
-# @brief Generate figures from quantities and gauges for each sww file.
-# @param plot_quantity  ??
-# @param file_loc ??
-# @param report ??
-# @param reportname ??
-# @param surface ??
-# @param leg_label ??
-# @param f_list ??
-# @param gauges ??
-# @param locations ??
-# @param elev ??
-# @param gauge_index ??
-# @param production_dirs ??
-# @param time_min ??
-# @param time_max ??
-# @param time_unit ??
-# @param title_on ??
-# @param label_id ??
-# @param generate_fig ??
-# @param verbose??
-# @return (texfile2, elev_output)
-def generate_figures(plot_quantity, file_loc, report, reportname, surface,
-                     leg_label, f_list, gauges, locations, elev, gauge_index,
-                     production_dirs, time_min, time_max, time_unit,
-                     title_on, label_id, generate_fig, verbose):
-    """ Generate figures based on required quantities and gauges for
-    each sww file
-    """
-    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
-
-    if generate_fig is True:
-        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
-             xlabel, ylabel, title, close, subplot
-    
-        if surface is True:
-            import pylab as p1
-            import mpl3d.mplot3d as p3
-        
-    if report == True:    
-        texdir = getcwd()+sep+'report'+sep
-        if access(texdir,F_OK) == 0:
-            mkdir (texdir)
-        if len(label_id) == 1:
-            label_id1 = label_id[0].replace(sep,'')
-            label_id2 = label_id1.replace('_','')
-            texfile = texdir + reportname + '%s' % label_id2
-            texfile2 = reportname + '%s' % label_id2
-            texfilename = texfile + '.tex'
-            fid = open(texfilename, 'w')
-
-            if verbose: log.critical('Latex output printed to %s' % texfilename)
-        else:
-            texfile = texdir+reportname 
-            texfile2 = reportname
-            texfilename = texfile + '.tex' 
-            fid = open(texfilename, 'w')
-
-            if verbose: log.critical('Latex output printed to %s' % texfilename)
-    else:
-        texfile = ''
-        texfile2 = ''
-
-    p = len(f_list)
-    n = []
-    n0 = 0
-    for i in range(len(f_list)):
-        n.append(len(f_list[i].get_time()))
-        if n[i] > n0: n0 = n[i]  
-    n0 = int(n0)
-    m = len(locations)
-    model_time = num.zeros((n0, m, p), num.float) 
-    stages = num.zeros((n0, m, p), num.float)
-    elevations = num.zeros((n0, m, p), num.float) 
-    momenta = num.zeros((n0, m, p), num.float)
-    xmom = num.zeros((n0, m, p), num.float)
-    ymom = num.zeros((n0, m, p), num.float)
-    speed = num.zeros((n0, m, p), num.float)
-    bearings = num.zeros((n0, m, p), num.float)
-    due_east = 90.0*num.ones((n0, 1), num.float)
-    due_west = 270.0*num.ones((n0, 1), num.float)
-    depths = num.zeros((n0, m, p), num.float)
-    eastings = num.zeros((n0, m, p), num.float)
-    min_stages = []
-    max_stages = []
-    min_momentums = []    
-    max_momentums = []
-    max_xmomentums = []
-    max_ymomentums = []
-    min_xmomentums = []
-    min_ymomentums = []
-    max_speeds = []
-    min_speeds = []    
-    max_depths = []
-    model_time_plot3d = num.zeros((n0, m), num.float)
-    stages_plot3d = num.zeros((n0, m), num.float)
-    eastings_plot3d = num.zeros((n0, m),num.float)
-    if time_unit is 'mins': scale = 60.0
-    if time_unit is 'hours': scale = 3600.0
-
-    ##### loop over each swwfile #####
-    for j, f in enumerate(f_list):
-        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
-
-        starttime = f.starttime
-        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
-        fid_compare = open(comparefile, 'w')
-        file0 = file_loc[j] + 'gauges_t0.csv'
-        fid_0 = open(file0, 'w')
-
-        ##### loop over each gauge #####
-        for k in gauge_index:
-            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
-
-            g = gauges[k]
-            min_stage = 10
-            max_stage = 0
-            max_momentum = max_xmomentum = max_ymomentum = 0
-            min_momentum = min_xmomentum = min_ymomentum = 100
-            max_speed = 0
-            min_speed = 0            
-            max_depth = 0            
-            gaugeloc = str(locations[k])
-            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
-                       + gaugeloc + '.csv'
-            if j == 0:
-                fid_out = open(thisfile, 'w')
-                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
-                fid_out.write(s)            
-
-            #### generate quantities #######
-            for i, t in enumerate(f.get_time()):
-                if time_min <= t <= time_max:
-                    w = f(t, point_id = k)[0]
-                    z = f(t, point_id = k)[1]
-                    uh = f(t, point_id = k)[2]
-                    vh = f(t, point_id = k)[3]
-                    depth = w-z      
-                    m = sqrt(uh*uh + vh*vh)
-                    if depth < 0.001:
-                        vel = 0.0
-                    else:
-                        vel = m / (depth + 1.e-6/depth) 
-                    bearing = calc_bearing(uh, vh)                    
-                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
-                    stages[i,k,j] = w
-                    elevations[i,k,j] = z 
-                    xmom[i,k,j] = uh 
-                    ymom[i,k,j] = vh 
-                    momenta[i,k,j] = m 
-                    speed[i,k,j] = vel 
-                    bearings[i,k,j] = bearing 
-                    depths[i,k,j] = depth
-                    thisgauge = gauges[k]
-                    eastings[i,k,j] = thisgauge[0]
-                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
-                            % (t, w, m, vel, z, uh, vh, bearing)
-                    fid_out.write(s)
-                    if t == 0:
-                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
-                        fid_0.write(s)
-                    if t/60.0 <= 13920: tindex = i
-                    if w > max_stage: max_stage = w
-                    if w < min_stage: min_stage = w
-                    if m > max_momentum: max_momentum = m
-                    if m < min_momentum: min_momentum = m                    
-                    if uh > max_xmomentum: max_xmomentum = uh
-                    if vh > max_ymomentum: max_ymomentum = vh
-                    if uh < min_xmomentum: min_xmomentum = uh
-                    if vh < min_ymomentum: min_ymomentum = vh
-                    if vel > max_speed: max_speed = vel
-                    if vel < min_speed: min_speed = vel                    
-                    if z > 0 and depth > max_depth: max_depth = depth
-                    
-                    
-            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
-                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
-            fid_compare.write(s)
-            max_stages.append(max_stage)
-            min_stages.append(min_stage)
-            max_momentums.append(max_momentum)
-            max_xmomentums.append(max_xmomentum)
-            max_ymomentums.append(max_ymomentum)
-            min_xmomentums.append(min_xmomentum)
-            min_ymomentums.append(min_ymomentum)
-            min_momentums.append(min_momentum)            
-            max_depths.append(max_depth)
-            max_speeds.append(max_speed)
-            min_speeds.append(min_speed)            
-            #### finished generating quantities for each swwfile #####
-        
-        model_time_plot3d[:,:] = model_time[:,:,j]
-        stages_plot3d[:,:] = stages[:,:,j]
-        eastings_plot3d[:,] = eastings[:,:,j]
-            
-        if surface is True:
-            log.critical('Printing surface figure')
-            for i in range(2):
-                fig = p1.figure(10)
-                ax = p3.Axes3D(fig)
-                if len(gauges) > 80:
-                    ax.plot_surface(model_time[:,:,j],
-                                    eastings[:,:,j],
-                                    stages[:,:,j])
-                else:
-                    ax.plot3D(num.ravel(eastings[:,:,j]),
-                              num.ravel(model_time[:,:,j]),
-                              num.ravel(stages[:,:,j]))
-                ax.set_xlabel('time')
-                ax.set_ylabel('x')
-                ax.set_zlabel('stage')
-                fig.add_axes(ax)
-                p1.show()
-                surfacefig = 'solution_surface%s' % leg_label[j]
-                p1.savefig(surfacefig)
-                p1.close()
-            
-    #### finished generating quantities for all swwfiles #####
-
-    # x profile for given time
-    if surface is True:
-        figure(11)
-        plot(eastings[tindex,:,j], stages[tindex,:,j])
-        xlabel('x')
-        ylabel('stage')
-        profilefig = 'solution_xprofile' 
-        savefig('profilefig')
-
-    elev_output = []
-    if generate_fig is True:
-        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
-                           max(max_depths)*1.1])
-        stage_axis = axis([starttime/scale, time_max/scale,
-                           min(min_stages), max(max_stages)*1.1])
-        vel_axis = axis([starttime/scale, time_max/scale,
-                         min(min_speeds), max(max_speeds)*1.1])
-        mom_axis = axis([starttime/scale, time_max/scale,
-                         min(min_momentums), max(max_momentums)*1.1])
-        xmom_axis = axis([starttime/scale, time_max/scale,
-                          min(min_xmomentums), max(max_xmomentums)*1.1])
-        ymom_axis = axis([starttime/scale, time_max/scale,
-                          min(min_ymomentums), max(max_ymomentums)*1.1])
-        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
-        nn = len(plot_quantity)
-        no_cols = 2
-        
-        if len(label_id) > 1: graphname_report = []
-        pp = 1
-        div = 11.
-        cc = 0
-        for k in gauge_index:
-            g = gauges[k]
-            count1 = 0
-            if report == True and len(label_id) > 1:
-                s = '\\begin{figure}[ht] \n' \
-                    '\\centering \n' \
-                    '\\begin{tabular}{cc} \n'
-                fid.write(s)
-            if len(label_id) > 1: graphname_report = []
-
-            #### generate figures for each gauge ####
-            for j, f in enumerate(f_list):
-                ion()
-                hold(True)
-                count = 0
-                where1 = 0
-                where2 = 0
-                word_quantity = ''
-                if report == True and len(label_id) == 1:
-                    s = '\\begin{figure}[hbt] \n' \
-                        '\\centering \n' \
-                        '\\begin{tabular}{cc} \n'
-                    fid.write(s)
-                    
-                for which_quantity in plot_quantity:
-                    count += 1
-                    where1 += 1
-                    figure(count, frameon = False)
-                    if which_quantity == 'depth':
-                        plot(model_time[0:n[j]-1,k,j],
-                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
-                        units = 'm'
-                        axis(depth_axis)
-                    if which_quantity == 'stage':
-                        if elevations[0,k,j] <= 0:
-                            plot(model_time[0:n[j]-1,k,j],
-                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
-                            axis(stage_axis)
-                        else:
-                            plot(model_time[0:n[j]-1,k,j],
-                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
-                            #axis(depth_axis)                 
-                        units = 'm'
-                    if which_quantity == 'momentum':
-                        plot(model_time[0:n[j]-1,k,j],
-                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
-                        axis(mom_axis)
-                        units = 'm^2 / sec'
-                    if which_quantity == 'xmomentum':
-                        plot(model_time[0:n[j]-1,k,j],
-                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
-                        axis(xmom_axis)
-                        units = 'm^2 / sec'
-                    if which_quantity == 'ymomentum':
-                        plot(model_time[0:n[j]-1,k,j],
-                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
-                        axis(ymom_axis)
-                        units = 'm^2 / sec'
-                    if which_quantity == 'speed':
-                        plot(model_time[0:n[j]-1,k,j],
-                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
-                        axis(vel_axis)
-                        units = 'm / sec'
-                    if which_quantity == 'bearing':
-                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
-                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
-                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
-                        units = 'degrees from North'
-                        #ax = axis([time_min, time_max, 0.0, 360.0])
-                        legend(('Bearing','West','East'))
-
-                    if time_unit is 'mins': xlabel('time (mins)')
-                    if time_unit is 'hours': xlabel('time (hours)')
-                    #if which_quantity == 'stage' \
-                    #   and elevations[0:n[j]-1,k,j] > 0:
-                    #    ylabel('%s (%s)' %('depth', units))
-                    #else:
-                    #    ylabel('%s (%s)' %(which_quantity, units))
-                        #ylabel('%s (%s)' %('wave height', units))
-                    ylabel('%s (%s)' %(which_quantity, units))
-                    if len(label_id) > 1: legend((leg_label),loc='upper right')
-
-                    #gaugeloc1 = gaugeloc.replace(' ','')
-                    #gaugeloc2 = gaugeloc1.replace('_','')
-                    gaugeloc2 = str(locations[k]).replace(' ','')
-                    graphname = '%sgauge%s_%s' %(file_loc[j],
-                                                 gaugeloc2,
-                                                 which_quantity)
-
-                    if report == True and len(label_id) > 1:
-                        figdir = getcwd()+sep+'report_figures'+sep
-                        if access(figdir,F_OK) == 0 :
-                            mkdir (figdir)
-                        latex_file_loc = figdir.replace(sep,altsep) 
-                        # storing files in production directory    
-                        graphname_latex = '%sgauge%s%s' \
-                                          % (latex_file_loc, gaugeloc2,
-                                             which_quantity)
-                        # giving location in latex output file
-                        graphname_report_input = '%sgauge%s%s' % \
-                                                 ('..' + altsep + 
-                                                      'report_figures' + altsep,
-                                                  gaugeloc2, which_quantity)
-                        graphname_report.append(graphname_report_input)
-                        
-                        # save figures in production directory for report
-                        savefig(graphname_latex)
-
-                    if report == True:
-                        figdir = getcwd() + sep + 'report_figures' + sep
-                        if access(figdir,F_OK) == 0:
-                            mkdir(figdir)
-                        latex_file_loc = figdir.replace(sep,altsep)    
-
-                        if len(label_id) == 1: 
-                            # storing files in production directory  
-                            graphname_latex = '%sgauge%s%s%s' % \
-                                              (latex_file_loc, gaugeloc2,
-                                               which_quantity, label_id2)
-                            # giving location in latex output file
-                            graphname_report = '%sgauge%s%s%s' % \
-                                               ('..' + altsep +
-                                                    'report_figures' + altsep,
-                                                gaugeloc2, which_quantity,
-                                                label_id2)
-                            s = '\includegraphics' \
-                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
-                                (graphname_report, '.png')
-                            fid.write(s)
-                            if where1 % 2 == 0:
-                                s = '\\\\ \n'
-                                where1 = 0
-                            else:
-                                s = '& \n'
-                            fid.write(s)
-                            savefig(graphname_latex)
-                    
-                    if title_on == True:
-                        title('%s scenario: %s at %s gauge' % \
-                              (label_id, which_quantity, gaugeloc2))
-                        #title('Gauge %s (MOST elevation %.2f, ' \
-                        #      'ANUGA elevation %.2f)' % \
-                        #      (gaugeloc2, elevations[10,k,0],
-                        #       elevations[10,k,1]))
-
-                    savefig(graphname) # save figures with sww file
-
-                if report == True and len(label_id) == 1:
-                    for i in range(nn-1):
-                        if nn > 2:
-                            if plot_quantity[i] == 'stage' \
-                               and elevations[0,k,j] > 0:
-                                word_quantity += 'depth' + ', '
-                            else:
-                                word_quantity += plot_quantity[i] + ', '
-                        else:
-                            if plot_quantity[i] == 'stage' \
-                               and elevations[0,k,j] > 0:
-                                word_quantity += 'depth' + ', '
-                            else:
-                                word_quantity += plot_quantity[i]
-                        
-                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
-                        word_quantity += ' and ' + 'depth'
-                    else:
-                        word_quantity += ' and ' + plot_quantity[nn-1]
-                    caption = 'Time series for %s at %s location ' \
-                              '(elevation %.2fm)' % \
-                              (word_quantity, locations[k], elev[k])
-                    if elev[k] == 0.0:
-                        caption = 'Time series for %s at %s location ' \
-                                  '(elevation %.2fm)' % \
-                                  (word_quantity, locations[k],
-                                   elevations[0,k,j])
-                        east = gauges[0]
-                        north = gauges[1]
-                        elev_output.append([locations[k], east, north,
-                                            elevations[0,k,j]])
-                    label = '%sgauge%s' % (label_id2, gaugeloc2)
-                    s = '\end{tabular} \n' \
-                        '\\caption{%s} \n' \
-                        '\label{fig:%s} \n' \
-                        '\end{figure} \n \n' % (caption, label)
-                    fid.write(s)
-                    cc += 1
-                    if cc % 6 == 0: fid.write('\\clearpage \n')
-                    savefig(graphname_latex)               
-                    
-            if report == True and len(label_id) > 1:
-                for i in range(nn-1):
-                    if nn > 2:
-                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
-                            word_quantity += 'depth' + ','
-                        else:
-                            word_quantity += plot_quantity[i] + ', '
-                    else:
-                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
-                            word_quantity += 'depth'
-                        else:
-                            word_quantity += plot_quantity[i]
-                    where1 = 0
-                    count1 += 1
-                    index = j*len(plot_quantity)
-                    for which_quantity in plot_quantity:
-                        where1 += 1
-                        s = '\includegraphics' \
-                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
-                            (graphname_report[index], '.png')
-                        index += 1
-                        fid.write(s)
-                        if where1 % 2 == 0:
-                            s = '\\\\ \n'
-                            where1 = 0
-                        else:
-                            s = '& \n'
-                        fid.write(s)
-                word_quantity += ' and ' + plot_quantity[nn-1]            
-                label = 'gauge%s' %(gaugeloc2) 
-                caption = 'Time series for %s at %s location ' \
-                          '(elevation %.2fm)' % \
-                          (word_quantity, locations[k], elev[k])
-                if elev[k] == 0.0:
-                        caption = 'Time series for %s at %s location ' \
-                                  '(elevation %.2fm)' % \
-                                  (word_quantity, locations[k],
-                                   elevations[0,k,j])
-                        thisgauge = gauges[k]
-                        east = thisgauge[0]
-                        north = thisgauge[1]
-                        elev_output.append([locations[k], east, north,
-                                            elevations[0,k,j]])
-                        
-                s = '\end{tabular} \n' \
-                    '\\caption{%s} \n' \
-                    '\label{fig:%s} \n' \
-                    '\end{figure} \n \n' % (caption, label)
-                fid.write(s)
-                if float((k+1)/div - pp) == 0.:
-                    fid.write('\\clearpage \n')
-                    pp += 1
-                #### finished generating figures ###
-
-            close('all')
-        
-    return texfile2, elev_output
 
 

anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

@@ -67 +67 @@
 
           vertex_coordinates: List of coordinate pairs [xi, eta] of
-              points constituting a mesh (or an m x 2 numeric array or
+          points constituting a mesh (or an m x 2 numeric array or
               a geospatial object)
               Points may appear multiple times

anuga_core/source/anuga/utilities/mainland_only.csv

@@ -2 +2 @@
 447203.12,7699910.99
 447515.37,7699565.56
-447802.51,7699594.13
-447897.04,7699314.42
-448096.87,7699309.5
-448049.53,7699137.81
-448241.01,7699138.4
-447998.25,7699230.62
-448202.22,7699229.04
-448076.58,7699475.45
-447944.63,7699405.32
-447823.17,7699641.78
-447701.51,7699610.41
-447805.71,7699904.03
-447882,7700136.68
-448239.36,7700008.3
-448380.84,7700024.23
-448136.17,7700058.89
-447721.58,7700183.77
-447669.07,7699667.86
 447500.56,7699975.01
 447108.86,7700098.85
@@ -36 +18 @@
 449882.5,7702139.37
 449798.8,7701938.79
-449950.91,7701889.44
-449858.65,7701764.11
-450041.48,7701881.96
-449808.01,7701988.62
-449874.32,7702091.75
-450046.73,7702217.33
-450427.05,7702076.79
-450713.91,7701864.04
-450914.68,7701542.56
-450865.57,7701253.56
-450830.65,7701097.4
-450841.67,7700886.05
-450609.91,7700776.91
-450814.8,7700822.89
-450892.05,7701098.69
-450959.46,7701180.79
-450946.17,7701453.01
-451146.1,7701420.39
-451556.06,7701458.1
-451694.32,7701513.84
-451782.48,7701260.65
-452236.82,7701064.95
-452539.29,7700837.82
-452541.07,7700944.07
-452332.44,7701111.71
-451953.31,7701211.34
-451784.01,7701451.02
-451712.63,7701663.3
-451519.51,7701500.05
-451204.36,7701429.42
-451096.7,7701584.05
-450949.31,7701800.54
-450701.11,7701969.14
-450562.26,7702109.29
-450268.03,7702351.9
-449866.18,7702377.26
-449649.11,7702566.98
-449368.93,7702974.52
-449318.54,7703114.93
-449401.98,7703402.93
-449536,7703490.76
-449636.2,7703747.83
-449878.9,7704041.85
-450026.94,7704315.65
-450229.4,7704493.33
-450389.76,7704479.42
-450662.87,7704351.85
-450856.13,7704121.11
-450964.33,7703783.86
-451102.63,7703477.7
-451378.82,7703357.87
-451788.47,7703518.42
-451947.33,7703661.65
-452211.68,7703676.79
-452215.58,7703404.54
-451836.88,7702974.04
-452106.41,7702632.83
-452265.58,7702295.73
 452234.95,7702449.48
 452124.96,7702695.97
@@ -145 +69 @@
 455206.33,7705662.8
 455432.49,7705568.22
-455328.81,7705799.25
-455542.92,7705931.52
-455994.78,7705918.32
-456256.1,7705910.15
-456387.59,7705789.86
-456538.72,7705720.53
-456515.56,7705419.44
-456530.5,7705278.92
-456578.82,7705513.68
-456527.4,7705670.7
-456655,7705842.57
-456363.24,7705945.85
-456639.63,7706150.2
-456878.51,7706373.27
-457188.46,7706488.06
-457472.53,7706139.05
-457782.93,7706071.22
-457529.4,7705878.01
-457412.68,7705515.81
-457074.97,7705273.68
-457279.32,7705151.36
-457463.68,7705110.88
-457604.41,7705033.76
-457437.51,7705169.47
-457233.96,7705384.76
-457552.19,7705512.84
-457676.08,7705923.75
-458180.98,7705929.45
-458540.56,7705758.8
-458546.42,7705494.31
-458514.81,7705228.61
-458816.11,7705046.74
-458783.48,7705192.75
-458484.17,7705410.04
-458712.21,7705389.57
-458976.32,7705515.28
-459128.1,7705600.87
-458840.13,7705439.69
-458713.57,7705683.97
-458529.7,7705942.49
-458272.31,7706041.45
-458040.76,7706214.64
-457889.15,7706477.66
-457710.6,7706681.95
-457623.71,7706874.3
-457755.5,7707041.75
-458081.25,7707092.37
-458067.52,7707170.92
-457663.81,7707071.4
-457565.14,7706977.08
-457478.27,7707157.26
-457775.78,7707262.04
-458245.06,7707370.57
-458501.69,7707588.12
-458731.37,7707758.02
-458848.75,7707870.09
-459056.24,7708172.73
-459281.45,7708474.31
-459417.46,7708638.43
-459761.72,7708796.41
-459800.63,7708634.93
-459653.3,7708419.87
-459920.41,7708619.71
-460181.3,7708814.01
-460300.74,7708941.57
-460473.14,7709130.12
-460650.4,7709095.42
-460376.97,7708408.31
-460330,7708016.42
-460051.17,7707926.11
-459976.52,7707795.34
-460031.85,7707733.5
-460445.72,7707952.5
-460662.87,7707704
-460841.78,7707329.24
-460579.31,7707372.89
-460519.52,7707562
-460347.01,7707424.36
-460582.11,7707510.13
-460583.74,7707261.12
-460815.17,7707132.18
-460931.27,7707350.48
-460817.97,7707715.43
-461029.59,7708051.26
-461072.53,7708395.55
-460945.07,7708581.18
-460970.87,7709126.85
-461133.74,7709389.53
-461269.97,7709467.31
-461394.48,7709657.95
 461705.68,7709721.75
 462681.38,7710184.34
@@ -496 +330 @@
 478255.43,7714831.11
 478127.84,7714613.18
-478268.97,7714032.09
-478416.29,7713672.6
-478259.24,7713504.18
-478105.32,7713332.45
-477593.89,7713369.4
-477803.58,7713107.39
-477759.28,7712741.01
-478042.92,7712468.03
-478124.29,7712364.11
-477873.96,7712651.52
-477877.52,7713114.13
-478139.79,7713255.02
-478421.16,7713121.47
-478507.01,7712777.4
-478477.48,7713061.78
-478315.61,7713405.76
-478521.27,7713840.96
-478291.77,7714115.13
-478290.13,7714582.15
-478434.75,7714706.29
-478887.94,7714624.97
-479148.9,7714167.12
-479227.67,7713637.11
-479325.14,7713144.75
-479127.16,7713221.97
-479404.41,7713050.78
-479509.43,7713188.14
-479662.75,7713004.61
-479850.6,7712682.79
-480148.8,7712394.29
-480355.23,7712206.39
-480446.42,7711972.48
-480574.6,7711576.89
-480263.09,7711542.43
-479801.34,7711470.75
-480151.44,7711255.73
-480675.22,7711463.86
-480935.74,7711377.02
-481274.96,7711141.7
-481648.11,7710828.91
-481866.13,7710472.79
 482021.51,7710263.79
 482078.8,7710247.25
@@ -625 +418 @@
 496055.54,7712910.29
 496465.85,7713002.24
-496605.35,7713285.59
-496814.63,7713583.33
-497167.66,7713809.16
-497289.51,7713874.48
-497134.28,7714116.81
-497252.99,7714279.52
-497506.09,7714258.53
-497594.65,7714015.08
-497996.72,7713721.86
-497856.15,7713478.37
-497844.78,7712887.39
-497806.32,7712316.33
-497615.77,7712118.2
-497856.29,7712414.84
-497913.53,7712783.37
-497913.46,7713306.83
-498066.51,7713686.45
-497809.22,7713948.71
-497824.79,7714289.57
-498351.81,7714378.17
-498580.96,7714271.95
-498520.59,7713890.13
-498933.04,7713824.87
-498678.93,7713583.6
-498747.66,7713600.2
-499010.11,7713746.3
-498782,7713983.12
-498544.53,7714024.05
-498621.57,7714410.29
-498869.47,7714234.34
-499184.03,7713996.43
-499113.22,7713789.47
-499548.58,7713736.37
-499555.86,7713822.69
-499192.37,7713813.82
-499265.27,7714007.5
-498936.13,7714264.23
-498637.18,7714516.53
-498190.35,7714502.1
-497623.76,7714388.04
-497548.74,7714554.03
-497514.33,7714811.88
-497892.43,7714815.25
-498068.43,7714974.64
-498256.93,7715211.5
-498531.91,7715324.41
-498686.05,7715515.88
-498900.61,7715666.4
-499294.35,7715710.7
-499576.63,7715688.57
-499948.49,7715861.23
-500144.32,7716010.63
-500342.23,7715994.03
-500548.48,7716033.87
-500731.82,7716156.71
-500929.73,7716177.73
-501103.7,7716323.8
-501317.24,7716344.82
-501469.34,7716462.11
-501611.01,7716540.67
-501820.39,7716489.75
-502037.06,7716486.4
-502243.31,7716486.38
-502468.32,7716540.58
-502691.25,7716596.99
-502913.14,7716624.62
 503095.45,7716705.38
 503302.75,7716750.71
@@ -845 +572 @@
 517248.64,7721707.83
 517357.62,7721981.01
-454913.9,7687033.4
+530000,7687033.4
+424913.9,7687033.4
 425539.85,7691615.14
 426064.25,7691355.15
@@ -854 +582 @@
 426898.56,7691103.16
 426763.28,7690643.24
-427055.91,7690088.91
-426796.88,7690332.37
-426857.75,7690921.46
-427646.12,7691424.07
-428763.81,7691794.12
-429123.61,7691823.32
-430084.18,7692218.07
-430563.65,7692990.4
-430763.67,7693423.99
 430997.89,7693894.23
 431118.95,7694305.35

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -76 +76 @@
     if v2 is None:
         v2 = [1.0, 0.0] # Unit vector along the x-axis
-        
+    
     v1 = ensure_numeric(v1, num.float)
     v2 = ensure_numeric(v2, num.float)    

anuga_core/source/anuga/utilities/polygon.py

@@ -322 +322 @@
     return False
 
+    
 def is_complex(polygon, verbose=False):
-    """Check if a polygon is complex (self-intersecting)
-    """
-    
+    """Check if a polygon is complex (self-intersecting).
+       Uses a sweep algorithm that is O(n^2) in the worst case, but
+       for most normal looking polygons it'll be O(n log n). 
+
+       polygon is a list of points that define a closed polygon.
+       verbose will print a list of the intersection points if true
+       
+       Return True if polygon is complex.
+    """            
+            
+    def key_xpos(item):
+        return (item[0][0])
+    
+    def segments_joined(seg0, seg1):
+        for i in seg0:
+            for j in seg1:    
+                if i == j: return True
+        return False
+        
     polygon = ensure_numeric(polygon, num.float)
 
+    # build a list of discrete segments from the polygon
+    unsorted_segs = []
     for i in range(0, len(polygon)-1):
-        for j in range(i+1, len(polygon)-1):    
-                (type, point) = intersection([polygon[i], polygon[i+1]], [polygon[j], polygon[j+1]])
-
-                if (abs(i-j) > 1 and type == 1) or (type == 2 and list(point[0]) != list(point[1])) or (type == 3) and type != 4:
+        unsorted_segs.append([list(polygon[i]), list(polygon[i+1])])
+    unsorted_segs.append([list(polygon[0]), list(polygon[-1])])
+    
+    # all segments must point in same direction
+    for val in unsorted_segs:
+        if val[0][0] > val[1][0]:
+            val[0], val[1] = val[1], val[0]    
+            
+    l_x = sorted(unsorted_segs, key=key_xpos)
+
+    comparisons = 0
+    
+    # loop through, only comparing lines that partially overlap in x
+    for index, leftmost in enumerate(l_x):
+        cmp = index+1
+        while cmp < len(l_x) and leftmost[1][0] > l_x[cmp][0][0]:
+            if not segments_joined(leftmost, l_x[cmp]):
+                (type, point) = intersection(leftmost, l_x[cmp])
+                comparisons += 1
+                if type != 0 and type != 4 or (type == 2 and list(point[0]) != list(point[1])):
                     if verbose:
-                        print 'Self-intersecting polygon found, type ', type, ' point', point, 'vertex indices ', i, j                
-                    return True
+                        print 'Self-intersecting polygon found, type ', type, ' point', point,
+                        print 'vertices: ', leftmost, ' - ', l_x[cmp]                
+                    return True            
+            cmp += 1
         
     return False
+    
     
 def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
@@ -996 +1034 @@
         polygon.append([float(fields[0]), float(fields[1])])
     
-    # check this is a valid polygon    
-    # if is_complex(polygon):    
-            # msg = 'ERROR: Self-intersecting polygon detected in file' + filename +'. '
-            # msg += 'Please fix.'
-            # log.critical(msg)
-# #            raise Exception, msg
+    # check this is a valid polygon. 
+    if is_complex(polygon, verbose=True):    
+        msg = 'ERROR: Self-intersecting polygon detected in file ' + filename +'. '
+        msg += 'A complex polygon will not necessarily break the algorithms within ANUGA, '
+        msg += 'but it usually signifies pathological data. Please fix this file.'
+        raise Exception, msg
     
     return polygon

anuga_core/source/anuga/utilities/test_polygon.py

@@ -1808 +1808 @@
         
     def test_is_polygon_complex(self):
+        """Test a concave and a complex poly with is_complex, to make sure it can detect
+           self-intersection.
+        """
         concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
         complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], [0, 10]]
 
-        not_complex = is_complex(concave_poly)        
-        complex = is_complex(complex_poly)        
-
-        assert not not_complex
-        assert complex
+        assert not is_complex(concave_poly)
+        assert is_complex(complex_poly)
 
     def test_is_polygon_complex2(self):
+        """Test a concave and a complex poly with is_complex, to make sure it can detect
+           self-intersection. This test uses more complicated polygons.
+        """    
         concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7,6], [10, 10], [1,5], [0, 10]]
-        complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]
-
-        not_complex = is_complex(concave_poly)        
-        complex = is_complex(complex_poly)        
-
-        assert not not_complex
-        assert complex
+        complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]       
+
+        assert not is_complex(concave_poly)
+        assert is_complex(complex_poly)
         
 ################################################################################