"""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. """ 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 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? 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 #Choose format #FIXME: Maybe these can be merged later on if line[:3] == 'CDF': return File_function_NetCDF(filename, domain, quantities, interpolation_points, verbose = verbose) else: return File_function_ASCII(filename, domain, quantities, interpolation_points) 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 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: File must contain at least one independent value' 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 #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 ' Reference:' print ' Lower left corner: [%f, %f]'\ %(self.xllcorner, self.yllcorner) 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 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: #print 'DOMAIN IS NONE!!!!!!!!!' tau = t #print 'D start', self.domain.starttime #print 'S start', self.starttime #print 't', t #print 'tau', tau 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 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 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 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 #Redefine these 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 ############################################## #Initialise module import compile if compile.can_use_C_extension('util_ext.c'): from util_ext import gradient, point_on_line separate_points_by_polygon = separate_points_by_polygon_c else: gradient = gradient_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