Changeset 3944


Ignore:
Timestamp:
Nov 8, 2006, 5:16:44 PM (16 years ago)
Author:
sexton
Message:

(i) update sww2timeseries so can handle gauges which don't fall within sww domain (ii) script for generating report for dampier based on sww files from parallel setup

Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/documentation/user_manual/anuga_user_manual.tex

    r3925 r3944  
    32193219http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
    32203220
    3221 
     3221\bibitem[grid250]{grid250}
     3222Australian Bathymetry and Topography Grid, June 2005.
     3223Webster, M.A. and Petkovic, P.
     3224Geoscience Australia Record 2005/12. ISBN: 1-920871-46-2.\\
     3225http://www.ga.gov.au/meta/ANZCW0703008022.html
    32223226
    32233227
  • anuga_core/documentation/user_manual/examples/project.py

    r3869 r3944  
    4545# demo poly
    4646j0 = [385000, 6280000]
    47 j1 = [360000, 6272500]
    48 j2 = [335000, 6272500]
     47j1 = [360000, 6273000]
     48j2 = [335000, 6273000]
    4949j3 = [330000, 6265000]
    5050j31 = [325000, 6260000]
    5151j4 = [316000, 6260000]
    52 j5 = [316000, 6247000]
    53 j6 = [350000, 6247000]
     52j5 = [316000, 6246750]
     53j6 = [350000, 6246750]
    5454j7 = [385000, 6238000]
    5555
    5656demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
    5757
    58 from utilities.polygon import read_polygon
    59 polygonptsfile4 = polygondir + 'poly1'
    60 polygonptsfile0 = polygondir + 'poly2'
    61 polygonptsfile1 = polygondir + 'poly3'
    62 polygonptsfile2 = polygondir + 'poly4'
    63 polygonptsfile3 = polygondir + 'poly5'
    64 northern_polygon = read_polygon(polygonptsfile0 + '.csv')
    65 manly_polygon = read_polygon(polygonptsfile1 + '.csv')
    66 harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
    67 southern_polygon = read_polygon(polygonptsfile3 + '.csv')
    68 top_polygon = read_polygon(polygonptsfile4 + '.csv')
     58from anuga.utilities.polygon import read_polygon, plot_polygons
     59#polygonptsfile4 = polygondir + 'poly1'
     60#polygonptsfile0 = polygondir + 'poly2'
     61#polygonptsfile1 = polygondir + 'poly3'
     62#polygonptsfile2 = polygondir + 'poly4'
     63#polygonptsfile3 = polygondir + 'poly5'
     64#northern_polygon = read_polygon(polygonptsfile0 + '.csv')
     65#manly_polygon = read_polygon(polygonptsfile1 + '.csv')
     66#harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
     67#southern_polygon = read_polygon(polygonptsfile3 + '.csv')
     68#top_polygon = read_polygon(polygonptsfile4 + '.csv')
     69coastal_polygon = read_polygon(polygondir+'coastal'+'.csv')
     70shallow_polygon = read_polygon(polygondir+'shallow'+'.csv')
     71
     72#plot_polygons([demopoly,northern_polygon,manly_polygon,harbour_polygon,southern_polygon,top_polygon],'model_setup',verbose=False)
     73plot_polygons([demopoly,coastal_polygon,shallow_polygon],'new_model_setup',verbose=False)
    6974
    7075slump_origin = [372500.0, 6255000.0] #Absolute UTM
  • anuga_core/documentation/user_manual/examples/runsydney.py

    r3869 r3944  
    5959# Interior regions and mesh resolutions
    6060interior_res = 5000
    61 interior_regions = [[project.northern_polygon, interior_res],
    62                     [project.harbour_polygon, interior_res],
    63                     [project.southern_polygon, interior_res],
    64                     [project.manly_polygon, interior_res],
    65                     [project.top_polygon, interior_res]]
    66 
     61shallow_res = 15000
     62##interior_regions = [[project.northern_polygon, interior_res],
     63##                    [project.harbour_polygon, interior_res],
     64##                    [project.southern_polygon, interior_res],
     65##                    [project.manly_polygon, interior_res],
     66##                    [project.top_polygon, interior_res]]
     67interior_regions = [[project.coastal_polygon, interior_res],
     68                    [project.shallow_polygon, shallow_res]]
    6769
    6870create_mesh_from_regions(project.demopoly,
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r3939 r3944  
    681681       
    682682    assert type(gauge_filename) == type(''),\
    683                'Gauge filename must be a string'
     683           'Gauge filename must be a string'
    684684   
    685685    try:
     
    693693        report = False
    694694       
    695    
    696695    if plot_quantity is None:
    697696        plot_quantity = ['depth', 'speed', 'bearing']
     
    733732                          verbose = True,
    734733                          use_cache = True)
    735  
     734
     735        # determine which gauges are contained in sww file
     736        count = 0
     737        gauge_index = []
     738        print 'swwfile', swwfile
     739        for k, g in enumerate(gauges):
     740            if f(0.0, point_id = k)[2] > 1.0e6:
     741                count += 1
     742                if count == 1: print 'Gauges not contained here:'
     743                print locations[k]
     744            else:
     745                gauge_index.append(k)
     746
     747        if len(gauge_index) > 0:
     748            print 'Gauges contained here: \n',
     749        else:
     750            print 'No gauges contained here. \n'
     751        for i in range(len(gauge_index)):
     752             print locations[gauge_index[i]]
     753             
    736754        index = swwfile.rfind(sep)
    737755        file_loc.append(swwfile[:index+1])
     
    759777            raise Exception, msg
    760778
    761     if verbose: print 'Inputs OK - going to generate figures'
    762 
    763     return generate_figures(plot_quantity, file_loc, report, reportname, surface,
    764                             leg_label, f_list, gauges, locations, elev, production_dirs,
    765                             time_min, time_max, title_on, label_id, verbose)
     779    if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures'
     780
     781    if len(gauge_index) <> 0:
     782        texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface,
     783                                                leg_label, f_list, gauges, locations, elev, gauge_index,
     784                                                production_dirs, time_min, time_max, title_on, label_id, verbose)
     785    else:
     786        texfile = ''
     787        elev_output = []
     788
     789    return texfile, elev_output
    766790                         
    767791#Fixme - Use geospatial to read this file - it's an xya file
    768792#Need to include other information into this filename, so xya + Name - required for report
    769793def get_gauges_from_file(filename):
     794    """ Read in gauge information from file
     795    """
    770796    from os import sep, getcwd, access, F_OK, mkdir
    771797    fid = open(filename)
     
    803829
    804830def check_list(quantity):
    805 
     831    """ Check that input quantities in quantity list are possible
     832    """
    806833    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
    807834                    'ymomentum', 'speed', 'bearing', 'elevation']
     
    817844
    818845def calc_bearing(uh, vh):
    819 
     846    """ Calculate velocity bearing from North
     847    """
    820848    from math import atan, degrees
    821849   
     
    836864
    837865def generate_figures(plot_quantity, file_loc, report, reportname, surface,
    838                      leg_label, f_list,
    839                      gauges, locations, elev, production_dirs, time_min, time_max,
    840                      title_on, label_id, verbose):
    841 
     866                     leg_label, f_list, gauges, locations, elev, gauge_index,
     867                     production_dirs, time_min, time_max, title_on, label_id,
     868                     verbose):
     869    """ Generate figures based on required quantities and gauges for each sww file
     870    """
    842871    from math import sqrt, atan, degrees
    843872    from Numeric import ones, allclose, zeros, Float, ravel
     
    849878    if surface is True:
    850879        import mpl3d.mplot3d as p3
    851 
     880       
    852881    if report == True:   
    853882        texdir = getcwd()+sep+'report'+sep
     
    904933        fid_compare = open(comparefile, 'w')
    905934        ##### loop over each gauge #####
    906         for k, g in enumerate(gauges):
     935        for k in gauge_index:
     936            g = gauges[k]
    907937            if verbose: print 'Gauge %d of %d' %(k, len(gauges))
    908938            min_stage = 10
     
    9921022        savefig('profilefig')
    9931023               
    994     #stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
    995     stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])   
     1024    stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
     1025    #stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])   
    9961026    vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
    9971027    mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1]) 
     
    10021032    elev_output = []
    10031033    if len(label_id) > 1: graphname_report = []
    1004     for k, g in enumerate(gauges):
     1034    for k in gauge_index:
     1035        g = gauges[k]
    10051036        count1 = 0
    10061037        if report == True and len(label_id) > 1:
     
    10081039            fid.write(s)
    10091040        if len(label_id) > 1: graphname_report = []
    1010         #### generate figures for each gauge ###
     1041        #### generate figures for each gauge ####
    10111042        for j, f in enumerate(f_list):
    10121043            ion()
     
    11091140                        word_quantity += plot_quantity[i]
    11101141                   
    1111                 word_quantity += ' and ' + plot_quantity[nn-1]               
     1142                word_quantity += ' and ' + plot_quantity[nn-1]
    11121143                caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' '))
    11131144                if elev[k] == 0.0:
     
    11641195    return texfile2, elev_output
    11651196
     1197
     1198def gauge_in_sww(swwfile,
     1199                 gauge_filename,
     1200                 verbose = False):
     1201   
     1202    """ Determine which gauges fall in sww file.
     1203
     1204    Input variables:
     1205
     1206    swwfiles        - sww files
     1207   
     1208    gauge_filename  - name of file containing gauge data
     1209                        - easting, northing, name , elevation
     1210                    - OR (this is not yet done)
     1211                        - structure which can be converted to a Numeric array,
     1212                          such as a geospatial data object
     1213                     
     1214    Output:
     1215   
     1216    - Modified gauge_filename containing gauges which fall is given sww file.
     1217      Name = gauges_filename_mod
     1218    Other important information:
     1219                         
     1220    """
     1221
     1222   
     1223    k = _gauge_in_sww(swwfile,
     1224                      gauge_filename,
     1225                      verbose)
     1226
     1227    return k
     1228
     1229def _gauge_in_sww(swwfile,
     1230                  gauge_filename,
     1231                  verbose = False):   
     1232
     1233    try:
     1234        fid = open(swwfile)
     1235    except Exception, e:
     1236        msg = 'File "%s" could not be opened: Error="%s"'\
     1237              %(swwfile, e)
     1238        raise msg
     1239       
     1240    assert type(gauge_filename) == type(''),\
     1241           'Gauge filename must be a string'
     1242   
     1243    try:
     1244        fid = open(gauge_filename)
     1245    except Exception, e:
     1246        msg = 'File "%s" could not be opened: Error="%s"'\
     1247                  %(gauge_filename, e)
     1248        raise msg
     1249
     1250    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
     1251    f = file_function(swwfile,
     1252                      quantities = sww_quantity,
     1253                      #interpolation_points = gauges,
     1254                      verbose = True,
     1255                      use_cache = True)
     1256    print dir(f)
     1257   
     1258    #from bstract_2d_finite_volumes.neighbour_mesh import get_boundary_polygon
     1259    from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     1260    #from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
     1261    #from anuga.fit_interpolate.interpolate import Interpolation_function
     1262    from anuga.fit_interpolate.interpolate import Interpolate
     1263
     1264    test = Mesh(f.vertex_coordinates)
     1265    print test.mesh.get_boundary_polygon()
     1266   
     1267    return check_gauge(gauge_filename,bounding_polygon)
     1268
     1269def check_gauge(gauge_filename,bounding_polygon):
     1270
     1271    from anuga.utilities.polygon import is_inside_polygon
     1272   
     1273    filename_out = gauge_filename + '_mod'
     1274    fid_out = open(filename_out, 'w')
     1275    fid_out.write(s)
     1276
     1277    gauges, locations, elev = get_gauges_from_file(gauge_filename)
     1278   
     1279    for i, gauge in enumerate(gauges):
     1280        v = is_inside_polygon(gauge,bounding_polygon,verbose=False)
     1281        if v == True:
     1282            s = '%.2f, %.2f, %.2f, %s\n' %(gauges[0], gauges[1], elev[i], locations[i])
     1283   
     1284    return fid_out
     1285
    11661286from os.path import basename
    11671287
  • anuga_work/development/cairns_2006/run_cairns.py

    r3846 r3944  
    99#-------------------------------
    1010import sys, os
    11 from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
     11from anuga.shallow_water import Domain, Reflective_boundary,\
    1212     File_boundary, Transmissive_Momentum_Set_Stage_boundary,\
    1313     Transmissive_boundary, Dirichlet_boundary
    14 from anuga.pyvolution.mesh_factory import rectangular_cross
    15 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
     14#from anuga.pyvolution.mesh_factory import rectangular_cross
     15#from anuga.pmesh2domain import pmesh_to_domain_instance
     16from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1617from Numeric import array, zeros, Float, allclose
    1718import cairns_project as project
     
    4142"""
    4243
    43 domain = cache(pmesh_to_domain_instance,
    44                (project.mesh_filename, Domain),
    45                dependencies = [project.mesh_filename])
     44#cache(pmesh_to_domain_instance,
     45#      (project.mesh_filename, Domain),
     46#      dependencies = [project.mesh_filename])
     47domain = Domain(project.mesh_filename, use_cache = True, verbose = True)
    4648
    4749domain.check_integrity()
  • anuga_work/production/broome_2006/export_results.py

    r3881 r3944  
    66from os import sep
    77
    8 time_dir = '20061026_220716'
     8time_dir = '20061101_024119'
    99
    1010directory = project.outputdir
     
    2222
    2323if which_var == 2:  # Depth
    24     outname = name + '_depth_bruny'
     24    outname = name + '_depth'
    2525    quantityname = 'stage-elevation'  #Depth
    2626
     
    3737sww2dem(name, basename_out = outname,
    3838            quantity = quantityname,
    39             cellsize = 25,      # would prefer this at 25
     39            cellsize = 100,      # would prefer this at 25
    4040            # define region for viz purposes
    4141            easting_min = project.eastingmin,
  • anuga_work/production/broome_2006/project.py

    r3908 r3944  
    138138
    139139# exporting asc grid
    140 eastingmin = 406215.87
    141 eastingmax = 440208.78
    142 northingmin = 7983427.73
    143 northingmax = 8032834.52
     140eastingmin = 412036.43
     141eastingmax = 427184.37
     142northingmin = 8007984
     143northingmax = 8029830
    144144
    145145       
  • anuga_work/production/hobart_2006/export_results.py

    r3852 r3944  
    4949##            format = 'asc')
    5050
    51 # for bruny
     51### for bruny
     52##sww2dem(name, basename_out = outname,
     53##            quantity = quantityname,
     54##            cellsize = 25,      # would prefer this at 25
     55##            # define region for viz purposes
     56##            easting_min = project.eastingmin25_2,
     57##            easting_max = project.eastingmax25_2,
     58##            northing_min = project.northingmin25_2,
     59##            northing_max = project.northingmax25_2,       
     60##            reduction = max, #this is because we want max quantityname
     61##            verbose = True,
     62##            format = 'asc')
     63
     64# for fixed timestep
    5265sww2dem(name, basename_out = outname,
    5366            quantity = quantityname,
    54             cellsize = 25,      # would prefer this at 25
     67            cellsize = 10,      # would prefer this at 25
    5568            # define region for viz purposes
    5669            easting_min = project.eastingmin25_2,
     
    5871            northing_min = project.northingmin25_2,
    5972            northing_max = project.northingmax25_2,       
    60             reduction = max, #this is because we want max quantityname
     73            #reduction = max, #this is because we want max quantityname
    6174            verbose = True,
    6275            format = 'asc')
     76
  • anuga_work/production/hobart_2006/make_report.py

    r3886 r3944  
    9999                                      #reportname = 'latexoutput',
    100100                                      plot_quantity = ['stage', 'speed'],
    101                                       surface = False,
     101                                      surface = True,
    102102                                      time_min = None,
    103103                                      time_max = None,
  • anuga_work/production/hobart_2006/report/hobart_2006_report_MRT.tex

    r3867 r3944  
    172172     \input{metadata}
    173173
    174 \pagebreak
     174\clearpage
    175175
    176176   \section{Time series}
Note: See TracChangeset for help on using the changeset viewer.