Changeset 7690


Ignore:
Timestamp:
Apr 21, 2010, 1:48:09 PM (10 years ago)
Author:
hudson
Message:

Check for complex polygons raises an exception.

Location:
anuga_core/source/anuga
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/file_function.py

    r7679 r7690  
    1717from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1818from anuga.utilities.numerical_tools import ensure_numeric
     19
     20import anuga.utilities.log as log
     21
    1922
    2023
     
    192195                   verbose=False,
    193196                   boundary_polygon=None,
    194                                    output_centroids=False):
     197                   output_centroids=False):
    195198    """Internal function
    196199   
  • anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

    r7685 r7690  
    1212
    1313from anuga.geospatial_data.geospatial_data import ensure_absolute
    14 from util import check_list, generate_figures
     14from util import check_list
    1515from file_function import file_function
    1616
     
    548548    if len(gauge_index) <> 0:
    549549        texfile, elev_output = \
    550             generate_figures(plot_quantity, file_loc, report, reportname,
     550            _generate_figures(plot_quantity, file_loc, report, reportname,
    551551                             surface, leg_label, f_list, gauges, locations,
    552552                             elev, gauge_index, production_dirs, time_min,
     
    638638
    639639    return gauges, gaugelocation, elev
     640
     641   
     642##
     643# @brief Generate figures from quantities and gauges for each sww file.
     644#        This is a legacy function, used only by sww2timeseries
     645# @param plot_quantity  ??
     646# @param file_loc ??
     647# @param report ??
     648# @param reportname ??
     649# @param surface ??
     650# @param leg_label ??
     651# @param f_list ??
     652# @param gauges ??
     653# @param locations ??
     654# @param elev ??
     655# @param gauge_index ??
     656# @param production_dirs ??
     657# @param time_min ??
     658# @param time_max ??
     659# @param time_unit ??
     660# @param title_on ??
     661# @param label_id ??
     662# @param generate_fig ??
     663# @param verbose??
     664# @return (texfile2, elev_output)
     665def _generate_figures(plot_quantity, file_loc, report, reportname, surface,
     666                     leg_label, f_list, gauges, locations, elev, gauge_index,
     667                     production_dirs, time_min, time_max, time_unit,
     668                     title_on, label_id, generate_fig, verbose):
     669    """ Generate figures based on required quantities and gauges for
     670    each sww file
     671    """
     672    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
     673
     674    if generate_fig is True:
     675        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
     676             xlabel, ylabel, title, close, subplot
     677   
     678        if surface is True:
     679            import pylab as p1
     680            import mpl3d.mplot3d as p3
     681       
     682    if report == True:   
     683        texdir = getcwd()+sep+'report'+sep
     684        if access(texdir,F_OK) == 0:
     685            mkdir (texdir)
     686        if len(label_id) == 1:
     687            label_id1 = label_id[0].replace(sep,'')
     688            label_id2 = label_id1.replace('_','')
     689            texfile = texdir + reportname + '%s' % label_id2
     690            texfile2 = reportname + '%s' % label_id2
     691            texfilename = texfile + '.tex'
     692            fid = open(texfilename, 'w')
     693
     694            if verbose: log.critical('Latex output printed to %s' % texfilename)
     695        else:
     696            texfile = texdir+reportname
     697            texfile2 = reportname
     698            texfilename = texfile + '.tex'
     699            fid = open(texfilename, 'w')
     700
     701            if verbose: log.critical('Latex output printed to %s' % texfilename)
     702    else:
     703        texfile = ''
     704        texfile2 = ''
     705
     706    p = len(f_list)
     707    n = []
     708    n0 = 0
     709    for i in range(len(f_list)):
     710        n.append(len(f_list[i].get_time()))
     711        if n[i] > n0: n0 = n[i] 
     712    n0 = int(n0)
     713    m = len(locations)
     714    model_time = num.zeros((n0, m, p), num.float)
     715    stages = num.zeros((n0, m, p), num.float)
     716    elevations = num.zeros((n0, m, p), num.float)
     717    momenta = num.zeros((n0, m, p), num.float)
     718    xmom = num.zeros((n0, m, p), num.float)
     719    ymom = num.zeros((n0, m, p), num.float)
     720    speed = num.zeros((n0, m, p), num.float)
     721    bearings = num.zeros((n0, m, p), num.float)
     722    due_east = 90.0*num.ones((n0, 1), num.float)
     723    due_west = 270.0*num.ones((n0, 1), num.float)
     724    depths = num.zeros((n0, m, p), num.float)
     725    eastings = num.zeros((n0, m, p), num.float)
     726    min_stages = []
     727    max_stages = []
     728    min_momentums = []   
     729    max_momentums = []
     730    max_xmomentums = []
     731    max_ymomentums = []
     732    min_xmomentums = []
     733    min_ymomentums = []
     734    max_speeds = []
     735    min_speeds = []   
     736    max_depths = []
     737    model_time_plot3d = num.zeros((n0, m), num.float)
     738    stages_plot3d = num.zeros((n0, m), num.float)
     739    eastings_plot3d = num.zeros((n0, m),num.float)
     740    if time_unit is 'mins': scale = 60.0
     741    if time_unit is 'hours': scale = 3600.0
     742
     743    ##### loop over each swwfile #####
     744    for j, f in enumerate(f_list):
     745        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
     746
     747        starttime = f.starttime
     748        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
     749        fid_compare = open(comparefile, 'w')
     750        file0 = file_loc[j] + 'gauges_t0.csv'
     751        fid_0 = open(file0, 'w')
     752
     753        ##### loop over each gauge #####
     754        for k in gauge_index:
     755            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
     756
     757            g = gauges[k]
     758            min_stage = 10
     759            max_stage = 0
     760            max_momentum = max_xmomentum = max_ymomentum = 0
     761            min_momentum = min_xmomentum = min_ymomentum = 100
     762            max_speed = 0
     763            min_speed = 0           
     764            max_depth = 0           
     765            gaugeloc = str(locations[k])
     766            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
     767                       + gaugeloc + '.csv'
     768            if j == 0:
     769                fid_out = open(thisfile, 'w')
     770                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
     771                fid_out.write(s)           
     772
     773            #### generate quantities #######
     774            for i, t in enumerate(f.get_time()):
     775                if time_min <= t <= time_max:
     776                    w = f(t, point_id = k)[0]
     777                    z = f(t, point_id = k)[1]
     778                    uh = f(t, point_id = k)[2]
     779                    vh = f(t, point_id = k)[3]
     780                    depth = w-z     
     781                    m = sqrt(uh*uh + vh*vh)
     782                    if depth < 0.001:
     783                        vel = 0.0
     784                    else:
     785                        vel = m / (depth + 1.e-6/depth)
     786                    bearing = calc_bearing(uh, vh)                   
     787                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
     788                    stages[i,k,j] = w
     789                    elevations[i,k,j] = z
     790                    xmom[i,k,j] = uh
     791                    ymom[i,k,j] = vh
     792                    momenta[i,k,j] = m
     793                    speed[i,k,j] = vel
     794                    bearings[i,k,j] = bearing
     795                    depths[i,k,j] = depth
     796                    thisgauge = gauges[k]
     797                    eastings[i,k,j] = thisgauge[0]
     798                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
     799                            % (t, w, m, vel, z, uh, vh, bearing)
     800                    fid_out.write(s)
     801                    if t == 0:
     802                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
     803                        fid_0.write(s)
     804                    if t/60.0 <= 13920: tindex = i
     805                    if w > max_stage: max_stage = w
     806                    if w < min_stage: min_stage = w
     807                    if m > max_momentum: max_momentum = m
     808                    if m < min_momentum: min_momentum = m                   
     809                    if uh > max_xmomentum: max_xmomentum = uh
     810                    if vh > max_ymomentum: max_ymomentum = vh
     811                    if uh < min_xmomentum: min_xmomentum = uh
     812                    if vh < min_ymomentum: min_ymomentum = vh
     813                    if vel > max_speed: max_speed = vel
     814                    if vel < min_speed: min_speed = vel                   
     815                    if z > 0 and depth > max_depth: max_depth = depth
     816                   
     817                   
     818            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
     819                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
     820            fid_compare.write(s)
     821            max_stages.append(max_stage)
     822            min_stages.append(min_stage)
     823            max_momentums.append(max_momentum)
     824            max_xmomentums.append(max_xmomentum)
     825            max_ymomentums.append(max_ymomentum)
     826            min_xmomentums.append(min_xmomentum)
     827            min_ymomentums.append(min_ymomentum)
     828            min_momentums.append(min_momentum)           
     829            max_depths.append(max_depth)
     830            max_speeds.append(max_speed)
     831            min_speeds.append(min_speed)           
     832            #### finished generating quantities for each swwfile #####
     833       
     834        model_time_plot3d[:,:] = model_time[:,:,j]
     835        stages_plot3d[:,:] = stages[:,:,j]
     836        eastings_plot3d[:,] = eastings[:,:,j]
     837           
     838        if surface is True:
     839            log.critical('Printing surface figure')
     840            for i in range(2):
     841                fig = p1.figure(10)
     842                ax = p3.Axes3D(fig)
     843                if len(gauges) > 80:
     844                    ax.plot_surface(model_time[:,:,j],
     845                                    eastings[:,:,j],
     846                                    stages[:,:,j])
     847                else:
     848                    ax.plot3D(num.ravel(eastings[:,:,j]),
     849                              num.ravel(model_time[:,:,j]),
     850                              num.ravel(stages[:,:,j]))
     851                ax.set_xlabel('time')
     852                ax.set_ylabel('x')
     853                ax.set_zlabel('stage')
     854                fig.add_axes(ax)
     855                p1.show()
     856                surfacefig = 'solution_surface%s' % leg_label[j]
     857                p1.savefig(surfacefig)
     858                p1.close()
     859           
     860    #### finished generating quantities for all swwfiles #####
     861
     862    # x profile for given time
     863    if surface is True:
     864        figure(11)
     865        plot(eastings[tindex,:,j], stages[tindex,:,j])
     866        xlabel('x')
     867        ylabel('stage')
     868        profilefig = 'solution_xprofile'
     869        savefig('profilefig')
     870
     871    elev_output = []
     872    if generate_fig is True:
     873        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
     874                           max(max_depths)*1.1])
     875        stage_axis = axis([starttime/scale, time_max/scale,
     876                           min(min_stages), max(max_stages)*1.1])
     877        vel_axis = axis([starttime/scale, time_max/scale,
     878                         min(min_speeds), max(max_speeds)*1.1])
     879        mom_axis = axis([starttime/scale, time_max/scale,
     880                         min(min_momentums), max(max_momentums)*1.1])
     881        xmom_axis = axis([starttime/scale, time_max/scale,
     882                          min(min_xmomentums), max(max_xmomentums)*1.1])
     883        ymom_axis = axis([starttime/scale, time_max/scale,
     884                          min(min_ymomentums), max(max_ymomentums)*1.1])
     885        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
     886        nn = len(plot_quantity)
     887        no_cols = 2
     888       
     889        if len(label_id) > 1: graphname_report = []
     890        pp = 1
     891        div = 11.
     892        cc = 0
     893        for k in gauge_index:
     894            g = gauges[k]
     895            count1 = 0
     896            if report == True and len(label_id) > 1:
     897                s = '\\begin{figure}[ht] \n' \
     898                    '\\centering \n' \
     899                    '\\begin{tabular}{cc} \n'
     900                fid.write(s)
     901            if len(label_id) > 1: graphname_report = []
     902
     903            #### generate figures for each gauge ####
     904            for j, f in enumerate(f_list):
     905                ion()
     906                hold(True)
     907                count = 0
     908                where1 = 0
     909                where2 = 0
     910                word_quantity = ''
     911                if report == True and len(label_id) == 1:
     912                    s = '\\begin{figure}[hbt] \n' \
     913                        '\\centering \n' \
     914                        '\\begin{tabular}{cc} \n'
     915                    fid.write(s)
     916                   
     917                for which_quantity in plot_quantity:
     918                    count += 1
     919                    where1 += 1
     920                    figure(count, frameon = False)
     921                    if which_quantity == 'depth':
     922                        plot(model_time[0:n[j]-1,k,j],
     923                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
     924                        units = 'm'
     925                        axis(depth_axis)
     926                    if which_quantity == 'stage':
     927                        if elevations[0,k,j] <= 0:
     928                            plot(model_time[0:n[j]-1,k,j],
     929                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
     930                            axis(stage_axis)
     931                        else:
     932                            plot(model_time[0:n[j]-1,k,j],
     933                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
     934                            #axis(depth_axis)                 
     935                        units = 'm'
     936                    if which_quantity == 'momentum':
     937                        plot(model_time[0:n[j]-1,k,j],
     938                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
     939                        axis(mom_axis)
     940                        units = 'm^2 / sec'
     941                    if which_quantity == 'xmomentum':
     942                        plot(model_time[0:n[j]-1,k,j],
     943                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
     944                        axis(xmom_axis)
     945                        units = 'm^2 / sec'
     946                    if which_quantity == 'ymomentum':
     947                        plot(model_time[0:n[j]-1,k,j],
     948                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
     949                        axis(ymom_axis)
     950                        units = 'm^2 / sec'
     951                    if which_quantity == 'speed':
     952                        plot(model_time[0:n[j]-1,k,j],
     953                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
     954                        axis(vel_axis)
     955                        units = 'm / sec'
     956                    if which_quantity == 'bearing':
     957                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
     958                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.',
     959                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
     960                        units = 'degrees from North'
     961                        #ax = axis([time_min, time_max, 0.0, 360.0])
     962                        legend(('Bearing','West','East'))
     963
     964                    if time_unit is 'mins': xlabel('time (mins)')
     965                    if time_unit is 'hours': xlabel('time (hours)')
     966                    #if which_quantity == 'stage' \
     967                    #   and elevations[0:n[j]-1,k,j] > 0:
     968                    #    ylabel('%s (%s)' %('depth', units))
     969                    #else:
     970                    #    ylabel('%s (%s)' %(which_quantity, units))
     971                        #ylabel('%s (%s)' %('wave height', units))
     972                    ylabel('%s (%s)' %(which_quantity, units))
     973                    if len(label_id) > 1: legend((leg_label),loc='upper right')
     974
     975                    #gaugeloc1 = gaugeloc.replace(' ','')
     976                    #gaugeloc2 = gaugeloc1.replace('_','')
     977                    gaugeloc2 = str(locations[k]).replace(' ','')
     978                    graphname = '%sgauge%s_%s' %(file_loc[j],
     979                                                 gaugeloc2,
     980                                                 which_quantity)
     981
     982                    if report == True and len(label_id) > 1:
     983                        figdir = getcwd()+sep+'report_figures'+sep
     984                        if access(figdir,F_OK) == 0 :
     985                            mkdir (figdir)
     986                        latex_file_loc = figdir.replace(sep,altsep)
     987                        # storing files in production directory   
     988                        graphname_latex = '%sgauge%s%s' \
     989                                          % (latex_file_loc, gaugeloc2,
     990                                             which_quantity)
     991                        # giving location in latex output file
     992                        graphname_report_input = '%sgauge%s%s' % \
     993                                                 ('..' + altsep +
     994                                                      'report_figures' + altsep,
     995                                                  gaugeloc2, which_quantity)
     996                        graphname_report.append(graphname_report_input)
     997                       
     998                        # save figures in production directory for report
     999                        savefig(graphname_latex)
     1000
     1001                    if report == True:
     1002                        figdir = getcwd() + sep + 'report_figures' + sep
     1003                        if access(figdir,F_OK) == 0:
     1004                            mkdir(figdir)
     1005                        latex_file_loc = figdir.replace(sep,altsep)   
     1006
     1007                        if len(label_id) == 1:
     1008                            # storing files in production directory 
     1009                            graphname_latex = '%sgauge%s%s%s' % \
     1010                                              (latex_file_loc, gaugeloc2,
     1011                                               which_quantity, label_id2)
     1012                            # giving location in latex output file
     1013                            graphname_report = '%sgauge%s%s%s' % \
     1014                                               ('..' + altsep +
     1015                                                    'report_figures' + altsep,
     1016                                                gaugeloc2, which_quantity,
     1017                                                label_id2)
     1018                            s = '\includegraphics' \
     1019                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
     1020                                (graphname_report, '.png')
     1021                            fid.write(s)
     1022                            if where1 % 2 == 0:
     1023                                s = '\\\\ \n'
     1024                                where1 = 0
     1025                            else:
     1026                                s = '& \n'
     1027                            fid.write(s)
     1028                            savefig(graphname_latex)
     1029                   
     1030                    if title_on == True:
     1031                        title('%s scenario: %s at %s gauge' % \
     1032                              (label_id, which_quantity, gaugeloc2))
     1033                        #title('Gauge %s (MOST elevation %.2f, ' \
     1034                        #      'ANUGA elevation %.2f)' % \
     1035                        #      (gaugeloc2, elevations[10,k,0],
     1036                        #       elevations[10,k,1]))
     1037
     1038                    savefig(graphname) # save figures with sww file
     1039
     1040                if report == True and len(label_id) == 1:
     1041                    for i in range(nn-1):
     1042                        if nn > 2:
     1043                            if plot_quantity[i] == 'stage' \
     1044                               and elevations[0,k,j] > 0:
     1045                                word_quantity += 'depth' + ', '
     1046                            else:
     1047                                word_quantity += plot_quantity[i] + ', '
     1048                        else:
     1049                            if plot_quantity[i] == 'stage' \
     1050                               and elevations[0,k,j] > 0:
     1051                                word_quantity += 'depth' + ', '
     1052                            else:
     1053                                word_quantity += plot_quantity[i]
     1054                       
     1055                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
     1056                        word_quantity += ' and ' + 'depth'
     1057                    else:
     1058                        word_quantity += ' and ' + plot_quantity[nn-1]
     1059                    caption = 'Time series for %s at %s location ' \
     1060                              '(elevation %.2fm)' % \
     1061                              (word_quantity, locations[k], elev[k])
     1062                    if elev[k] == 0.0:
     1063                        caption = 'Time series for %s at %s location ' \
     1064                                  '(elevation %.2fm)' % \
     1065                                  (word_quantity, locations[k],
     1066                                   elevations[0,k,j])
     1067                        east = gauges[0]
     1068                        north = gauges[1]
     1069                        elev_output.append([locations[k], east, north,
     1070                                            elevations[0,k,j]])
     1071                    label = '%sgauge%s' % (label_id2, gaugeloc2)
     1072                    s = '\end{tabular} \n' \
     1073                        '\\caption{%s} \n' \
     1074                        '\label{fig:%s} \n' \
     1075                        '\end{figure} \n \n' % (caption, label)
     1076                    fid.write(s)
     1077                    cc += 1
     1078                    if cc % 6 == 0: fid.write('\\clearpage \n')
     1079                    savefig(graphname_latex)               
     1080                   
     1081            if report == True and len(label_id) > 1:
     1082                for i in range(nn-1):
     1083                    if nn > 2:
     1084                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
     1085                            word_quantity += 'depth' + ','
     1086                        else:
     1087                            word_quantity += plot_quantity[i] + ', '
     1088                    else:
     1089                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
     1090                            word_quantity += 'depth'
     1091                        else:
     1092                            word_quantity += plot_quantity[i]
     1093                    where1 = 0
     1094                    count1 += 1
     1095                    index = j*len(plot_quantity)
     1096                    for which_quantity in plot_quantity:
     1097                        where1 += 1
     1098                        s = '\includegraphics' \
     1099                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
     1100                            (graphname_report[index], '.png')
     1101                        index += 1
     1102                        fid.write(s)
     1103                        if where1 % 2 == 0:
     1104                            s = '\\\\ \n'
     1105                            where1 = 0
     1106                        else:
     1107                            s = '& \n'
     1108                        fid.write(s)
     1109                word_quantity += ' and ' + plot_quantity[nn-1]           
     1110                label = 'gauge%s' %(gaugeloc2)
     1111                caption = 'Time series for %s at %s location ' \
     1112                          '(elevation %.2fm)' % \
     1113                          (word_quantity, locations[k], elev[k])
     1114                if elev[k] == 0.0:
     1115                        caption = 'Time series for %s at %s location ' \
     1116                                  '(elevation %.2fm)' % \
     1117                                  (word_quantity, locations[k],
     1118                                   elevations[0,k,j])
     1119                        thisgauge = gauges[k]
     1120                        east = thisgauge[0]
     1121                        north = thisgauge[1]
     1122                        elev_output.append([locations[k], east, north,
     1123                                            elevations[0,k,j]])
     1124                       
     1125                s = '\end{tabular} \n' \
     1126                    '\\caption{%s} \n' \
     1127                    '\label{fig:%s} \n' \
     1128                    '\end{figure} \n \n' % (caption, label)
     1129                fid.write(s)
     1130                if float((k+1)/div - pp) == 0.:
     1131                    fid.write('\\clearpage \n')
     1132                    pp += 1
     1133                #### finished generating figures ###
     1134
     1135            close('all')
     1136       
     1137    return texfile2, elev_output
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

    r7685 r7690  
    2929    return x+y
    3030
    31 class Test_Util(unittest.TestCase):
     31class Test_Gauge(unittest.TestCase):
    3232    def setUp(self):
    33         pass
    34 
    35     def tearDown(self):
    36         pass
    37 
    38     def test_sww2csv(self):
    3933
    4034        def elevation_function(x, y):
    4135            return -x
    4236       
    43         """Most of this test was copied from test_interpolate
    44         test_interpole_sww2csv
    45        
    46         This is testing the sww2csv_gauges function, by creating a sww file and
    47         then exporting the gauges and checking the results.
    48         """
    49        
    50         # Create mesh
     37        """ Setup for all tests. """
     38       
    5139        mesh_file = tempfile.mktemp(".tsh")   
    5240        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
     
    6048        domain = Domain(mesh_file)
    6149        os.remove(mesh_file)
    62         
    63         domain.default_order=2
    64        
     50 
     51        domain.default_order = 2
     52
    6553        # This test was made before tight_slope_limiters were introduced
    6654        # Since were are testing interpolation values this is OK
    67         domain.tight_slope_limiters = 0
    68        
    69 
     55        domain.tight_slope_limiters = 0
     56               
    7057        # Set some field values
    7158        domain.set_quantity('elevation', elevation_function)
     
    8370        domain.set_quantity('stage', 1.0)
    8471
    85 
    8672        domain.set_name('datatest' + str(time.time()))
    8773        domain.smooth = True
    8874        domain.reduction = mean
    89 
    90 
    91         sww = SWW_file(domain)
    92         sww.store_connectivity()
    93         sww.store_timestep()
    94         domain.set_quantity('stage', 10.0) # This is automatically limited
     75       
     76        self.domain = domain
     77       
     78       
     79    def tearDown(self):
     80        """Called at end of each test."""
     81        if self.sww:
     82            os.remove(self.sww.filename)
     83
     84    def _create_sww(self):
     85        self.sww = SWW_file(self.domain)
     86        self.sww.store_connectivity()
     87        self.sww.store_timestep()
     88        self.domain.set_quantity('stage', 10.0) # This is automatically limited
    9589        # so it will not be less than the elevation
    96         domain.time = 2.
    97         sww.store_timestep()
    98 
     90        self.domain.time = 2.
     91        self.sww.store_timestep()
     92       
     93       
     94    def test_sww2csv(self):
     95
     96        """Most of this test was copied from test_interpolate
     97        test_interpole_sww2csv
     98       
     99        This is testing the sww2csv_gauges function, by creating a sww file and
     100        then exporting the gauges and checking the results.
     101        """
     102       
     103        domain = self.domain
     104        self._create_sww()
     105       
    99106        # test the function
    100107        points = [[5.0,1.],[0.5,2.]]
     
    109116
    110117       
    111         sww2csv_gauges(sww.filename,
     118        sww2csv_gauges(self.sww.filename,
    112119                       points_file,
    113120                       verbose=False,
     
    146153        point1_handle.close()
    147154        point2_handle.close()
    148         #print "sww.filename",sww.filename
    149         os.remove(sww.filename)
    150155        os.remove(points_file)
    151156        os.remove(point1_filename)
     
    159164        import time
    160165        import string
    161 
    162         def elevation_function(x, y):
    163             return -x
    164166       
    165167        """Most of this test was copied from test_interpolate
     
    173175        """
    174176       
    175         # Create mesh
    176         mesh_file = tempfile.mktemp(".tsh")   
    177         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    178         m = Mesh()
    179         m.add_vertices(points)
    180         m.auto_segment()
    181         m.generate_mesh(verbose=False)
    182         m.export_mesh_file(mesh_file)
    183        
    184         # Create shallow water domain
    185         domain = Domain(mesh_file)
    186         os.remove(mesh_file)
    187        
    188         domain.default_order=2
    189 
    190         # Set some field values
    191         domain.set_quantity('elevation', elevation_function)
    192         domain.set_quantity('friction', 0.03)
    193         domain.set_quantity('xmomentum', 3.0)
    194         domain.set_quantity('ymomentum', 4.0)
    195 
    196         ######################
    197         # Boundary conditions
    198         B = Transmissive_boundary(domain)
    199         domain.set_boundary( {'exterior': B})
    200 
    201         # This call mangles the stage values.
    202         domain.distribute_to_vertices_and_edges()
    203         domain.set_quantity('stage', 1.0)
    204 
    205 
    206         domain.set_name('datatest' + str(time.time()))
    207         domain.smooth = True
    208         domain.reduction = mean
    209 
    210         sww = SWW_file(domain)
    211         sww.store_connectivity()
    212         sww.store_timestep()
    213         domain.set_quantity('stage', 10.0) # This is automatically limited
    214         # so it will not be less than the elevation
    215         domain.time = 2.
    216         sww.store_timestep()
    217 
     177        domain = self.domain
     178        self._create_sww()
     179       
    218180        # test the function
    219181        points = [[5.0,1.],[0.5,2.]]
     
    227189        file_id.close()
    228190
    229         sww2csv_gauges(sww.filename,
     191        sww2csv_gauges(self.sww.filename,
    230192                            points_file,
    231193                            quantities=['stage', 'elevation'],
     
    263225        # clean up
    264226        point1_handle.close()
    265         point2_handle.close()
    266         #print "sww.filename",sww.filename
    267         os.remove(sww.filename)
     227        point2_handle.close()
    268228        os.remove(points_file)
    269229        os.remove(point1_filename)
    270         os.remove(point2_filename)
    271 
     230        os.remove(point2_filename)       
     231       
    272232
    273233    def test_sww2csv_gauges2(self):
    274 
    275         def elevation_function(x, y):
    276             return -x
    277234       
    278235        """Most of this test was copied from test_interpolate
     
    286243        """
    287244       
    288         # Create mesh
    289         mesh_file = tempfile.mktemp(".tsh")   
    290         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    291         m = Mesh()
    292         m.add_vertices(points)
    293         m.auto_segment()
    294         m.generate_mesh(verbose=False)
    295         m.export_mesh_file(mesh_file)
    296        
    297         # Create shallow water domain
    298         domain = Domain(mesh_file)
    299         os.remove(mesh_file)
    300        
    301         domain.default_order=2
    302 
    303         # This test was made before tight_slope_limiters were introduced
    304         # Since were are testing interpolation values this is OK
    305         domain.tight_slope_limiters = 0         
    306 
    307         # Set some field values
    308         domain.set_quantity('elevation', elevation_function)
    309         domain.set_quantity('friction', 0.03)
    310         domain.set_quantity('xmomentum', 3.0)
    311         domain.set_quantity('ymomentum', 4.0)
     245        domain = self.domain
    312246        domain.set_starttime(5)
    313 
    314         ######################
    315         # Boundary conditions
    316         B = Transmissive_boundary(domain)
    317         domain.set_boundary( {'exterior': B})
    318 
    319         # This call mangles the stage values.
    320         domain.distribute_to_vertices_and_edges()
    321         domain.set_quantity('stage', 1.0)
    322        
    323 
    324 
    325         domain.set_name('datatest' + str(time.time()))
    326         domain.smooth = True
    327         domain.reduction = mean
    328 
    329         sww = SWW_file(domain)
    330         sww.store_connectivity()
    331         sww.store_timestep()
    332         domain.set_quantity('stage', 10.0) # This is automatically limited
    333         # so it will not be less than the elevation
    334         domain.time = 2.
    335         sww.store_timestep()
    336 
     247        self._create_sww()
     248       
    337249        # test the function
    338250        points = [[5.0,1.],[0.5,2.]]
     
    346258        file_id.close()
    347259
    348        
    349         sww2csv_gauges(sww.filename,
     260        sww2csv_gauges(self.sww.filename,
    350261                            points_file,
    351262                            verbose=False,
     
    364275            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    365276                         float(row[4]), float(row[5]), float(row[6])])
    366             #print 'assert line',line[i],'point1',point1_answers_array[i]
     277            #print 'assert line',line[i],'answer',point1_answers_array[i]
    367278            assert num.allclose(line[i], point1_answers_array[i])
    368279
     
    384295        point1_handle.close()
    385296        point2_handle.close()
    386         #print "sww.filename",sww.filename
    387         os.remove(sww.filename)
    388297        os.remove(points_file)
    389298        os.remove(point1_filename)
    390299        os.remove(point2_filename)
    391        
    392 
     300
     301
     302       
     303    # def test_sww2csv_gauge_point_off_mesh(self):
     304        # from anuga.pmesh.mesh import Mesh
     305        # from anuga.shallow_water import Domain, Transmissive_boundary
     306        # from csv import reader,writer
     307        # import time
     308        # import string
     309       
     310        # """Most of this test was copied from test_interpolate
     311        # test_interpole_sww2csv
     312       
     313        # This is testing the sww2csv_gauges function with one gauge off the mesh, by creating a sww file and
     314        # then exporting the gauges and checking the results.
     315       
     316        # This tests the correct values for when a gauge is off the mesh, which is important for parallel.
     317        # """
     318       
     319        # domain = self.domain
     320        # self._create_sww()     
     321     
     322        # # test the function
     323        # points = [[50.0,1.],[50.5,-20.25]]
     324
     325# #        points_file = tempfile.mktemp(".csv")
     326        # points_file = 'test_point.csv'
     327        # file_id = open(points_file,"w")
     328        # file_id.write("name,easting,northing \n\
     329# point1, 50.0, 1.0\n\
     330# point2, 50.5, 20.25\n")
     331        # file_id.close()
     332
     333        # sww2csv_gauges(sww.filename,
     334                            # points_file,
     335                            # quantities=['stage', 'elevation'],
     336                            # use_cache=False,
     337                            # verbose=False)
     338
     339        # point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
     340        # point1_filename = 'gauge_point1.csv'
     341        # point1_handle = file(point1_filename)
     342        # point1_reader = reader(point1_handle)
     343        # point1_reader.next()
     344
     345        # line=[]
     346        # for i,row in enumerate(point1_reader):
     347# #            print 'i',i,'row',row
     348            # # note the 'hole' (element 1) below - skip the new 'hours' field
     349            # line.append([float(row[0]),float(row[2]),float(row[3])])
     350            # #print 'line',line[i],'point1',point1_answers_array[i]
     351            # assert num.allclose(line[i], point1_answers_array[i])
     352
     353        # point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
     354        # point2_filename = 'gauge_point2.csv'
     355        # point2_handle = file(point2_filename)
     356        # point2_reader = reader(point2_handle)
     357        # point2_reader.next()
     358                       
     359        # line=[]
     360        # for i,row in enumerate(point2_reader):
     361# #            print 'i',i,'row',row
     362            # # note the 'hole' (element 1) below - skip the new 'hours' field
     363            # line.append([float(row[0]),float(row[2]),float(row[3])])
     364# #            print 'line',line[i],'point1',point1_answers_array[i]
     365            # assert num.allclose(line[i], point2_answers_array[i])
     366                         
     367        # # clean up
     368        # point1_handle.close()
     369        # point2_handle.close()
     370        # os.remove(points_file)
     371# #        os.remove(point1_filename)
     372# #        os.remove(point2_filename)
     373       
     374       
    393375    def test_sww2csv_centroid(self):
    394 
    395         def elevation_function(x, y):
    396             return -x
    397376       
    398377        """Check sww2csv timeseries at centroid.
     
    402381        """
    403382       
    404         # Create rectangular mesh
    405         mesh_file = tempfile.mktemp(".tsh")   
    406         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    407         m = Mesh()
    408         m.add_vertices(points)
    409         m.auto_segment()
    410         m.generate_mesh(verbose=False)
    411         m.export_mesh_file(mesh_file)
    412        
    413         # Create shallow water domain
    414         domain = Domain(mesh_file)
    415        
    416         domain.default_order=2
    417        
    418         # This test was made before tight_slope_limiters were introduced
    419         # Since were are testing interpolation values this is OK
    420         domain.tight_slope_limiters = 0
    421        
    422 
    423         # Set some field values
    424         domain.set_quantity('elevation', elevation_function)
    425         domain.set_quantity('friction', 0.03)
    426         domain.set_quantity('xmomentum', 3.0)
    427         domain.set_quantity('ymomentum', 4.0)
    428 
    429         ######################
    430         # Boundary conditions
    431         B = Transmissive_boundary(domain)
    432         domain.set_boundary( {'exterior': B})
    433 
    434         # This call mangles the stage values.
    435         domain.distribute_to_vertices_and_edges()
    436         domain.set_quantity('stage', 1.0)
    437 
    438 
    439         domain.set_name('datatest' + str(time.time()))
    440         domain.smooth = True
    441         domain.reduction = mean
    442 
    443 
    444         sww = SWW_file(domain)
    445         sww.store_connectivity()
    446         sww.store_timestep()
    447         domain.set_quantity('stage', 10.0) # This is automatically limited
    448         # so it will not be less than the elevation
    449         domain.time = 2.
    450         sww.store_timestep()
    451 
     383        domain = self.domain
     384        sww = self._create_sww()
     385       
    452386        # create a csv file containing our gauge points
    453387        points_file = tempfile.mktemp(".csv")
     
    466400        file_id.close()
    467401
    468         sww2csv_gauges(sww.filename,
     402        sww2csv_gauges(self.sww.filename,
    469403                       points_file,
    470404                       verbose=False,
     
    501435        point1_handle.close()
    502436        point2_handle.close()
    503         os.remove(mesh_file)       
    504         os.remove(sww.filename)
    505437        os.remove(points_file)
    506438        os.remove(point1_filename)
     
    509441
    510442    def test_sww2csv_output_centroid_attribute(self):
    511 
    512         def elevation_function(x, y):
    513             return -x
    514443       
    515444        """Check sww2csv timeseries at centroid, then output the centroid coordinates.
     
    519448        """
    520449       
    521         # Create rectangular mesh
    522         mesh_file = tempfile.mktemp(".tsh")   
    523         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    524         m = Mesh()
    525         m.add_vertices(points)
    526         m.auto_segment()
    527         m.generate_mesh(verbose=False)
    528         m.export_mesh_file(mesh_file)
    529        
    530         # Create shallow water domain
    531         domain = Domain(mesh_file)
    532        
    533         domain.default_order=2
    534        
    535         # This test was made before tight_slope_limiters were introduced
    536         # Since were are testing interpolation values this is OK
    537         domain.tight_slope_limiters = 0
    538        
    539 
    540         # Set some field values
    541         domain.set_quantity('elevation', elevation_function)
    542         domain.set_quantity('friction', 0.03)
    543 
    544         ######################
    545         # Boundary conditions
    546         B = Transmissive_boundary(domain)
    547         domain.set_boundary( {'exterior': B})
    548 
    549         # This call mangles the stage values.
    550         domain.distribute_to_vertices_and_edges()
    551         domain.set_quantity('stage', 1.0)
    552 
    553 
    554         domain.set_name('datatest' + str(time.time()))
    555         domain.smooth = True
    556         domain.reduction = mean
    557 
    558 
    559         sww = SWW_file(domain)
    560         sww.store_connectivity()
    561         sww.store_timestep()
    562         domain.set_quantity('stage', 10.0) # This is automatically limited
    563         # so it will not be less than the elevation
    564         domain.time = 2.
    565         sww.store_timestep()
    566 
     450        domain = self.domain       
     451        self._create_sww()
     452       
    567453        # create a csv file containing our gauge points
    568454        points_file = tempfile.mktemp(".csv")
     
    575461        file_id.close()
    576462
    577         sww2csv_gauges(sww.filename,
     463        sww2csv_gauges(self.sww.filename,
    578464                       points_file,
    579465                       quantities=['stage', 'xcentroid', 'ycentroid'],
     
    595481
    596482        # clean up
    597         point1_handle.close()
    598         os.remove(mesh_file)       
    599         os.remove(sww.filename)
     483        point1_handle.close()       
    600484        os.remove(points_file)
    601485        os.remove(point1_filename)
     
    606490
    607491if __name__ == "__main__":
    608     suite = unittest.makeSuite(Test_Util, 'test')
     492    suite = unittest.makeSuite(Test_Gauge, 'test')
    609493#    runner = unittest.TextTestRunner(verbosity=2)
    610494    runner = unittest.TextTestRunner(verbosity=1)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r7685 r7690  
    317317
    318318    return bearing
    319 
    320 
    321 ##
    322 # @brief Generate figures from quantities and gauges for each sww file.
    323 # @param plot_quantity  ??
    324 # @param file_loc ??
    325 # @param report ??
    326 # @param reportname ??
    327 # @param surface ??
    328 # @param leg_label ??
    329 # @param f_list ??
    330 # @param gauges ??
    331 # @param locations ??
    332 # @param elev ??
    333 # @param gauge_index ??
    334 # @param production_dirs ??
    335 # @param time_min ??
    336 # @param time_max ??
    337 # @param time_unit ??
    338 # @param title_on ??
    339 # @param label_id ??
    340 # @param generate_fig ??
    341 # @param verbose??
    342 # @return (texfile2, elev_output)
    343 def generate_figures(plot_quantity, file_loc, report, reportname, surface,
    344                      leg_label, f_list, gauges, locations, elev, gauge_index,
    345                      production_dirs, time_min, time_max, time_unit,
    346                      title_on, label_id, generate_fig, verbose):
    347     """ Generate figures based on required quantities and gauges for
    348     each sww file
    349     """
    350     from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
    351 
    352     if generate_fig is True:
    353         from pylab import ion, hold, plot, axis, figure, legend, savefig, \
    354              xlabel, ylabel, title, close, subplot
    355    
    356         if surface is True:
    357             import pylab as p1
    358             import mpl3d.mplot3d as p3
    359        
    360     if report == True:   
    361         texdir = getcwd()+sep+'report'+sep
    362         if access(texdir,F_OK) == 0:
    363             mkdir (texdir)
    364         if len(label_id) == 1:
    365             label_id1 = label_id[0].replace(sep,'')
    366             label_id2 = label_id1.replace('_','')
    367             texfile = texdir + reportname + '%s' % label_id2
    368             texfile2 = reportname + '%s' % label_id2
    369             texfilename = texfile + '.tex'
    370             fid = open(texfilename, 'w')
    371 
    372             if verbose: log.critical('Latex output printed to %s' % texfilename)
    373         else:
    374             texfile = texdir+reportname
    375             texfile2 = reportname
    376             texfilename = texfile + '.tex'
    377             fid = open(texfilename, 'w')
    378 
    379             if verbose: log.critical('Latex output printed to %s' % texfilename)
    380     else:
    381         texfile = ''
    382         texfile2 = ''
    383 
    384     p = len(f_list)
    385     n = []
    386     n0 = 0
    387     for i in range(len(f_list)):
    388         n.append(len(f_list[i].get_time()))
    389         if n[i] > n0: n0 = n[i] 
    390     n0 = int(n0)
    391     m = len(locations)
    392     model_time = num.zeros((n0, m, p), num.float)
    393     stages = num.zeros((n0, m, p), num.float)
    394     elevations = num.zeros((n0, m, p), num.float)
    395     momenta = num.zeros((n0, m, p), num.float)
    396     xmom = num.zeros((n0, m, p), num.float)
    397     ymom = num.zeros((n0, m, p), num.float)
    398     speed = num.zeros((n0, m, p), num.float)
    399     bearings = num.zeros((n0, m, p), num.float)
    400     due_east = 90.0*num.ones((n0, 1), num.float)
    401     due_west = 270.0*num.ones((n0, 1), num.float)
    402     depths = num.zeros((n0, m, p), num.float)
    403     eastings = num.zeros((n0, m, p), num.float)
    404     min_stages = []
    405     max_stages = []
    406     min_momentums = []   
    407     max_momentums = []
    408     max_xmomentums = []
    409     max_ymomentums = []
    410     min_xmomentums = []
    411     min_ymomentums = []
    412     max_speeds = []
    413     min_speeds = []   
    414     max_depths = []
    415     model_time_plot3d = num.zeros((n0, m), num.float)
    416     stages_plot3d = num.zeros((n0, m), num.float)
    417     eastings_plot3d = num.zeros((n0, m),num.float)
    418     if time_unit is 'mins': scale = 60.0
    419     if time_unit is 'hours': scale = 3600.0
    420 
    421     ##### loop over each swwfile #####
    422     for j, f in enumerate(f_list):
    423         if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
    424 
    425         starttime = f.starttime
    426         comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
    427         fid_compare = open(comparefile, 'w')
    428         file0 = file_loc[j] + 'gauges_t0.csv'
    429         fid_0 = open(file0, 'w')
    430 
    431         ##### loop over each gauge #####
    432         for k in gauge_index:
    433             if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
    434 
    435             g = gauges[k]
    436             min_stage = 10
    437             max_stage = 0
    438             max_momentum = max_xmomentum = max_ymomentum = 0
    439             min_momentum = min_xmomentum = min_ymomentum = 100
    440             max_speed = 0
    441             min_speed = 0           
    442             max_depth = 0           
    443             gaugeloc = str(locations[k])
    444             thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
    445                        + gaugeloc + '.csv'
    446             if j == 0:
    447                 fid_out = open(thisfile, 'w')
    448                 s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
    449                 fid_out.write(s)           
    450 
    451             #### generate quantities #######
    452             for i, t in enumerate(f.get_time()):
    453                 if time_min <= t <= time_max:
    454                     w = f(t, point_id = k)[0]
    455                     z = f(t, point_id = k)[1]
    456                     uh = f(t, point_id = k)[2]
    457                     vh = f(t, point_id = k)[3]
    458                     depth = w-z     
    459                     m = sqrt(uh*uh + vh*vh)
    460                     if depth < 0.001:
    461                         vel = 0.0
    462                     else:
    463                         vel = m / (depth + 1.e-6/depth)
    464                     bearing = calc_bearing(uh, vh)                   
    465                     model_time[i,k,j] = (t + starttime)/scale #t/60.0
    466                     stages[i,k,j] = w
    467                     elevations[i,k,j] = z
    468                     xmom[i,k,j] = uh
    469                     ymom[i,k,j] = vh
    470                     momenta[i,k,j] = m
    471                     speed[i,k,j] = vel
    472                     bearings[i,k,j] = bearing
    473                     depths[i,k,j] = depth
    474                     thisgauge = gauges[k]
    475                     eastings[i,k,j] = thisgauge[0]
    476                     s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
    477                             % (t, w, m, vel, z, uh, vh, bearing)
    478                     fid_out.write(s)
    479                     if t == 0:
    480                         s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
    481                         fid_0.write(s)
    482                     if t/60.0 <= 13920: tindex = i
    483                     if w > max_stage: max_stage = w
    484                     if w < min_stage: min_stage = w
    485                     if m > max_momentum: max_momentum = m
    486                     if m < min_momentum: min_momentum = m                   
    487                     if uh > max_xmomentum: max_xmomentum = uh
    488                     if vh > max_ymomentum: max_ymomentum = vh
    489                     if uh < min_xmomentum: min_xmomentum = uh
    490                     if vh < min_ymomentum: min_ymomentum = vh
    491                     if vel > max_speed: max_speed = vel
    492                     if vel < min_speed: min_speed = vel                   
    493                     if z > 0 and depth > max_depth: max_depth = depth
    494                    
    495                    
    496             s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
    497                     % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
    498             fid_compare.write(s)
    499             max_stages.append(max_stage)
    500             min_stages.append(min_stage)
    501             max_momentums.append(max_momentum)
    502             max_xmomentums.append(max_xmomentum)
    503             max_ymomentums.append(max_ymomentum)
    504             min_xmomentums.append(min_xmomentum)
    505             min_ymomentums.append(min_ymomentum)
    506             min_momentums.append(min_momentum)           
    507             max_depths.append(max_depth)
    508             max_speeds.append(max_speed)
    509             min_speeds.append(min_speed)           
    510             #### finished generating quantities for each swwfile #####
    511        
    512         model_time_plot3d[:,:] = model_time[:,:,j]
    513         stages_plot3d[:,:] = stages[:,:,j]
    514         eastings_plot3d[:,] = eastings[:,:,j]
    515            
    516         if surface is True:
    517             log.critical('Printing surface figure')
    518             for i in range(2):
    519                 fig = p1.figure(10)
    520                 ax = p3.Axes3D(fig)
    521                 if len(gauges) > 80:
    522                     ax.plot_surface(model_time[:,:,j],
    523                                     eastings[:,:,j],
    524                                     stages[:,:,j])
    525                 else:
    526                     ax.plot3D(num.ravel(eastings[:,:,j]),
    527                               num.ravel(model_time[:,:,j]),
    528                               num.ravel(stages[:,:,j]))
    529                 ax.set_xlabel('time')
    530                 ax.set_ylabel('x')
    531                 ax.set_zlabel('stage')
    532                 fig.add_axes(ax)
    533                 p1.show()
    534                 surfacefig = 'solution_surface%s' % leg_label[j]
    535                 p1.savefig(surfacefig)
    536                 p1.close()
    537            
    538     #### finished generating quantities for all swwfiles #####
    539 
    540     # x profile for given time
    541     if surface is True:
    542         figure(11)
    543         plot(eastings[tindex,:,j], stages[tindex,:,j])
    544         xlabel('x')
    545         ylabel('stage')
    546         profilefig = 'solution_xprofile'
    547         savefig('profilefig')
    548 
    549     elev_output = []
    550     if generate_fig is True:
    551         depth_axis = axis([starttime/scale, time_max/scale, -0.1,
    552                            max(max_depths)*1.1])
    553         stage_axis = axis([starttime/scale, time_max/scale,
    554                            min(min_stages), max(max_stages)*1.1])
    555         vel_axis = axis([starttime/scale, time_max/scale,
    556                          min(min_speeds), max(max_speeds)*1.1])
    557         mom_axis = axis([starttime/scale, time_max/scale,
    558                          min(min_momentums), max(max_momentums)*1.1])
    559         xmom_axis = axis([starttime/scale, time_max/scale,
    560                           min(min_xmomentums), max(max_xmomentums)*1.1])
    561         ymom_axis = axis([starttime/scale, time_max/scale,
    562                           min(min_ymomentums), max(max_ymomentums)*1.1])
    563         cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
    564         nn = len(plot_quantity)
    565         no_cols = 2
    566        
    567         if len(label_id) > 1: graphname_report = []
    568         pp = 1
    569         div = 11.
    570         cc = 0
    571         for k in gauge_index:
    572             g = gauges[k]
    573             count1 = 0
    574             if report == True and len(label_id) > 1:
    575                 s = '\\begin{figure}[ht] \n' \
    576                     '\\centering \n' \
    577                     '\\begin{tabular}{cc} \n'
    578                 fid.write(s)
    579             if len(label_id) > 1: graphname_report = []
    580 
    581             #### generate figures for each gauge ####
    582             for j, f in enumerate(f_list):
    583                 ion()
    584                 hold(True)
    585                 count = 0
    586                 where1 = 0
    587                 where2 = 0
    588                 word_quantity = ''
    589                 if report == True and len(label_id) == 1:
    590                     s = '\\begin{figure}[hbt] \n' \
    591                         '\\centering \n' \
    592                         '\\begin{tabular}{cc} \n'
    593                     fid.write(s)
    594                    
    595                 for which_quantity in plot_quantity:
    596                     count += 1
    597                     where1 += 1
    598                     figure(count, frameon = False)
    599                     if which_quantity == 'depth':
    600                         plot(model_time[0:n[j]-1,k,j],
    601                              depths[0:n[j]-1,k,j], '-', c = cstr[j])
    602                         units = 'm'
    603                         axis(depth_axis)
    604                     if which_quantity == 'stage':
    605                         if elevations[0,k,j] <= 0:
    606                             plot(model_time[0:n[j]-1,k,j],
    607                                  stages[0:n[j]-1,k,j], '-', c = cstr[j])
    608                             axis(stage_axis)
    609                         else:
    610                             plot(model_time[0:n[j]-1,k,j],
    611                                  depths[0:n[j]-1,k,j], '-', c = cstr[j])
    612                             #axis(depth_axis)                 
    613                         units = 'm'
    614                     if which_quantity == 'momentum':
    615                         plot(model_time[0:n[j]-1,k,j],
    616                              momenta[0:n[j]-1,k,j], '-', c = cstr[j])
    617                         axis(mom_axis)
    618                         units = 'm^2 / sec'
    619                     if which_quantity == 'xmomentum':
    620                         plot(model_time[0:n[j]-1,k,j],
    621                              xmom[0:n[j]-1,k,j], '-', c = cstr[j])
    622                         axis(xmom_axis)
    623                         units = 'm^2 / sec'
    624                     if which_quantity == 'ymomentum':
    625                         plot(model_time[0:n[j]-1,k,j],
    626                              ymom[0:n[j]-1,k,j], '-', c = cstr[j])
    627                         axis(ymom_axis)
    628                         units = 'm^2 / sec'
    629                     if which_quantity == 'speed':
    630                         plot(model_time[0:n[j]-1,k,j],
    631                              speed[0:n[j]-1,k,j], '-', c = cstr[j])
    632                         axis(vel_axis)
    633                         units = 'm / sec'
    634                     if which_quantity == 'bearing':
    635                         plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
    636                              model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.',
    637                              model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
    638                         units = 'degrees from North'
    639                         #ax = axis([time_min, time_max, 0.0, 360.0])
    640                         legend(('Bearing','West','East'))
    641 
    642                     if time_unit is 'mins': xlabel('time (mins)')
    643                     if time_unit is 'hours': xlabel('time (hours)')
    644                     #if which_quantity == 'stage' \
    645                     #   and elevations[0:n[j]-1,k,j] > 0:
    646                     #    ylabel('%s (%s)' %('depth', units))
    647                     #else:
    648                     #    ylabel('%s (%s)' %(which_quantity, units))
    649                         #ylabel('%s (%s)' %('wave height', units))
    650                     ylabel('%s (%s)' %(which_quantity, units))
    651                     if len(label_id) > 1: legend((leg_label),loc='upper right')
    652 
    653                     #gaugeloc1 = gaugeloc.replace(' ','')
    654                     #gaugeloc2 = gaugeloc1.replace('_','')
    655                     gaugeloc2 = str(locations[k]).replace(' ','')
    656                     graphname = '%sgauge%s_%s' %(file_loc[j],
    657                                                  gaugeloc2,
    658                                                  which_quantity)
    659 
    660                     if report == True and len(label_id) > 1:
    661                         figdir = getcwd()+sep+'report_figures'+sep
    662                         if access(figdir,F_OK) == 0 :
    663                             mkdir (figdir)
    664                         latex_file_loc = figdir.replace(sep,altsep)
    665                         # storing files in production directory   
    666                         graphname_latex = '%sgauge%s%s' \
    667                                           % (latex_file_loc, gaugeloc2,
    668                                              which_quantity)
    669                         # giving location in latex output file
    670                         graphname_report_input = '%sgauge%s%s' % \
    671                                                  ('..' + altsep +
    672                                                       'report_figures' + altsep,
    673                                                   gaugeloc2, which_quantity)
    674                         graphname_report.append(graphname_report_input)
    675                        
    676                         # save figures in production directory for report
    677                         savefig(graphname_latex)
    678 
    679                     if report == True:
    680                         figdir = getcwd() + sep + 'report_figures' + sep
    681                         if access(figdir,F_OK) == 0:
    682                             mkdir(figdir)
    683                         latex_file_loc = figdir.replace(sep,altsep)   
    684 
    685                         if len(label_id) == 1:
    686                             # storing files in production directory 
    687                             graphname_latex = '%sgauge%s%s%s' % \
    688                                               (latex_file_loc, gaugeloc2,
    689                                                which_quantity, label_id2)
    690                             # giving location in latex output file
    691                             graphname_report = '%sgauge%s%s%s' % \
    692                                                ('..' + altsep +
    693                                                     'report_figures' + altsep,
    694                                                 gaugeloc2, which_quantity,
    695                                                 label_id2)
    696                             s = '\includegraphics' \
    697                                 '[width=0.49\linewidth, height=50mm]{%s%s}' % \
    698                                 (graphname_report, '.png')
    699                             fid.write(s)
    700                             if where1 % 2 == 0:
    701                                 s = '\\\\ \n'
    702                                 where1 = 0
    703                             else:
    704                                 s = '& \n'
    705                             fid.write(s)
    706                             savefig(graphname_latex)
    707                    
    708                     if title_on == True:
    709                         title('%s scenario: %s at %s gauge' % \
    710                               (label_id, which_quantity, gaugeloc2))
    711                         #title('Gauge %s (MOST elevation %.2f, ' \
    712                         #      'ANUGA elevation %.2f)' % \
    713                         #      (gaugeloc2, elevations[10,k,0],
    714                         #       elevations[10,k,1]))
    715 
    716                     savefig(graphname) # save figures with sww file
    717 
    718                 if report == True and len(label_id) == 1:
    719                     for i in range(nn-1):
    720                         if nn > 2:
    721                             if plot_quantity[i] == 'stage' \
    722                                and elevations[0,k,j] > 0:
    723                                 word_quantity += 'depth' + ', '
    724                             else:
    725                                 word_quantity += plot_quantity[i] + ', '
    726                         else:
    727                             if plot_quantity[i] == 'stage' \
    728                                and elevations[0,k,j] > 0:
    729                                 word_quantity += 'depth' + ', '
    730                             else:
    731                                 word_quantity += plot_quantity[i]
    732                        
    733                     if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
    734                         word_quantity += ' and ' + 'depth'
    735                     else:
    736                         word_quantity += ' and ' + plot_quantity[nn-1]
    737                     caption = 'Time series for %s at %s location ' \
    738                               '(elevation %.2fm)' % \
    739                               (word_quantity, locations[k], elev[k])
    740                     if elev[k] == 0.0:
    741                         caption = 'Time series for %s at %s location ' \
    742                                   '(elevation %.2fm)' % \
    743                                   (word_quantity, locations[k],
    744                                    elevations[0,k,j])
    745                         east = gauges[0]
    746                         north = gauges[1]
    747                         elev_output.append([locations[k], east, north,
    748                                             elevations[0,k,j]])
    749                     label = '%sgauge%s' % (label_id2, gaugeloc2)
    750                     s = '\end{tabular} \n' \
    751                         '\\caption{%s} \n' \
    752                         '\label{fig:%s} \n' \
    753                         '\end{figure} \n \n' % (caption, label)
    754                     fid.write(s)
    755                     cc += 1
    756                     if cc % 6 == 0: fid.write('\\clearpage \n')
    757                     savefig(graphname_latex)               
    758                    
    759             if report == True and len(label_id) > 1:
    760                 for i in range(nn-1):
    761                     if nn > 2:
    762                         if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
    763                             word_quantity += 'depth' + ','
    764                         else:
    765                             word_quantity += plot_quantity[i] + ', '
    766                     else:
    767                         if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
    768                             word_quantity += 'depth'
    769                         else:
    770                             word_quantity += plot_quantity[i]
    771                     where1 = 0
    772                     count1 += 1
    773                     index = j*len(plot_quantity)
    774                     for which_quantity in plot_quantity:
    775                         where1 += 1
    776                         s = '\includegraphics' \
    777                             '[width=0.49\linewidth, height=50mm]{%s%s}' % \
    778                             (graphname_report[index], '.png')
    779                         index += 1
    780                         fid.write(s)
    781                         if where1 % 2 == 0:
    782                             s = '\\\\ \n'
    783                             where1 = 0
    784                         else:
    785                             s = '& \n'
    786                         fid.write(s)
    787                 word_quantity += ' and ' + plot_quantity[nn-1]           
    788                 label = 'gauge%s' %(gaugeloc2)
    789                 caption = 'Time series for %s at %s location ' \
    790                           '(elevation %.2fm)' % \
    791                           (word_quantity, locations[k], elev[k])
    792                 if elev[k] == 0.0:
    793                         caption = 'Time series for %s at %s location ' \
    794                                   '(elevation %.2fm)' % \
    795                                   (word_quantity, locations[k],
    796                                    elevations[0,k,j])
    797                         thisgauge = gauges[k]
    798                         east = thisgauge[0]
    799                         north = thisgauge[1]
    800                         elev_output.append([locations[k], east, north,
    801                                             elevations[0,k,j]])
    802                        
    803                 s = '\end{tabular} \n' \
    804                     '\\caption{%s} \n' \
    805                     '\label{fig:%s} \n' \
    806                     '\end{figure} \n \n' % (caption, label)
    807                 fid.write(s)
    808                 if float((k+1)/div - pp) == 0.:
    809                     fid.write('\\clearpage \n')
    810                     pp += 1
    811                 #### finished generating figures ###
    812 
    813             close('all')
    814        
    815     return texfile2, elev_output
    816319
    817320
  • anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

    r7685 r7690  
    6767
    6868          vertex_coordinates: List of coordinate pairs [xi, eta] of
    69               points constituting a mesh (or an m x 2 numeric array or
     69          points constituting a mesh (or an m x 2 numeric array or
    7070              a geospatial object)
    7171              Points may appear multiple times
  • anuga_core/source/anuga/utilities/mainland_only.csv

    r5086 r7690  
    22447203.12,7699910.99
    33447515.37,7699565.56
    4 447802.51,7699594.13
    5 447897.04,7699314.42
    6 448096.87,7699309.5
    7 448049.53,7699137.81
    8 448241.01,7699138.4
    9 447998.25,7699230.62
    10 448202.22,7699229.04
    11 448076.58,7699475.45
    12 447944.63,7699405.32
    13 447823.17,7699641.78
    14 447701.51,7699610.41
    15 447805.71,7699904.03
    16 447882,7700136.68
    17 448239.36,7700008.3
    18 448380.84,7700024.23
    19 448136.17,7700058.89
    20 447721.58,7700183.77
    21 447669.07,7699667.86
    224447500.56,7699975.01
    235447108.86,7700098.85
     
    3618449882.5,7702139.37
    3719449798.8,7701938.79
    38 449950.91,7701889.44
    39 449858.65,7701764.11
    40 450041.48,7701881.96
    41 449808.01,7701988.62
    42 449874.32,7702091.75
    43 450046.73,7702217.33
    44 450427.05,7702076.79
    45 450713.91,7701864.04
    46 450914.68,7701542.56
    47 450865.57,7701253.56
    48 450830.65,7701097.4
    49 450841.67,7700886.05
    50 450609.91,7700776.91
    51 450814.8,7700822.89
    52 450892.05,7701098.69
    53 450959.46,7701180.79
    54 450946.17,7701453.01
    55 451146.1,7701420.39
    56 451556.06,7701458.1
    57 451694.32,7701513.84
    58 451782.48,7701260.65
    59 452236.82,7701064.95
    60 452539.29,7700837.82
    61 452541.07,7700944.07
    62 452332.44,7701111.71
    63 451953.31,7701211.34
    64 451784.01,7701451.02
    65 451712.63,7701663.3
    66 451519.51,7701500.05
    67 451204.36,7701429.42
    68 451096.7,7701584.05
    69 450949.31,7701800.54
    70 450701.11,7701969.14
    71 450562.26,7702109.29
    72 450268.03,7702351.9
    73 449866.18,7702377.26
    74 449649.11,7702566.98
    75 449368.93,7702974.52
    76 449318.54,7703114.93
    77 449401.98,7703402.93
    78 449536,7703490.76
    79 449636.2,7703747.83
    80 449878.9,7704041.85
    81 450026.94,7704315.65
    82 450229.4,7704493.33
    83 450389.76,7704479.42
    84 450662.87,7704351.85
    85 450856.13,7704121.11
    86 450964.33,7703783.86
    87 451102.63,7703477.7
    88 451378.82,7703357.87
    89 451788.47,7703518.42
    90 451947.33,7703661.65
    91 452211.68,7703676.79
    92 452215.58,7703404.54
    93 451836.88,7702974.04
    94 452106.41,7702632.83
    95 452265.58,7702295.73
    9620452234.95,7702449.48
    9721452124.96,7702695.97
     
    14569455206.33,7705662.8
    14670455432.49,7705568.22
    147 455328.81,7705799.25
    148 455542.92,7705931.52
    149 455994.78,7705918.32
    150 456256.1,7705910.15
    151 456387.59,7705789.86
    152 456538.72,7705720.53
    153 456515.56,7705419.44
    154 456530.5,7705278.92
    155 456578.82,7705513.68
    156 456527.4,7705670.7
    157 456655,7705842.57
    158 456363.24,7705945.85
    159 456639.63,7706150.2
    160 456878.51,7706373.27
    161 457188.46,7706488.06
    162 457472.53,7706139.05
    163 457782.93,7706071.22
    164 457529.4,7705878.01
    165 457412.68,7705515.81
    166 457074.97,7705273.68
    167 457279.32,7705151.36
    168 457463.68,7705110.88
    169 457604.41,7705033.76
    170 457437.51,7705169.47
    171 457233.96,7705384.76
    172 457552.19,7705512.84
    173 457676.08,7705923.75
    174 458180.98,7705929.45
    175 458540.56,7705758.8
    176 458546.42,7705494.31
    177 458514.81,7705228.61
    178 458816.11,7705046.74
    179 458783.48,7705192.75
    180 458484.17,7705410.04
    181 458712.21,7705389.57
    182 458976.32,7705515.28
    183 459128.1,7705600.87
    184 458840.13,7705439.69
    185 458713.57,7705683.97
    186 458529.7,7705942.49
    187 458272.31,7706041.45
    188 458040.76,7706214.64
    189 457889.15,7706477.66
    190 457710.6,7706681.95
    191 457623.71,7706874.3
    192 457755.5,7707041.75
    193 458081.25,7707092.37
    194 458067.52,7707170.92
    195 457663.81,7707071.4
    196 457565.14,7706977.08
    197 457478.27,7707157.26
    198 457775.78,7707262.04
    199 458245.06,7707370.57
    200 458501.69,7707588.12
    201 458731.37,7707758.02
    202 458848.75,7707870.09
    203 459056.24,7708172.73
    204 459281.45,7708474.31
    205 459417.46,7708638.43
    206 459761.72,7708796.41
    207 459800.63,7708634.93
    208 459653.3,7708419.87
    209 459920.41,7708619.71
    210 460181.3,7708814.01
    211 460300.74,7708941.57
    212 460473.14,7709130.12
    213 460650.4,7709095.42
    214 460376.97,7708408.31
    215 460330,7708016.42
    216 460051.17,7707926.11
    217 459976.52,7707795.34
    218 460031.85,7707733.5
    219 460445.72,7707952.5
    220 460662.87,7707704
    221 460841.78,7707329.24
    222 460579.31,7707372.89
    223 460519.52,7707562
    224 460347.01,7707424.36
    225 460582.11,7707510.13
    226 460583.74,7707261.12
    227 460815.17,7707132.18
    228 460931.27,7707350.48
    229 460817.97,7707715.43
    230 461029.59,7708051.26
    231 461072.53,7708395.55
    232 460945.07,7708581.18
    233 460970.87,7709126.85
    234 461133.74,7709389.53
    235 461269.97,7709467.31
    236 461394.48,7709657.95
    23771461705.68,7709721.75
    23872462681.38,7710184.34
     
    496330478255.43,7714831.11
    497331478127.84,7714613.18
    498 478268.97,7714032.09
    499 478416.29,7713672.6
    500 478259.24,7713504.18
    501 478105.32,7713332.45
    502 477593.89,7713369.4
    503 477803.58,7713107.39
    504 477759.28,7712741.01
    505 478042.92,7712468.03
    506 478124.29,7712364.11
    507 477873.96,7712651.52
    508 477877.52,7713114.13
    509 478139.79,7713255.02
    510 478421.16,7713121.47
    511 478507.01,7712777.4
    512 478477.48,7713061.78
    513 478315.61,7713405.76
    514 478521.27,7713840.96
    515 478291.77,7714115.13
    516 478290.13,7714582.15
    517 478434.75,7714706.29
    518 478887.94,7714624.97
    519 479148.9,7714167.12
    520 479227.67,7713637.11
    521 479325.14,7713144.75
    522 479127.16,7713221.97
    523 479404.41,7713050.78
    524 479509.43,7713188.14
    525 479662.75,7713004.61
    526 479850.6,7712682.79
    527 480148.8,7712394.29
    528 480355.23,7712206.39
    529 480446.42,7711972.48
    530 480574.6,7711576.89
    531 480263.09,7711542.43
    532 479801.34,7711470.75
    533 480151.44,7711255.73
    534 480675.22,7711463.86
    535 480935.74,7711377.02
    536 481274.96,7711141.7
    537 481648.11,7710828.91
    538 481866.13,7710472.79
    539332482021.51,7710263.79
    540333482078.8,7710247.25
     
    625418496055.54,7712910.29
    626419496465.85,7713002.24
    627 496605.35,7713285.59
    628 496814.63,7713583.33
    629 497167.66,7713809.16
    630 497289.51,7713874.48
    631 497134.28,7714116.81
    632 497252.99,7714279.52
    633 497506.09,7714258.53
    634 497594.65,7714015.08
    635 497996.72,7713721.86
    636 497856.15,7713478.37
    637 497844.78,7712887.39
    638 497806.32,7712316.33
    639 497615.77,7712118.2
    640 497856.29,7712414.84
    641 497913.53,7712783.37
    642 497913.46,7713306.83
    643 498066.51,7713686.45
    644 497809.22,7713948.71
    645 497824.79,7714289.57
    646 498351.81,7714378.17
    647 498580.96,7714271.95
    648 498520.59,7713890.13
    649 498933.04,7713824.87
    650 498678.93,7713583.6
    651 498747.66,7713600.2
    652 499010.11,7713746.3
    653 498782,7713983.12
    654 498544.53,7714024.05
    655 498621.57,7714410.29
    656 498869.47,7714234.34
    657 499184.03,7713996.43
    658 499113.22,7713789.47
    659 499548.58,7713736.37
    660 499555.86,7713822.69
    661 499192.37,7713813.82
    662 499265.27,7714007.5
    663 498936.13,7714264.23
    664 498637.18,7714516.53
    665 498190.35,7714502.1
    666 497623.76,7714388.04
    667 497548.74,7714554.03
    668 497514.33,7714811.88
    669 497892.43,7714815.25
    670 498068.43,7714974.64
    671 498256.93,7715211.5
    672 498531.91,7715324.41
    673 498686.05,7715515.88
    674 498900.61,7715666.4
    675 499294.35,7715710.7
    676 499576.63,7715688.57
    677 499948.49,7715861.23
    678 500144.32,7716010.63
    679 500342.23,7715994.03
    680 500548.48,7716033.87
    681 500731.82,7716156.71
    682 500929.73,7716177.73
    683 501103.7,7716323.8
    684 501317.24,7716344.82
    685 501469.34,7716462.11
    686 501611.01,7716540.67
    687 501820.39,7716489.75
    688 502037.06,7716486.4
    689 502243.31,7716486.38
    690 502468.32,7716540.58
    691 502691.25,7716596.99
    692 502913.14,7716624.62
    693420503095.45,7716705.38
    694421503302.75,7716750.71
     
    845572517248.64,7721707.83
    846573517357.62,7721981.01
    847 454913.9,7687033.4
     574530000,7687033.4
     575424913.9,7687033.4
    848576425539.85,7691615.14
    849577426064.25,7691355.15
     
    854582426898.56,7691103.16
    855583426763.28,7690643.24
    856 427055.91,7690088.91
    857 426796.88,7690332.37
    858 426857.75,7690921.46
    859 427646.12,7691424.07
    860 428763.81,7691794.12
    861 429123.61,7691823.32
    862 430084.18,7692218.07
    863 430563.65,7692990.4
    864 430763.67,7693423.99
    865584430997.89,7693894.23
    866585431118.95,7694305.35
  • anuga_core/source/anuga/utilities/numerical_tools.py

    r7317 r7690  
    7676    if v2 is None:
    7777        v2 = [1.0, 0.0] # Unit vector along the x-axis
    78        
     78   
    7979    v1 = ensure_numeric(v1, num.float)
    8080    v2 = ensure_numeric(v2, num.float)   
  • anuga_core/source/anuga/utilities/polygon.py

    r7687 r7690  
    322322    return False
    323323
     324   
    324325def is_complex(polygon, verbose=False):
    325     """Check if a polygon is complex (self-intersecting)
    326     """
    327    
     326    """Check if a polygon is complex (self-intersecting).
     327       Uses a sweep algorithm that is O(n^2) in the worst case, but
     328       for most normal looking polygons it'll be O(n log n).
     329
     330       polygon is a list of points that define a closed polygon.
     331       verbose will print a list of the intersection points if true
     332       
     333       Return True if polygon is complex.
     334    """           
     335           
     336    def key_xpos(item):
     337        return (item[0][0])
     338   
     339    def segments_joined(seg0, seg1):
     340        for i in seg0:
     341            for j in seg1:   
     342                if i == j: return True
     343        return False
     344       
    328345    polygon = ensure_numeric(polygon, num.float)
    329346
     347    # build a list of discrete segments from the polygon
     348    unsorted_segs = []
    330349    for i in range(0, len(polygon)-1):
    331         for j in range(i+1, len(polygon)-1):   
    332                 (type, point) = intersection([polygon[i], polygon[i+1]], [polygon[j], polygon[j+1]])
    333 
    334                 if (abs(i-j) > 1 and type == 1) or (type == 2 and list(point[0]) != list(point[1])) or (type == 3) and type != 4:
     350        unsorted_segs.append([list(polygon[i]), list(polygon[i+1])])
     351    unsorted_segs.append([list(polygon[0]), list(polygon[-1])])
     352   
     353    # all segments must point in same direction
     354    for val in unsorted_segs:
     355        if val[0][0] > val[1][0]:
     356            val[0], val[1] = val[1], val[0]   
     357           
     358    l_x = sorted(unsorted_segs, key=key_xpos)
     359
     360    comparisons = 0
     361   
     362    # loop through, only comparing lines that partially overlap in x
     363    for index, leftmost in enumerate(l_x):
     364        cmp = index+1
     365        while cmp < len(l_x) and leftmost[1][0] > l_x[cmp][0][0]:
     366            if not segments_joined(leftmost, l_x[cmp]):
     367                (type, point) = intersection(leftmost, l_x[cmp])
     368                comparisons += 1
     369                if type != 0 and type != 4 or (type == 2 and list(point[0]) != list(point[1])):
    335370                    if verbose:
    336                         print 'Self-intersecting polygon found, type ', type, ' point', point, 'vertex indices ', i, j               
    337                     return True
     371                        print 'Self-intersecting polygon found, type ', type, ' point', point,
     372                        print 'vertices: ', leftmost, ' - ', l_x[cmp]               
     373                    return True           
     374            cmp += 1
    338375       
    339376    return False
     377   
    340378   
    341379def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
     
    9961034        polygon.append([float(fields[0]), float(fields[1])])
    9971035   
    998     # check this is a valid polygon   
    999     # if is_complex(polygon):   
    1000             # msg = 'ERROR: Self-intersecting polygon detected in file' + filename +'. '
    1001             # msg += 'Please fix.'
    1002             # log.critical(msg)
    1003 # #            raise Exception, msg
     1036    # check this is a valid polygon.
     1037    if is_complex(polygon, verbose=True):   
     1038        msg = 'ERROR: Self-intersecting polygon detected in file ' + filename +'. '
     1039        msg += 'A complex polygon will not necessarily break the algorithms within ANUGA, '
     1040        msg += 'but it usually signifies pathological data. Please fix this file.'
     1041        raise Exception, msg
    10041042   
    10051043    return polygon
  • anuga_core/source/anuga/utilities/test_polygon.py

    r7687 r7690  
    18081808       
    18091809    def test_is_polygon_complex(self):
     1810        """Test a concave and a complex poly with is_complex, to make sure it can detect
     1811           self-intersection.
     1812        """
    18101813        concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
    18111814        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], [0, 10]]
    18121815
    1813         not_complex = is_complex(concave_poly)       
    1814         complex = is_complex(complex_poly)       
    1815 
    1816         assert not not_complex
    1817         assert complex
     1816        assert not is_complex(concave_poly)
     1817        assert is_complex(complex_poly)
    18181818
    18191819    def test_is_polygon_complex2(self):
     1820        """Test a concave and a complex poly with is_complex, to make sure it can detect
     1821           self-intersection. This test uses more complicated polygons.
     1822        """   
    18201823        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7,6], [10, 10], [1,5], [0, 10]]
    1821         complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]
    1822 
    1823         not_complex = is_complex(concave_poly)       
    1824         complex = is_complex(complex_poly)       
    1825 
    1826         assert not not_complex
    1827         assert complex
     1824        complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]       
     1825
     1826        assert not is_complex(concave_poly)
     1827        assert is_complex(complex_poly)
    18281828       
    18291829################################################################################
Note: See TracChangeset for help on using the changeset viewer.