source: inundation/pyvolution/util.py @ 1740

Last change on this file since 1740 was 1681, checked in by ole, 19 years ago

Played with set_quantity

File size: 24.8 KB
Line 
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7#FIXME: Move polygon stuff out
8
9def angle(v):
10    """Compute angle between e1 (the unit vector in the x-direction)
11    and the specified vector
12    """
13
14    from math import acos, pi, sqrt
15    from Numeric import sum, array
16
17    l = sqrt( sum (array(v)**2))
18    v1 = v[0]/l
19    v2 = v[1]/l
20
21
22    theta = acos(v1)
23   
24    #try:
25    #   theta = acos(v1)
26    #except ValueError, e:
27    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
28    #
29    #    #FIXME, hack to avoid acos(1.0) Value error
30    #    # why is it happening?
31    #    # is it catching something we should avoid?
32    #    #FIXME (Ole): When did this happen? We need a unit test to
33    #    #reveal this condition or otherwise remove the hack
34    #   
35    #    s = 1e-6
36    #    if (v1+s >  1.0) and (v1-s < 1.0) :
37    #        theta = 0.0
38    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
39    #        theta = 3.1415926535897931
40    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
41    #          %(v1, theta)
42
43    if v2 < 0:
44        #Quadrant 3 or 4
45        theta = 2*pi-theta
46
47    return theta
48
49
50def anglediff(v0, v1):
51    """Compute difference between angle of vector x0, y0 and x1, y1.
52    This is used for determining the ordering of vertices,
53    e.g. for checking if they are counter clockwise.
54
55    Always return a positive value
56    """
57
58    from math import pi
59
60    a0 = angle(v0)
61    a1 = angle(v1)
62
63    #Ensure that difference will be positive
64    if a0 < a1:
65        a0 += 2*pi
66
67    return a0-a1
68
69
70def mean(x):
71    from Numeric import sum
72    return sum(x)/len(x)
73
74
75def point_on_line(x, y, x0, y0, x1, y1):
76    """Determine whether a point is on a line segment
77
78    Input: x, y, x0, x0, x1, y1: where
79        point is given by x, y
80        line is given by (x0, y0) and (x1, y1)
81
82    """
83
84    from Numeric import array, dot, allclose
85    from math import sqrt
86
87    a = array([x - x0, y - y0])
88    a_normal = array([a[1], -a[0]])
89
90    b = array([x1 - x0, y1 - y0])
91
92    if dot(a_normal, b) == 0:
93        #Point is somewhere on the infinite extension of the line
94
95        len_a = sqrt(sum(a**2))
96        len_b = sqrt(sum(b**2))
97        if dot(a, b) >= 0 and len_a <= len_b:
98           return True
99        else:
100           return False
101    else:
102      return False
103
104def ensure_numeric(A, typecode = None):
105    """Ensure that sequence is a Numeric array.
106    Inputs:
107        A: Sequance. If A is already a Numeric array it will be returned
108                     unaltered
109                     If not, an attempt is made to convert it to a Numeric
110                     array
111        typecode: Numeric type. If specified, use this in the conversion.
112                                If not, let Numeric decide
113
114    This function is necessary as array(A) can cause memory overflow.
115    """
116
117    from Numeric import ArrayType, array
118
119    if typecode is None:
120        if type(A) == ArrayType:
121            return A
122        else:
123            return array(A)
124    else:
125        if type(A) == ArrayType:
126            if A.typecode == typecode:
127                return array(A)
128            else:
129                return A.astype(typecode)
130        else:
131            return array(A).astype(typecode)
132
133
134
135def file_function(filename,
136                  domain = None,
137                  quantities = None,
138                  interpolation_points = None, verbose = False):
139    """If quantities is not specified, derive them from domain
140    (if that is specified)
141
142
143    see get_netcdf_file_function
144    """
145
146
147    #FIXME (OLE): Should check origin of domain against that of file
148    #In fact, this is where origin should be converted to that of domain
149    #Also, check that file covers domain fully.
150
151    assert type(filename) == type(''),\
152               'First argument to File_function must be a string'
153
154    try:
155        fid = open(filename)
156    except Exception, e:
157        msg = 'File "%s" could not be opened: Error="%s"'\
158                  %(filename, e)
159        raise msg
160
161    line = fid.readline()
162    fid.close()
163
164    if quantities is None: 
165        if domain is not None:
166            quantities = domain.conserved_quantities
167
168
169
170    if line[:3] == 'CDF':
171        return get_netcdf_file_function(filename, domain, quantities,
172                                        interpolation_points, verbose = verbose)
173    else:
174        raise 'Must be a NetCDF File'
175
176
177
178def get_netcdf_file_function(filename, domain=None, quantity_names=None,
179                             interpolation_points=None, verbose = False):
180    """Read time history of spatial data from NetCDF sww file and
181    return a callable object f(t,x,y)
182    which will return interpolated values based on the input file.
183
184    If domain is specified, model time (domain.starttime)
185    will be checked and possibly modified
186   
187    All times are assumed to be in UTC
188
189    See Interpolation function for further documetation
190   
191
192    """
193   
194    #verbose = True
195   
196    #FIXME: Check that model origin is the same as file's origin
197    #(both in UTM coordinates)
198    #If not - modify those from file to match domain
199    #Take this code from e.g. dem2pts in data_manager.py
200    #FIXME: Use geo_reference to read and write xllcorner...
201       
202
203    import time, calendar, types
204    from config import time_format
205    from Scientific.IO.NetCDF import NetCDFFile
206    from Numeric import array, zeros, Float, alltrue, concatenate, reshape       
207
208    #Open NetCDF file
209    if verbose: print 'Reading', filename
210    fid = NetCDFFile(filename, 'r')
211
212    if type(quantity_names) == types.StringType:
213        quantity_names = [quantity_names]       
214   
215    if quantity_names is None or len(quantity_names) < 1:
216        #Get quantities from file
217        quantity_names = []
218        for name in fid.variables.keys():
219            if name not in ['x', 'y', 'time']:
220                quantity_names.append(name)
221               
222        #msg = 'ERROR: At least one independent value must be specified'
223        #raise msg
224
225    missing = []
226       
227    for quantity in ['time'] + quantity_names:
228        if not fid.variables.has_key(quantity):
229            missing.append(quantity)
230
231    if len(missing) > 0:
232        msg = 'Quantities %s could not be found in file %s'\
233              %(str(missing), filename)
234        fid.close()
235        raise msg
236
237    #Decide whether this data has a spatial dimension
238    #FIXME: Let us split this up into two functions: One that reads
239    #sww and one that reads tms (time series only)
240    spatial = True
241    for quantity in ['x', 'y']:
242        if not fid.variables.has_key(quantity):
243            spatial = False
244
245
246
247    #Get first timestep
248    try:
249        starttime = fid.starttime[0]
250    except ValueError:
251        msg = 'Could not read starttime from file %s' %filename
252        raise msg
253
254    #Get variables
255    if verbose: print 'Get variables'   
256    time = fid.variables['time'][:]   
257
258    if spatial:
259        #Get origin
260        xllcorner = fid.xllcorner[0]
261        yllcorner = fid.yllcorner[0]
262
263        x = fid.variables['x'][:]
264        y = fid.variables['y'][:]
265        triangles = fid.variables['volumes'][:]
266
267        x = reshape(x, (len(x),1))
268        y = reshape(y, (len(y),1))
269        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
270
271
272    if domain is not None:
273        #Update domain.startime if it is *earlier* than starttime
274        if starttime > domain.starttime:
275            msg = 'WARNING: Start time as specified in domain (%f)'\
276                  %domain.starttime
277            msg += ' is earlier than the starttime of file %s (%f).'\
278                     %(filename, starttime)
279            msg += ' Modifying domain starttime accordingly.'
280           
281            if verbose: print msg
282            domain.starttime = starttime #Modifying model time
283            if verbose: print 'Domain starttime is now set to %f'\
284               %domain.starttime
285
286
287        #If domain.startime is *later* than starttime,
288        #move time back - relative to domain's time
289        if domain.starttime > starttime:
290            time = time - domain.starttime + starttime
291
292
293    if verbose:
294        print 'File_function data obtained from: %s' %filename
295        print '  References:'
296        #print '    Datum ....' #FIXME
297        if spatial:
298            print '    Lower left corner: [%f, %f]'\
299                  %(xllcorner, yllcorner)
300        print '    Start time:   %f' %starttime               
301       
302   
303    #Produce values for desired data points at
304    #each timestep
305
306    quantities = {}
307    for i, name in enumerate(quantity_names):
308        quantities[name] = fid.variables[name][:]
309    fid.close()
310   
311
312    from least_squares import Interpolation_function
313
314    if not spatial:
315        vertex_coordinates = triangles = interpolation_points = None         
316       
317    return Interpolation_function(time,
318                                  quantities,
319                                  quantity_names,
320                                  vertex_coordinates,
321                                  triangles,
322                                  interpolation_points,
323                                  verbose = verbose)
324   
325
326
327
328
329def read_xya(filename, format = 'netcdf'):
330    """Read simple xya file, possibly with a header in the first line, just
331    x y [attributes]
332    separated by whitespace
333
334    Format can be either ASCII or NetCDF
335
336    Return list of points, list of attributes and
337    attribute names if present otherwise None
338    """
339    #FIXME: Probably obsoleted by similar function in load_ASCII
340    #FIXME: Name read_data_points (or see 'load' in pylab)
341
342    from Scientific.IO.NetCDF import NetCDFFile
343
344    if format.lower() == 'netcdf':
345        #Get NetCDF
346
347        fid = NetCDFFile(filename, 'r')
348
349        # Get the variables
350        points = fid.variables['points']
351        keys = fid.variables.keys()
352        attributes = {}
353        for key in keys:
354            if key != 'points':
355                attributes[key] = fid.variables[key]
356        #Don't close - arrays are needed outside this function,
357        #alternatively take a copy here
358    else:
359        #Read as ASCII file assuming that it is separated by whitespace
360        fid = open(filename)
361        lines = fid.readlines()
362        fid.close()
363
364        #Check if there is a header line
365        fields = lines[0].strip().split()
366        try:
367            float(fields[0])
368        except:
369            #This must be a header line
370            attribute_names = fields
371            lines = lines[1:]
372        else:
373            attribute_names = ['elevation']  #HACK
374
375        attributes = {}
376        for key in attribute_names:
377            attributes[key] = []
378
379        points = []
380        for line in lines:
381            fields = line.strip().split()
382            points.append( (float(fields[0]), float(fields[1])) )
383            for i, key in enumerate(attribute_names):
384                attributes[key].append( float(fields[2+i]) )
385
386    return points, attributes
387
388
389#####################################
390#POLYGON STUFF
391#
392#FIXME: All these should be put into new module polygon.py
393
394
395
396#These have been redefined to use separate_by_polygon which will
397#put all inside indices in the first part of the indices array and
398#outside indices in the last
399
400def inside_polygon(points, polygon, closed = True, verbose = False):
401    """See separate_points_by_polygon for documentation
402    """
403
404    from Numeric import array, Float, reshape
405
406    if verbose: print 'Checking input to inside_polygon'
407    try:
408        points = ensure_numeric(points, Float)
409    except:
410        msg = 'Points could not be converted to Numeric array'
411        raise msg
412
413    try:
414        polygon = ensure_numeric(polygon, Float)
415    except:
416        msg = 'Polygon could not be converted to Numeric array'
417        raise msg
418
419
420
421    if len(points.shape) == 1:
422        one_point = True
423        points = reshape(points, (1,2))
424    else:
425        one_point = False
426
427    indices, count = separate_points_by_polygon(points, polygon,
428                                                closed, verbose)
429
430    if one_point:
431        return count == 1
432    else:
433        return indices[:count]
434
435def outside_polygon(points, polygon, closed = True, verbose = False):
436    """See separate_points_by_polygon for documentation
437    """
438
439    from Numeric import array, Float, reshape
440
441    if verbose: print 'Checking input to outside_polygon'
442    try:
443        points = ensure_numeric(points, Float)
444    except:
445        msg = 'Points could not be converted to Numeric array'
446        raise msg
447
448    try:
449        polygon = ensure_numeric(polygon, Float)
450    except:
451        msg = 'Polygon could not be converted to Numeric array'
452        raise msg
453
454
455
456    if len(points.shape) == 1:
457        one_point = True
458        points = reshape(points, (1,2))
459    else:
460        one_point = False
461
462    indices, count = separate_points_by_polygon(points, polygon,
463                                                closed, verbose)
464
465    if one_point:
466        return count != 1
467    else:
468        return indices[count:][::-1]  #return reversed
469
470
471def separate_points_by_polygon(points, polygon,
472                               closed = True, verbose = False):
473    """Determine whether points are inside or outside a polygon
474
475    Input:
476       points - Tuple of (x, y) coordinates, or list of tuples
477       polygon - list of vertices of polygon
478       closed - (optional) determine whether points on boundary should be
479       regarded as belonging to the polygon (closed = True)
480       or not (closed = False)
481
482    Outputs:
483       indices: array of same length as points with indices of points falling
484       inside the polygon listed from the beginning and indices of points
485       falling outside listed from the end.
486
487       count: count of points falling inside the polygon
488
489       The indices of points inside are obtained as indices[:count]
490       The indices of points outside are obtained as indices[count:]
491
492
493    Examples:
494       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
495
496       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
497       will return the indices [0, 2, 1] and count == 2 as only the first
498       and the last point are inside the unit square
499
500    Remarks:
501       The vertices may be listed clockwise or counterclockwise and
502       the first point may optionally be repeated.
503       Polygons do not need to be convex.
504       Polygons can have holes in them and points inside a hole is
505       regarded as being outside the polygon.
506
507    Algorithm is based on work by Darel Finley,
508    http://www.alienryderflex.com/polygon/
509
510    """
511
512    from Numeric import array, Float, reshape, Int, zeros
513
514
515    #Input checks
516    try:
517        points = ensure_numeric(points, Float)
518    except:
519        msg = 'Points could not be converted to Numeric array'
520        raise msg
521
522    try:
523        polygon = ensure_numeric(polygon, Float)
524    except:
525        msg = 'Polygon could not be converted to Numeric array'
526        raise msg
527
528    assert len(polygon.shape) == 2,\
529       'Polygon array must be a 2d array of vertices'
530
531    assert polygon.shape[1] == 2,\
532       'Polygon array must have two columns'
533
534    assert len(points.shape) == 2,\
535       'Points array must be a 2d array'
536
537    assert points.shape[1] == 2,\
538       'Points array must have two columns'
539
540    N = polygon.shape[0] #Number of vertices in polygon
541    M = points.shape[0]  #Number of points
542
543    px = polygon[:,0]
544    py = polygon[:,1]
545
546    #Used for an optimisation when points are far away from polygon
547    minpx = min(px); maxpx = max(px)
548    minpy = min(py); maxpy = max(py)
549
550
551    #Begin main loop
552    indices = zeros(M, Int)
553
554    inside_index = 0    #Keep track of points inside
555    outside_index = M-1 #Keep track of points outside (starting from end)
556
557    for k in range(M):
558
559        if verbose:
560            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
561
562        #for each point
563        x = points[k, 0]
564        y = points[k, 1]
565
566        inside = False
567
568        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
569            #Check polygon
570            for i in range(N):
571                j = (i+1)%N
572
573                #Check for case where point is contained in line segment
574                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
575                    if closed:
576                        inside = True
577                    else:
578                        inside = False
579                    break
580                else:
581                    #Check if truly inside polygon
582                    if py[i] < y and py[j] >= y or\
583                       py[j] < y and py[i] >= y:
584                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
585                            inside = not inside
586
587        if inside:
588            indices[inside_index] = k
589            inside_index += 1
590        else:
591            indices[outside_index] = k
592            outside_index -= 1
593
594    return indices, inside_index
595
596
597def separate_points_by_polygon_c(points, polygon,
598                                 closed = True, verbose = False):
599    """Determine whether points are inside or outside a polygon
600
601    C-wrapper
602    """
603
604    from Numeric import array, Float, reshape, zeros, Int
605
606
607    if verbose: print 'Checking input to separate_points_by_polygon'
608    #Input checks
609    try:
610        points = ensure_numeric(points, Float)
611    except:
612        msg = 'Points could not be converted to Numeric array'
613        raise msg
614
615    #if verbose: print 'Checking input to separate_points_by_polygon 2'
616    try:
617        polygon = ensure_numeric(polygon, Float)
618    except:
619        msg = 'Polygon could not be converted to Numeric array'
620        raise msg
621
622    if verbose: print 'check'
623
624    assert len(polygon.shape) == 2,\
625       'Polygon array must be a 2d array of vertices'
626
627    assert polygon.shape[1] == 2,\
628       'Polygon array must have two columns'
629
630    assert len(points.shape) == 2,\
631       'Points array must be a 2d array'
632
633    assert points.shape[1] == 2,\
634       'Points array must have two columns'
635
636    N = polygon.shape[0] #Number of vertices in polygon
637    M = points.shape[0]  #Number of points
638
639    from util_ext import separate_points_by_polygon
640
641    if verbose: print 'Allocating array for indices'
642
643    indices = zeros( M, Int )
644
645    if verbose: print 'Calling C-version of inside poly'
646    count = separate_points_by_polygon(points, polygon, indices,
647                                       int(closed), int(verbose))
648
649    return indices, count
650
651
652
653class Polygon_function:
654    """Create callable object f: x,y -> z, where a,y,z are vectors and
655    where f will return different values depending on whether x,y belongs
656    to specified polygons.
657
658    To instantiate:
659
660       Polygon_function(polygons)
661
662    where polygons is a list of tuples of the form
663
664      [ (P0, v0), (P1, v1), ...]
665
666      with Pi being lists of vertices defining polygons and vi either
667      constants or functions of x,y to be applied to points with the polygon.
668
669    The function takes an optional argument, default which is the value
670    (or function) to used for points not belonging to any polygon.
671    For example:
672
673       Polygon_function(polygons, default = 0.03)
674
675    If omitted the default value will be 0.0
676
677    Note: If two polygons overlap, the one last in the list takes precedence
678
679    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???
680
681    """
682
683    def __init__(self, regions, default = 0.0):
684
685        try:
686            len(regions)
687        except:
688            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
689            raise msg
690
691
692        T = regions[0]
693        try:
694            a = len(T)
695        except:
696            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
697            raise msg
698
699        assert a == 2, 'Must have two component each: %s' %T
700
701        self.regions = regions
702        self.default = default
703
704
705    def __call__(self, x, y):
706        from util import inside_polygon
707        from Numeric import ones, Float, concatenate, array, reshape, choose
708
709        x = array(x).astype(Float)
710        y = array(y).astype(Float)
711
712        N = len(x)
713        assert len(y) == N
714
715        points = concatenate( (reshape(x, (N, 1)),
716                               reshape(y, (N, 1))), axis=1 )
717
718        if callable(self.default):
719            z = self.default(x,y)
720        else:
721            z = ones(N, Float) * self.default
722
723        for polygon, value in self.regions:
724            indices = inside_polygon(points, polygon)
725
726            #FIXME: This needs to be vectorised
727            if callable(value):
728                for i in indices:
729                    xx = array([x[i]])
730                    yy = array([y[i]])
731                    z[i] = value(xx, yy)[0]
732            else:
733                for i in indices:
734                    z[i] = value
735
736        return z
737
738def read_polygon(filename):
739    """Read points assumed to form a polygon
740       There must be exactly two numbers in each line
741    """
742
743    #Get polygon
744    fid = open(filename)
745    lines = fid.readlines()
746    fid.close()
747    polygon = []
748    for line in lines:
749        fields = line.split(',')
750        polygon.append( [float(fields[0]), float(fields[1])] )
751
752    return polygon
753
754def populate_polygon(polygon, number_of_points, seed = None):
755    """Populate given polygon with uniformly distributed points.
756
757    Input:
758       polygon - list of vertices of polygon
759       number_of_points - (optional) number of points
760       seed - seed for random number generator (default=None)
761
762    Output:
763       points - list of points inside polygon
764
765    Examples:
766       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
767       will return five randomly selected points inside the unit square
768    """
769
770    from random import uniform, seed
771
772    seed(seed)
773
774    points = []
775
776    #Find outer extent of polygon
777    max_x = min_x = polygon[0][0]
778    max_y = min_y = polygon[0][1]
779    for point in polygon[1:]:
780        x = point[0]
781        if x > max_x: max_x = x
782        if x < min_x: min_x = x
783        y = point[1]
784        if y > max_y: max_y = y
785        if y < min_y: min_y = y
786
787
788    while len(points) < number_of_points:
789        x = uniform(min_x, max_x)
790        y = uniform(min_y, max_y)
791
792        if inside_polygon( [x,y], polygon ):
793            points.append([x,y])
794
795    return points
796
797def tsh2sww(filename, verbose=True):
798    """
799    to check if a tsh/msh file 'looks' good.
800    """
801
802    #FIXME Move to data_manager
803    from shallow_water import Domain
804    from pmesh2domain import pmesh_to_domain_instance
805    import time, os
806    from data_manager import get_dataobject
807    from os import sep, path
808    #from util import mean
809
810    if verbose == True:print 'Creating domain from', filename
811    domain = pmesh_to_domain_instance(filename, Domain)
812    if verbose == True:print "Number of triangles = ", len(domain)
813
814    domain.smooth = True
815    domain.format = 'sww'   #Native netcdf visualisation format
816    file_path, filename = path.split(filename)
817    filename, ext = path.splitext(filename)
818    domain.filename = filename
819    domain.reduction = mean
820    if verbose == True:print "file_path",file_path
821    if file_path == "":file_path = "."
822    domain.set_datadir(file_path)
823
824    if verbose == True:
825        print "Output written to " + domain.get_datadir() + sep + \
826              domain.filename + "." + domain.format
827    sww = get_dataobject(domain)
828    sww.store_connectivity()
829    sww.store_timestep('stage')
830
831####################################################################
832#Python versions of function that are also implemented in util_gateway.c
833#
834
835def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
836    """
837    """
838
839    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
840    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
841    a /= det
842
843    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
844    b /= det
845
846    return a, b
847
848
849def gradient2_python(x0, y0, x1, y1, q0, q1):
850    """Compute radient based on two points and enforce zero gradient
851    in the direction orthogonal to (x1-x0), (y1-y0)
852    """
853
854    #Old code
855    #det = x0*y1 - x1*y0
856    #if det != 0.0:
857    #    a = (y1*q0 - y0*q1)/det
858    #    b = (x0*q1 - x1*q0)/det
859
860    #Correct code (ON)
861    det = (x1-x0)**2 + (y1-y0)**2
862    if det != 0.0:
863        a = (x1-x0)*(q1-q0)/det
864        b = (y1-y0)*(q1-q0)/det
865       
866    return a, b       
867
868
869##############################################
870#Initialise module
871
872import compile
873if compile.can_use_C_extension('util_ext.c'):
874    from util_ext import gradient, gradient2, point_on_line
875    separate_points_by_polygon = separate_points_by_polygon_c
876else:
877    gradient = gradient_python
878    gradient2 = gradient2_python   
879
880
881if __name__ == "__main__":
882    pass
883
884
885
886
Note: See TracBrowser for help on using the repository browser.