Context Navigation

← Previous Change
Next Change →

util.py

Timestamp:

Aug 2, 2005, 4:48:56 PM (20 years ago)

Author:

ole

Message:

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

File:

: 1 edited

inundation/ga/storm_surge/pyvolution/util.py (modified) (13 diffs)

Legend:

: Unmodified
: Added
: Removed

inundation/ga/storm_surge/pyvolution/util.py

-                      r1670
+                      r1671
 #FIXME: Move polygon stuff out
-#FIXME: Finish Interpolation function and get rid of File_function_ASCII
 def angle(v):
 …
+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)
     """
 …
     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
 …
                                         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'
 …
     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)
 …
         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
 …
         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:
 …
         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
 …
     from least_squares import Interpolation_function
+    if not spatial:
+        vertex_coordinates = triangles = interpolation_points = None
     return Interpolation_function(time,
                                   quantities,
 …
                                   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
+/08/04 00:00:00, 1.328223 0 0
+/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
+/08/04 00:00:00, 1.328223 0 0
+/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
 …
     """
     #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
 …
     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
 …
-#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

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 1671 for inundation/ga/storm_surge/pyvolution/util.py

Legend:

inundation/ga/storm_surge/pyvolution/util.py

Download in other formats: