source: inundation/ga/storm_surge/pyvolution/util.py @ 1668

"""This module contains various auxiliary function used by pyvolution.

It is also a clearing house for functions that may later earn a module
of their own.
"""

#FIXME: Move polygon stuff out
#FIXME: Finish Interpolation function and get rid of File_function_ASCII


def angle(v):
    """Compute angle between e1 (the unit vector in the x-direction)
    and the specified vector
    """

    from math import acos, pi, sqrt
    from Numeric import sum, array

    l = sqrt( sum (array(v)**2))
    v1 = v[0]/l
    v2 = v[1]/l


    theta = acos(v1)
    
    #try:
    #   theta = acos(v1)
    #except ValueError, e:
    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
    #
    #    #FIXME, hack to avoid acos(1.0) Value error
    #    # why is it happening?
    #    # is it catching something we should avoid?
    #    #FIXME (Ole): When did this happen? We need a unit test to
    #    #reveal this condition or otherwise remove the hack
    #    
    #    s = 1e-6
    #    if (v1+s >  1.0) and (v1-s < 1.0) :
    #        theta = 0.0
    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
    #        theta = 3.1415926535897931
    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
    #          %(v1, theta)

    if v2 < 0:
        #Quadrant 3 or 4
        theta = 2*pi-theta

    return theta


def anglediff(v0, v1):
    """Compute difference between angle of vector x0, y0 and x1, y1.
    This is used for determining the ordering of vertices,
    e.g. for checking if they are counter clockwise.

    Always return a positive value
    """

    from math import pi

    a0 = angle(v0)
    a1 = angle(v1)

    #Ensure that difference will be positive
    if a0 < a1:
        a0 += 2*pi

    return a0-a1


def mean(x):
    from Numeric import sum
    return sum(x)/len(x)


def point_on_line(x, y, x0, y0, x1, y1):
    """Determine whether a point is on a line segment

    Input: x, y, x0, x0, x1, y1: where
        point is given by x, y
        line is given by (x0, y0) and (x1, y1)

    """

    from Numeric import array, dot, allclose
    from math import sqrt

    a = array([x - x0, y - y0])
    a_normal = array([a[1], -a[0]])

    b = array([x1 - x0, y1 - y0])

    if dot(a_normal, b) == 0:
        #Point is somewhere on the infinite extension of the line

        len_a = sqrt(sum(a**2))
        len_b = sqrt(sum(b**2))
        if dot(a, b) >= 0 and len_a <= len_b:
           return True
        else:
           return False
    else:
      return False

def ensure_numeric(A, typecode = None):
    """Ensure that sequence is a Numeric array.
    Inputs:
        A: Sequance. If A is already a Numeric array it will be returned
                     unaltered
                     If not, an attempt is made to convert it to a Numeric
                     array
        typecode: Numeric type. If specified, use this in the conversion.
                                If not, let Numeric decide

    This function is necessary as array(A) can cause memory overflow.
    """

    from Numeric import ArrayType, array

    if typecode is None:
        if type(A) == ArrayType:
            return A
        else:
            return array(A)
    else:
        if type(A) == ArrayType:
            if A.typecode == typecode:
                return array(A)
            else:
                return A.astype(typecode)
        else:
            return array(A).astype(typecode)



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
    """


    #FIXME (OLE): Should check origin of domain against that of file
    #In fact, this is where origin should be converted to that of domain
    #Also, check that file covers domain fully.

    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()

    if domain is not None:
        quantities = domain.conserved_quantities
    else:
        quantities = None


    if line[:3] == 'CDF':
        return get_netcdf_file_function(filename, domain, quantities,
                                        interpolation_points, verbose = verbose)
    else:
        #FIXME Should be phased out
        return File_function_ASCII(filename, domain, quantities,
                                   interpolation_points)
    
    








def get_netcdf_file_function(filename, domain=None, quantity_names=None,
                             interpolation_points=None, verbose = False):
    """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.

    If domain is specified, model time (domain.starttime)
    will be checked and possibly modified
    
    All times are assumed to be in UTC

    See Interpolation function for further documetation
    

    """
    
    #verbose = True
    
    #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
    from Numeric import array, zeros, Float, alltrue, concatenate, reshape        

    #Open NetCDF file
    if verbose: print 'Reading', filename
    fid = NetCDFFile(filename, 'r')

    if quantity_names is None or len(quantity_names) < 1:
        msg = 'ERROR: At least one independent value must be specified'
        raise msg

    missing = []
    for quantity in ['x', 'y', 'time'] + quantity_names:
        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:
        starttime = fid.starttime[0]
    except ValueError:
        msg = 'Could not read starttime from file %s' %filename
        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 domain is not None:
        #Update domain.startime if it is *earlier* than starttime
        if starttime > domain.starttime:
            msg = 'WARNING: Start time as specified in domain (%f)'\
                  %domain.starttime
            msg += ' is earlier than the starttime of file %s (%f).'\
                     %(filename, starttime)
            msg += ' Modifying domain starttime accordingly.'
            
            if verbose: print msg
            domain.starttime = starttime #Modifying model time
            if verbose: print 'Domain starttime is now set to %f'\
               %domain.starttime


        #If domain.startime is *later* than starttime,
        #move time back - relative to domain's time
        if domain.starttime > starttime:
            time = time - domain.starttime + starttime


    if verbose:
        print 'File_function data obtained from: %s' %filename
        print '  References:'
        #print '    Datum ....' #FIXME
        print '    Lower left corner: [%f, %f]'\
              %(xllcorner, yllcorner)
        print '    Start time:   %f' %starttime                
        
    
    #Produce values for desired data points at
    #each timestep

    quantities = {}
    for i, name in enumerate(quantity_names):
        quantities[name] = fid.variables[name][:]
    fid.close()
    

    return Interpolation_function(vertex_coordinates,
                                  triangles,
                                  time,
                                  quantities,
                                  quantity_names,  
                                  interpolation_points,
                                  alpha = 0,
                                  verbose = verbose)





class Interpolation_function:
    """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
    which is interpolated from time series defined at vertices of
    triangular mesh 

    Input
        vertex_coordinates: mx2 array of coordinates (Float)
        triangles:          nx3 array of indices into vertex_coordinates (Int)
        time:               px1 array of times (Float)               
        quantities:         Dictionary of pxm values
        quantity_names:
        interpolation_points:
        alpha:
        verbose:
    

    
    The quantities returned by the callable object 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']

    The parameter interpolation_points decides at which points interpolated
    quantities are to be computed whenever object is called.
    If None, return average value
    """

    

    
    def __init__(self,
                 vertex_coordinates,
                 triangles,
                 time,
                 quantities,
                 quantity_names = None,  
                 interpolation_points = None,
                 alpha = 0,
                 verbose = False):
        """
        """

        from Numeric import array, zeros, Float, alltrue, concatenate,\
             reshape, ArrayType

        from config import time_format
        from least_squares import Interpolation

        #Check
        msg = 'Time must be a monotonuosly '
        msg += 'increasing sequence %s' %time
        assert alltrue(time[1:] - time[:-1] > 0 ), msg


        #Check if quantities is a single array only
        if type(quantities) == ArrayType:
            quantity_names = 'Attribute'
            quantities = {quantity_names: quantities}            

        #Use keys if no names specified
        if quantity_names is not None:
            self.quantity_names = quantity_names
        else:
            self.quantity_names = quantities.keys()


        self.interpolation_points = interpolation_points
        self.T = time[:]  #Time assumed to be relative to starttime
        self.index = 0    #Initial time index
        self.values = {}
            

        if interpolation_points is not None:
            #Spatial interpolation

            try:
                interpolation_points = 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


            for name in quantity_names:
                self.values[name] = zeros( (len(self.T),
                                            len(interpolation_points)),
                                            Float)

            #Build interpolator
            interpol = Interpolation(vertex_coordinates,
                                     triangles,
                                     point_coordinates = interpolation_points,
                                     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 quantity_names:
                    self.values[name][i, :] =\
                    interpol.interpolate(quantities[name][i,:])

            #Report
            if verbose:
                x = vertex_coordinates[:,0]
                y = vertex_coordinates[:,1]                
                
                print '------------------------------------------------'
                print 'Interpolation_function statistics:'
                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 quantity_names:
                    q = quantities[name][:].flat
                    print '    %s in [%f, %f]' %(name, min(q), max(q))


                print '  Interpolation points (xi, eta):'\
                      'number of points == %d ' %interpolation_points.shape[0]
                print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
                                             max(interpolation_points[:,0]))
                print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
                                              max(interpolation_points[:,1]))
                print '  Interpolated quantities (over all timesteps):'
                for name in quantity_names:
                    q = self.values[name][:].flat
                    print '    %s at interpolation points in [%f, %f]'\
                          %(name, min(q), max(q))
                print '------------------------------------------------'
        else:
            #Return an average, making this a time series
            #FIXME: Needs a test
            for name in quantity_names:
                self.values[name] = avg(quantities[name][i,:])


    def __repr__(self):
        return 'Interpolation function'

    def __call__(self, t, point_id = None, x = None, y = None):
        """Evaluate f(t), f(t, point_id) or f(t, x, y)

        Inputs:
          t: time - Model time. Must lie within existing timesteps
          point_id: index of one of the preprocessed points.
                    If point_id is None all preprocessed points are computed
                    
          If no x,y info is present, point_id and x,y arguments are ignored
          making f a function of time only.
          
          FIXME: point_id could also be a slice
          FIXME: One could allow arbitrary x, y coordinates
          (requires computation of a new interpolator)
          FIXME: What if x and y are vectors?
          FIXME: Does point_id or x,y take precedence?
          FIXME: What about f(x,y) without t
          
        """

        from math import pi, cos, sin, sqrt
        from Numeric import zeros, Float


        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
        msg += ' does not match model time: %s\n' %t
        if t < self.T[0]: raise msg
        if t > self.T[-1]: raise msg

        oldindex = self.index #Time index

        #Find current time slot
        while t > self.T[self.index]: self.index += 1
        while t < self.T[self.index]: self.index -= 1

        if t == self.T[self.index]:
            #Protect against case where t == T[-1] (last time)
            # - also works in general when t == T[i]
            ratio = 0
        else:
            #t is now between index and index+1
            ratio = (t - self.T[self.index])/\
                    (self.T[self.index+1] - self.T[self.index])

        #Compute interpolated values
        q = zeros(len(self.quantity_names), Float)

        for i, name in enumerate(self.quantity_names):
            Q = self.values[name]

            if point_id is None:
                Q0 = Q[self.index]

                try:
                    Q1 = Q[self.index+1]
                except:
                    pass                
            else:    
                Q0 = Q[self.index, point_id]

                try:
                    Q1 = Q[self.index+1, point_id]
                except:
                    pass
                

            #Linear temporal interpolation    
            if ratio > 0:
                q[i] = Q0 + ratio*(Q1 - Q0)
            else:
                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
        return q










#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


def read_xya(filename, format = 'netcdf'):
    """Read simple xya file, possibly with a header in the first line, just
    x y [attributes]
    separated by whitespace

    Format can be either ASCII or NetCDF

    Return list of points, list of attributes and
    attribute names if present otherwise None
    """
    #FIXME: Probably obsoleted by similar function in load_ASCII
    #FIXE: Name read_data_points (or see 'load' in pylab)

    from Scientific.IO.NetCDF import NetCDFFile

    if format.lower() == 'netcdf':
        #Get NetCDF

        fid = NetCDFFile(filename, 'r')

        # Get the variables
        points = fid.variables['points']
        keys = fid.variables.keys()
        attributes = {}
        for key in keys:
            attributes[key] = fid.variables[key]
        #Don't close - arrays are needed outside this function,
        #alternatively take a copy here
    else:
        #Read as ASCII file assuming that it is separated by whitespace
        fid = open(filename)
        lines = fid.readlines()
        fid.close()

        #Check if there is a header line
        fields = lines[0].strip().split()
        try:
            float(fields[0])
        except:
            #This must be a header line
            attribute_names = fields
            lines = lines[1:]
        else:
            attribute_names = ['elevation']  #HACK

        attributes = {}
        for key in attribute_names:
            attributes[key] = []

        points = []
        for line in lines:
            fields = line.strip().split()
            points.append( (float(fields[0]), float(fields[1])) )
            for i, key in enumerate(attribute_names):
                attributes[key].append( float(fields[2+i]) )

    return points, attributes


#####################################
#POLYGON STUFF
#
#FIXME: All these should be put into new module polygon.py



#These have been redefined to use separate_by_polygon which will
#put all inside indices in the first part of the indices array and
#outside indices in the last

def inside_polygon(points, polygon, closed = True, verbose = False):
    """See separate_points_by_polygon for documentation
    """

    from Numeric import array, Float, reshape

    if verbose: print 'Checking input to inside_polygon'
    try:
        points = ensure_numeric(points, Float)
    except:
        msg = 'Points could not be converted to Numeric array'
        raise msg

    try:
        polygon = ensure_numeric(polygon, Float)
    except:
        msg = 'Polygon could not be converted to Numeric array'
        raise msg



    if len(points.shape) == 1:
        one_point = True
        points = reshape(points, (1,2))
    else:
        one_point = False

    indices, count = separate_points_by_polygon(points, polygon,
                                                closed, verbose)

    if one_point:
        return count == 1
    else:
        return indices[:count]

def outside_polygon(points, polygon, closed = True, verbose = False):
    """See separate_points_by_polygon for documentation
    """

    from Numeric import array, Float, reshape

    if verbose: print 'Checking input to outside_polygon'
    try:
        points = ensure_numeric(points, Float)
    except:
        msg = 'Points could not be converted to Numeric array'
        raise msg

    try:
        polygon = ensure_numeric(polygon, Float)
    except:
        msg = 'Polygon could not be converted to Numeric array'
        raise msg



    if len(points.shape) == 1:
        one_point = True
        points = reshape(points, (1,2))
    else:
        one_point = False

    indices, count = separate_points_by_polygon(points, polygon,
                                                closed, verbose)

    if one_point:
        return count != 1
    else:
        return indices[count:][::-1]  #return reversed


def separate_points_by_polygon(points, polygon,
                               closed = True, verbose = False):
    """Determine whether points are inside or outside a polygon

    Input:
       points - 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)

    Outputs:
       indices: array of same length as points with indices of points falling
       inside the polygon listed from the beginning and indices of points
       falling outside listed from the end.

       count: count of points falling inside the polygon

       The indices of points inside are obtained as indices[:count]
       The indices of points outside are obtained as indices[count:]


    Examples:
       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square

       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
       will return the indices [0, 2, 1] and count == 2 as only the first
       and the last point are 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, Int, zeros


    #Input checks
    try:
        points = ensure_numeric(points, Float)
    except:
        msg = 'Points could not be converted to Numeric array'
        raise msg

    try:
        polygon = ensure_numeric(polygon, Float)
    except:
        msg = 'Polygon could not be converted to Numeric array'
        raise msg

    assert len(polygon.shape) == 2,\
       'Polygon array must be a 2d array of vertices'

    assert polygon.shape[1] == 2,\
       'Polygon array must have two columns'

    assert len(points.shape) == 2,\
       'Points array must be a 2d array'

    assert points.shape[1] == 2,\
       'Points array must have two columns'

    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
    indices = zeros(M, Int)

    inside_index = 0    #Keep track of points inside
    outside_index = M-1 #Keep track of points outside (starting from end)

    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

        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
            #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[inside_index] = k
            inside_index += 1
        else:
            indices[outside_index] = k
            outside_index -= 1

    return indices, inside_index


def separate_points_by_polygon_c(points, polygon,
                                 closed = True, verbose = False):
    """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 separate_points_by_polygon'
    #Input checks
    try:
        points = ensure_numeric(points, Float)
    except:
        msg = 'Points could not be converted to Numeric array'
        raise msg

    #if verbose: print 'Checking input to separate_points_by_polygon 2'
    try:
        polygon = ensure_numeric(polygon, Float)
    except:
        msg = 'Polygon could not be converted to Numeric array'
        raise msg

    if verbose: print 'check'

    assert len(polygon.shape) == 2,\
       'Polygon array must be a 2d array of vertices'

    assert polygon.shape[1] == 2,\
       'Polygon array must have two columns'

    assert len(points.shape) == 2,\
       'Points array must be a 2d array'

    assert points.shape[1] == 2,\
       'Points array must have two columns'

    N = polygon.shape[0] #Number of vertices in polygon
    M = points.shape[0]  #Number of points

    from util_ext import separate_points_by_polygon

    if verbose: print 'Allocating array for indices'

    indices = zeros( M, Int )

    if verbose: print 'Calling C-version of inside poly'
    count = separate_points_by_polygon(points, polygon, indices,
                                       int(closed), int(verbose))

    return indices, count



class Polygon_function:
    """Create callable object f: x,y -> z, where a,y,z are vectors and
    where f will return different values depending on whether x,y belongs
    to specified polygons.

    To instantiate:

       Polygon_function(polygons)

    where polygons is a list of tuples of the form

      [ (P0, v0), (P1, v1), ...]

      with Pi being lists of vertices defining polygons and vi either
      constants or functions of x,y to be applied to points with the polygon.

    The function takes an optional argument, default which is the value
    (or function) to used for points not belonging to any polygon.
    For example:

       Polygon_function(polygons, default = 0.03)

    If omitted the default value will be 0.0

    Note: If two polygons overlap, the one last in the list takes precedence

    FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference???

    """

    def __init__(self, regions, default = 0.0):

        try:
            len(regions)
        except:
            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
            raise msg


        T = regions[0]
        try:
            a = len(T)
        except:
            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
            raise msg

        assert a == 2, 'Must have two component each: %s' %T

        self.regions = regions
        self.default = default


    def __call__(self, x, y):
        from util import inside_polygon
        from Numeric import ones, Float, concatenate, array, reshape, choose

        x = array(x).astype(Float)
        y = array(y).astype(Float)

        N = len(x)
        assert len(y) == N

        points = concatenate( (reshape(x, (N, 1)),
                               reshape(y, (N, 1))), axis=1 )

        if callable(self.default):
            z = self.default(x,y)
        else:
            z = ones(N, Float) * self.default

        for polygon, value in self.regions:
            indices = inside_polygon(points, polygon)

            #FIXME: This needs to be vectorised
            if callable(value):
                for i in indices:
                    xx = array([x[i]])
                    yy = array([y[i]])
                    z[i] = value(xx, yy)[0]
            else:
                for i in indices:
                    z[i] = value

        return z

def read_polygon(filename):
    """Read points assumed to form a polygon
       There must be exactly two numbers in each line
    """

    #Get polygon
    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    polygon = []
    for line in lines:
        fields = line.split(',')
        polygon.append( [float(fields[0]), float(fields[1])] )

    return polygon

def populate_polygon(polygon, number_of_points, seed = None):
    """Populate given polygon with uniformly distributed points.

    Input:
       polygon - list of vertices of polygon
       number_of_points - (optional) number of points
       seed - seed for random number generator (default=None)

    Output:
       points - list of points inside polygon

    Examples:
       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
       will return five randomly selected points inside the unit square
    """

    from random import uniform, seed

    seed(seed)

    points = []

    #Find outer extent of polygon
    max_x = min_x = polygon[0][0]
    max_y = min_y = polygon[0][1]
    for point in polygon[1:]:
        x = point[0]
        if x > max_x: max_x = x
        if x < min_x: min_x = x
        y = point[1]
        if y > max_y: max_y = y
        if y < min_y: min_y = y


    while len(points) < number_of_points:
        x = uniform(min_x, max_x)
        y = uniform(min_y, max_y)

        if inside_polygon( [x,y], polygon ):
            points.append([x,y])

    return points

def tsh2sww(filename, verbose=True):
    """
    to check if a tsh/msh file 'looks' good.
    """
    from shallow_water import Domain
    from pmesh2domain import pmesh_to_domain_instance
    import time, os
    from data_manager import get_dataobject
    from os import sep, path
    #from util import mean

    if verbose == True:print 'Creating domain from', filename
    domain = pmesh_to_domain_instance(filename, Domain)
    if verbose == True:print "Number of triangles = ", len(domain)

    domain.smooth = True
    domain.format = 'sww'   #Native netcdf visualisation format
    file_path, filename = path.split(filename)
    filename, ext = path.splitext(filename)
    domain.filename = filename
    domain.reduction = mean
    if verbose == True:print "file_path",file_path
    if file_path == "":file_path = "."
    domain.set_datadir(file_path)

    if verbose == True:
        print "Output written to " + domain.get_datadir() + sep + \
              domain.filename + "." + domain.format
    sww = get_dataobject(domain)
    sww.store_connectivity()
    sww.store_timestep('stage')

####################################################################
#Python versions of function that are also implemented in util_gateway.c
#

def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
    """
    """

    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
    a /= det

    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
    b /= det

    return a, b


def gradient2_python(x0, y0, x1, y1, q0, q1):
    """Compute radient based on two points and enforce zero gradient
    in the direction orthogonal to (x1-x0), (y1-y0) 
    """

    #Old code
    #det = x0*y1 - x1*y0
    #if det != 0.0:
    #    a = (y1*q0 - y0*q1)/det
    #    b = (x0*q1 - x1*q0)/det

    #Correct code (ON)
    det = (x1-x0)**2 + (y1-y0)**2
    if det != 0.0:
        a = (x1-x0)*(q1-q0)/det
        b = (y1-y0)*(q1-q0)/det
        
    return a, b        


##############################################
#Initialise module

import compile
if compile.can_use_C_extension('util_ext.c'):
    from util_ext import gradient, gradient2, point_on_line
    separate_points_by_polygon = separate_points_by_polygon_c
else:
    gradient = gradient_python
    gradient2 = gradient2_python    


if __name__ == "__main__":
    pass






#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