Changeset 1280


Ignore:
Timestamp:
May 3, 2005, 5:33:04 PM (20 years ago)
Author:
steve
Message:

Consolidated files in to parallel directory

Location:
inundation/ga/storm_surge
Files:
6 added
7 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/domain.py

    r1187 r1280  
    8181        #Checkpointing and storage
    8282        from config import default_datadir
    83         self.datadir = default_datadir
    84         self.filename = 'domain'
    85         self.checkpoint = False
     83        self.datadir = default_datadir
     84        self.filename = 'domain'
     85        self.checkpoint = False
    8686
    8787
    … …  
    322322        """Output statistics about boundary forcing
    323323
    324        
    325         """
    326        
     324
     325        """
     326
    327327        if quantities is None:
    328328            quantities = self.conserved_quantities
    … …  
    344344                        v = q.boundary_values[i]
    345345                        if minval is None or v < minval: minval = v
    346                         if maxval is None or v > maxval: maxval = v                       
     346                        if maxval is None or v > maxval: maxval = v
    347347
    348348                if minval is None or maxval is None:
    349349                    print 'Sorry no information about tag %s' %tag
    350                 else:   
     350                else:
    351351                    print '    Quantity %s, tag %s: min = %12.8f, max = %12.8f'\
    352352                          %(name, tag, minval, maxval)
    353                        
    354 
    355 
    356                
    357 
    358        
    359        
    360            
    361        
     353
     354
     355
     356
     357
     358
     359
     360
     361
    362362
    363363    def get_name(self):
    … …  
    412412        if yieldstep is None:
    413413            yieldstep = max_timestep
    414         else:
     414        else:
    415415            yieldstep = float(yieldstep)
    416416
    … …  
    463463            #Yield results
    464464            if finaltime is not None and abs(self.time - finaltime) < epsilon:
    465            
    466                 #FIXME: There is a rare situation where the 
     465
     466                #FIXME: There is a rare situation where the
    467467                #final time step is stored twice. Can we make a test?
    468468
    … …  
    637637    except:
    638638        import os
    639         if os.name == 'posix' and os.uname()[4] == 'x86_64':       
     639        if os.name == 'posix' and os.uname()[4] == 'x86_64':
    640640            pass
    641641            #Psyco isn't supported on 64 bit systems, but it doesn't matter

  • inundation/ga/storm_surge/pyvolution/island.py

    r912 r1280  
    66
    77######################
    8 # Module imports 
     8# Module imports
    99#
    1010from os import sep, path
    … …  
    1818
    1919#Create basic mesh
    20 N = 20
     20N = 40
    2121points, vertices, boundary = rectangular(N, N, 100, 100)
    2222
    2323#Create shallow water domain
    2424domain = Domain(points, vertices, boundary)
    25 domain.store = True
     25domain.store = False
    2626domain.smooth = False
    27 domain.visualise = False
     27domain.visualise = True
    2828domain.set_name('island')
    2929print 'Output being written to ' + domain.get_datadir() + sep + \
    3030              domain.get_name() + '.' + domain.format
    31        
     31
    3232
    3333domain.default_order=2
    … …  
    3939#
    4040# beta_h == 1.0 means that the largest gradients (on h) are allowed
    41 # beta_h == 0.0 means that constant (1st order) gradients are introduced 
     41# beta_h == 0.0 means that constant (1st order) gradients are introduced
    4242# on h. This is equivalent to the constant depth used previously.
    43 domain.beta_h = 0.2
     43domain.beta_h = 0.9
    4444
    4545
    … …  
    4949    z = 0*x
    5050    for i in range(len(x)):
    51         z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )   
    52     return z   
     51        z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
     52    return z
    5353
    5454domain.set_quantity('friction', 0.0)
    5555domain.set_quantity('elevation', island)
    5656domain.set_quantity('stage', 1)
    57        
     57
    5858
    5959
    … …  
    6969#Evolution
    7070import time
    71 for t in domain.evolve(yieldstep = 1, finaltime = 20):
     71for t in domain.evolve(yieldstep = 1, finaltime = 40):
    7272    domain.write_time()
    7373
    74 print 'Done'   
    75 
     74print 'Done'

  • inundation/ga/storm_surge/pyvolution/netherlands.py

    r1276 r1280  
    7979N = 130 #size = 33800
    8080N = 600 #Size = 720000
    81 N = 60
     81N = 30
    8282
    8383#N = 15

  • inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

    r1278 r1280  
    22
    33class Surface:
    4     def __init__(self,domain,scale_z=1.0):
     4    def __init__(self,domain,scale_z=1.0,range_z=None):
    55        """Create visualisation of domain
    66        """
    77
    88        self.frame  = frame()
    9         autoscale=0
     9
    1010        self.bed_model   = faces(frame=self.frame)
    1111        self.stage_model = faces(frame=self.frame)
    … …  
    1414
    1515        self.vertices = domain.vertex_coordinates
    16         self.bed   = domain.quantities['elevation'].vertex_values
     16
    1717        self.stage = domain.quantities['stage'].vertex_values
    18         self.xmomentum = domain.quantities['xmomentum'].vertex_values
    19         self.ymomentum = domain.quantities['ymomentum'].vertex_values
     18
     19        try:
     20            self.bed   = domain.quantities['elevation'].vertex_values
     21
     22            self.xmomentum = domain.quantities['xmomentum'].vertex_values
     23            self.ymomentum = domain.quantities['ymomentum'].vertex_values
     24        except:
     25            self.bed = None
     26            self.xmomentum = None
     27            self.ymomentum = None
     28
    2029
    2130        self.max_x = max(max(self.vertices[:,0],self.vertices[:,2],self.vertices[:,4]))
    … …  
    2837        self.range_xy = max(self.range_x, self.range_y)
    2938
    30         self.max_bed = max(max(self.bed))
    31         self.min_bed = min(min(self.bed))
    32         self.range_bed = max(self.max_bed - self.min_bed, 1e-10)*2.0
    33 
    34         print 'min_bed=',self.min_bed
    35         print 'max_bed=',self.max_bed
    36         print 'range_bed=',self.range_bed
     39        if self.bed != None and range_z == None:
     40            self.max_z = max(max(self.bed))
     41            self.min_z = min(min(self.bed))
     42            self.range_z = max(self.max_z - self.min_z, 1.0e-10)
     43
     44            print 'min_bed=',self.min_z
     45            print 'max_bed=',self.max_z
     46            print 'range_z=',self.range_z
     47        else:
     48            self.max_z = range_z/2.0
     49            self.min_z = -range_z/2.0
     50            self.range_z = max(self.max_z - self.min_z, 1.0e-10)
     51
    3752
    3853        #print self.max_x, self.min_x, self.max_y, self.min_y,
    … …  
    4863        #print 'shape of stage',shape(self.stage)
    4964
    50         self.border_model = curve(frame = self.frame, pos=[(0,0),(0,1),(1,1),(1,0),(0,0)])
     65        self.border_model = curve(frame = self.frame,
     66                                  pos=[(0,0),(0,1),(1,1),(1,0),(0,0)])
    5167        self.timer=label(pos=(0.75,0.5,0.5),text='Time=%10.5e'%self.domain.time)
    5268
     69        #scene.autoscale=0
    5370        #self.update_all()
    5471
    … …  
    6784        c1 = 0.3
    6885        c2 = 0.3
    69         self.update_arrays(self.bed,  (c0,c1,c2) )
     86        try:
     87            self.update_arrays(self.bed,  (c0,c1,c2) )
     88        except:
     89            pass
    7090
    7191        #print 'update bed image'
    … …  
    118138        min_y = self.min_y
    119139        range_xy = self.range_xy
    120         range_bed = self.range_bed
     140        range_z = self.range_z
    121141
    122142        vertices = self.vertices
    … …  
    127147        try:
    128148            update_arrays_weave(vertices,quantity,col,pos,normals,colour,N,
    129                                 min_x,min_y,range_xy,range_bed,scale_z)
     149                                min_x,min_y,range_xy,range_z,scale_z)
    130150        except:
    131151            update_arrays_python(vertices,quantity,col,pos,normals,colour,N,
    132                                  min_x,min_y,range_xy,range_bed,scale_z)
     152                                 min_x,min_y,range_xy,range_z,scale_z)
    133153
    134154
    135155
    136156def  update_arrays_python(vertices,quantity,col,pos,normals,colour,N,
    137                           min_x,min_y,range_xy,range_bed,scale_z):
     157                          min_x,min_y,range_xy,range_z,scale_z):
    138158
    139159    from math import sqrt
    … …  
    147167            v[j,0] = (vertices[i,2*j  ]-min_x)/range_xy
    148168            v[j,1] = (vertices[i,2*j+1]-min_y)/range_xy
    149             v[j,2] = quantity[i,j]/range_bed*scale_z
     169            v[j,2] = quantity[i,j]/range_z*scale_z*0.5
    150170
    151171        v10 = v[1,:]-v[0,:]
    … …  
    190210
    191211def  update_arrays_weave(vertices,quantity,col,pos,normals,colour,
    192                          N,min_x,min_y,range_xy,range_bed,scale_z):
     212                         N,min_x,min_y,range_xy,range_z,scale_z):
    193213
    194214    import weave
    … …  
    208228                v(j,0) = (vertices(i,2*j  )-min_x)/range_xy;
    209229                v(j,1) = (vertices(i,2*j+1)-min_y)/range_xy;
    210                 v(j,2) = quantity(i,j)/range_bed*scale_z;
     230                v(j,2) = quantity(i,j)/range_z*scale_z*0.5;
    211231            }
    212232
    … …  
    256276
    257277    weave.inline(code1, ['vertices','quantity','col','pos','normals','colour','N',
    258                          'min_x','min_y','range_xy','range_bed','scale_z','v','normal','v10','v20'],
     278                         'min_x','min_y','range_xy','range_z','scale_z','v','normal','v10','v20'],
    259279                 type_converters = converters.blitz, compiler='gcc');
    260280
    … …  
    286306
    287307
    288 def create_surface(domain):
    289 
    290     surface = Surface(domain)
     308def create_surface(domain,range_z=None):
     309
     310    surface = Surface(domain,range_z=range_z)
    291311    domain.surface = surface
    292312
    293     surface.update_bed()
     313    #surface.update_bed()
    294314    surface.update_stage()
    295315    #surface.update_xmomentum()

  • inundation/ga/storm_surge/pyvolution/util.py

    r1207 r1280  
    22
    33It is also a clearing house for functions that may later earn a module
    4 of their own. 
     4of their own.
    55"""
    66
    … …  
    1717    v2 = v[1]/l
    1818
    19     try:       
     19    try:
    2020        theta = acos(v1)
    2121    except ValueError, e:
    … …  
    3131            theta = 3.1415926535897931
    3232        print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
    33               %(v1, theta)   
    34              
     33              %(v1, theta)
     34
    3535    if v2 < 0:
    3636        #Quadrant 3 or 4
    … …  
    4343    """Compute difference between angle of vector x0, y0 and x1, y1.
    4444    This is used for determining the ordering of vertices,
    45     e.g. for checking if they are counter clockwise. 
    46    
     45    e.g. for checking if they are counter clockwise.
     46
    4747    Always return a positive value
    4848    """
    49        
     49
    5050    from math import pi
    51    
     51
    5252    a0 = angle(v0)
    5353    a1 = angle(v1)
    … …  
    5656    if a0 < a1:
    5757        a0 += 2*pi
    58            
     58
    5959    return a0-a1
    6060
    … …  
    6565
    6666
    67 def point_on_line(x, y, x0, y0, x1, y1): 
    68     """Determine whether a point is on a line segment 
    69    
     67def point_on_line(x, y, x0, y0, x1, y1):
     68    """Determine whether a point is on a line segment
     69
    7070    Input: x, y, x0, x0, x1, y1: where
    7171        point is given by x, y
    7272        line is given by (x0, y0) and (x1, y1)
    73        
    74     """
    75        
     73
     74    """
     75
    7676    from Numeric import array, dot, allclose
    7777    from math import sqrt
    7878
    79     a = array([x - x0, y - y0]) 
     79    a = array([x - x0, y - y0])
    8080    a_normal = array([a[1], -a[0]])
    81                        
     81
    8282    b = array([x1 - x0, y1 - y0])
    83    
     83
    8484    if dot(a_normal, b) == 0:
    8585        #Point is somewhere on the infinite extension of the line
    8686
    87         len_a = sqrt(sum(a**2))               
    88         len_b = sqrt(sum(b**2))                                 
     87        len_a = sqrt(sum(a**2))
     88        len_b = sqrt(sum(b**2))
    8989        if dot(a, b) >= 0 and len_a <= len_b:
    9090           return True
    91         else:   
     91        else:
    9292           return False
    9393    else:
    94       return False 
    95    
     94      return False
     95
    9696def ensure_numeric(A, typecode = None):
    9797    """Ensure that sequence is a Numeric array.
    … …  
    117117        if type(A) == ArrayType:
    118118            if A.typecode == typecode:
    119                 return array(A)           
     119                return array(A)
    120120            else:
    121121                return A.astype(typecode)
    122122        else:
    123123            return array(A).astype(typecode)
    124        
    125    
     124
     125
    126126
    127127def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False):
    … …  
    133133    #In fact, this is where origin should be converted to that of domain
    134134    #Also, check that file covers domain fully.
    135    
     135
    136136    assert type(filename) == type(''),\
    137137               'First argument to File_function must be a string'
    … …  
    148148
    149149    if domain is not None:
    150         quantities = domain.conserved_quantities   
     150        quantities = domain.conserved_quantities
    151151    else:
    152152        quantities = None
    153        
    154     #Choose format 
     153
     154    #Choose format
    155155    #FIXME: Maybe these can be merged later on
    156156    if line[:3] == 'CDF':
    157157        return File_function_NetCDF(filename, domain, quantities, interpolation_points, verbose = verbose)
    158158    else:
    159         return File_function_ASCII(filename, domain, quantities, interpolation_points)   
    160            
    161              
     159        return File_function_ASCII(filename, domain, quantities, interpolation_points)
     160
     161
    162162class File_function_NetCDF:
    163     """Read time history of spatial data from NetCDF sww file and 
    164     return a callable object f(t,x,y) 
    165     which will return interpolated values based on the input file.   
    166    
     163    """Read time history of spatial data from NetCDF sww file and
     164    return a callable object f(t,x,y)
     165    which will return interpolated values based on the input file.
     166
    167167    x, y may be either scalars or vectors
    168    
     168
    169169    #FIXME: More about format, interpolation  and order of quantities
    170    
     170
    171171    The quantities returned by the callable objects are specified by
    172     the list quantities which must contain the names of the 
     172    the list quantities which must contain the names of the
    173173    quantities to be returned and also reflect the order, e.g. for
    174     the shallow water wave equation, on would have 
     174    the shallow water wave equation, on would have
    175175    quantities = ['stage', 'xmomentum', 'ymomentum']
    176    
    177     interpolation_points decides at which points interpolated 
     176
     177    interpolation_points decides at which points interpolated
    178178    quantities are to be computed whenever object is called.
    179    
    180    
    181    
     179
     180
     181
    182182    If None, return average value
    183183    """
    184    
     184
    185185    def  __init__(self, filename, domain=None, quantities=None,
    186186                  interpolation_points=None, verbose = False):
    … …  
    189189        If domain is specified, model time (domain.starttime)
    190190        will be checked and possibly modified
    191        
     191
    192192        All times are assumed to be in UTC
    193193        """
    … …  
    196196        #(both in UTM coordinates)
    197197        #If not - modify those from file to match domain
    198        
     198
    199199        import time, calendar
    200200        from config import time_format
    201         from Scientific.IO.NetCDF import NetCDFFile     
     201        from Scientific.IO.NetCDF import NetCDFFile
    202202
    203203        #Open NetCDF file
    … …  
    213213             if not fid.variables.has_key(quantity):
    214214                 missing.append(quantity)
    215                  
     215
    216216        if len(missing) > 0:
    217217            msg = 'Quantities %s could not be found in file %s'\
    218218                  %(str(missing), filename)
    219219            raise msg
    220                  
     220
    221221
    222222        #Get first timestep
    223223        try:
    224224            self.starttime = fid.starttime[0]
    225         except ValueError:   
     225        except ValueError:
    226226            msg = 'Could not read starttime from file %s' %filename
    227227            raise msg
    … …  
    229229        #Get origin
    230230        self.xllcorner = fid.xllcorner[0]
    231         self.yllcorner = fid.yllcorner[0]           
    232    
     231        self.yllcorner = fid.yllcorner[0]
     232
    233233        self.number_of_values = len(quantities)
    234234        self.fid = fid
    235235        self.filename = filename
    236         self.quantities = quantities
    237         self.domain = domain
    238         self.interpolation_points = interpolation_points
     236        self.quantities = quantities
     237        self.domain = domain
     238        self.interpolation_points = interpolation_points
    239239
    240240        if domain is not None:
    … …  
    249249                if verbose: print 'Domain starttime is now set to %f' %domain.starttime
    250250
    251         #Read all data in and produce values for desired data points at each timestep
    252         self.spatial_interpolation(interpolation_points, verbose = verbose) 
     251        #Read all data in and produce values for desired data points at each timestep
     252        self.spatial_interpolation(interpolation_points, verbose = verbose)
    253253        fid.close()
    254254
    255255    def spatial_interpolation(self, interpolation_points, verbose = False):
    256         """For each timestep in fid: read surface, interpolate to desired points 
    257         and store values for use when object is called.
     256        """For each timestep in fid: read surface, interpolate to desired points
     257           and store values for use when object is called.
    258258        """
    259        
     259
    260260        from Numeric import array, zeros, Float, alltrue, concatenate, reshape
    261261
    262262        from config import time_format
    263         from least_squares import Interpolation
     263        from least_squares import Interpolation
    264264        import time, calendar
    265        
     265
    266266        fid = self.fid
    267267
    … …  
    273273        triangles = fid.variables['volumes'][:]
    274274        time = fid.variables['time'][:]
    275        
    276         #Check
     275
     276        #Check
    277277        msg = 'File %s must list time as a monotonuosly ' %self.filename
    278278        msg += 'increasing sequence'
    279         assert alltrue(time[1:] - time[:-1] > 0 ), msg 
    280        
    281 
    282         if interpolation_points is not None: 
    283        
     279        assert alltrue(time[1:] - time[:-1] > 0 ), msg
     280
     281
     282        if interpolation_points is not None:
     283
    284284            try:
    285                 P = ensure_numeric(interpolation_points) 
     285                P = ensure_numeric(interpolation_points)
    286286            except:
    287287                msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
    288288                msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
    289289                raise msg
    290            
    291            
     290
     291
    292292            self.T = time[:]  #Time assumed to be relative to starttime
    293             self.index = 0    #Initial time index       
    294 
    295             self.values = {}               
     293            self.index = 0    #Initial time index
     294
     295            self.values = {}
    296296            for name in self.quantities:
    297                 self.values[name] = zeros( (len(self.T), 
    298                                             len(interpolation_points)), 
    299                                             Float)                 
     297                self.values[name] = zeros( (len(self.T),
     298                                            len(interpolation_points)),
     299                                            Float)
    300300
    301301            #Build interpolator
    302302            if verbose: print 'Build interpolation matrix'
    303303            x = reshape(x, (len(x),1))
    304             y = reshape(y, (len(y),1))               
     304            y = reshape(y, (len(y),1))
    305305            vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
    306306
    307             interpol = Interpolation(vertex_coordinates, 
     307            interpol = Interpolation(vertex_coordinates,
    308308                                     triangles,
    309309                                     point_coordinates = P,
    310                                      alpha = 0, 
     310                                     alpha = 0,
    311311                                     verbose = verbose)
    312312
    313313
    314314            if verbose: print 'Interpolate'
    315             for i, t in enumerate(self.T):
     315            for i, t in enumerate(self.T):
    316316                #Interpolate quantities at this timestep
    317317                #print ' time step %d of %d' %(i, len(self.T))
    318318                for name in self.quantities:
    319                     self.values[name][i, :] =\
    320                     interpol.interpolate(fid.variables[name][i,:])   
     319                    self.values[name][i, :] =\
     320                                   interpol.interpolate(fid.variables[name][i,:])
    321321
    322322            #Report
    … …  
    327327                print '  Reference:'
    328328                print '    Lower left corner: [%f, %f]'\
    329                       %(self.xllcorner, self.yllcorner) 
     329                      %(self.xllcorner, self.yllcorner)
    330330                print '    Start time (file):   %f' %self.starttime
    331331                #print '    Start time (domain): %f' %domain.starttime
    … …  
    342342                    print '    %s in [%f, %f]' %(name, min(q), max(q))
    343343
    344                    
     344
    345345                print '  Interpolation points (xi, eta):'\
    346346                      'number of points == %d ' %P.shape[0]
    … …  
    349349                print '  Interpolated quantities (over all timesteps):'
    350350                for name in self.quantities:
    351                     q = self.values[name][:].flat                   
     351                    q = self.values[name][:].flat
    352352                    print '    %s at interpolation points in [%f, %f]'\
    353                           %(name, min(q), max(q))               
    354                
    355                
    356                
     353                          %(name, min(q), max(q))
    357354
    358355                print '------------------------------------------------'
    … …  
    367364    def __call__(self, t, x=None, y=None, point_id = None):
    368365        """Evaluate f(t, point_id)
    369        
    370         Inputs:
    371           t: time - Model time (tau = domain.starttime-self.starttime+t) 
    372                     must lie within existing timesteps   
    373           point_id: index of one of the preprocessed points.
    374                     If point_id is None all preprocessed points are computed
    375 
    376           FIXME: point_id could also be a slice     
    377           FIXME: One could allow arbitrary x, y coordinates
    378           (requires computation of a new interpolator)
    379           Maybe not,.,.
    380           FIXME: What if x and y are vectors?
     366
     367           Inputs:
     368           t: time - Model time (tau = domain.starttime-self.starttime+t)
     369                     must lie within existing timesteps
     370           point_id: index of one of the preprocessed points.
     371                  If point_id is None all preprocessed points are computed
     372
     373           FIXME: point_id could also be a slice
     374           FIXME: One could allow arbitrary x, y coordinates
     375           (requires computation of a new interpolator)
     376           Maybe not,.,.
     377           FIXME: What if x and y are vectors?
    381378        """
    382379
    … …  
    391388
    392389        #Find time tau such that
    393         #
    394         # domain.starttime + t == self.starttime + tau
    395        
    396         if self.domain is not None:
     390        #
     391        # domain.starttime + t == self.starttime + tau
     392
     393        if self.domain is not None:
    397394            tau = self.domain.starttime-self.starttime+t
    398         else:
     395        else:
    399396            #print 'DOMAIN IS NONE!!!!!!!!!'
    400397            tau = t
    … …  
    403400        #print 'S start', self.starttime
    404401        #print 't', t
    405         #print 'tau', tau             
    406                
     402        #print 'tau', tau
     403
    407404        msg = 'Time interval derived from file %s [%s:%s]'\
    408405              %(self.filename, self.T[0], self.T[1])
    … …  
    411408
    412409        if tau < self.T[0]: raise msg
    413         if tau > self.T[-1]: raise msg       
     410        if tau > self.T[-1]: raise msg
    414411
    415412
    … …  
    424421            #Protect against case where tau == T[-1] (last time)
    425422            # - also works in general when tau == T[i]
    426             ratio = 0 
     423            ratio = 0
    427424        else:
    428             #t is now between index and index+1           
     425            #t is now between index and index+1
    429426            ratio = (tau - self.T[self.index])/\
    430427                    (self.T[self.index+1] - self.T[self.index])
    … …  
    432429
    433430
    434                
     431
    435432        #Compute interpolated values
    436433        q = zeros( len(self.quantities), Float)
    437434
    438435        for i, name in enumerate(self.quantities):
    439             Q = self.values[name]   
     436            Q = self.values[name]
    440437
    441438            if ratio > 0:
    … …  
    450447        #Return vector of interpolated values
    451448        return q
    452                
    453              
     449
     450
    454451class File_function_ASCII:
    455452    """Read time series from file and return a callable object:
    456     f(t,x,y)
    457     which will return interpolated values based on the input file.
    458 
    459     The file format is assumed to be either two fields separated by a comma:
    460 
    461         time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
    462 
    463     or four comma separated fields
    464 
    465         time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
    466 
    467     In either case, the callable object will return a tuple of
    468     interpolated values, one each value stated in the file.
    469 
    470    
    471     E.g
    472 
    473       31/08/04 00:00:00, 1.328223 0 0
    474       31/08/04 00:15:00, 1.292912 0 0
    475 
    476     will provide a time dependent function f(t,x=None,y=None),
    477     ignoring x and y, which returns three values per call.
    478 
    479    
    480     NOTE: At this stage the function is assumed to depend on
    481     time only, i.e  no spatial dependency!!!!!
    482     When that is needed we can use the least_squares interpolation.
    483    
    484     #FIXME: This should work with netcdf (e.g. sww) and thus render the
    485     #spatio-temporal boundary condition in shallow water fully general
    486    
    487     #FIXME: Specified quantites not used here -
    488     #always return whatever is in the file
    489     """
    490 
    491    
     453       f(t,x,y)
     454       which will return interpolated values based on the input file.
     455
     456       The file format is assumed to be either two fields separated by a comma:
     457
     458       time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
     459
     460       or four comma separated fields
     461
     462       time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
     463
     464       In either case, the callable object will return a tuple of
     465       interpolated values, one each value stated in the file.
     466
     467
     468       E.g
     469
     470       31/08/04 00:00:00, 1.328223 0 0
     471       31/08/04 00:15:00, 1.292912 0 0
     472
     473       will provide a time dependent function f(t,x=None,y=None),
     474       ignoring x and y, which returns three values per call.
     475
     476
     477       NOTE: At this stage the function is assumed to depend on
     478       time only, i.e  no spatial dependency!!!!!
     479       When that is needed we can use the least_squares interpolation.
     480
     481       #FIXME: This should work with netcdf (e.g. sww) and thus render the
     482       #spatio-temporal boundary condition in shallow water fully general
     483
     484       #FIXME: Specified quantites not used here -
     485       #always return whatever is in the file
     486    """
     487
     488
    492489    def __init__(self, filename, domain=None, quantities = None, interpolation_points=None):
    493490        """Initialise callable object from a file.
    … …  
    507504        line = fid.readline()
    508505        fid.close()
    509    
     506
    510507        fields = line.split(',')
    511508        msg = 'File %s must have the format date, value0 value1 value2 ...'
    … …  
    516513        try:
    517514            starttime = calendar.timegm(time.strptime(fields[0], time_format))
    518         except ValueError:   
     515        except ValueError:
    519516            msg = 'First field in file %s must be' %filename
    520517            msg += ' date-time with format %s.\n' %time_format
    521             msg += 'I got %s instead.' %fields[0] 
     518            msg += 'I got %s instead.' %fields[0]
    522519            raise msg
    523520
    … …  
    529526
    530527        q = ensure_numeric(values)
    531        
     528
    532529        msg = 'ERROR: File must contain at least one independent value'
    533530        assert len(q.shape) == 1, msg
    534            
     531
    535532        self.number_of_values = len(q)
    536533
    537534        self.filename = filename
    538535        self.starttime = starttime
    539         self.domain = domain
     536        self.domain = domain
    540537
    541538        if domain is not None:
    … …  
    549546                ##if verbose: print msg
    550547                domain.starttime = self.starttime #Modifying model time
    551        
     548
    552549            #if domain.starttime is None:
    553550            #    domain.starttime = self.starttime
    … …  
    563560            #        domain.starttime = self.starttime #Modifying model time
    564561
    565         if mode == 2:       
     562        if mode == 2:
    566563            self.read_times() #Now read all times in.
    567564        else:
    568565            raise 'x,y dependency not yet implemented'
    569566
    570        
     567
    571568    def read_times(self):
    572         """Read time and values 
     569        """Read time and values
    573570        """
    574571        from Numeric import zeros, Float, alltrue
    575572        from config import time_format
    576573        import time, calendar
    577        
    578         fid = open(self.filename)       
     574
     575        fid = open(self.filename)
    579576        lines = fid.readlines()
    580577        fid.close()
    581        
     578
    582579        N = len(lines)
    583580        d = self.number_of_values
    … …  
    585582        T = zeros(N, Float)       #Time
    586583        Q = zeros((N, d), Float)  #Values
    587        
     584
    588585        for i, line in enumerate(lines):
    589586            fields = line.split(',')
    … …  
    603600        self.index = 0 #Initial index
    604601
    605        
     602
    606603    def __repr__(self):
    607604        return 'File function'
    … …  
    613610        result is a vector of same length as x and y
    614611        FIXME: Naaaa
    615        
    616         FIXME: Rethink semantics when x,y are vectors.
     612
     613        FIXME: Rethink semantics when x,y are vectors.
    617614        """
    618615
    619616        from math import pi, cos, sin, sqrt
    620        
     617
    621618
    622619        #Find time tau such that
    623         #
    624         # domain.starttime + t == self.starttime + tau
    625        
    626         if self.domain is not None:
     620        #
     621        # domain.starttime + t == self.starttime + tau
     622
     623        if self.domain is not None:
    627624            tau = self.domain.starttime-self.starttime+t
    628         else:
     625        else:
    629626            tau = t
    630        
    631                
     627
     628
    632629        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
    633630              %(self.filename, self.T[0], self.T[1], tau)
    634631        if tau < self.T[0]: raise msg
    635         if tau > self.T[-1]: raise msg       
     632        if tau > self.T[-1]: raise msg
    636633
    637634        oldindex = self.index
    … …  
    647644        #Compute interpolated values
    648645        q = self.Q[self.index,:] +\
    649             ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
     646            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
    650647
    651648        #Return vector of interpolated values
    … …  
    662659                N = len(x)
    663660                assert len(y) == N, 'x and y must have same length'
    664              
     661
    665662                res = []
    666663                for col in q:
    667664                    res.append(col*ones(N, Float))
    668            
     665
    669666                return res
    670            
    671            
    672 #FIXME: TEMP       
     667
     668
     669#FIXME: TEMP
    673670class File_function_copy:
    674671    """Read time series from file and return a callable object:
    675     f(t,x,y)
    676     which will return interpolated values based on the input file.
    677 
    678     The file format is assumed to be either two fields separated by a comma:
    679 
    680         time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
    681 
    682     or four comma separated fields
    683 
    684         time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
    685 
    686     In either case, the callable object will return a tuple of
    687     interpolated values, one each value stated in the file.
    688 
    689    
    690     E.g
    691 
    692       31/08/04 00:00:00, 1.328223 0 0
    693       31/08/04 00:15:00, 1.292912 0 0
    694 
    695     will provide a time dependent function f(t,x=None,y=None),
    696     ignoring x and y, which returns three values per call.
    697 
    698    
    699     NOTE: At this stage the function is assumed to depend on
    700     time only, i.e  no spatial dependency!!!!!
    701     When that is needed we can use the least_squares interpolation.
    702    
    703     #FIXME: This should work with netcdf (e.g. sww) and thus render the
    704     #spatio-temporal boundary condition in shallow water fully general
    705     """
    706 
    707    
     672       f(t,x,y)
     673       which will return interpolated values based on the input file.
     674
     675       The file format is assumed to be either two fields separated by a comma:
     676
     677           time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
     678
     679       or four comma separated fields
     680
     681           time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
     682
     683       In either case, the callable object will return a tuple of
     684       interpolated values, one each value stated in the file.
     685
     686
     687       E.g
     688
     689         31/08/04 00:00:00, 1.328223 0 0
     690         31/08/04 00:15:00, 1.292912 0 0
     691
     692       will provide a time dependent function f(t,x=None,y=None),
     693       ignoring x and y, which returns three values per call.
     694
     695
     696       NOTE: At this stage the function is assumed to depend on
     697       time only, i.e  no spatial dependency!!!!!
     698       When that is needed we can use the least_squares interpolation.
     699
     700       #FIXME: This should work with netcdf (e.g. sww) and thus render the
     701       #spatio-temporal boundary condition in shallow water fully general
     702    """
     703
     704
    708705    def __init__(self, filename, domain=None):
    709706        """Initialise callable object from a file.
    … …  
    742739        try:
    743740            starttime = calendar.timegm(time.strptime(fields[0], time_format))
    744         except ValueError:   
     741        except ValueError:
    745742            msg = 'First field in file %s must be' %filename
    746743            msg += ' date-time with format %s.\n' %time_format
    747             msg += 'I got %s instead.' %fields[0] 
     744            msg += 'I got %s instead.' %fields[0]
    748745            raise msg
    749746
    … …  
    755752
    756753        q = ensure_numeric(values)
    757        
     754
    758755        msg = 'ERROR: File must contain at least one independent value'
    759756        assert len(q.shape) == 1, msg
    760            
     757
    761758        self.number_of_values = len(q)
    762759
    763760        self.filename = filename
    764761        self.starttime = starttime
    765         self.domain = domain
     762        self.domain = domain
    766763
    767764        if domain is not None:
    … …  
    779776                    domain.starttime = self.starttime #Modifying model time
    780777
    781         if mode == 2:       
     778        if mode == 2:
    782779            self.read_times() #Now read all times in.
    783780        else:
    784781            raise 'x,y dependency not yet implemented'
    785782
    786        
     783
    787784    def read_times(self):
    788         """Read time and values 
     785        """Read time and values
    789786        """
    790787        from Numeric import zeros, Float, alltrue
    791788        from config import time_format
    792789        import time, calendar
    793        
    794         fid = open(self.filename)       
     790
     791        fid = open(self.filename)
    795792        lines = fid.readlines()
    796793        fid.close()
    797        
     794
    798795        N = len(lines)
    799796        d = self.number_of_values
    … …  
    801798        T = zeros(N, Float)       #Time
    802799        Q = zeros((N, d), Float)  #Values
    803        
     800
    804801        for i, line in enumerate(lines):
    805802            fields = line.split(',')
    … …  
    819816        self.index = 0 #Initial index
    820817
    821        
     818
    822819    def __repr__(self):
    823820        return 'File function'
    824821
    825        
     822
    826823
    827824    def __call__(self, t, x=None, y=None):
    … …  
    834831
    835832        from math import pi, cos, sin, sqrt
    836        
     833
    837834
    838835        #Find time tau such that
    839         #
    840         # domain.starttime + t == self.starttime + tau
    841        
    842         if self.domain is not None:
     836        #
     837        # domain.starttime + t == self.starttime + tau
     838
     839        if self.domain is not None:
    843840            tau = self.domain.starttime-self.starttime+t
    844         else:
     841        else:
    845842            tau = t
    846        
    847                
     843
     844
    848845        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
    849846              %(self.filename, self.T[0], self.T[1], tau)
    850847        if tau < self.T[0]: raise msg
    851         if tau > self.T[-1]: raise msg       
     848        if tau > self.T[-1]: raise msg
    852849
    853850        oldindex = self.index
    … …  
    863860        #Compute interpolated values
    864861        q = self.Q[self.index,:] +\
    865             ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
     862            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
    866863
    867864        #Return vector of interpolated values
    … …  
    878875                N = len(x)
    879876                assert len(y) == N, 'x and y must have same length'
    880              
     877
    881878                res = []
    882879                for col in q:
    883880                    res.append(col*ones(N, Float))
    884            
     881
    885882                return res
    886            
     883
    887884
    888885def read_xya(filename, format = 'netcdf'):
    … …  
    892889
    893890    Format can be either ASCII or NetCDF
    894    
     891
    895892    Return list of points, list of attributes and
    896893    attribute names if present otherwise None
    897894    """
    898895    #FIXME: Probably obsoleted by similar function in load_ASCII
    899    
     896
    900897    from Scientific.IO.NetCDF import NetCDFFile
    901898
    902     if format.lower() == 'netcdf': 
     899    if format.lower() == 'netcdf':
    903900        #Get NetCDF
    904        
    905         fid = NetCDFFile(filename, 'r') 
    906      
     901
     902        fid = NetCDFFile(filename, 'r')
     903
    907904        # Get the variables
    908905        points = fid.variables['points']
    … …  
    913910        #Don't close - arrays are needed outside this function,
    914911        #alternatively take a copy here
    915     else:   
     912    else:
    916913        #Read as ASCII file assuming that it is separated by whitespace
    917914        fid = open(filename)
    … …  
    930927            attribute_names = ['elevation']  #HACK
    931928
    932         attributes = {}   
     929        attributes = {}
    933930        for key in attribute_names:
    934931            attributes[key] = []
    … …  
    940937            for i, key in enumerate(attribute_names):
    941938                attributes[key].append( float(fields[2+i]) )
    942        
     939
    943940    return points, attributes
    944    
     941
    945942
    946943#####################################
    … …  
    950947
    951948
    952              
    953 #Redefine these to use separate_by_polygon which will put all inside indices 
     949
     950#Redefine these to use separate_by_polygon which will put all inside indices
    954951#in the first part of the indices array and outside indices in the last
    955952
    … …  
    957954    """See separate_points_by_polygon for documentation
    958955    """
    959    
    960     from Numeric import array, Float, reshape   
    961    
    962     if verbose: print 'Checking input to inside_polygon'   
     956
     957    from Numeric import array, Float, reshape
     958
     959    if verbose: print 'Checking input to inside_polygon'
    963960    try:
    964         points = ensure_numeric(points, Float)             
     961        points = ensure_numeric(points, Float)
    965962    except:
    966963        msg = 'Points could not be converted to Numeric array'
    967         raise msg       
    968        
     964        raise msg
     965
    969966    try:
    970         polygon = ensure_numeric(polygon, Float)           
     967        polygon = ensure_numeric(polygon, Float)
    971968    except:
    972969        msg = 'Polygon could not be converted to Numeric array'
    973         raise msg               
    974    
    975    
    976    
     970        raise msg
     971
     972
     973
    977974    if len(points.shape) == 1:
    978975        one_point = True
    979976        points = reshape(points, (1,2))
    980     else:       
     977    else:
    981978        one_point = False
    982979
    983     indices, count = separate_points_by_polygon(points, polygon, 
    984                                                 closed, verbose)     
    985    
     980    indices, count = separate_points_by_polygon(points, polygon,
     981                                                closed, verbose)
     982
    986983    if one_point:
    987984        return count == 1
    988985    else:
    989         return indices[:count]   
    990        
     986        return indices[:count]
     987
    991988def outside_polygon(points, polygon, closed = True, verbose = False):
    992989    """See separate_points_by_polygon for documentation
    993990    """
    994    
    995     from Numeric import array, Float, reshape   
    996    
    997     if verbose: print 'Checking input to outside_polygon'   
     991
     992    from Numeric import array, Float, reshape
     993
     994    if verbose: print 'Checking input to outside_polygon'
    998995    try:
    999         points = ensure_numeric(points, Float)             
     996        points = ensure_numeric(points, Float)
    1000997    except:
    1001998        msg = 'Points could not be converted to Numeric array'
    1002         raise msg       
    1003        
     999        raise msg
     1000
    10041001    try:
    1005         polygon = ensure_numeric(polygon, Float)           
     1002        polygon = ensure_numeric(polygon, Float)
    10061003    except:
    10071004        msg = 'Polygon could not be converted to Numeric array'
    1008         raise msg               
    1009    
    1010    
    1011    
     1005        raise msg
     1006
     1007
     1008
    10121009    if len(points.shape) == 1:
    10131010        one_point = True
    10141011        points = reshape(points, (1,2))
    1015     else:       
     1012    else:
    10161013        one_point = False
    10171014
    1018     indices, count = separate_points_by_polygon(points, polygon, 
    1019                                                 closed, verbose)     
    1020    
     1015    indices, count = separate_points_by_polygon(points, polygon,
     1016                                                closed, verbose)
     1017
    10211018    if one_point:
    10221019        return count != 1
    10231020    else:
    1024         return indices[count:][::-1]  #return reversed   
    1025        
    1026 
    1027 def separate_points_by_polygon(points, polygon, 
     1021        return indices[count:][::-1]  #return reversed
     1022
     1023
     1024def separate_points_by_polygon(points, polygon,
    10281025                               closed = True, verbose = False):
    10291026    """Determine whether points are inside or outside a polygon
    1030    
     1027
    10311028    Input:
    10321029       points - Tuple of (x, y) coordinates, or list of tuples
    10331030       polygon - list of vertices of polygon
    1034        closed - (optional) determine whether points on boundary should be 
    1035        regarded as belonging to the polygon (closed = True) 
    1036        or not (closed = False) 
    1037    
     1031       closed - (optional) determine whether points on boundary should be
     1032       regarded as belonging to the polygon (closed = True)
     1033       or not (closed = False)
     1034
    10381035    Outputs:
    1039        indices: array of same length as points with indices of points falling 
    1040        inside the polygon listed from the beginning and indices of points 
     1036       indices: array of same length as points with indices of points falling
     1037       inside the polygon listed from the beginning and indices of points
    10411038       falling outside listed from the end.
    1042        
     1039
    10431040       count: count of points falling inside the polygon
    1044        
     1041
    10451042       The indices of points inside are obtained as indices[:count]
    1046        The indices of points outside are obtained as indices[count:]       
    1047        
    1048        
     1043       The indices of points outside are obtained as indices[count:]
     1044
     1045
    10491046    Examples:
    10501047       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
    1051        
     1048
    10521049       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
    1053        will return the indices [0, 2, 1] and count == 2 as only the first 
     1050       will return the indices [0, 2, 1] and count == 2 as only the first
    10541051       and the last point are inside the unit square
    1055        
     1052
    10561053    Remarks:
    10571054       The vertices may be listed clockwise or counterclockwise and
    1058        the first point may optionally be repeated.   
     1055       the first point may optionally be repeated.
    10591056       Polygons do not need to be convex.
    1060        Polygons can have holes in them and points inside a hole is 
    1061        regarded as being outside the polygon.         
    1062                
    1063     Algorithm is based on work by Darel Finley, 
    1064     http://www.alienryderflex.com/polygon/ 
    1065    
    1066     """   
     1057       Polygons can have holes in them and points inside a hole is
     1058       regarded as being outside the polygon.
     1059
     1060    Algorithm is based on work by Darel Finley,
     1061    http://www.alienryderflex.com/polygon/
     1062
     1063    """
    10671064
    10681065    from Numeric import array, Float, reshape, Int, zeros
    1069    
    1070    
    1071     #Input checks
    1072     try:
    1073         points = ensure_numeric(points, Float)             
    1074     except:
    1075         msg = 'Points could not be converted to Numeric array'
    1076         raise msg       
    1077        
    1078     try:
    1079         polygon = ensure_numeric(polygon, Float)           
    1080     except:
    1081         msg = 'Polygon could not be converted to Numeric array'
    1082         raise msg               
    1083 
    1084     assert len(polygon.shape) == 2,\
    1085        'Polygon array must be a 2d array of vertices'
    1086 
    1087     assert polygon.shape[1] == 2,\
    1088        'Polygon array must have two columns'
    1089 
    1090     assert len(points.shape) == 2,\
    1091        'Points array must be a 2d array'
    1092        
    1093     assert points.shape[1] == 2,\
    1094        'Points array must have two columns'
    1095        
    1096     N = polygon.shape[0] #Number of vertices in polygon
    1097     M = points.shape[0]  #Number of points
    1098    
    1099     px = polygon[:,0]
    1100     py = polygon[:,1]   
    1101 
    1102     #Used for an optimisation when points are far away from polygon
    1103     minpx = min(px); maxpx = max(px)
    1104     minpy = min(py); maxpy = max(py)   
    1105 
    1106 
    1107     #Begin main loop
    1108     indices = zeros(M, Int)
    1109    
    1110     inside_index = 0    #Keep track of points inside
    1111     outside_index = M-1 #Keep track of points outside (starting from end)
    1112    
    1113     for k in range(M):
    1114 
    1115         if verbose:
    1116             if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
    1117        
    1118         #for each point   
    1119         x = points[k, 0]
    1120         y = points[k, 1]       
    1121 
    1122         inside = False
    1123 
    1124         if not x > maxpx or x < minpx or y > maxpy or y < minpy:
    1125             #Check polygon
    1126             for i in range(N):
    1127                 j = (i+1)%N
    1128                
    1129                 #Check for case where point is contained in line segment       
    1130                 if point_on_line(x, y, px[i], py[i], px[j], py[j]):
    1131                     if closed:
    1132                         inside = True
    1133                     else:       
    1134                         inside = False
    1135                     break
    1136                 else:   
    1137                     #Check if truly inside polygon                 
    1138                     if py[i] < y and py[j] >= y or\
    1139                        py[j] < y and py[i] >= y:
    1140                         if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
    1141                             inside = not inside
    1142        
    1143         if inside:
    1144             indices[inside_index] = k
    1145             inside_index += 1
    1146         else:
    1147             indices[outside_index] = k
    1148             outside_index -= 1   
    1149            
    1150     return indices, inside_index
    1151 
    1152 
    1153 def separate_points_by_polygon_c(points, polygon,
    1154                                  closed = True, verbose = False):
    1155     """Determine whether points are inside or outside a polygon
    1156 
    1157     C-wrapper
    1158     """   
    1159 
    1160     from Numeric import array, Float, reshape, zeros, Int
    1161    
    1162 
    1163     if verbose: print 'Checking input to separate_points_by_polygon'   
     1066
     1067
    11641068    #Input checks
    11651069    try:
    … …  
    11691073        raise msg
    11701074
    1171     #if verbose: print 'Checking input to separate_points_by_polygon 2'
    11721075    try:
    11731076        polygon = ensure_numeric(polygon, Float)
    11741077    except:
    11751078        msg = 'Polygon could not be converted to Numeric array'
    1176         raise msg               
    1177 
    1178     if verbose: print 'check'
    1179    
     1079        raise msg
     1080
    11801081    assert len(polygon.shape) == 2,\
    11811082       'Polygon array must be a 2d array of vertices'
    … …  
    11861087    assert len(points.shape) == 2,\
    11871088       'Points array must be a 2d array'
    1188        
     1089
    11891090    assert points.shape[1] == 2,\
    11901091       'Points array must have two columns'
    1191        
    1192     N = polygon.shape[0] #Number of vertices in polygon 
     1092
     1093    N = polygon.shape[0] #Number of vertices in polygon
    11931094    M = points.shape[0]  #Number of points
    11941095
     1096    px = polygon[:,0]
     1097    py = polygon[:,1]
     1098
     1099    #Used for an optimisation when points are far away from polygon
     1100    minpx = min(px); maxpx = max(px)
     1101    minpy = min(py); maxpy = max(py)
     1102
     1103
     1104    #Begin main loop
     1105    indices = zeros(M, Int)
     1106
     1107    inside_index = 0    #Keep track of points inside
     1108    outside_index = M-1 #Keep track of points outside (starting from end)
     1109
     1110    for k in range(M):
     1111
     1112        if verbose:
     1113            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
     1114
     1115        #for each point
     1116        x = points[k, 0]
     1117        y = points[k, 1]
     1118
     1119        inside = False
     1120
     1121        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
     1122            #Check polygon
     1123            for i in range(N):
     1124                j = (i+1)%N
     1125
     1126                #Check for case where point is contained in line segment
     1127                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
     1128                    if closed:
     1129                        inside = True
     1130                    else:
     1131                        inside = False
     1132                    break
     1133                else:
     1134                    #Check if truly inside polygon
     1135                    if py[i] < y and py[j] >= y or\
     1136                       py[j] < y and py[i] >= y:
     1137                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
     1138                            inside = not inside
     1139
     1140        if inside:
     1141            indices[inside_index] = k
     1142            inside_index += 1
     1143        else:
     1144            indices[outside_index] = k
     1145            outside_index -= 1
     1146
     1147    return indices, inside_index
     1148
     1149
     1150def separate_points_by_polygon_c(points, polygon,
     1151                                 closed = True, verbose = False):
     1152    """Determine whether points are inside or outside a polygon
     1153
     1154    C-wrapper
     1155    """
     1156
     1157    from Numeric import array, Float, reshape, zeros, Int
     1158
     1159
     1160    if verbose: print 'Checking input to separate_points_by_polygon'
     1161    #Input checks
     1162    try:
     1163        points = ensure_numeric(points, Float)
     1164    except:
     1165        msg = 'Points could not be converted to Numeric array'
     1166        raise msg
     1167
     1168    #if verbose: print 'Checking input to separate_points_by_polygon 2'
     1169    try:
     1170        polygon = ensure_numeric(polygon, Float)
     1171    except:
     1172        msg = 'Polygon could not be converted to Numeric array'
     1173        raise msg
     1174
     1175    if verbose: print 'check'
     1176
     1177    assert len(polygon.shape) == 2,\
     1178       'Polygon array must be a 2d array of vertices'
     1179
     1180    assert polygon.shape[1] == 2,\
     1181       'Polygon array must have two columns'
     1182
     1183    assert len(points.shape) == 2,\
     1184       'Points array must be a 2d array'
     1185
     1186    assert points.shape[1] == 2,\
     1187       'Points array must have two columns'
     1188
     1189    N = polygon.shape[0] #Number of vertices in polygon
     1190    M = points.shape[0]  #Number of points
     1191
    11951192    from util_ext import separate_points_by_polygon
    11961193
    … …  
    11991196    indices = zeros( M, Int )
    12001197
    1201     if verbose: print 'Calling C-version of inside poly' 
     1198    if verbose: print 'Calling C-version of inside poly'
    12021199    count = separate_points_by_polygon(points, polygon, indices,
    12031200                                       int(closed), int(verbose))
    12041201
    1205     return indices, count 
     1202    return indices, count
    12061203
    12071204
    12081205
    12091206class Polygon_function:
    1210     """Create callable object f: x,y -> z, where a,y,z are vectors and 
    1211     where f will return different values depending on whether x,y belongs 
     1207    """Create callable object f: x,y -> z, where a,y,z are vectors and
     1208    where f will return different values depending on whether x,y belongs
    12121209    to specified polygons.
    1213    
     1210
    12141211    To instantiate:
    1215      
     1212
    12161213       Polygon_function(polygons)
    1217        
     1214
    12181215    where polygons is a list of tuples of the form
    1219    
     1216
    12201217      [ (P0, v0), (P1, v1), ...]
    1221      
    1222       with Pi being lists of vertices defining polygons and vi either 
     1218
     1219      with Pi being lists of vertices defining polygons and vi either
    12231220      constants or functions of x,y to be applied to points with the polygon.
    1224      
    1225     The function takes an optional argument, default which is the value 
     1221
     1222    The function takes an optional argument, default which is the value
    12261223    (or function) to used for points not belonging to any polygon.
    12271224    For example:
    1228    
    1229        Polygon_function(polygons, default = 0.03)       
    1230      
    1231     If omitted the default value will be 0.0 
    1232    
     1225
     1226       Polygon_function(polygons, default = 0.03)
     1227
     1228    If omitted the default value will be 0.0
     1229
    12331230    Note: If two polygons overlap, the one last in the list takes precedence
    12341231
    12351232    """
    1236    
     1233
    12371234    def __init__(self, regions, default = 0.0):
    1238        
     1235
    12391236        try:
    12401237            len(regions)
    12411238        except:
    1242             msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
     1239            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
    12431240            raise msg
    12441241
    12451242
    1246         T = regions[0]                                   
     1243        T = regions[0]
    12471244        try:
    12481245            a = len(T)
    1249         except:   
    1250             msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
    1251             raise msg   
    1252                    
    1253         assert a == 2, 'Must have two component each: %s' %T       
    1254 
    1255         self.regions = regions             
     1246        except:
     1247            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
     1248            raise msg
     1249
     1250        assert a == 2, 'Must have two component each: %s' %T
     1251
     1252        self.regions = regions
    12561253        self.default = default
    1257          
    1258          
     1254
     1255
    12591256    def __call__(self, x, y):
    12601257        from util import inside_polygon
    12611258        from Numeric import ones, Float, concatenate, array, reshape, choose
    1262        
     1259
    12631260        x = array(x).astype(Float)
    1264         y = array(y).astype(Float)     
    1265        
     1261        y = array(y).astype(Float)
     1262
    12661263        N = len(x)
    12671264        assert len(y) == N
    1268        
    1269         points = concatenate( (reshape(x, (N, 1)), 
     1265
     1266        points = concatenate( (reshape(x, (N, 1)),
    12701267                               reshape(y, (N, 1))), axis=1 )
    12711268
    12721269        if callable(self.default):
    12731270            z = self.default(x,y)
    1274         else:                   
     1271        else:
    12751272            z = ones(N, Float) * self.default
    1276            
     1273
    12771274        for polygon, value in self.regions:
    12781275            indices = inside_polygon(points, polygon)
    1279            
    1280             #FIXME: This needs to be vectorised 
     1276
     1277            #FIXME: This needs to be vectorised
    12811278            if callable(value):
    12821279                for i in indices:
    … …  
    12841281                    yy = array([y[i]])
    12851282                    z[i] = value(xx, yy)[0]
    1286             else:   
     1283            else:
    12871284                for i in indices:
    1288                     z[i] = value 
    1289 
    1290         return z                 
     1285                    z[i] = value
     1286
     1287        return z
    12911288
    12921289def read_polygon(filename):
    12931290    """Read points assumed to form a polygon
    12941291       There must be exactly two numbers in each line
    1295     """             
    1296    
     1292    """
     1293
    12971294    #Get polygon
    12981295    fid = open(filename)
    … …  
    13041301        polygon.append( [float(fields[0]), float(fields[1])] )
    13051302
    1306     return polygon     
    1307    
     1303    return polygon
     1304
    13081305def populate_polygon(polygon, number_of_points, seed = None):
    13091306    """Populate given polygon with uniformly distributed points.
    … …  
    13111308    Input:
    13121309       polygon - list of vertices of polygon
    1313        number_of_points - (optional) number of points 
     1310       number_of_points - (optional) number of points
    13141311       seed - seed for random number generator (default=None)
    1315        
     1312
    13161313    Output:
    13171314       points - list of points inside polygon
    1318        
     1315
    13191316    Examples:
    13201317       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
    1321        will return five randomly selected points inside the unit square 
     1318       will return five randomly selected points inside the unit square
    13221319    """
    13231320
    … …  
    13321329    max_y = min_y = polygon[0][1]
    13331330    for point in polygon[1:]:
    1334         x = point[0] 
     1331        x = point[0]
    13351332        if x > max_x: max_x = x
    1336         if x < min_x: min_x = x           
    1337         y = point[1] 
     1333        if x < min_x: min_x = x
     1334        y = point[1]
    13381335        if y > max_y: max_y = y
    1339         if y < min_y: min_y = y           
    1340 
    1341 
    1342     while len(points) < number_of_points:     
     1336        if y < min_y: min_y = y
     1337
     1338
     1339    while len(points) < number_of_points:
    13431340        x = uniform(min_x, max_x)
    1344         y = uniform(min_y, max_y)       
     1341        y = uniform(min_y, max_y)
    13451342
    13461343        if inside_polygon( [x,y], polygon ):
    … …  
    13511348def tsh2sww(filename, verbose=True):
    13521349    """
    1353     to check if a tsh/msh file 'looks' good. 
     1350    to check if a tsh/msh file 'looks' good.
    13541351    """
    13551352    from shallow_water import Domain
    13561353    from pmesh2domain import pmesh_to_domain_instance
    1357     import time, os 
    1358     from data_manager import get_dataobject 
    1359     from os import sep, path 
     1354    import time, os
     1355    from data_manager import get_dataobject
     1356    from os import sep, path
    13601357    #from util import mean
    1361    
     1358
    13621359    if verbose == True:print 'Creating domain from', filename
    13631360    domain = pmesh_to_domain_instance(filename, Domain)
    … …  
    13761373    if verbose == True:
    13771374        print "Output written to " + domain.get_datadir() + sep + \
    1378               domain.filename + "." + domain.format 
     1375              domain.filename + "." + domain.format
    13791376    sww = get_dataobject(domain)
    13801377    sww.store_connectivity()
    … …  
    13881385    """
    13891386    """
    1390    
    1391     det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)           
     1387
     1388    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
    13921389    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
    13931390    a /= det
    13941391
    13951392    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
    1396     b /= det           
     1393    b /= det
    13971394
    13981395    return a, b
    … …  
    14141411    pass
    14151412
    1416    
    1417    
    1418    
    1419    
    1420 
    1421 #OBSOLETED STUFF   
     1413
     1414
     1415
     1416
     1417
     1418#OBSOLETED STUFF
    14221419def inside_polygon_old(point, polygon, closed = True, verbose = False):
    14231420    #FIXME Obsoleted
    14241421    """Determine whether points are inside or outside a polygon
    1425    
     1422
    14261423    Input:
    14271424       point - Tuple of (x, y) coordinates, or list of tuples
    14281425       polygon - list of vertices of polygon
    1429        closed - (optional) determine whether points on boundary should be 
    1430        regarded as belonging to the polygon (closed = True) 
    1431        or not (closed = False) 
    1432    
     1426       closed - (optional) determine whether points on boundary should be
     1427       regarded as belonging to the polygon (closed = True)
     1428       or not (closed = False)
     1429
    14331430    Output:
    14341431       If one point is considered, True or False is returned.
    1435        If multiple points are passed in, the function returns indices 
    1436        of those points that fall inside the polygon     
    1437        
     1432       If multiple points are passed in, the function returns indices
     1433       of those points that fall inside the polygon
     1434
    14381435    Examples:
    14391436       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
    14401437       inside_polygon( [0.5, 0.5], U)
    1441        will evaluate to True as the point 0.5, 0.5 is inside the unit square 
    1442        
     1438       will evaluate to True as the point 0.5, 0.5 is inside the unit square
     1439
    14431440       inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
    1444        will return the indices [0, 2] as only the first and the last point 
     1441       will return the indices [0, 2] as only the first and the last point
    14451442       is inside the unit square
    1446        
     1443
    14471444    Remarks:
    14481445       The vertices may be listed clockwise or counterclockwise and
    1449        the first point may optionally be repeated.   
     1446       the first point may optionally be repeated.
    14501447       Polygons do not need to be convex.
    1451        Polygons can have holes in them and points inside a hole is 
    1452        regarded as being outside the polygon.         
    1453                
    1454        
    1455     Algorithm is based on work by Darel Finley, 
    1456     http://www.alienryderflex.com/polygon/ 
    1457    
    1458     """   
     1448       Polygons can have holes in them and points inside a hole is
     1449       regarded as being outside the polygon.
     1450
     1451
     1452    Algorithm is based on work by Darel Finley,
     1453    http://www.alienryderflex.com/polygon/
     1454
     1455    """
    14591456
    14601457    from Numeric import array, Float, reshape
    1461    
    1462    
     1458
     1459
    14631460    #Input checks
    14641461    try:
    1465         point = array(point).astype(Float)         
     1462        point = array(point).astype(Float)
    14661463    except:
    14671464        msg = 'Point %s could not be converted to Numeric array' %point
    1468         raise msg       
    1469        
     1465        raise msg
     1466
    14701467    try:
    1471         polygon = array(polygon).astype(Float)             
     1468        polygon = array(polygon).astype(Float)
    14721469    except:
    14731470        msg = 'Polygon %s could not be converted to Numeric array' %polygon
    1474         raise msg               
    1475    
     1471        raise msg
     1472
    14761473    assert len(polygon.shape) == 2,\
    14771474       'Polygon array must be a 2d array of vertices: %s' %polygon
    14781475
    1479        
     1476
    14801477    assert polygon.shape[1] == 2,\
    1481        'Polygon array must have two columns: %s' %polygon   
     1478       'Polygon array must have two columns: %s' %polygon
    14821479
    14831480
    … …  
    14851482        one_point = True
    14861483        points = reshape(point, (1,2))
    1487     else:       
     1484    else:
    14881485        one_point = False
    14891486        points = point
    14901487
    1491     N = polygon.shape[0] #Number of vertices in polygon 
     1488    N = polygon.shape[0] #Number of vertices in polygon
    14921489    M = points.shape[0]  #Number of points
    1493    
     1490
    14941491    px = polygon[:,0]
    1495     py = polygon[:,1]   
     1492    py = polygon[:,1]
    14961493
    14971494    #Used for an optimisation when points are far away from polygon
    14981495    minpx = min(px); maxpx = max(px)
    1499     minpy = min(py); maxpy = max(py)   
     1496    minpy = min(py); maxpy = max(py)
    15001497
    15011498
    … …  
    15061503        if verbose:
    15071504            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
    1508        
    1509         #for each point   
     1505
     1506        #for each point
    15101507        x = points[k, 0]
    1511         y = points[k, 1]       
     1508        y = points[k, 1]
    15121509
    15131510        inside = False
    … …  
    15151512        #Optimisation
    15161513        if x > maxpx or x < minpx: continue
    1517         if y > maxpy or y < minpy: continue       
    1518 
    1519         #Check polygon 
     1514        if y > maxpy or y < minpy: continue
     1515
     1516        #Check polygon
    15201517        for i in range(N):
    15211518            j = (i+1)%N
    1522            
    1523             #Check for case where point is contained in line segment   
     1519
     1520            #Check for case where point is contained in line segment
    15241521            if point_on_line(x, y, px[i], py[i], px[j], py[j]):
    15251522                if closed:
    15261523                    inside = True
    1527                 else:       
     1524                else:
    15281525                    inside = False
    15291526                break
    1530             else:       
    1531                 #Check if truly inside polygon             
     1527            else:
     1528                #Check if truly inside polygon
    15321529                if py[i] < y and py[j] >= y or\
    15331530                   py[j] < y and py[i] >= y:
    15341531                    if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
    15351532                        inside = not inside
    1536                    
     1533
    15371534        if inside: indices.append(k)
    15381535
    1539     if one_point: 
     1536    if one_point:
    15401537        return inside
    15411538    else:
    15421539        return indices
    1543              
     1540
    15441541
    15451542#def remove_from(A, B):
    … …  
    15501547#    ind = search_sorted(A, B)
    15511548
    1552    
     1549
    15531550
    15541551def outside_polygon_old(point, polygon, closed = True, verbose = False):
    15551552    #OBSOLETED
    15561553    """Determine whether points are outside a polygon
    1557    
     1554
    15581555    Input:
    15591556       point - Tuple of (x, y) coordinates, or list of tuples
    15601557       polygon - list of vertices of polygon
    1561        closed - (optional) determine whether points on boundary should be 
    1562        regarded as belonging to the polygon (closed = True) 
    1563        or not (closed = False) 
    1564    
     1558       closed - (optional) determine whether points on boundary should be
     1559       regarded as belonging to the polygon (closed = True)
     1560       or not (closed = False)
     1561
    15651562    Output:
    15661563       If one point is considered, True or False is returned.
    1567        If multiple points are passed in, the function returns indices 
    1568        of those points that fall outside the polygon     
    1569        
     1564       If multiple points are passed in, the function returns indices
     1565       of those points that fall outside the polygon
     1566
    15701567    Examples:
    1571        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
     1568       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
    15721569       outside_polygon( [0.5, 0.5], U )
    15731570       will evaluate to False as the point 0.5, 0.5 is inside the unit square
    … …  
    15751572       ouside_polygon( [1.5, 0.5], U )
    15761573       will evaluate to True as the point 1.5, 0.5 is outside the unit square
    1577        
     1574
    15781575       outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
    1579        will return the indices [1] as only the first and the last point 
     1576       will return the indices [1] as only the first and the last point
    15801577       is inside the unit square
    15811578    """
    15821579
    15831580    #FIXME: This is too slow
    1584    
     1581
    15851582    res  = inside_polygon(point, polygon, closed, verbose)
    15861583
    … …  
    16001597
    16011598    C-wrapper
    1602     """   
     1599    """
    16031600
    16041601    from Numeric import array, Float, reshape, zeros, Int
    1605    
    1606 
    1607     if verbose: print 'Checking input to inside_polygon'   
     1602
     1603
     1604    if verbose: print 'Checking input to inside_polygon'
    16081605    #Input checks
    16091606    try:
    1610         point = array(point).astype(Float)         
     1607        point = array(point).astype(Float)
    16111608    except:
    16121609        msg = 'Point %s could not be converted to Numeric array' %point
    1613         raise msg       
    1614        
     1610        raise msg
     1611
    16151612    try:
    1616         polygon = array(polygon).astype(Float)             
     1613        polygon = array(polygon).astype(Float)
    16171614    except:
    16181615        msg = 'Polygon %s could not be converted to Numeric array' %polygon
    1619         raise msg               
    1620    
     1616        raise msg
     1617
    16211618    assert len(polygon.shape) == 2,\
    16221619       'Polygon array must be a 2d array of vertices: %s' %polygon
    16231620
    1624        
     1621
    16251622    assert polygon.shape[1] == 2,\
    1626        'Polygon array must have two columns: %s' %polygon   
     1623       'Polygon array must have two columns: %s' %polygon
    16271624
    16281625
    … …  
    16301627        one_point = True
    16311628        points = reshape(point, (1,2))
    1632     else:       
     1629    else:
    16331630        one_point = False
    16341631        points = point
    … …  
    16401637    indices = zeros( points.shape[0], Int )
    16411638
    1642     if verbose: print 'Calling C-version of inside poly' 
     1639    if verbose: print 'Calling C-version of inside poly'
    16431640    count = inside_polygon(points, polygon, indices,
    16441641                           int(closed), int(verbose))
    … …  
    16481645    else:
    16491646        if verbose: print 'Got %d points' %count
    1650        
     1647
    16511648        return indices[:count]   #Return those indices that were inside
    1652 

  • inundation/ga/storm_surge/zeus/anuga-workspace.zwi

    r1271 r1280  
    22<workspace name="anuga-workspace">
    33    <mode>Debug</mode>
    4     <active>steve</active>
     4    <active>parallel</active>
     5    <project name="parallel">parallel.zpi</project>
    56    <project name="pyvolution">pyvolution.zpi</project>
    67    <project name="steve">steve.zpi</project>
    7     <project name="pyvolution-parallel">pyvolution-parallel.zpi</project>
    88</workspace>
Note: See TracChangeset for help on using the changeset viewer.