Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1671)
@@ -734,6 +734,8 @@
 
 
-
+#########################################################
 #Conversion routines
+########################################################
+
 def sww2obj(basefilename, size):
     """Convert netcdf based data output to obj
@@ -1943,4 +1945,133 @@
     outfile.close()
 
+
+
+
+def timefile2swww(filename, quantity_names = None):
+    """Template for converting typical text files with time series to
+    NetCDF sww file.
+
+
+    The file format is assumed to be either two fields separated by a comma:
+
+        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
+
+    E.g
+
+      31/08/04 00:00:00, 1.328223 0 0
+      31/08/04 00:15:00, 1.292912 0 0
+
+    will provide a time dependent function f(t) with three attributes
+
+    filename is assumed to be the rootname with extenisons .txt and .sww
+    """
+
+    import time, calendar
+    from Numeric import array
+    from config import time_format
+    from util import ensure_numeric
+
+
+    fid = open(filename + '.txt')
+    line = fid.readline()
+    fid.close()
+
+    fields = line.split(',')
+    msg = 'File %s must have the format date, value0 value1 value2 ...'
+    assert len(fields) == 2, msg
+
+    try:
+        starttime = calendar.timegm(time.strptime(fields[0], time_format))
+    except ValueError:
+        msg = 'First field in file %s must be' %filename
+        msg += ' date-time with format %s.\n' %time_format
+        msg += 'I got %s instead.' %fields[0]
+        raise msg
+
+
+    #Split values
+    values = []
+    for value in fields[1].split():
+        values.append(float(value))
+
+    q = ensure_numeric(values)
+
+    msg = 'ERROR: File must contain at least one independent value'
+    assert len(q.shape) == 1, msg
+
+
+
+    #Read times proper
+    from Numeric import zeros, Float, alltrue
+    from config import time_format
+    import time, calendar
+
+    fid = open(filename + '.txt')
+    lines = fid.readlines()
+    fid.close()
+
+    N = len(lines)
+    d = len(q)
+
+    T = zeros(N, Float)       #Time
+    Q = zeros((N, d), Float)  #Values
+
+    for i, line in enumerate(lines):
+        fields = line.split(',')
+        realtime = calendar.timegm(time.strptime(fields[0], time_format))
+        
+        T[i] = realtime - starttime
+
+        for j, value in enumerate(fields[1].split()):
+            Q[i, j] = float(value)
+
+    msg = 'File %s must list time as a monotonuosly ' %filename
+    msg += 'increasing sequence'
+    assert alltrue( T[1:] - T[:-1] > 0 ), msg
+
+
+    #Create sww file
+    from Scientific.IO.NetCDF import NetCDFFile
+
+        
+
+    fid = NetCDFFile(filename + '.sww', 'w')
+
+
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'Time series'
+
+
+    #Reference point
+    #Start time in seconds since the epoch (midnight 1/1/1970)
+    #FIXME: Use Georef
+    fid.starttime = starttime
+    
+    
+    # dimension definitions
+    #fid.createDimension('number_of_volumes', self.number_of_volumes)
+    #fid.createDimension('number_of_vertices', 3)
+
+
+    fid.createDimension('number_of_timesteps', len(T))
+
+    #if geo_reference is not None:
+    #    geo_reference.write_NetCDF(fid)
+		
+
+    fid.createVariable('time', Float, ('number_of_timesteps',))
+
+    fid.variables['time'][:] = T
+    
+    for i in range(Q.shape[1]):
+        try:
+            name = quantity_names[i]
+        except:    
+            name = 'Attribute%d'%i
+            
+        fid.createVariable(name, Float, ('number_of_timesteps',))
+        fid.variables[name][:] = Q[:,i]        
+
+    fid.close()
 
 
Index: /inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1671)
@@ -1156,17 +1156,36 @@
                 q[i] = Q0
 
-            #if ratio > 0:
-            #    q[i] = Q[self.index, point_id] +\
-            #           ratio*(Q[self.index+1, point_id] -\
-            #                  Q[self.index, point_id])
-            #else:
-            #    q[i] = Q[self.index, point_id]
 
         #Return vector of interpolated values
-        if len(q) == 1:
-            return q[0]
+        #if len(q) == 1:
+        #    return q[0]
+        #else:
+        #    return q
+
+
+        #Return vector of interpolated values
+        #FIXME:
+        if self.spatial is True:
+            return q
         else:
-            return q
-
+            #Replicate q according to x and y
+            #This is e.g used for Wind_stress
+            if x == None or y == None: 
+                return q
+            else:
+                try:
+                    N = len(x)
+                except:
+                    return q
+                else:
+                    from Numeric import ones, Float
+                    #x is a vector - Create one constant column for each value
+                    N = len(x)
+                    assert len(y) == N, 'x and y must have same length'
+                    res = []
+                    for col in q:
+                        res.append(col*ones(N, Float))
+                        
+                return res
 
 
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1671)
@@ -1141,5 +1141,5 @@
         y = ones(3, Float)
         try:
-            q = f(1.0, x, y)
+            q = f(1.0, x=x, y=y)
         except Exception, e:
             msg = 'Function %s could not be executed:\n%s' %(f, e)
@@ -1155,5 +1155,6 @@
             raise msg
 
-        msg = 'Return vector from function %s' %f
+        #Is this really what we want?
+        msg = 'Return vector from function %s ' %f
         msg += 'must have same lenght as input vectors'
         assert len(q) == N, msg
@@ -1223,6 +1224,6 @@
             #Assume vector function returning (s, phi)(t,x,y)
             vector_function = args[0]
-            s = lambda t,x,y: vector_function(t,x,y)[0]
-            phi = lambda t,x,y: vector_function(t,x,y)[1]
+            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
         else:
            #Assume info is in 2 keyword arguments
Index: /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 1671)
@@ -138,5 +138,5 @@
         #Write file
         filename = 'boundarytest' + str(time.time())
-        fid = open(filename, 'w')
+        fid = open(filename + '.txt', 'w')
         start = time.mktime(time.strptime('2000', '%Y'))
         dt = 5*60  #Five minute intervals
@@ -149,6 +149,16 @@
 
 
-        F = File_boundary(filename, domain)
-        os.remove(filename)
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        
+        from data_manager import timefile2swww                
+        timefile2swww(filename, quantity_names = ['stage', 'ymomentum'])
+        
+
+
+        F = File_boundary(filename + '.sww', domain)
+
+        
+        os.remove(filename + '.txt')
+        os.remove(filename + '.sww')        
 
 
@@ -182,5 +192,5 @@
 
     def test_fileboundary_exception(self):
-        """Test that boundary object complians if number of
+        """Test that boundary object complains if number of
         conserved quantities are wrong
         """
@@ -221,5 +231,5 @@
         #Write file (with only two values)
         filename = 'boundarytest' + str(time.time())
-        fid = open(filename, 'w')
+        fid = open(filename + '.txt', 'w')
         start = time.mktime(time.strptime('2000', '%Y'))
         dt = 5*60  #Five minute intervals
@@ -232,12 +242,20 @@
 
 
-        try:
-            F = File_boundary(filename, domain)
-        except AssertionError:
-            pass
-        else:
-            raise 'Should have raised an exception'
-
-        os.remove(filename)
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        from data_manager import timefile2swww                
+        timefile2swww(filename, quantity_names = ['stage', 'xmomentum'])
+
+        
+        try:
+            F = File_boundary(filename + '.sww', domain)            
+        except:
+            pass
+        else:
+            raise 'Should have raised an exception'        
+        
+        os.remove(filename + '.txt')
+        os.remove(filename + '.sww')        
+
+
 
 
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1671)
@@ -994,7 +994,7 @@
         #Write wind stress file (ensure that domain.time is covered)
         #Take x=1 and y=0
-        filename = 'test_windstress_from_file.txt'
+        filename = 'test_windstress_from_file'
         start = time.mktime(time.strptime('2000', '%Y'))
-        fid = open(filename, 'w')
+        fid = open(filename + '.txt', 'w')
         dt = 1  #One second interval
         t = 0.0
@@ -1009,7 +1009,31 @@
         fid.close()
 
+
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        from data_manager import timefile2swww        
+        timefile2swww(filename)
+
+
+        
         #Setup wind stress
-        F = file_function(filename)
+        F = file_function(filename + '.sww', quantities = ['Attribute0',
+                                                           'Attribute1'])
+
+        #print 'F(5)', F(5)
+
+        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
+        
+        #print dir(F)
+        #print F.T
+        #print F.precomputed_values
+        #
+        #F = file_function(filename + '.txt')        
+        #
+        #print dir(F)
+        #print F.T        
+        #print F.Q
+        
         W = Wind_stress(F)
+
         domain.forcing_terms = []
         domain.forcing_terms.append(W)
@@ -1043,5 +1067,5 @@
             assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
 
-        os.remove(filename)
+        #os.remove(filename + '.txt')
 
     def test_wind_stress_error_condition(self):
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 1671)
@@ -8,4 +8,5 @@
 from util import *
 from config import epsilon
+from data_manager import timefile2swww
 
 
@@ -171,8 +172,54 @@
         assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
 
-
-
-
-    def test_file_function_time(self):
+    def test_ensure_numeric(self):
+        from util import ensure_numeric
+        from Numeric import ArrayType, Float, array
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
+
+
+        A = [1,2,3.14,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'        
+        assert A == B    
+        assert A is B   #Same object
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'        
+        assert A == B    
+        assert A is not B   #Not the same object        
+
+
+
+    def test_file_function_time1(self):
         """Test that File function interpolates correctly
         between given times. No x,y dependency here.
@@ -184,7 +231,8 @@
         from math import sin, pi
 
+        #Typical ASCII file
         finaltime = 1200
-        filename = 'test_file_function.txt'
-        fid = open(filename, 'w')
+        filename = 'test_file_function'
+        fid = open(filename + '.txt', 'w')
         start = time.mktime(time.strptime('2000', '%Y'))
         dt = 60  #One minute intervals
@@ -197,6 +245,13 @@
         fid.close()
 
-        F = file_function(filename)
-
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        timefile2swww(filename)
+
+
+        #Create file function
+        F = file_function(filename + '.sww', quantities = ['Attribute0',
+                                                           'Attribute1',
+                                                           'Attribute2'])
+        
         #Now try interpolation
         for i in range(20):
@@ -223,53 +278,6 @@
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
-        os.remove(filename)
-
-
-    def test_ensure_numeric(self):
-        from util import ensure_numeric
-        from Numeric import ArrayType, Float, array
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'
-        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
-
-
-        A = [1,2,3.14,4]
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
-
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
-
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
-
-
-        A = array([1,2,3,4])
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'        
-        assert A == B    
-        assert A is B   #Same object
-
-
-        A = array([1,2,3,4])
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'        
-        assert A == B    
-        assert A is not B   #Not the same object        
+        os.remove(filename + '.txt')
+        os.remove(filename + '.sww')        
 
 
@@ -649,6 +657,6 @@
 
         finaltime = 1200
-        filename = 'test_file_function.txt'
-        fid = open(filename, 'w')
+        filename = 'test_file_function'
+        fid = open(filename + '.txt', 'w')
         start = time.mktime(time.strptime('2000', '%Y'))
         dt = 60  #One minute intervals
@@ -661,4 +669,10 @@
         fid.close()
 
+
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        timefile2swww(filename)
+
+
+
         a = [0.0, 0.0]
         b = [4.0, 0.0]
@@ -670,20 +684,27 @@
 
         #Check that domain.starttime is updated if non-existing
-        F = file_function(filename, domain)
+        F = file_function(filename + '.sww', domain)
+        
         assert allclose(domain.starttime, start)
 
         #Check that domain.starttime is updated if too early
         domain.starttime = start - 1
-        F = file_function(filename, domain)
+        F = file_function(filename + '.sww', domain)
         assert allclose(domain.starttime, start)
 
         #Check that domain.starttime isn't updated if later
         domain.starttime = start + 1
-        F = file_function(filename, domain)
+        F = file_function(filename + '.sww', domain)        
         assert allclose(domain.starttime, start+1)
 
         domain.starttime = start
-
-
+        F = file_function(filename + '.sww', domain,
+                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
+        
+
+        #print F.T
+        #print F.precomputed_values
+        #print 'F(60)', F(60)
+        
         #Now try interpolation
         for i in range(20):
@@ -710,6 +731,6 @@
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
-        os.remove(filename)
-
+        os.remove(filename + '.sww')
+        os.remove(filename + '.txt')        
 
     def test_file_function_time_with_domain_different_start(self):
@@ -728,6 +749,6 @@
 
         finaltime = 1200
-        filename = 'test_file_function.txt'
-        fid = open(filename, 'w')
+        filename = 'test_file_function'
+        fid = open(filename + '.txt', 'w')
         start = time.mktime(time.strptime('2000', '%Y'))
         dt = 60  #One minute intervals
@@ -740,4 +761,7 @@
         fid.close()
 
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        timefile2swww(filename)        
+
         a = [0.0, 0.0]
         b = [4.0, 0.0]
@@ -752,5 +776,7 @@
 	delta = 23
         domain.starttime = start + delta
-        F = file_function(filename, domain)
+        #F = file_function(filename + '.sww', domain)
+        F = file_function(filename + '.sww', domain,
+                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])        
         assert allclose(domain.starttime, start+delta)
 
@@ -782,5 +808,7 @@
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
-        os.remove(filename)
+
+        os.remove(filename + '.sww')
+        os.remove(filename + '.txt')                
 
 
@@ -1088,6 +1116,6 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Util,'test_spatio_temporal_file_function_time')
     suite = unittest.makeSuite(Test_Util,'test')
+    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 1670)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 1671)
@@ -6,6 +6,4 @@
 
 #FIXME: Move polygon stuff out
-#FIXME: Finish Interpolation function and get rid of File_function_ASCII
-
 
 def angle(v):
@@ -135,6 +133,10 @@
 
 
-def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False):
-    """If domain is specified, don't specify quantites as they are automatically derived
+def file_function(filename,
+                  domain = None,
+                  quantities = None,
+                  interpolation_points = None, verbose = False):
+    """If quantities is not specified, derive them from domain
+    (if that is specified) 
     """
 
@@ -157,8 +159,8 @@
     fid.close()
 
-    if domain is not None:
-        quantities = domain.conserved_quantities
-    else:
-        quantities = None
+    if quantities is None: 
+        if domain is not None:
+            quantities = domain.conserved_quantities
+
 
 
@@ -167,14 +169,5 @@
                                         interpolation_points, verbose = verbose)
     else:
-        #FIXME Should be phased out
-        return File_function_ASCII(filename, domain, quantities,
-                                   interpolation_points)
-    
-    
-
-
-
-
-
+        raise 'Must be a NetCDF File'
 
 
@@ -215,9 +208,16 @@
 
     if quantity_names is None or len(quantity_names) < 1:
-        msg = 'ERROR: At least one independent value must be specified'
-        raise msg
+        #Get quantities from file
+        quantity_names = []
+        for name in fid.variables.keys():
+            if name not in ['x', 'y', 'time']:
+                quantity_names.append(name)
+                
+        #msg = 'ERROR: At least one independent value must be specified'
+        #raise msg
 
     missing = []
-    for quantity in ['x', 'y', 'time'] + quantity_names:
+        
+    for quantity in ['time'] + quantity_names:
         if not fid.variables.has_key(quantity):
             missing.append(quantity)
@@ -226,5 +226,13 @@
         msg = 'Quantities %s could not be found in file %s'\
               %(str(missing), filename)
+        fid.close()
         raise msg
+
+    #Decide whether this data has a spatial dimension
+    spatial = True
+    for quantity in ['x', 'y']:
+        if not fid.variables.has_key(quantity):
+            spatial = False
+
 
 
@@ -236,19 +244,21 @@
         raise msg
 
-    #Get origin
-    xllcorner = fid.xllcorner[0]
-    yllcorner = fid.yllcorner[0]
-
-
     #Get variables
-    if verbose: print 'Get variables'
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    triangles = fid.variables['volumes'][:]
-    time = fid.variables['time'][:]
-
-    x = reshape(x, (len(x),1))
-    y = reshape(y, (len(y),1))
-    vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+    if verbose: print 'Get variables'    
+    time = fid.variables['time'][:]    
+
+    if spatial:
+        #Get origin
+        xllcorner = fid.xllcorner[0]
+        yllcorner = fid.yllcorner[0]
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        triangles = fid.variables['volumes'][:]
+
+        x = reshape(x, (len(x),1))
+        y = reshape(y, (len(y),1))
+        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
 
     if domain is not None:
@@ -277,6 +287,7 @@
         print '  References:'
         #print '    Datum ....' #FIXME
-        print '    Lower left corner: [%f, %f]'\
-              %(xllcorner, yllcorner)
+        if spatial:
+            print '    Lower left corner: [%f, %f]'\
+                  %(xllcorner, yllcorner)
         print '    Start time:   %f' %starttime                
         
@@ -292,4 +303,8 @@
 
     from least_squares import Interpolation_function
+
+    if not spatial:
+        vertex_coordinates = triangles = interpolation_points = None         
+        
     return Interpolation_function(time,
                                   quantities,
@@ -299,744 +314,7 @@
                                   interpolation_points,
                                   verbose = verbose)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#FIXME Obsolete
-class File_function_NetCDF:
-    """Read time history of spatial data from NetCDF sww file and
-    return a callable object f(t,x,y)
-    which will return interpolated values based on the input file.
-
-    x, y may be either scalars or vectors
-
-    #FIXME: More about format, interpolation and order of quantities
-
-    The quantities returned by the callable objects are specified by
-    the list quantities which must contain the names of the
-    quantities to be returned and also reflect the order, e.g. for
-    the shallow water wave equation, on would have
-    quantities = ['stage', 'xmomentum', 'ymomentum']
-
-    interpolation_points decides at which points interpolated
-    quantities are to be computed whenever object is called.
-    If None, return average value
-    """
-
-    def  __init__(self, filename, domain=None, quantities=None,
-                  interpolation_points=None, verbose = False):
-        """Initialise callable object from a file.
-
-        If domain is specified, model time (domain.starttime)
-        will be checked and possibly modified
-
-        All times are assumed to be in UTC
-        """
-
-        #FIXME: Check that model origin is the same as file's origin
-        #(both in UTM coordinates)
-        #If not - modify those from file to match domain
-        #Take this code from e.g. dem2pts in data_manager.py
-        #FIXME: Use geo_reference to read and write xllcorner...
-        
-
-        import time, calendar
-        from config import time_format
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Open NetCDF file
-        if verbose: print 'Reading', filename
-        fid = NetCDFFile(filename, 'r')
-
-        if quantities is None or len(quantities) < 1:
-            msg = 'ERROR: At least one independent value must be specified'
-            raise msg
-
-        missing = []
-        for quantity in ['x', 'y', 'time'] + quantities:
-             if not fid.variables.has_key(quantity):
-                 missing.append(quantity)
-
-        if len(missing) > 0:
-            msg = 'Quantities %s could not be found in file %s'\
-	          %(str(missing), filename)
-            raise msg
-
-
-        #Get first timestep
-        try:
-            self.starttime = fid.starttime[0]
-        except ValueError:
-            msg = 'Could not read starttime from file %s' %filename
-            raise msg
-
-        #Get origin
-        self.xllcorner = fid.xllcorner[0]
-        self.yllcorner = fid.yllcorner[0]
-
-        self.number_of_values = len(quantities)
-        self.fid = fid
-        self.filename = filename
-	self.quantities = quantities
-	self.domain = domain
-	self.interpolation_points = interpolation_points
-
-        if domain is not None:
-            msg = 'WARNING: Start time as specified in domain (%f)'\
-                  %domain.starttime
-            msg += ' is earlier than the starttime of file %s (%f).'\
-                     %(self.filename, self.starttime)
-            msg += ' Modifying domain starttime accordingly.'
-            if self.starttime > domain.starttime:
-                if verbose: print msg
-                domain.starttime = self.starttime #Modifying model time
-                if verbose: print 'Domain starttime is now set to %f'\
-                   %domain.starttime
-
-	#Read all data in and produce values for desired data points at
-        #each timestep
-        self.spatial_interpolation(interpolation_points, verbose = verbose)
-        fid.close()
-
-    def spatial_interpolation(self, interpolation_points, verbose = False):
-        """For each timestep in fid: read surface,
-        interpolate to desired points and store values for use when
-        object is called.
-        """
-
-        from Numeric import array, zeros, Float, alltrue, concatenate, reshape
-
-        from config import time_format
-	from least_squares import Interpolation
-        import time, calendar
-
-        fid = self.fid
-
-        #Get variables
-        if verbose: print 'Get variables'
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        #z = fid.variables['elevation'][:]
-        triangles = fid.variables['volumes'][:]
-        time = fid.variables['time'][:]
-
-	#Check
-        msg = 'File %s must list time as a monotonuosly ' %self.filename
-        msg += 'increasing sequence'
-        assert alltrue(time[1:] - time[:-1] > 0 ), msg
-
-
-        if interpolation_points is not None:
-
-            try:
-	        P = ensure_numeric(interpolation_points)
-            except:
-	        msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
-		msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
-                raise msg
-
-
-	    self.T = time[:]  #Time assumed to be relative to starttime
-            self.index = 0    #Initial time index
-
-            self.values = {}
-	    for name in self.quantities:
-                self.values[name] = zeros( (len(self.T),
-		                            len(interpolation_points)),
-					    Float)
-
-            #Build interpolator
-            if verbose: print 'Build interpolation matrix'
-            x = reshape(x, (len(x),1))
-            y = reshape(y, (len(y),1))
-            vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
-
-            interpol = Interpolation(vertex_coordinates,
-                                     triangles,
-                                     point_coordinates = P,
-                                     alpha = 0,
-                                     verbose = verbose)
-
-
-            if verbose: print 'Interpolate'
-	    for i, t in enumerate(self.T):
-                #Interpolate quantities at this timestep
-                if verbose: print ' time step %d of %d' %(i, len(self.T))
-                for name in self.quantities:
-		    self.values[name][i, :] =\
-		    interpol.interpolate(fid.variables[name][i,:])
-
-            #Report
-            if verbose:
-                print '------------------------------------------------'
-                print 'File_function statistics:'
-                print '  Name: %s' %self.filename
-                print '  References:'
-                #print '    Datum ....' #FIXME
-                print '    Lower left corner: [%f, %f]'\
-                      %(self.xllcorner, self.yllcorner)
-                print '    Start time:   %f' %self.starttime                
-                #print '    Start time (file):   %f' %self.starttime
-                #print '    Start time (domain): %f' %domain.starttime
-                print '  Extent:'
-                print '    x in [%f, %f], len(x) == %d'\
-                      %(min(x.flat), max(x.flat), len(x.flat))
-                print '    y in [%f, %f], len(y) == %d'\
-                      %(min(y.flat), max(y.flat), len(y.flat))
-                print '    t in [%f, %f], len(t) == %d'\
-                      %(min(self.T), max(self.T), len(self.T))
-                print '  Quantities:'
-                for name in self.quantities:
-                    q = fid.variables[name][:].flat
-                    print '    %s in [%f, %f]' %(name, min(q), max(q))
-
-
-                print '  Interpolation points (xi, eta):'\
-                      'number of points == %d ' %P.shape[0]
-                print '    xi in [%f, %f]' %(min(P[:,0]), max(P[:,0]))
-                print '    eta in [%f, %f]' %(min(P[:,1]), max(P[:,1]))
-                print '  Interpolated quantities (over all timesteps):'
-                for name in self.quantities:
-                    q = self.values[name][:].flat
-                    print '    %s at interpolation points in [%f, %f]'\
-                          %(name, min(q), max(q))
-                print '------------------------------------------------'
-        else:
-            msg = 'File_function_NetCDF must be invoked with '
-            msg += 'a list of interpolation points'
-            raise msg
-            #FIXME: Should not be an error.
-            #__call__ could return an average or in case of time series only
-            #no point s are needed
-
-
-    def __repr__(self):
-        return 'File function'
-
-    def __call__(self, t, x=None, y=None, point_id = None):
-        """Evaluate f(t, point_id)
-
-	Inputs:
-	  t: time - Model time (tau = domain.starttime-self.starttime+t)
-	            must lie within existing timesteps
-	  point_id: index of one of the preprocessed points.
-	            If point_id is None all preprocessed points are computed
-
-	  FIXME: point_id could also be a slice
-	  FIXME: One could allow arbitrary x, y coordinates
-	  (requires computation of a new interpolator)
-          Maybe not,.,.
-	  FIXME: What if x and y are vectors?
-        """
-
-        from math import pi, cos, sin, sqrt
-        from Numeric import zeros, Float
-
-
-        if point_id is None:
-            msg = 'NetCDF File function needs a point_id when invoked'
-            raise msg
-
-
-        #Find time tau such that
-	#
-	# domain.starttime + t == self.starttime + tau
-
-	if self.domain is not None:
-            tau = self.domain.starttime-self.starttime+t
-	else:
-            tau = t
-
-
-        msg = 'Time interval derived from file %s [%s:%s]'\
-              %(self.filename, self.T[0], self.T[1])
-        msg += ' does not match model time: %s\n' %tau
-        msg += 'Domain says its starttime == %f' %(self.domain.starttime)
-
-        if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg
-
-
-        oldindex = self.index #Time index
-
-        #Find time slot
-        while tau > self.T[self.index]: self.index += 1
-        while tau < self.T[self.index]: self.index -= 1
-
-
-        if tau == self.T[self.index]:
-            #Protect against case where tau == T[-1] (last time)
-            # - also works in general when tau == T[i]
-            ratio = 0
-        else:
-            #t is now between index and index+1
-            ratio = (tau - self.T[self.index])/\
-                    (self.T[self.index+1] - self.T[self.index])
-
-
-
-
-        #Compute interpolated values
-        q = zeros( len(self.quantities), Float)
-
-	for i, name in enumerate(self.quantities):
-            Q = self.values[name]
-
-            if ratio > 0:
-                q[i] = Q[self.index, point_id] +\
-                       ratio*(Q[self.index+1, point_id] -\
-                              Q[self.index, point_id])
-            else:
-                q[i] = Q[self.index, point_id]
-
-
-
-        #Return vector of interpolated values
-        return q
-
-#FIXME Obsolete
-class File_function_ASCII:
-    """Read time series from file and return a callable object:
-    f(t,x,y)
-    which will return interpolated values based on the input file.
-
-    The file format is assumed to be either two fields separated by a comma:
-
-        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
-
-    or four comma separated fields
-
-        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
-
-    In either case, the callable object will return a tuple of
-    interpolated values, one each value stated in the file.
-
-
-    E.g
-
-      31/08/04 00:00:00, 1.328223 0 0
-      31/08/04 00:15:00, 1.292912 0 0
-
-    will provide a time dependent function f(t,x=None,y=None),
-    ignoring x and y, which returns three values per call.
-
-
-    NOTE: At this stage the function is assumed to depend on
-    time only, i.e  no spatial dependency!!!!!
-    When that is needed we can use the least_squares interpolation.
-
-    #FIXME: This should work with netcdf (e.g. sww) and thus render the
-    #spatio-temporal boundary condition in shallow water fully general
-
-    #FIXME: Specified quantites not used here -
-    #always return whatever is in the file
-    """
-
-
-    def __init__(self, filename, domain=None, quantities = None, interpolation_points=None):
-        """Initialise callable object from a file.
-        See docstring for class File_function
-
-        If domain is specified, model time (domain,starttime)
-        will be checked and possibly modified
-
-        All times are assumed to be in UTC
-        """
-
-        import time, calendar
-        from Numeric import array
-        from config import time_format
-
-        fid = open(filename)
-        line = fid.readline()
-        fid.close()
-
-        fields = line.split(',')
-        msg = 'File %s must have the format date, value0 value1 value2 ...'
-        msg += ' or date, x, y, value0 value1 value2 ...'
-        mode = len(fields)
-        assert mode in [2,4], msg
-
-        #time.strptime(fields[0], time_format)
-        try:
-            starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:
-            msg = 'First field in file %s must be' %filename
-            msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0]
-            raise msg
-
-
-        #Split values
-        values = []
-        for value in fields[mode-1].split():
-            values.append(float(value))
-
-        q = ensure_numeric(values)
-
-        msg = 'ERROR: File must contain at least one independent value'
-        assert len(q.shape) == 1, msg
-
-        self.number_of_values = len(q)
-
-        self.filename = filename
-        self.starttime = starttime
-	self.domain = domain
-
-        if domain is not None:
-            msg = 'WARNING: Start time as specified in domain (%s)'\
-                  %domain.starttime
-            msg += ' is earlier than the starttime of file %s: %s.'\
-                     %(self.filename, self.starttime)
-            msg += 'Modifying starttime accordingly.'
-            if self.starttime > domain.starttime:
-                #FIXME: Print depending on some verbosity setting
-                ##if verbose: print msg
-                domain.starttime = self.starttime #Modifying model time
-
-            #if domain.starttime is None:
-            #    domain.starttime = self.starttime
-            #else:
-            #    msg = 'WARNING: Start time as specified in domain (%s)'\
-            #          %domain.starttime
-            #    msg += ' is earlier than the starttime of file %s: %s.'\
-            #             %(self.filename, self.starttime)
-            #    msg += 'Modifying starttime accordingly.'
-            #    if self.starttime > domain.starttime:
-            #        #FIXME: Print depending on some verbosity setting
-            #        #print msg
-            #        domain.starttime = self.starttime #Modifying model time
-
-        if mode == 2:
-            self.read_times() #Now read all times in.
-        else:
-            raise 'x,y dependency not yet implemented'
-
-
-    def read_times(self):
-        """Read time and values
-        """
-        from Numeric import zeros, Float, alltrue
-        from config import time_format
-        import time, calendar
-
-        fid = open(self.filename)
-        lines = fid.readlines()
-        fid.close()
-
-        N = len(lines)
-        d = self.number_of_values
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, d), Float)  #Values
-
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-
-            T[i] = realtime - self.starttime
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-        msg = 'File %s must list time as a monotonuosly ' %self.filename
-        msg += 'increasing sequence'
-        assert alltrue( T[1:] - T[:-1] > 0 ), msg
-
-        self.T = T     #Time
-        self.Q = Q     #Values
-        self.index = 0 #Initial index
-
-
-    def __repr__(self):
-        return 'File function'
-
-    def __call__(self, t, x=None, y=None, point_id=None):
-        """Evaluate f(t,x,y)
-
-        FIXME: x, y dependency not yet implemented except that
-        result is a vector of same length as x and y
-        FIXME: Naaaa
-
-	FIXME: Rethink semantics when x,y are vectors.
-        """
-
-        from math import pi, cos, sin, sqrt
-
-
-        #Find time tau such that
-	#
-	# domain.starttime + t == self.starttime + tau
-
-	if self.domain is not None:
-            tau = self.domain.starttime-self.starttime+t
-	else:
-            tau = t
-
-
-        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
-              %(self.filename, self.T[0], self.T[1], tau)
-        if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg
-
-        oldindex = self.index
-
-        #Find slot
-        while tau > self.T[self.index]: self.index += 1
-        while tau < self.T[self.index]: self.index -= 1
-
-        #t is now between index and index+1
-        ratio = (tau - self.T[self.index])/\
-                (self.T[self.index+1] - self.T[self.index])
-
-        #Compute interpolated values
-        q = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
-
-        #Return vector of interpolated values
-        if x == None and y == None:
-            return q
-        else:
-            try:
-                N = len(x)
-            except:
-                return q
-            else:
-                from Numeric import ones, Float
-                #x is a vector - Create one constant column for each value
-                N = len(x)
-                assert len(y) == N, 'x and y must have same length'
-
-                res = []
-                for col in q:
-                    res.append(col*ones(N, Float))
-
-                return res
-
-
-#FIXME: TEMP
-class File_function_copy:
-    """Read time series from file and return a callable object:
-    f(t,x,y)
-    which will return interpolated values based on the input file.
-
-    The file format is assumed to be either two fields separated by a comma:
-
-        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
-
-    or four comma separated fields
-
-        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
-
-    In either case, the callable object will return a tuple of
-    interpolated values, one each value stated in the file.
-
-
-    E.g
-
-      31/08/04 00:00:00, 1.328223 0 0
-      31/08/04 00:15:00, 1.292912 0 0
-
-    will provide a time dependent function f(t,x=None,y=None),
-    ignoring x and y, which returns three values per call.
-
-
-    NOTE: At this stage the function is assumed to depend on
-    time only, i.e  no spatial dependency!!!!!
-    When that is needed we can use the least_squares interpolation.
-
-    #FIXME: This should work with netcdf (e.g. sww) and thus render the
-    #spatio-temporal boundary condition in shallow water fully general
-    """
-
-
-    def __init__(self, filename, domain=None):
-        """Initialise callable object from a file.
-        See docstring for class File_function
-
-        If domain is specified, model time (domain,starttime)
-        will be checked and possibly modified
-
-        All times are assumed to be in UTC
-        """
-
-        import time, calendar
-        from Numeric import array
-        from config import time_format
-
-        assert type(filename) == type(''),\
-               'First argument to File_function must be a string'
-
-
-        try:
-            fid = open(filename)
-        except Exception, e:
-            msg = 'File "%s" could not be opened: Error="%s"'\
-                  %(filename, e)
-            raise msg
-
-
-        line = fid.readline()
-        fid.close()
-        fields = line.split(',')
-        msg = 'File %s must have the format date, value0 value1 value2 ...'
-        msg += ' or date, x, y, value0 value1 value2 ...'
-        mode = len(fields)
-        assert mode in [2,4], msg
-
-        try:
-            starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:
-            msg = 'First field in file %s must be' %filename
-            msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0]
-            raise msg
-
-
-        #Split values
-        values = []
-        for value in fields[mode-1].split():
-            values.append(float(value))
-
-        q = ensure_numeric(values)
-
-        msg = 'ERROR: File must contain at least one independent value'
-        assert len(q.shape) == 1, msg
-
-        self.number_of_values = len(q)
-
-        self.filename = filename
-        self.starttime = starttime
-	self.domain = domain
-
-        if domain is not None:
-            if domain.starttime is None:
-                domain.starttime = self.starttime
-            else:
-                msg = 'WARNING: Start time as specified in domain (%s)'\
-                      %domain.starttime
-                msg += ' is earlier than the starttime of file %s: %s.'\
-                         %(self.filename, self.starttime)
-                msg += 'Modifying starttime accordingly.'
-                if self.starttime > domain.starttime:
-                    #FIXME: Print depending on some verbosity setting
-                    #print msg
-                    domain.starttime = self.starttime #Modifying model time
-
-        if mode == 2:
-            self.read_times() #Now read all times in.
-        else:
-            raise 'x,y dependency not yet implemented'
-
-
-    def read_times(self):
-        """Read time and values
-        """
-        from Numeric import zeros, Float, alltrue
-        from config import time_format
-        import time, calendar
-
-        fid = open(self.filename)
-        lines = fid.readlines()
-        fid.close()
-
-        N = len(lines)
-        d = self.number_of_values
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, d), Float)  #Values
-
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-
-            T[i] = realtime - self.starttime
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-        msg = 'File %s must list time as a monotonuosly ' %self.filename
-        msg += 'increasing sequence'
-        assert alltrue( T[1:] - T[:-1] > 0 ), msg
-
-        self.T = T     #Time
-        self.Q = Q     #Values
-        self.index = 0 #Initial index
-
-
-    def __repr__(self):
-        return 'File function'
-
-
-
-    def __call__(self, t, x=None, y=None):
-        """Evaluate f(t,x,y)
-
-        FIXME: x, y dependency not yet implemented except that
-        result is a vector of same length as x and y
-        FIXME: Naaaa
-        """
-
-        from math import pi, cos, sin, sqrt
-
-
-        #Find time tau such that
-	#
-	# domain.starttime + t == self.starttime + tau
-
-	if self.domain is not None:
-            tau = self.domain.starttime-self.starttime+t
-	else:
-            tau = t
-
-
-        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
-              %(self.filename, self.T[0], self.T[1], tau)
-        if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg
-
-        oldindex = self.index
-
-        #Find slot
-        while tau > self.T[self.index]: self.index += 1
-        while tau < self.T[self.index]: self.index -= 1
-
-        #t is now between index and index+1
-        ratio = (tau - self.T[self.index])/\
-                (self.T[self.index+1] - self.T[self.index])
-
-        #Compute interpolated values
-        q = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
-
-        #Return vector of interpolated values
-        if x == None and y == None:
-            return q
-        else:
-            try:
-                N = len(x)
-            except:
-                return q
-            else:
-                from Numeric import ones, Float
-                #x is a vector - Create one constant column for each value
-                N = len(x)
-                assert len(y) == N, 'x and y must have same length'
-
-                res = []
-                for col in q:
-                    res.append(col*ones(N, Float))
-
-                return res
+    
+
+
 
 
@@ -1052,5 +330,5 @@
     """
     #FIXME: Probably obsoleted by similar function in load_ASCII
-    #FIXE: Name read_data_points (or see 'load' in pylab)
+    #FIXME: Name read_data_points (or see 'load' in pylab)
 
     from Scientific.IO.NetCDF import NetCDFFile
@@ -1512,4 +790,6 @@
     to check if a tsh/msh file 'looks' good.
     """
+
+    #FIXME Move to data_manager
     from shallow_water import Domain
     from pmesh2domain import pmesh_to_domain_instance
@@ -1596,235 +876,2 @@
 
 
-
-
-#OBSOLETED STUFF
-def inside_polygon_old(point, polygon, closed = True, verbose = False):
-    #FIXME Obsoleted
-    """Determine whether points are inside or outside a polygon
-
-    Input:
-       point - Tuple of (x, y) coordinates, or list of tuples
-       polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be
-       regarded as belonging to the polygon (closed = True)
-       or not (closed = False)
-
-    Output:
-       If one point is considered, True or False is returned.
-       If multiple points are passed in, the function returns indices
-       of those points that fall inside the polygon
-
-    Examples:
-       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
-       inside_polygon( [0.5, 0.5], U)
-       will evaluate to True as the point 0.5, 0.5 is inside the unit square
-
-       inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
-       will return the indices [0, 2] as only the first and the last point
-       is inside the unit square
-
-    Remarks:
-       The vertices may be listed clockwise or counterclockwise and
-       the first point may optionally be repeated.
-       Polygons do not need to be convex.
-       Polygons can have holes in them and points inside a hole is
-       regarded as being outside the polygon.
-
-
-    Algorithm is based on work by Darel Finley,
-    http://www.alienryderflex.com/polygon/
-
-    """
-
-    from Numeric import array, Float, reshape
-
-
-    #Input checks
-    try:
-        point = array(point).astype(Float)
-    except:
-        msg = 'Point %s could not be converted to Numeric array' %point
-	raise msg
-
-    try:
-        polygon = array(polygon).astype(Float)
-    except:
-        msg = 'Polygon %s could not be converted to Numeric array' %polygon
-	raise msg
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices: %s' %polygon
-
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns: %s' %polygon
-
-
-    if len(point.shape) == 1:
-        one_point = True
-	points = reshape(point, (1,2))
-    else:
-        one_point = False
-	points = point
-
-    N = polygon.shape[0] #Number of vertices in polygon
-    M = points.shape[0]  #Number of points
-
-    px = polygon[:,0]
-    py = polygon[:,1]
-
-    #Used for an optimisation when points are far away from polygon
-    minpx = min(px); maxpx = max(px)
-    minpy = min(py); maxpy = max(py)
-
-
-    #Begin main loop (FIXME: It'd be crunchy to have this written in C:-)
-    indices = []
-    for k in range(M):
-
-        if verbose:
-            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
-
-        #for each point
-        x = points[k, 0]
-        y = points[k, 1]
-
-        inside = False
-
-        #Optimisation
-        if x > maxpx or x < minpx: continue
-        if y > maxpy or y < minpy: continue
-
-        #Check polygon
-        for i in range(N):
-	    j = (i+1)%N
-
-	    #Check for case where point is contained in line segment
-	    if point_on_line(x, y, px[i], py[i], px[j], py[j]):
-	        if closed:
-		    inside = True
-		else:
-	            inside = False
-		break
-	    else:
- 	        #Check if truly inside polygon
-                if py[i] < y and py[j] >= y or\
-                   py[j] < y and py[i] >= y:
-		    if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
-		        inside = not inside
-
-        if inside: indices.append(k)
-
-    if one_point:
-        return inside
-    else:
-        return indices
-
-
-#def remove_from(A, B):
-#    """Assume that A
-#    """
-#    from Numeric import search_sorted##
-#
-#    ind = search_sorted(A, B)
-
-
-
-def outside_polygon_old(point, polygon, closed = True, verbose = False):
-    #OBSOLETED
-    """Determine whether points are outside a polygon
-
-    Input:
-       point - Tuple of (x, y) coordinates, or list of tuples
-       polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be
-       regarded as belonging to the polygon (closed = True)
-       or not (closed = False)
-
-    Output:
-       If one point is considered, True or False is returned.
-       If multiple points are passed in, the function returns indices
-       of those points that fall outside the polygon
-
-    Examples:
-       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
-       outside_polygon( [0.5, 0.5], U )
-       will evaluate to False as the point 0.5, 0.5 is inside the unit square
-
-       ouside_polygon( [1.5, 0.5], U )
-       will evaluate to True as the point 1.5, 0.5 is outside the unit square
-
-       outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-       will return the indices [1] as only the first and the last point
-       is inside the unit square
-    """
-
-    #FIXME: This is too slow
-
-    res  = inside_polygon(point, polygon, closed, verbose)
-
-    if res is True or res is False:
-        return not res
-
-    #Now invert indices
-    from Numeric import arrayrange, compress
-    outside_indices = arrayrange(len(point))
-    for i in res:
-        outside_indices = compress(outside_indices != i, outside_indices)
-    return outside_indices
-
-def inside_polygon_c(point, polygon, closed = True, verbose = False):
-    #FIXME: Obsolete
-    """Determine whether points are inside or outside a polygon
-
-    C-wrapper
-    """
-
-    from Numeric import array, Float, reshape, zeros, Int
-
-
-    if verbose: print 'Checking input to inside_polygon'
-    #Input checks
-    try:
-        point = array(point).astype(Float)
-    except:
-        msg = 'Point %s could not be converted to Numeric array' %point
-	raise msg
-
-    try:
-        polygon = array(polygon).astype(Float)
-    except:
-        msg = 'Polygon %s could not be converted to Numeric array' %polygon
-	raise msg
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices: %s' %polygon
-
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns: %s' %polygon
-
-
-    if len(point.shape) == 1:
-        one_point = True
-	points = reshape(point, (1,2))
-    else:
-        one_point = False
-	points = point
-
-    from util_ext import inside_polygon
-
-    if verbose: print 'Allocating array for indices'
-
-    indices = zeros( points.shape[0], Int )
-
-    if verbose: print 'Calling C-version of inside poly'
-    count = inside_polygon(points, polygon, indices,
-                           int(closed), int(verbose))
-
-    if one_point:
-        return count == 1 #Return True if the point was inside
-    else:
-        if verbose: print 'Got %d points' %count
-
-        return indices[:count]   #Return those indices that were inside
