Changeset 2622


Ignore:
Timestamp:
Mar 28, 2006, 2:36:03 PM (18 years ago)
Author:
steve
Message:

Updated merimbula viewer

Files:
51 added
8 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/view_tsh.py

    r1363 r2622  
    99from shallow_water import Domain
    1010from pmesh2domain import pmesh_to_domain_instance
    11 from util import file_function, Polygon_function, read_polygon
     11from pyvolution.util import file_function
     12from utilities.polygon import Polygon_function, read_polygon
    1213from Numeric import zeros, Float, maximum, minimum
    1314from realtime_visualisation_new import *
  • inundation/utilities/polygon.py

    r2531 r2622  
    1010    print 'Could not find scipy - using Numeric'
    1111    from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot
    12    
     12
    1313
    1414from math import sqrt
     
    4545def inside_polygon(points, polygon, closed = True, verbose = False):
    4646    """Determine points inside a polygon
    47    
     47
    4848       Functions inside_polygon and outside_polygon have been defined in
    49        terms af separate_by_polygon which will put all inside indices in 
     49       terms af separate_by_polygon which will put all inside indices in
    5050       the first part of the indices array and outside indices in the last
    51        
    52        See separate_points_by_polygon for documentation         
     51
     52       See separate_points_by_polygon for documentation
    5353    """
    5454
    5555    if verbose: print 'Checking input to inside_polygon'
    56    
     56
    5757    try:
    5858        points = ensure_numeric(points, Float)
     
    6666        polygon = ensure_numeric(polygon, Float)
    6767    except NameError, e:
    68         raise NameError, e       
     68        raise NameError, e
    6969    except:
    7070        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
     
    9090def outside_polygon(points, polygon, closed = True, verbose = False):
    9191    """Determine points outside a polygon
    92    
     92
    9393       Functions inside_polygon and outside_polygon have been defined in
    94        terms af separate_by_polygon which will put all inside indices in 
     94       terms af separate_by_polygon which will put all inside indices in
    9595       the first part of the indices array and outside indices in the last
    96        
    97        See separate_points_by_polygon for documentation         
     96
     97       See separate_points_by_polygon for documentation
    9898    """
    9999
     
    102102        points = ensure_numeric(points, Float)
    103103    except NameError, e:
    104         raise NameError, e       
     104        raise NameError, e
    105105    except:
    106106        msg = 'Points could not be converted to Numeric array'
     
    110110        polygon = ensure_numeric(polygon, Float)
    111111    except NameError, e:
    112         raise NameError, e       
     112        raise NameError, e
    113113    except:
    114114        msg = 'Polygon could not be converted to Numeric array'
     
    133133            # No points are outside
    134134            return []
    135         else:   
     135        else:
    136136            return indices[count:][::-1]  #return reversed
    137137
     
    180180    #Input checks
    181181
    182    
     182
    183183    try:
    184184        points = ensure_numeric(points, Float)
    185185    except NameError, e:
    186         raise NameError, e       
     186        raise NameError, e
    187187    except:
    188188        msg = 'Points could not be converted to Numeric array'
     
    192192        polygon = ensure_numeric(polygon, Float)
    193193    except NameError, e:
    194         raise NameError, e       
     194        raise NameError, e
    195195    except:
    196196        msg = 'Polygon could not be converted to Numeric array'
     
    276276    if verbose: print 'Checking input to separate_points_by_polygon'
    277277
    278    
     278
    279279    #Input checks
    280280
    281281    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
    282     assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean' 
    283 
    284    
     282    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
     283
     284
    285285    try:
    286286        points = ensure_numeric(points, Float)
    287287    except NameError, e:
    288         raise NameError, e       
     288        raise NameError, e
    289289    except:
    290290        msg = 'Points could not be converted to Numeric array'
     
    295295        polygon = ensure_numeric(polygon, Float)
    296296    except NameError, e:
    297         raise NameError, e       
     297        raise NameError, e
    298298    except:
    299299        msg = 'Polygon could not be converted to Numeric array'
     
    391391
    392392        #Make points in polygons relative to geo_reference
    393         self.regions = []           
     393        self.regions = []
    394394        for polygon, value in regions:
    395             P = geo_reference.change_points_geo_ref(polygon)           
     395            P = geo_reference.change_points_geo_ref(polygon)
    396396            self.regions.append( (P, value) )
    397397
     
    429429        return z
    430430
    431 def read_polygon(filename):
     431def read_polygon(filename,split=','):
    432432    """Read points assumed to form a polygon.
    433433       There must be exactly two numbers in each line separated by a comma.
     
    441441    polygon = []
    442442    for line in lines:
    443         fields = line.split(',')
     443        fields = line.split(split)
    444444        polygon.append( [float(fields[0]), float(fields[1])] )
    445445
     
    453453       number_of_points - (optional) number of points
    454454       seed - seed for random number generator (default=None)
    455        exclude - list of polygons (inside main polygon) from where points should be excluded 
     455       exclude - list of polygons (inside main polygon) from where points should be excluded
    456456
    457457    Output:
     
    489489
    490490            append = True
    491            
     491
    492492            #Check exclusions
    493493            if exclude is not None:
     
    495495                    if inside_polygon( [x,y], ex_poly ):
    496496                        append = False
    497                        
    498                    
    499         if append is True:               
    500             points.append([x,y])               
     497
     498
     499        if append is True:
     500            points.append([x,y])
    501501
    502502    return points
     
    530530                        else:
    531531                            x_delta = x+x_mult*x*delta
    532                            
     532
    533533                        if y == 0:
    534534                            y_delta = y_mult*delta
    535535                        else:
    536536                            y_delta = y+y_mult*y*delta
    537                            
     537
    538538                        point = [x_delta, y_delta]
    539539                        #print "point",point
    540540                        if inside_polygon( point, polygon, closed = False ):
    541541                            raise Found
    542         except Found: 
     542        except Found:
    543543            point_in = True
    544544        else:
    545             delta = delta*0.1 
     545            delta = delta*0.1
    546546    return point
    547547
  • inundation/zeus/Merimbula.zpi

    r2255 r2622  
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
    75     <file>..\merimbula\prepare.py</file>
    76     <file>..\merimbula\project.py</file>
    77     <file>..\merimbula\run_cairns.py</file>
    78     <file>..\merimbula\run_least_squares_merimbula.py</file>
    79     <file>..\merimbula\run_merimbula_lake.py</file>
    80     <file>..\merimbula\run_merimbula_lake_checkpoint.py</file>
    81     <file>..\merimbula\run_new_meribula.py</file>
    82     <file>..\merimbula\view_sww.py</file>
     75    <file>..\..\production\merimbula_2005\prepare.py</file>
     76    <file>..\..\production\merimbula_2005\project.py</file>
     77    <file>..\..\production\merimbula_2005\run_merimbula_lake.py</file>
     78    <file>..\..\production\merimbula_2005\run_merimbula_lake_weed.py</file>
     79    <file>..\..\production\merimbula_2005\run_new_meribula.py</file>
     80    <file>..\..\production\merimbula_2005\run_new_meribula_weed.py</file>
    8381    <folder name="Header Files" />
    8482    <folder name="Resource Files" />
  • inundation/zeus/anuga-workspace.zwi

    r2386 r2622  
    22<workspace name="anuga-workspace">
    33    <mode>Debug</mode>
    4     <active>pyvolution</active>
     4    <active>Merimbula</active>
    55    <project name="analytic_solutions">analytic_solutions.zpi</project>
    66    <project name="euler">euler.zpi</project>
  • inundation/zeus/pyvolution.zpi

    r2255 r2622  
    7676    <file>..\pyvolution\bed_w_eden_boundary.py</file>
    7777    <file>..\pyvolution\bed_w_file_boundary.py</file>
    78     <file>..\pyvolution\cg_solve.py</file>
    7978    <file>..\pyvolution\check_sww_tsh.py</file>
    8079    <file>..\pyvolution\combine_pts.py</file>
     
    9089    <file>..\pyvolution\island.py</file>
    9190    <file>..\pyvolution\least_squares.py</file>
    92     <file>..\pyvolution\load_mesh\loadASCII.py</file>
    9391    <file>..\pyvolution\log.ini</file>
    9492    <file>..\pyvolution\mesh.py</file>
     
    112110    <file>..\pyvolution\shallow_water_vtk.py</file>
    113111    <file>..\pyvolution\show_balanced_limiters.py</file>
    114     <file>..\pyvolution\sparse.py</file>
    115     <file>..\pyvolution\sparse_ext.c</file>
    116112    <file>..\pyvolution\test_advection.py</file>
    117113    <file>..\pyvolution\test_all.py</file>
    118     <file>..\pyvolution\test_cg_solve.py</file>
    119114    <file>..\pyvolution\test_combine_pts.py</file>
    120115    <file>..\pyvolution\test_data_manager.py</file>
     
    131126    <file>..\pyvolution\test_region.py</file>
    132127    <file>..\pyvolution\test_shallow_water.py</file>
    133     <file>..\pyvolution\test_sparse.py</file>
    134128    <file>..\pyvolution\test_util.py</file>
    135129    <file>..\pyvolution\treenode.py</file>
     
    137131    <file>..\pyvolution\twolevels.py</file>
    138132    <file>..\pyvolution\util.py</file>
    139     <file>..\pyvolution\util_ext.c</file>
    140     <file>..\pyvolution\util_ext.h</file>
    141133    <file>..\pyvolution\view_tsh.py</file>
    142134    <file>..\pyvolution\vtk_realtime_visualiser.py</file>
  • production/merimbula_2005/prepare.py

    r2225 r2622  
    11
    22
    3 def prepare_timeboundary(filename):
    4     """Converting tide time series to NetCDF tms file.
     3def prepare_wind_stress(filename):
     4    """Converting wind timeseries to NetCDF tms file.
    55    This is a 'throw-away' code taylor made for files like
    66    'Benchmark_2_input.txt' from the LWRU2 benchmark
    77    """
    88
    9     print 'Preparing time boundary from %s' %filename
     9    print 'Preparing wind timeseries %s' %filename
    1010    from Numeric import array, zeros, Float, asarray
    1111
     
    1313
    1414    #Skip first line
    15     line = fid.readline()
     15    #line = fid.readline()
    1616
    1717    #Read remaining lines
     
    2222    N = len(lines)
    2323    T = zeros(N, Float)  #Time
    24     Q = zeros(N, Float)  #Values
     24    S = zeros(N, Float)  #Speed
     25    B = zeros(N, Float)  #Bearing
    2526
    2627    Told = 0.0
    2728    Sold = ' '
    2829    Lold = ' '
     30
    2931    for i, line in enumerate(lines):
    3032        fields = line.split()
     
    4648        T[i] = T[i] - Tstart
    4749        #this is specific to this data set. deals with daylight saving
    48         if i>3270:
    49             T[i] = T[i]+3600
    50 
     50#        if i>3270:
     51#            T[i] = T[i]+3600
     52#
    5153        if T[i]<Told :
    5254            print Lold
     
    5860            print i, T[i]-Told
    5961
    60         Q[i] = float(fields[2])
     62        S[i] = float(fields[2])
     63        B[i] = float(fields[3])
    6164
    6265        Told = T[i]
     
    6568
    6669
     70    #print T
     71    #Create tms file
     72    from Scientific.IO.NetCDF import NetCDFFile
     73
     74    outfile = filename[:-4] + '.tms'
     75    print 'Writing to', outfile
     76    fid = NetCDFFile(outfile, 'w')
     77
     78    fid.institution = 'Australian National University'
     79    fid.description = 'Input wind for Merimbula'
     80    fid.starttime = 0.0
     81    fid.createDimension('number_of_timesteps', len(T))
     82    fid.createVariable('time', Float, ('number_of_timesteps',))
     83    fid.variables['time'][:] = T
     84
     85    fid.createVariable('speed', Float, ('number_of_timesteps',))
     86    fid.variables['speed'][:] = S[:]
     87
     88    fid.createVariable('bearing', Float, ('number_of_timesteps',))
     89    fid.variables['bearing'][:] = B[:]
     90
     91
     92    fid.close()
     93
     94
     95def prepare_timeboundary(filename):
     96    """Converting tide time series to NetCDF tms file.
     97    This is a 'throw-away' code taylor made for files like
     98    'Benchmark_2_input.txt' from the LWRU2 benchmark
     99    """
     100
     101    print 'Preparing time boundary from %s' %filename
     102    from Numeric import array, zeros, Float, asarray
     103
     104    fid = open(filename)
     105
     106    #Skip first line
     107    line = fid.readline()
     108
     109    #Read remaining lines
     110    lines = fid.readlines()
     111    fid.close()
     112
     113
     114    N = len(lines)
     115    T = zeros(N, Float)  #Time
     116    Q = zeros(N, Float)  #Values
     117
     118    Told = 0.0
     119    Sold = ' '
     120    Lold = ' '
     121    for i, line in enumerate(lines):
     122        fields = line.split()
     123
     124        #print fields
     125
     126        l_time = (fields[0]+' '+fields[1])[0:-1]
     127        from time import strptime, mktime
     128
     129        s_time = strptime(l_time,'%d/%m/%y  %H:%M:%S')
     130
     131        #print s_time
     132
     133        T[i] = float(mktime(s_time))
     134
     135        if i==0:
     136            Tstart = T[0]
     137
     138        T[i] = T[i] - Tstart
     139        #this is specific to this data set. deals with daylight saving
     140        if i>3270:
     141            T[i] = T[i]+3600
     142
     143        if T[i]<Told :
     144            print Lold
     145            print l_time
     146            print Sold
     147            print s_time
     148            print Told
     149            print T[i]
     150            print i, T[i]-Told
     151
     152        Q[i] = float(fields[2])
     153
     154        Told = T[i]
     155        Sold = s_time
     156        Lold = l_time
     157
     158
     159    #print T
    67160    #Create tms file
    68161    from Scientific.IO.NetCDF import NetCDFFile
     
    90183    fid.close()
    91184
    92 
    93185#-------------------------------------------------------------
    94186if __name__ == "__main__":
     
    96188    print 'Prepare Open sea boundary condition from ',project.original_boundary_filename
    97189    prepare_timeboundary(project.original_boundary_filename )
     190
     191
     192    print 'Prepare wind from ',project.original_wind_filename
     193    prepare_wind_stress(project.original_wind_filename )
    98194
    99195    #Preparing points
  • production/merimbula_2005/project.py

    r2225 r2622  
    33
    44original_boundary_filename = 'Eden_tide_Sept03.dat'
     5original_wind_filename = 'merimbula_wind_sept_2003.dat'
    56boundary_filename = 'Eden_tide_Sept03.tms'
    67bathymetry_filename = 'merimbula_bathymetry.xya'
    78mesh_filename = 'merimbula_10785.tsh'
    8 
Note: See TracChangeset for help on using the changeset viewer.