Changeset 1671


Ignore:
Timestamp:
Aug 2, 2005, 4:48:56 PM (20 years ago)
Author:
ole
Message:

Moved old file function stuff out and re-tested.
File_function is now all in one function which is based on NetCDF

Location:
inundation/ga/storm_surge/pyvolution
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r1669 r1671  
    734734
    735735
    736 
     736#########################################################
    737737#Conversion routines
     738########################################################
     739
    738740def sww2obj(basefilename, size):
    739741    """Convert netcdf based data output to obj
    … …  
    19431945    outfile.close()
    19441946
     1947
     1948
     1949
     1950def timefile2swww(filename, quantity_names = None):
     1951    """Template for converting typical text files with time series to
     1952    NetCDF sww file.
     1953
     1954
     1955    The file format is assumed to be either two fields separated by a comma:
     1956
     1957        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
     1958
     1959    E.g
     1960
     1961      31/08/04 00:00:00, 1.328223 0 0
     1962      31/08/04 00:15:00, 1.292912 0 0
     1963
     1964    will provide a time dependent function f(t) with three attributes
     1965
     1966    filename is assumed to be the rootname with extenisons .txt and .sww
     1967    """
     1968
     1969    import time, calendar
     1970    from Numeric import array
     1971    from config import time_format
     1972    from util import ensure_numeric
     1973
     1974
     1975    fid = open(filename + '.txt')
     1976    line = fid.readline()
     1977    fid.close()
     1978
     1979    fields = line.split(',')
     1980    msg = 'File %s must have the format date, value0 value1 value2 ...'
     1981    assert len(fields) == 2, msg
     1982
     1983    try:
     1984        starttime = calendar.timegm(time.strptime(fields[0], time_format))
     1985    except ValueError:
     1986        msg = 'First field in file %s must be' %filename
     1987        msg += ' date-time with format %s.\n' %time_format
     1988        msg += 'I got %s instead.' %fields[0]
     1989        raise msg
     1990
     1991
     1992    #Split values
     1993    values = []
     1994    for value in fields[1].split():
     1995        values.append(float(value))
     1996
     1997    q = ensure_numeric(values)
     1998
     1999    msg = 'ERROR: File must contain at least one independent value'
     2000    assert len(q.shape) == 1, msg
     2001
     2002
     2003
     2004    #Read times proper
     2005    from Numeric import zeros, Float, alltrue
     2006    from config import time_format
     2007    import time, calendar
     2008
     2009    fid = open(filename + '.txt')
     2010    lines = fid.readlines()
     2011    fid.close()
     2012
     2013    N = len(lines)
     2014    d = len(q)
     2015
     2016    T = zeros(N, Float)       #Time
     2017    Q = zeros((N, d), Float)  #Values
     2018
     2019    for i, line in enumerate(lines):
     2020        fields = line.split(',')
     2021        realtime = calendar.timegm(time.strptime(fields[0], time_format))
     2022       
     2023        T[i] = realtime - starttime
     2024
     2025        for j, value in enumerate(fields[1].split()):
     2026            Q[i, j] = float(value)
     2027
     2028    msg = 'File %s must list time as a monotonuosly ' %filename
     2029    msg += 'increasing sequence'
     2030    assert alltrue( T[1:] - T[:-1] > 0 ), msg
     2031
     2032
     2033    #Create sww file
     2034    from Scientific.IO.NetCDF import NetCDFFile
     2035
     2036       
     2037
     2038    fid = NetCDFFile(filename + '.sww', 'w')
     2039
     2040
     2041    fid.institution = 'Geoscience Australia'
     2042    fid.description = 'Time series'
     2043
     2044
     2045    #Reference point
     2046    #Start time in seconds since the epoch (midnight 1/1/1970)
     2047    #FIXME: Use Georef
     2048    fid.starttime = starttime
     2049   
     2050   
     2051    # dimension definitions
     2052    #fid.createDimension('number_of_volumes', self.number_of_volumes)
     2053    #fid.createDimension('number_of_vertices', 3)
     2054
     2055
     2056    fid.createDimension('number_of_timesteps', len(T))
     2057
     2058    #if geo_reference is not None:
     2059    #    geo_reference.write_NetCDF(fid)
     2060               
     2061
     2062    fid.createVariable('time', Float, ('number_of_timesteps',))
     2063
     2064    fid.variables['time'][:] = T
     2065   
     2066    for i in range(Q.shape[1]):
     2067        try:
     2068            name = quantity_names[i]
     2069        except:   
     2070            name = 'Attribute%d'%i
     2071           
     2072        fid.createVariable(name, Float, ('number_of_timesteps',))
     2073        fid.variables[name][:] = Q[:,i]       
     2074
     2075    fid.close()
    19452076
    19462077

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

    r1670 r1671  
    11561156                q[i] = Q0
    11571157
    1158             #if ratio > 0:
    1159             #    q[i] = Q[self.index, point_id] +\
    1160             #           ratio*(Q[self.index+1, point_id] -\
    1161             #                  Q[self.index, point_id])
    1162             #else:
    1163             #    q[i] = Q[self.index, point_id]
    11641158
    11651159        #Return vector of interpolated values
    1166         if len(q) == 1:
    1167             return q[0]
     1160        #if len(q) == 1:
     1161        #    return q[0]
     1162        #else:
     1163        #    return q
     1164
     1165
     1166        #Return vector of interpolated values
     1167        #FIXME:
     1168        if self.spatial is True:
     1169            return q
    11681170        else:
    1169             return q
    1170 
     1171            #Replicate q according to x and y
     1172            #This is e.g used for Wind_stress
     1173            if x == None or y == None:
     1174                return q
     1175            else:
     1176                try:
     1177                    N = len(x)
     1178                except:
     1179                    return q
     1180                else:
     1181                    from Numeric import ones, Float
     1182                    #x is a vector - Create one constant column for each value
     1183                    N = len(x)
     1184                    assert len(y) == N, 'x and y must have same length'
     1185                    res = []
     1186                    for col in q:
     1187                        res.append(col*ones(N, Float))
     1188                       
     1189                return res
    11711190
    11721191

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

    r1657 r1671  
    11411141        y = ones(3, Float)
    11421142        try:
    1143             q = f(1.0, x, y)
     1143            q = f(1.0, x=x, y=y)
    11441144        except Exception, e:
    11451145            msg = 'Function %s could not be executed:\n%s' %(f, e)
    … …  
    11551155            raise msg
    11561156
    1157         msg = 'Return vector from function %s' %f
     1157        #Is this really what we want?
     1158        msg = 'Return vector from function %s ' %f
    11581159        msg += 'must have same lenght as input vectors'
    11591160        assert len(q) == N, msg
    … …  
    12231224            #Assume vector function returning (s, phi)(t,x,y)
    12241225            vector_function = args[0]
    1225             s = lambda t,x,y: vector_function(t,x,y)[0]
    1226             phi = lambda t,x,y: vector_function(t,x,y)[1]
     1226            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
     1227            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
    12271228        else:
    12281229           #Assume info is in 2 keyword arguments

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

    r1018 r1671  
    138138        #Write file
    139139        filename = 'boundarytest' + str(time.time())
    140         fid = open(filename, 'w')
     140        fid = open(filename + '.txt', 'w')
    141141        start = time.mktime(time.strptime('2000', '%Y'))
    142142        dt = 5*60  #Five minute intervals
    … …  
    149149
    150150
    151         F = File_boundary(filename, domain)
    152         os.remove(filename)
     151        #Convert ASCII file to NetCDF (Which is what we really like!)
     152       
     153        from data_manager import timefile2swww               
     154        timefile2swww(filename, quantity_names = ['stage', 'ymomentum'])
     155       
     156
     157
     158        F = File_boundary(filename + '.sww', domain)
     159
     160       
     161        os.remove(filename + '.txt')
     162        os.remove(filename + '.sww')       
    153163
    154164
    … …  
    182192
    183193    def test_fileboundary_exception(self):
    184         """Test that boundary object complians if number of
     194        """Test that boundary object complains if number of
    185195        conserved quantities are wrong
    186196        """
    … …  
    221231        #Write file (with only two values)
    222232        filename = 'boundarytest' + str(time.time())
    223         fid = open(filename, 'w')
     233        fid = open(filename + '.txt', 'w')
    224234        start = time.mktime(time.strptime('2000', '%Y'))
    225235        dt = 5*60  #Five minute intervals
    … …  
    232242
    233243
    234         try:
    235             F = File_boundary(filename, domain)
    236         except AssertionError:
    237             pass
    238         else:
    239             raise 'Should have raised an exception'
    240 
    241         os.remove(filename)
     244        #Convert ASCII file to NetCDF (Which is what we really like!)
     245        from data_manager import timefile2swww               
     246        timefile2swww(filename, quantity_names = ['stage', 'xmomentum'])
     247
     248       
     249        try:
     250            F = File_boundary(filename + '.sww', domain)           
     251        except:
     252            pass
     253        else:
     254            raise 'Should have raised an exception'       
     255       
     256        os.remove(filename + '.txt')
     257        os.remove(filename + '.sww')       
     258
     259
    242260
    243261

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

    r1636 r1671  
    994994        #Write wind stress file (ensure that domain.time is covered)
    995995        #Take x=1 and y=0
    996         filename = 'test_windstress_from_file.txt'
     996        filename = 'test_windstress_from_file'
    997997        start = time.mktime(time.strptime('2000', '%Y'))
    998         fid = open(filename, 'w')
     998        fid = open(filename + '.txt', 'w')
    999999        dt = 1  #One second interval
    10001000        t = 0.0
    … …  
    10091009        fid.close()
    10101010
     1011
     1012        #Convert ASCII file to NetCDF (Which is what we really like!)
     1013        from data_manager import timefile2swww       
     1014        timefile2swww(filename)
     1015
     1016
     1017       
    10111018        #Setup wind stress
    1012         F = file_function(filename)
     1019        F = file_function(filename + '.sww', quantities = ['Attribute0',
     1020                                                           'Attribute1'])
     1021
     1022        #print 'F(5)', F(5)
     1023
     1024        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
     1025       
     1026        #print dir(F)
     1027        #print F.T
     1028        #print F.precomputed_values
     1029        #
     1030        #F = file_function(filename + '.txt')       
     1031        #
     1032        #print dir(F)
     1033        #print F.T       
     1034        #print F.Q
     1035       
    10131036        W = Wind_stress(F)
     1037
    10141038        domain.forcing_terms = []
    10151039        domain.forcing_terms.append(W)
    … …  
    10431067            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
    10441068
    1045         os.remove(filename)
     1069        #os.remove(filename + '.txt')
    10461070
    10471071    def test_wind_stress_error_condition(self):

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

    r1670 r1671  
    88from util import *
    99from config import epsilon
     10from data_manager import timefile2swww
    1011
    1112
    … …  
    171172        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
    172173
    173 
    174 
    175 
    176     def test_file_function_time(self):
     174    def test_ensure_numeric(self):
     175        from util import ensure_numeric
     176        from Numeric import ArrayType, Float, array
     177
     178        A = [1,2,3,4]
     179        B = ensure_numeric(A)
     180        assert type(B) == ArrayType
     181        assert B.typecode() == 'l'
     182        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
     183
     184
     185        A = [1,2,3.14,4]
     186        B = ensure_numeric(A)
     187        assert type(B) == ArrayType
     188        assert B.typecode() == 'd'
     189        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
     190
     191
     192        A = [1,2,3,4]
     193        B = ensure_numeric(A, Float)
     194        assert type(B) == ArrayType
     195        assert B.typecode() == 'd'
     196        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
     197
     198
     199        A = [1,2,3,4]
     200        B = ensure_numeric(A, Float)
     201        assert type(B) == ArrayType
     202        assert B.typecode() == 'd'
     203        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
     204
     205
     206        A = array([1,2,3,4])
     207        B = ensure_numeric(A)
     208        assert type(B) == ArrayType
     209        assert B.typecode() == 'l'       
     210        assert A == B   
     211        assert A is B   #Same object
     212
     213
     214        A = array([1,2,3,4])
     215        B = ensure_numeric(A, Float)
     216        assert type(B) == ArrayType
     217        assert B.typecode() == 'd'       
     218        assert A == B   
     219        assert A is not B   #Not the same object       
     220
     221
     222
     223    def test_file_function_time1(self):
    177224        """Test that File function interpolates correctly
    178225        between given times. No x,y dependency here.
    … …  
    184231        from math import sin, pi
    185232
     233        #Typical ASCII file
    186234        finaltime = 1200
    187         filename = 'test_file_function.txt'
    188         fid = open(filename, 'w')
     235        filename = 'test_file_function'
     236        fid = open(filename + '.txt', 'w')
    189237        start = time.mktime(time.strptime('2000', '%Y'))
    190238        dt = 60  #One minute intervals
    … …  
    197245        fid.close()
    198246
    199         F = file_function(filename)
    200 
     247        #Convert ASCII file to NetCDF (Which is what we really like!)
     248        timefile2swww(filename)
     249
     250
     251        #Create file function
     252        F = file_function(filename + '.sww', quantities = ['Attribute0',
     253                                                           'Attribute1',
     254                                                           'Attribute2'])
     255       
    201256        #Now try interpolation
    202257        for i in range(20):
    … …  
    223278        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
    224279
    225         os.remove(filename)
    226 
    227 
    228     def test_ensure_numeric(self):
    229         from util import ensure_numeric
    230         from Numeric import ArrayType, Float, array
    231 
    232         A = [1,2,3,4]
    233         B = ensure_numeric(A)
    234         assert type(B) == ArrayType
    235         assert B.typecode() == 'l'
    236         assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
    237 
    238 
    239         A = [1,2,3.14,4]
    240         B = ensure_numeric(A)
    241         assert type(B) == ArrayType
    242         assert B.typecode() == 'd'
    243         assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
    244 
    245 
    246         A = [1,2,3,4]
    247         B = ensure_numeric(A, Float)
    248         assert type(B) == ArrayType
    249         assert B.typecode() == 'd'
    250         assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
    251 
    252 
    253         A = [1,2,3,4]
    254         B = ensure_numeric(A, Float)
    255         assert type(B) == ArrayType
    256         assert B.typecode() == 'd'
    257         assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
    258 
    259 
    260         A = array([1,2,3,4])
    261         B = ensure_numeric(A)
    262         assert type(B) == ArrayType
    263         assert B.typecode() == 'l'       
    264         assert A == B   
    265         assert A is B   #Same object
    266 
    267 
    268         A = array([1,2,3,4])
    269         B = ensure_numeric(A, Float)
    270         assert type(B) == ArrayType
    271         assert B.typecode() == 'd'       
    272         assert A == B   
    273         assert A is not B   #Not the same object       
     280        os.remove(filename + '.txt')
     281        os.remove(filename + '.sww')       
    274282
    275283
    … …  
    649657
    650658        finaltime = 1200
    651         filename = 'test_file_function.txt'
    652         fid = open(filename, 'w')
     659        filename = 'test_file_function'
     660        fid = open(filename + '.txt', 'w')
    653661        start = time.mktime(time.strptime('2000', '%Y'))
    654662        dt = 60  #One minute intervals
    … …  
    661669        fid.close()
    662670
     671
     672        #Convert ASCII file to NetCDF (Which is what we really like!)
     673        timefile2swww(filename)
     674
     675
     676
    663677        a = [0.0, 0.0]
    664678        b = [4.0, 0.0]
    … …  
    670684
    671685        #Check that domain.starttime is updated if non-existing
    672         F = file_function(filename, domain)
     686        F = file_function(filename + '.sww', domain)
     687       
    673688        assert allclose(domain.starttime, start)
    674689
    675690        #Check that domain.starttime is updated if too early
    676691        domain.starttime = start - 1
    677         F = file_function(filename, domain)
     692        F = file_function(filename + '.sww', domain)
    678693        assert allclose(domain.starttime, start)
    679694
    680695        #Check that domain.starttime isn't updated if later
    681696        domain.starttime = start + 1
    682         F = file_function(filename, domain)
     697        F = file_function(filename + '.sww', domain)       
    683698        assert allclose(domain.starttime, start+1)
    684699
    685700        domain.starttime = start
    686 
    687 
     701        F = file_function(filename + '.sww', domain,
     702                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
     703       
     704
     705        #print F.T
     706        #print F.precomputed_values
     707        #print 'F(60)', F(60)
     708       
    688709        #Now try interpolation
    689710        for i in range(20):
    … …  
    710731        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
    711732
    712         os.remove(filename)
    713 
     733        os.remove(filename + '.sww')
     734        os.remove(filename + '.txt')       
    714735
    715736    def test_file_function_time_with_domain_different_start(self):
    … …  
    728749
    729750        finaltime = 1200
    730         filename = 'test_file_function.txt'
    731         fid = open(filename, 'w')
     751        filename = 'test_file_function'
     752        fid = open(filename + '.txt', 'w')
    732753        start = time.mktime(time.strptime('2000', '%Y'))
    733754        dt = 60  #One minute intervals
    … …  
    740761        fid.close()
    741762
     763        #Convert ASCII file to NetCDF (Which is what we really like!)
     764        timefile2swww(filename)       
     765
    742766        a = [0.0, 0.0]
    743767        b = [4.0, 0.0]
    … …  
    752776        delta = 23
    753777        domain.starttime = start + delta
    754         F = file_function(filename, domain)
     778        #F = file_function(filename + '.sww', domain)
     779        F = file_function(filename + '.sww', domain,
     780                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])       
    755781        assert allclose(domain.starttime, start+delta)
    756782
    … …  
    782808        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
    783809
    784         os.remove(filename)
     810
     811        os.remove(filename + '.sww')
     812        os.remove(filename + '.txt')               
    785813
    786814
    … …  
    10881116#-------------------------------------------------------------
    10891117if __name__ == "__main__":
    1090     #suite = unittest.makeSuite(Test_Util,'test_spatio_temporal_file_function_time')
    10911118    suite = unittest.makeSuite(Test_Util,'test')
     1119    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
    10921120    runner = unittest.TextTestRunner()
    10931121    runner.run(suite)

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

    r1670 r1671  
    66
    77#FIXME: Move polygon stuff out
    8 #FIXME: Finish Interpolation function and get rid of File_function_ASCII
    9 
    108
    119def angle(v):
    … …  
    135133
    136134
    137 def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False):
    138     """If domain is specified, don't specify quantites as they are automatically derived
     135def file_function(filename,
     136                  domain = None,
     137                  quantities = None,
     138                  interpolation_points = None, verbose = False):
     139    """If quantities is not specified, derive them from domain
     140    (if that is specified)
    139141    """
    140142
    … …  
    157159    fid.close()
    158160
    159     if domain is not None:
    160         quantities = domain.conserved_quantities
    161     else:
    162         quantities = None
     161    if quantities is None:
     162        if domain is not None:
     163            quantities = domain.conserved_quantities
     164
    163165
    164166
    … …  
    167169                                        interpolation_points, verbose = verbose)
    168170    else:
    169         #FIXME Should be phased out
    170         return File_function_ASCII(filename, domain, quantities,
    171                                    interpolation_points)
    172    
    173    
    174 
    175 
    176 
    177 
    178 
     171        raise 'Must be a NetCDF File'
    179172
    180173
    … …  
    215208
    216209    if quantity_names is None or len(quantity_names) < 1:
    217         msg = 'ERROR: At least one independent value must be specified'
    218         raise msg
     210        #Get quantities from file
     211        quantity_names = []
     212        for name in fid.variables.keys():
     213            if name not in ['x', 'y', 'time']:
     214                quantity_names.append(name)
     215               
     216        #msg = 'ERROR: At least one independent value must be specified'
     217        #raise msg
    219218
    220219    missing = []
    221     for quantity in ['x', 'y', 'time'] + quantity_names:
     220       
     221    for quantity in ['time'] + quantity_names:
    222222        if not fid.variables.has_key(quantity):
    223223            missing.append(quantity)
    … …  
    226226        msg = 'Quantities %s could not be found in file %s'\
    227227              %(str(missing), filename)
     228        fid.close()
    228229        raise msg
     230
     231    #Decide whether this data has a spatial dimension
     232    spatial = True
     233    for quantity in ['x', 'y']:
     234        if not fid.variables.has_key(quantity):
     235            spatial = False
     236
    229237
    230238
    … …  
    236244        raise msg
    237245
    238     #Get origin
    239     xllcorner = fid.xllcorner[0]
    240     yllcorner = fid.yllcorner[0]
    241 
    242 
    243246    #Get variables
    244     if verbose: print 'Get variables'
    245     x = fid.variables['x'][:]
    246     y = fid.variables['y'][:]
    247     triangles = fid.variables['volumes'][:]
    248     time = fid.variables['time'][:]
    249 
    250     x = reshape(x, (len(x),1))
    251     y = reshape(y, (len(y),1))
    252     vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
     247    if verbose: print 'Get variables'   
     248    time = fid.variables['time'][:]   
     249
     250    if spatial:
     251        #Get origin
     252        xllcorner = fid.xllcorner[0]
     253        yllcorner = fid.yllcorner[0]
     254
     255        x = fid.variables['x'][:]
     256        y = fid.variables['y'][:]
     257        triangles = fid.variables['volumes'][:]
     258
     259        x = reshape(x, (len(x),1))
     260        y = reshape(y, (len(y),1))
     261        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
     262
    253263
    254264    if domain is not None:
    … …  
    277287        print '  References:'
    278288        #print '    Datum ....' #FIXME
    279         print '    Lower left corner: [%f, %f]'\
    280               %(xllcorner, yllcorner)
     289        if spatial:
     290            print '    Lower left corner: [%f, %f]'\
     291                  %(xllcorner, yllcorner)
    281292        print '    Start time:   %f' %starttime               
    282293       
    … …  
    292303
    293304    from least_squares import Interpolation_function
     305
     306    if not spatial:
     307        vertex_coordinates = triangles = interpolation_points = None         
     308       
    294309    return Interpolation_function(time,
    295310                                  quantities,
    … …  
    299314                                  interpolation_points,
    300315                                  verbose = verbose)
    301 
    302 
    303 
    304 
    305 
    306 
    307 
    308 
    309 
    310 
    311 
    312 
    313 
    314 
    315 #FIXME Obsolete
    316 class File_function_NetCDF:
    317     """Read time history of spatial data from NetCDF sww file and
    318     return a callable object f(t,x,y)
    319     which will return interpolated values based on the input file.
    320 
    321     x, y may be either scalars or vectors
    322 
    323     #FIXME: More about format, interpolation and order of quantities
    324 
    325     The quantities returned by the callable objects are specified by
    326     the list quantities which must contain the names of the
    327     quantities to be returned and also reflect the order, e.g. for
    328     the shallow water wave equation, on would have
    329     quantities = ['stage', 'xmomentum', 'ymomentum']
    330 
    331     interpolation_points decides at which points interpolated
    332     quantities are to be computed whenever object is called.
    333     If None, return average value
    334     """
    335 
    336     def  __init__(self, filename, domain=None, quantities=None,
    337                   interpolation_points=None, verbose = False):
    338         """Initialise callable object from a file.
    339 
    340         If domain is specified, model time (domain.starttime)
    341         will be checked and possibly modified
    342 
    343         All times are assumed to be in UTC
    344         """
    345 
    346         #FIXME: Check that model origin is the same as file's origin
    347         #(both in UTM coordinates)
    348         #If not - modify those from file to match domain
    349         #Take this code from e.g. dem2pts in data_manager.py
    350         #FIXME: Use geo_reference to read and write xllcorner...
    351        
    352 
    353         import time, calendar
    354         from config import time_format
    355         from Scientific.IO.NetCDF import NetCDFFile
    356 
    357         #Open NetCDF file
    358         if verbose: print 'Reading', filename
    359         fid = NetCDFFile(filename, 'r')
    360 
    361         if quantities is None or len(quantities) < 1:
    362             msg = 'ERROR: At least one independent value must be specified'
    363             raise msg
    364 
    365         missing = []
    366         for quantity in ['x', 'y', 'time'] + quantities:
    367              if not fid.variables.has_key(quantity):
    368                  missing.append(quantity)
    369 
    370         if len(missing) > 0:
    371             msg = 'Quantities %s could not be found in file %s'\
    372                   %(str(missing), filename)
    373             raise msg
    374 
    375 
    376         #Get first timestep
    377         try:
    378             self.starttime = fid.starttime[0]
    379         except ValueError:
    380             msg = 'Could not read starttime from file %s' %filename
    381             raise msg
    382 
    383         #Get origin
    384         self.xllcorner = fid.xllcorner[0]
    385         self.yllcorner = fid.yllcorner[0]
    386 
    387         self.number_of_values = len(quantities)
    388         self.fid = fid
    389         self.filename = filename
    390         self.quantities = quantities
    391         self.domain = domain
    392         self.interpolation_points = interpolation_points
    393 
    394         if domain is not None:
    395             msg = 'WARNING: Start time as specified in domain (%f)'\
    396                   %domain.starttime
    397             msg += ' is earlier than the starttime of file %s (%f).'\
    398                      %(self.filename, self.starttime)
    399             msg += ' Modifying domain starttime accordingly.'
    400             if self.starttime > domain.starttime:
    401                 if verbose: print msg
    402                 domain.starttime = self.starttime #Modifying model time
    403                 if verbose: print 'Domain starttime is now set to %f'\
    404                    %domain.starttime
    405 
    406         #Read all data in and produce values for desired data points at
    407         #each timestep
    408         self.spatial_interpolation(interpolation_points, verbose = verbose)
    409         fid.close()
    410 
    411     def spatial_interpolation(self, interpolation_points, verbose = False):
    412         """For each timestep in fid: read surface,
    413         interpolate to desired points and store values for use when
    414         object is called.
    415         """
    416 
    417         from Numeric import array, zeros, Float, alltrue, concatenate, reshape
    418 
    419         from config import time_format
    420         from least_squares import Interpolation
    421         import time, calendar
    422 
    423         fid = self.fid
    424 
    425         #Get variables
    426         if verbose: print 'Get variables'
    427         x = fid.variables['x'][:]
    428         y = fid.variables['y'][:]
    429         #z = fid.variables['elevation'][:]
    430         triangles = fid.variables['volumes'][:]
    431         time = fid.variables['time'][:]
    432 
    433         #Check
    434         msg = 'File %s must list time as a monotonuosly ' %self.filename
    435         msg += 'increasing sequence'
    436         assert alltrue(time[1:] - time[:-1] > 0 ), msg
    437 
    438 
    439         if interpolation_points is not None:
    440 
    441             try:
    442                 P = ensure_numeric(interpolation_points)
    443             except:
    444                 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
    445                 msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
    446                 raise msg
    447 
    448 
    449             self.T = time[:]  #Time assumed to be relative to starttime
    450             self.index = 0    #Initial time index
    451 
    452             self.values = {}
    453             for name in self.quantities:
    454                 self.values[name] = zeros( (len(self.T),
    455                                             len(interpolation_points)),
    456                                             Float)
    457 
    458             #Build interpolator
    459             if verbose: print 'Build interpolation matrix'
    460             x = reshape(x, (len(x),1))
    461             y = reshape(y, (len(y),1))
    462             vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
    463 
    464             interpol = Interpolation(vertex_coordinates,
    465                                      triangles,
    466                                      point_coordinates = P,
    467                                      alpha = 0,
    468                                      verbose = verbose)
    469 
    470 
    471             if verbose: print 'Interpolate'
    472             for i, t in enumerate(self.T):
    473                 #Interpolate quantities at this timestep
    474                 if verbose: print ' time step %d of %d' %(i, len(self.T))
    475                 for name in self.quantities:
    476                     self.values[name][i, :] =\
    477                     interpol.interpolate(fid.variables[name][i,:])
    478 
    479             #Report
    480             if verbose:
    481                 print '------------------------------------------------'
    482                 print 'File_function statistics:'
    483                 print '  Name: %s' %self.filename
    484                 print '  References:'
    485                 #print '    Datum ....' #FIXME
    486                 print '    Lower left corner: [%f, %f]'\
    487                       %(self.xllcorner, self.yllcorner)
    488                 print '    Start time:   %f' %self.starttime               
    489                 #print '    Start time (file):   %f' %self.starttime
    490                 #print '    Start time (domain): %f' %domain.starttime
    491                 print '  Extent:'
    492                 print '    x in [%f, %f], len(x) == %d'\
    493                       %(min(x.flat), max(x.flat), len(x.flat))
    494                 print '    y in [%f, %f], len(y) == %d'\
    495                       %(min(y.flat), max(y.flat), len(y.flat))
    496                 print '    t in [%f, %f], len(t) == %d'\
    497                       %(min(self.T), max(self.T), len(self.T))
    498                 print '  Quantities:'
    499                 for name in self.quantities:
    500                     q = fid.variables[name][:].flat
    501                     print '    %s in [%f, %f]' %(name, min(q), max(q))
    502 
    503 
    504                 print '  Interpolation points (xi, eta):'\
    505                       'number of points == %d ' %P.shape[0]
    506                 print '    xi in [%f, %f]' %(min(P[:,0]), max(P[:,0]))
    507                 print '    eta in [%f, %f]' %(min(P[:,1]), max(P[:,1]))
    508                 print '  Interpolated quantities (over all timesteps):'
    509                 for name in self.quantities:
    510                     q = self.values[name][:].flat
    511                     print '    %s at interpolation points in [%f, %f]'\
    512                           %(name, min(q), max(q))
    513                 print '------------------------------------------------'
    514         else:
    515             msg = 'File_function_NetCDF must be invoked with '
    516             msg += 'a list of interpolation points'
    517             raise msg
    518             #FIXME: Should not be an error.
    519             #__call__ could return an average or in case of time series only
    520             #no point s are needed
    521 
    522 
    523     def __repr__(self):
    524         return 'File function'
    525 
    526     def __call__(self, t, x=None, y=None, point_id = None):
    527         """Evaluate f(t, point_id)
    528 
    529         Inputs:
    530           t: time - Model time (tau = domain.starttime-self.starttime+t)
    531                     must lie within existing timesteps
    532           point_id: index of one of the preprocessed points.
    533                     If point_id is None all preprocessed points are computed
    534 
    535           FIXME: point_id could also be a slice
    536           FIXME: One could allow arbitrary x, y coordinates
    537           (requires computation of a new interpolator)
    538           Maybe not,.,.
    539           FIXME: What if x and y are vectors?
    540         """
    541 
    542         from math import pi, cos, sin, sqrt
    543         from Numeric import zeros, Float
    544 
    545 
    546         if point_id is None:
    547             msg = 'NetCDF File function needs a point_id when invoked'
    548             raise msg
    549 
    550 
    551         #Find time tau such that
    552         #
    553         # domain.starttime + t == self.starttime + tau
    554 
    555         if self.domain is not None:
    556             tau = self.domain.starttime-self.starttime+t
    557         else:
    558             tau = t
    559 
    560 
    561         msg = 'Time interval derived from file %s [%s:%s]'\
    562               %(self.filename, self.T[0], self.T[1])
    563         msg += ' does not match model time: %s\n' %tau
    564         msg += 'Domain says its starttime == %f' %(self.domain.starttime)
    565 
    566         if tau < self.T[0]: raise msg
    567         if tau > self.T[-1]: raise msg
    568 
    569 
    570         oldindex = self.index #Time index
    571 
    572         #Find time slot
    573         while tau > self.T[self.index]: self.index += 1
    574         while tau < self.T[self.index]: self.index -= 1
    575 
    576 
    577         if tau == self.T[self.index]:
    578             #Protect against case where tau == T[-1] (last time)
    579             # - also works in general when tau == T[i]
    580             ratio = 0
    581         else:
    582             #t is now between index and index+1
    583             ratio = (tau - self.T[self.index])/\
    584                     (self.T[self.index+1] - self.T[self.index])
    585 
    586 
    587 
    588 
    589         #Compute interpolated values
    590         q = zeros( len(self.quantities), Float)
    591 
    592         for i, name in enumerate(self.quantities):
    593             Q = self.values[name]
    594 
    595             if ratio > 0:
    596                 q[i] = Q[self.index, point_id] +\
    597                        ratio*(Q[self.index+1, point_id] -\
    598                               Q[self.index, point_id])
    599             else:
    600                 q[i] = Q[self.index, point_id]
    601 
    602 
    603 
    604         #Return vector of interpolated values
    605         return q
    606 
    607 #FIXME Obsolete
    608 class File_function_ASCII:
    609     """Read time series from file and return a callable object:
    610     f(t,x,y)
    611     which will return interpolated values based on the input file.
    612 
    613     The file format is assumed to be either two fields separated by a comma:
    614 
    615         time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
    616 
    617     or four comma separated fields
    618 
    619         time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
    620 
    621     In either case, the callable object will return a tuple of
    622     interpolated values, one each value stated in the file.
    623 
    624 
    625     E.g
    626 
    627       31/08/04 00:00:00, 1.328223 0 0
    628       31/08/04 00:15:00, 1.292912 0 0
    629 
    630     will provide a time dependent function f(t,x=None,y=None),
    631     ignoring x and y, which returns three values per call.
    632 
    633 
    634     NOTE: At this stage the function is assumed to depend on
    635     time only, i.e  no spatial dependency!!!!!
    636     When that is needed we can use the least_squares interpolation.
    637 
    638     #FIXME: This should work with netcdf (e.g. sww) and thus render the
    639     #spatio-temporal boundary condition in shallow water fully general
    640 
    641     #FIXME: Specified quantites not used here -
    642     #always return whatever is in the file
    643     """
    644 
    645 
    646     def __init__(self, filename, domain=None, quantities = None, interpolation_points=None):
    647         """Initialise callable object from a file.
    648         See docstring for class File_function
    649 
    650         If domain is specified, model time (domain,starttime)
    651         will be checked and possibly modified
    652 
    653         All times are assumed to be in UTC
    654         """
    655 
    656         import time, calendar
    657         from Numeric import array
    658         from config import time_format
    659 
    660         fid = open(filename)
    661         line = fid.readline()
    662         fid.close()
    663 
    664         fields = line.split(',')
    665         msg = 'File %s must have the format date, value0 value1 value2 ...'
    666         msg += ' or date, x, y, value0 value1 value2 ...'
    667         mode = len(fields)
    668         assert mode in [2,4], msg
    669 
    670         #time.strptime(fields[0], time_format)
    671         try:
    672             starttime = calendar.timegm(time.strptime(fields[0], time_format))
    673         except ValueError:
    674             msg = 'First field in file %s must be' %filename
    675             msg += ' date-time with format %s.\n' %time_format
    676             msg += 'I got %s instead.' %fields[0]
    677             raise msg
    678 
    679 
    680         #Split values
    681         values = []
    682         for value in fields[mode-1].split():
    683             values.append(float(value))
    684 
    685         q = ensure_numeric(values)
    686 
    687         msg = 'ERROR: File must contain at least one independent value'
    688         assert len(q.shape) == 1, msg
    689 
    690         self.number_of_values = len(q)
    691 
    692         self.filename = filename
    693         self.starttime = starttime
    694         self.domain = domain
    695 
    696         if domain is not None:
    697             msg = 'WARNING: Start time as specified in domain (%s)'\
    698                   %domain.starttime
    699             msg += ' is earlier than the starttime of file %s: %s.'\
    700                      %(self.filename, self.starttime)
    701             msg += 'Modifying starttime accordingly.'
    702             if self.starttime > domain.starttime:
    703                 #FIXME: Print depending on some verbosity setting
    704                 ##if verbose: print msg
    705                 domain.starttime = self.starttime #Modifying model time
    706 
    707             #if domain.starttime is None:
    708             #    domain.starttime = self.starttime
    709             #else:
    710             #    msg = 'WARNING: Start time as specified in domain (%s)'\
    711             #          %domain.starttime
    712             #    msg += ' is earlier than the starttime of file %s: %s.'\
    713             #             %(self.filename, self.starttime)
    714             #    msg += 'Modifying starttime accordingly.'
    715             #    if self.starttime > domain.starttime:
    716             #        #FIXME: Print depending on some verbosity setting
    717             #        #print msg
    718             #        domain.starttime = self.starttime #Modifying model time
    719 
    720         if mode == 2:
    721             self.read_times() #Now read all times in.
    722         else:
    723             raise 'x,y dependency not yet implemented'
    724 
    725 
    726     def read_times(self):
    727         """Read time and values
    728         """
    729         from Numeric import zeros, Float, alltrue
    730         from config import time_format
    731         import time, calendar
    732 
    733         fid = open(self.filename)
    734         lines = fid.readlines()
    735         fid.close()
    736 
    737         N = len(lines)
    738         d = self.number_of_values
    739 
    740         T = zeros(N, Float)       #Time
    741         Q = zeros((N, d), Float)  #Values
    742 
    743         for i, line in enumerate(lines):
    744             fields = line.split(',')
    745             realtime = calendar.timegm(time.strptime(fields[0], time_format))
    746 
    747             T[i] = realtime - self.starttime
    748 
    749             for j, value in enumerate(fields[1].split()):
    750                 Q[i, j] = float(value)
    751 
    752         msg = 'File %s must list time as a monotonuosly ' %self.filename
    753         msg += 'increasing sequence'
    754         assert alltrue( T[1:] - T[:-1] > 0 ), msg
    755 
    756         self.T = T     #Time
    757         self.Q = Q     #Values
    758         self.index = 0 #Initial index
    759 
    760 
    761     def __repr__(self):
    762         return 'File function'
    763 
    764     def __call__(self, t, x=None, y=None, point_id=None):
    765         """Evaluate f(t,x,y)
    766 
    767         FIXME: x, y dependency not yet implemented except that
    768         result is a vector of same length as x and y
    769         FIXME: Naaaa
    770 
    771         FIXME: Rethink semantics when x,y are vectors.
    772         """
    773 
    774         from math import pi, cos, sin, sqrt
    775 
    776 
    777         #Find time tau such that
    778         #
    779         # domain.starttime + t == self.starttime + tau
    780 
    781         if self.domain is not None:
    782             tau = self.domain.starttime-self.starttime+t
    783         else:
    784             tau = t
    785 
    786 
    787         msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
    788               %(self.filename, self.T[0], self.T[1], tau)
    789         if tau < self.T[0]: raise msg
    790         if tau > self.T[-1]: raise msg
    791 
    792         oldindex = self.index
    793 
    794         #Find slot
    795         while tau > self.T[self.index]: self.index += 1
    796         while tau < self.T[self.index]: self.index -= 1
    797 
    798         #t is now between index and index+1
    799         ratio = (tau - self.T[self.index])/\
    800                 (self.T[self.index+1] - self.T[self.index])
    801 
    802         #Compute interpolated values
    803         q = self.Q[self.index,:] +\
    804             ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
    805 
    806         #Return vector of interpolated values
    807         if x == None and y == None:
    808             return q
    809         else:
    810             try:
    811                 N = len(x)
    812             except:
    813                 return q
    814             else:
    815                 from Numeric import ones, Float
    816                 #x is a vector - Create one constant column for each value
    817                 N = len(x)
    818                 assert len(y) == N, 'x and y must have same length'
    819 
    820                 res = []
    821                 for col in q:
    822                     res.append(col*ones(N, Float))
    823 
    824                 return res
    825 
    826 
    827 #FIXME: TEMP
    828 class File_function_copy:
    829     """Read time series from file and return a callable object:
    830     f(t,x,y)
    831     which will return interpolated values based on the input file.
    832 
    833     The file format is assumed to be either two fields separated by a comma:
    834 
    835         time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
    836 
    837     or four comma separated fields
    838 
    839         time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
    840 
    841     In either case, the callable object will return a tuple of
    842     interpolated values, one each value stated in the file.
    843 
    844 
    845     E.g
    846 
    847       31/08/04 00:00:00, 1.328223 0 0
    848       31/08/04 00:15:00, 1.292912 0 0
    849 
    850     will provide a time dependent function f(t,x=None,y=None),
    851     ignoring x and y, which returns three values per call.
    852 
    853 
    854     NOTE: At this stage the function is assumed to depend on
    855     time only, i.e  no spatial dependency!!!!!
    856     When that is needed we can use the least_squares interpolation.
    857 
    858     #FIXME: This should work with netcdf (e.g. sww) and thus render the
    859     #spatio-temporal boundary condition in shallow water fully general
    860     """
    861 
    862 
    863     def __init__(self, filename, domain=None):
    864         """Initialise callable object from a file.
    865         See docstring for class File_function
    866 
    867         If domain is specified, model time (domain,starttime)
    868         will be checked and possibly modified
    869 
    870         All times are assumed to be in UTC
    871         """
    872 
    873         import time, calendar
    874         from Numeric import array
    875         from config import time_format
    876 
    877         assert type(filename) == type(''),\
    878                'First argument to File_function must be a string'
    879 
    880 
    881         try:
    882             fid = open(filename)
    883         except Exception, e:
    884             msg = 'File "%s" could not be opened: Error="%s"'\
    885                   %(filename, e)
    886             raise msg
    887 
    888 
    889         line = fid.readline()
    890         fid.close()
    891         fields = line.split(',')
    892         msg = 'File %s must have the format date, value0 value1 value2 ...'
    893         msg += ' or date, x, y, value0 value1 value2 ...'
    894         mode = len(fields)
    895         assert mode in [2,4], msg
    896 
    897         try:
    898             starttime = calendar.timegm(time.strptime(fields[0], time_format))
    899         except ValueError:
    900             msg = 'First field in file %s must be' %filename
    901             msg += ' date-time with format %s.\n' %time_format
    902             msg += 'I got %s instead.' %fields[0]
    903             raise msg
    904 
    905 
    906         #Split values
    907         values = []
    908         for value in fields[mode-1].split():
    909             values.append(float(value))
    910 
    911         q = ensure_numeric(values)
    912 
    913         msg = 'ERROR: File must contain at least one independent value'
    914         assert len(q.shape) == 1, msg
    915 
    916         self.number_of_values = len(q)
    917 
    918         self.filename = filename
    919         self.starttime = starttime
    920         self.domain = domain
    921 
    922         if domain is not None:
    923             if domain.starttime is None:
    924                 domain.starttime = self.starttime
    925             else:
    926                 msg = 'WARNING: Start time as specified in domain (%s)'\
    927                       %domain.starttime
    928                 msg += ' is earlier than the starttime of file %s: %s.'\
    929                          %(self.filename, self.starttime)
    930                 msg += 'Modifying starttime accordingly.'
    931                 if self.starttime > domain.starttime:
    932                     #FIXME: Print depending on some verbosity setting
    933                     #print msg
    934                     domain.starttime = self.starttime #Modifying model time
    935 
    936         if mode == 2:
    937             self.read_times() #Now read all times in.
    938         else:
    939             raise 'x,y dependency not yet implemented'
    940 
    941 
    942     def read_times(self):
    943         """Read time and values
    944         """
    945         from Numeric import zeros, Float, alltrue
    946         from config import time_format
    947         import time, calendar
    948 
    949         fid = open(self.filename)
    950         lines = fid.readlines()
    951         fid.close()
    952 
    953         N = len(lines)
    954         d = self.number_of_values
    955 
    956         T = zeros(N, Float)       #Time
    957         Q = zeros((N, d), Float)  #Values
    958 
    959         for i, line in enumerate(lines):
    960             fields = line.split(',')
    961             realtime = calendar.timegm(time.strptime(fields[0], time_format))
    962 
    963             T[i] = realtime - self.starttime
    964 
    965             for j, value in enumerate(fields[1].split()):
    966                 Q[i, j] = float(value)
    967 
    968         msg = 'File %s must list time as a monotonuosly ' %self.filename
    969         msg += 'increasing sequence'
    970         assert alltrue( T[1:] - T[:-1] > 0 ), msg
    971 
    972         self.T = T     #Time
    973         self.Q = Q     #Values
    974         self.index = 0 #Initial index
    975 
    976 
    977     def __repr__(self):
    978         return 'File function'
    979 
    980 
    981 
    982     def __call__(self, t, x=None, y=None):
    983         """Evaluate f(t,x,y)
    984 
    985         FIXME: x, y dependency not yet implemented except that
    986         result is a vector of same length as x and y
    987         FIXME: Naaaa
    988         """
    989 
    990         from math import pi, cos, sin, sqrt
    991 
    992 
    993         #Find time tau such that
    994         #
    995         # domain.starttime + t == self.starttime + tau
    996 
    997         if self.domain is not None:
    998             tau = self.domain.starttime-self.starttime+t
    999         else:
    1000             tau = t
    1001 
    1002 
    1003         msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
    1004               %(self.filename, self.T[0], self.T[1], tau)
    1005         if tau < self.T[0]: raise msg
    1006         if tau > self.T[-1]: raise msg
    1007 
    1008         oldindex = self.index
    1009 
    1010         #Find slot
    1011         while tau > self.T[self.index]: self.index += 1
    1012         while tau < self.T[self.index]: self.index -= 1
    1013 
    1014         #t is now between index and index+1
    1015         ratio = (tau - self.T[self.index])/\
    1016                 (self.T[self.index+1] - self.T[self.index])
    1017 
    1018         #Compute interpolated values
    1019         q = self.Q[self.index,:] +\
    1020             ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
    1021 
    1022         #Return vector of interpolated values
    1023         if x == None and y == None:
    1024             return q
    1025         else:
    1026             try:
    1027                 N = len(x)
    1028             except:
    1029                 return q
    1030             else:
    1031                 from Numeric import ones, Float
    1032                 #x is a vector - Create one constant column for each value
    1033                 N = len(x)
    1034                 assert len(y) == N, 'x and y must have same length'
    1035 
    1036                 res = []
    1037                 for col in q:
    1038                     res.append(col*ones(N, Float))
    1039 
    1040                 return res
     316   
     317
     318
    1041319
    1042320
    … …  
    1052330    """
    1053331    #FIXME: Probably obsoleted by similar function in load_ASCII
    1054     #FIXE: Name read_data_points (or see 'load' in pylab)
     332    #FIXME: Name read_data_points (or see 'load' in pylab)
    1055333
    1056334    from Scientific.IO.NetCDF import NetCDFFile
    … …  
    1512790    to check if a tsh/msh file 'looks' good.
    1513791    """
     792
     793    #FIXME Move to data_manager
    1514794    from shallow_water import Domain
    1515795    from pmesh2domain import pmesh_to_domain_instance
    … …  
    1596876
    1597877
    1598 
    1599 
    1600 #OBSOLETED STUFF
    1601 def inside_polygon_old(point, polygon, closed = True, verbose = False):
    1602     #FIXME Obsoleted
    1603     """Determine whether points are inside or outside a polygon
    1604 
    1605     Input:
    1606        point - Tuple of (x, y) coordinates, or list of tuples
    1607        polygon - list of vertices of polygon
    1608        closed - (optional) determine whether points on boundary should be
    1609        regarded as belonging to the polygon (closed = True)
    1610        or not (closed = False)
    1611 
    1612     Output:
    1613        If one point is considered, True or False is returned.
    1614        If multiple points are passed in, the function returns indices
    1615        of those points that fall inside the polygon
    1616 
    1617     Examples:
    1618        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
    1619        inside_polygon( [0.5, 0.5], U)
    1620        will evaluate to True as the point 0.5, 0.5 is inside the unit square
    1621 
    1622        inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
    1623        will return the indices [0, 2] as only the first and the last point
    1624        is inside the unit square
    1625 
    1626     Remarks:
    1627        The vertices may be listed clockwise or counterclockwise and
    1628        the first point may optionally be repeated.
    1629        Polygons do not need to be convex.
    1630        Polygons can have holes in them and points inside a hole is
    1631        regarded as being outside the polygon.
    1632 
    1633 
    1634     Algorithm is based on work by Darel Finley,
    1635     http://www.alienryderflex.com/polygon/
    1636 
    1637     """
    1638 
    1639     from Numeric import array, Float, reshape
    1640 
    1641 
    1642     #Input checks
    1643     try:
    1644         point = array(point).astype(Float)
    1645     except:
    1646         msg = 'Point %s could not be converted to Numeric array' %point
    1647         raise msg
    1648 
    1649     try:
    1650         polygon = array(polygon).astype(Float)
    1651     except:
    1652         msg = 'Polygon %s could not be converted to Numeric array' %polygon
    1653         raise msg
    1654 
    1655     assert len(polygon.shape) == 2,\
    1656        'Polygon array must be a 2d array of vertices: %s' %polygon
    1657 
    1658 
    1659     assert polygon.shape[1] == 2,\
    1660        'Polygon array must have two columns: %s' %polygon
    1661 
    1662 
    1663     if len(point.shape) == 1:
    1664         one_point = True
    1665         points = reshape(point, (1,2))
    1666     else:
    1667         one_point = False
    1668         points = point
    1669 
    1670     N = polygon.shape[0] #Number of vertices in polygon
    1671     M = points.shape[0]  #Number of points
    1672 
    1673     px = polygon[:,0]
    1674     py = polygon[:,1]
    1675 
    1676     #Used for an optimisation when points are far away from polygon
    1677     minpx = min(px); maxpx = max(px)
    1678     minpy = min(py); maxpy = max(py)
    1679 
    1680 
    1681     #Begin main loop (FIXME: It'd be crunchy to have this written in C:-)
    1682     indices = []
    1683     for k in range(M):
    1684 
    1685         if verbose:
    1686             if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
    1687 
    1688         #for each point
    1689         x = points[k, 0]
    1690         y = points[k, 1]
    1691 
    1692         inside = False
    1693 
    1694         #Optimisation
    1695         if x > maxpx or x < minpx: continue
    1696         if y > maxpy or y < minpy: continue
    1697 
    1698         #Check polygon
    1699         for i in range(N):
    1700             j = (i+1)%N
    1701 
    1702             #Check for case where point is contained in line segment
    1703             if point_on_line(x, y, px[i], py[i], px[j], py[j]):
    1704                 if closed:
    1705                     inside = True
    1706                 else:
    1707                     inside = False
    1708                 break
    1709             else:
    1710                 #Check if truly inside polygon
    1711                 if py[i] < y and py[j] >= y or\
    1712                    py[j] < y and py[i] >= y:
    1713                     if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
    1714                         inside = not inside
    1715 
    1716         if inside: indices.append(k)
    1717 
    1718     if one_point:
    1719         return inside
    1720     else:
    1721         return indices
    1722 
    1723 
    1724 #def remove_from(A, B):
    1725 #    """Assume that A
    1726 #    """
    1727 #    from Numeric import search_sorted##
    1728 #
    1729 #    ind = search_sorted(A, B)
    1730 
    1731 
    1732 
    1733 def outside_polygon_old(point, polygon, closed = True, verbose = False):
    1734     #OBSOLETED
    1735     """Determine whether points are outside a polygon
    1736 
    1737     Input:
    1738        point - Tuple of (x, y) coordinates, or list of tuples
    1739        polygon - list of vertices of polygon
    1740        closed - (optional) determine whether points on boundary should be
    1741        regarded as belonging to the polygon (closed = True)
    1742        or not (closed = False)
    1743 
    1744     Output:
    1745        If one point is considered, True or False is returned.
    1746        If multiple points are passed in, the function returns indices
    1747        of those points that fall outside the polygon
    1748 
    1749     Examples:
    1750        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
    1751        outside_polygon( [0.5, 0.5], U )
    1752        will evaluate to False as the point 0.5, 0.5 is inside the unit square
    1753 
    1754        ouside_polygon( [1.5, 0.5], U )
    1755        will evaluate to True as the point 1.5, 0.5 is outside the unit square
    1756 
    1757        outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
    1758        will return the indices [1] as only the first and the last point
    1759        is inside the unit square
    1760     """
    1761 
    1762     #FIXME: This is too slow
    1763 
    1764     res  = inside_polygon(point, polygon, closed, verbose)
    1765 
    1766     if res is True or res is False:
    1767         return not res
    1768 
    1769     #Now invert indices
    1770     from Numeric import arrayrange, compress
    1771     outside_indices = arrayrange(len(point))
    1772     for i in res:
    1773         outside_indices = compress(outside_indices != i, outside_indices)
    1774     return outside_indices
    1775 
    1776 def inside_polygon_c(point, polygon, closed = True, verbose = False):
    1777     #FIXME: Obsolete
    1778     """Determine whether points are inside or outside a polygon
    1779 
    1780     C-wrapper
    1781     """
    1782 
    1783     from Numeric import array, Float, reshape, zeros, Int
    1784 
    1785 
    1786     if verbose: print 'Checking input to inside_polygon'
    1787     #Input checks
    1788     try:
    1789         point = array(point).astype(Float)
    1790     except:
    1791         msg = 'Point %s could not be converted to Numeric array' %point
    1792         raise msg
    1793 
    1794     try:
    1795         polygon = array(polygon).astype(Float)
    1796     except:
    1797         msg = 'Polygon %s could not be converted to Numeric array' %polygon
    1798         raise msg
    1799 
    1800     assert len(polygon.shape) == 2,\
    1801        'Polygon array must be a 2d array of vertices: %s' %polygon
    1802 
    1803 
    1804     assert polygon.shape[1] == 2,\
    1805        'Polygon array must have two columns: %s' %polygon
    1806 
    1807 
    1808     if len(point.shape) == 1:
    1809         one_point = True
    1810         points = reshape(point, (1,2))
    1811     else:
    1812         one_point = False
    1813         points = point
    1814 
    1815     from util_ext import inside_polygon
    1816 
    1817     if verbose: print 'Allocating array for indices'
    1818 
    1819     indices = zeros( points.shape[0], Int )
    1820 
    1821     if verbose: print 'Calling C-version of inside poly'
    1822     count = inside_polygon(points, polygon, indices,
    1823                            int(closed), int(verbose))
    1824 
    1825     if one_point:
    1826         return count == 1 #Return True if the point was inside
    1827     else:
    1828         if verbose: print 'Got %d points' %count
    1829 
    1830         return indices[:count]   #Return those indices that were inside
Note: See TracChangeset for help on using the changeset viewer.