source: inundation/pyvolution/test_util.py @ 1903

Last change on this file since 1903 was 1903, checked in by duncan, 19 years ago

removing / gutting duplication of methods reading xya files
File size: 42.0 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7
8from util import *
9from config import epsilon
10from data_manager import timefile2netcdf
11
12
13def test_function(x, y):
14    return x+y
15
16class Test_Util(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23
24    def test_gradient(self):
25        x0 = 0.0; y0 = 0.0; z0 = 0.0
26        x1 = 1.0; y1 = 0.0; z1 = -1.0
27        x2 = 0.0; y2 = 1.0; z2 = 0.0
28
29        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
30
31        assert zx == -1.0
32        assert zy == 0.0
33
34    def test_gradient_more(self):
35        x0 = 2.0/3; y0 = 2.0/3
36        x1=  8.0/3; y1 = 2.0/3
37        x2 = 2.0/3; y2 = 8.0/3
38
39        q0 = 2.0+2.0/3
40        q1 = 8.0+2.0/3
41        q2 = 2.0+8.0/3
42
43        #Gradient of fitted pwl surface
44        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
45
46        assert abs(a - 3.0) < epsilon
47        assert abs(b - 1.0) < epsilon
48
49
50    def test_gradient2(self):
51        """Test two-point gradient
52        """
53       
54        x0 = 5.0; y0 = 5.0; z0 = 10.0
55        x1 = 8.0; y1 = 2.0; z1 = 1.0
56        x2 = 8.0; y2 = 8.0; z2 = 10.0
57
58        #Reference
59        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
60        a, b = gradient2(x0, y0, x1, y1, z0, z1)
61
62        assert zx == a
63        assert zy == b
64
65        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
66        assert z2_computed == z2
67       
68    def test_gradient2_more(self):
69        """Test two-point gradient more
70        """
71        x0 = 2.0; y0 = 2.0
72        x1 = 8.0; y1 = 3.0
73        x2 = 1.0; y2 = 8.0
74
75        q0 = 2.0
76        q1 = 8.0
77        q2 = q0
78
79        #Gradient of fitted pwl surface
80        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
81        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
82
83        assert a == a_ref
84        assert b == b_ref
85
86
87    def test_that_C_extension_compiles(self):
88        FN = 'util_ext.c'
89        try:
90            import util_ext
91        except:
92            from compile import compile
93
94            try:
95                compile(FN)
96            except:
97                raise 'Could not compile %s' %FN
98            else:
99                import util_ext
100
101
102    def test_gradient_C_extension(self):
103        from util_ext import gradient as gradient_c
104
105        x0 = 2.0/3; y0 = 2.0/3
106        x1=  8.0/3; y1 = 2.0/3
107        x2 = 2.0/3; y2 = 8.0/3
108
109        q0 = 2.0+2.0/3
110        q1 = 8.0+2.0/3
111        q2 = 2.0+8.0/3
112
113        #Gradient of fitted pwl surface
114        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
115
116        assert abs(a - 3.0) < epsilon
117        assert abs(b - 1.0) < epsilon
118
119
120    def test_gradient_C_extension3(self):
121        from util_ext import gradient as gradient_c
122
123        from RandomArray import uniform, seed
124        seed(17, 53)
125
126        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
127
128        q0 = uniform(0.0, 10.0, 4)
129        q1 = uniform(1.0, 3.0, 4)
130        q2 = uniform(7.0, 20.0, 4)
131
132
133        for i in range(4):
134            #Gradient of fitted pwl surface
135            from util import gradient_python
136            a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
137                                    q0[i], q1[i], q2[i])
138
139            #print a_ref, b_ref
140            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
141                              q0[i], q1[i], q2[i])
142
143            #print a, a_ref, b, b_ref
144            assert abs(a - a_ref) < epsilon
145            assert abs(b - b_ref) < epsilon
146
147
148
149    #Geometric
150    #def test_distance(self):
151    #    from util import distance#
152    #
153    #    self.failUnless( distance([4,2],[7,6]) == 5.0,
154    #                     'Distance is wrong!')
155    #    self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
156    #                    'distance is wrong!')
157    #    self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
158    #                    'distance is wrong!')
159    #
160    #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
161    #                    'distance is wrong!')
162
163    def test_angle(self):
164        from util import angle
165
166        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
167
168
169    def test_anglediff(self):
170        from util import anglediff
171
172        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
173
174    def test_ensure_numeric(self):
175        from util import ensure_numeric
176        from Numeric import ArrayType, Float, array
177
178        A = [1,2,3,4]
179        B = ensure_numeric(A)
180        assert type(B) == ArrayType
181        assert B.typecode() == 'l'
182        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
183
184
185        A = [1,2,3.14,4]
186        B = ensure_numeric(A)
187        assert type(B) == ArrayType
188        assert B.typecode() == 'd'
189        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
190
191
192        A = [1,2,3,4]
193        B = ensure_numeric(A, Float)
194        assert type(B) == ArrayType
195        assert B.typecode() == 'd'
196        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
197
198
199        A = [1,2,3,4]
200        B = ensure_numeric(A, Float)
201        assert type(B) == ArrayType
202        assert B.typecode() == 'd'
203        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
204
205
206        A = array([1,2,3,4])
207        B = ensure_numeric(A)
208        assert type(B) == ArrayType
209        assert B.typecode() == 'l'       
210        assert A == B   
211        assert A is B   #Same object
212
213
214        A = array([1,2,3,4])
215        B = ensure_numeric(A, Float)
216        assert type(B) == ArrayType
217        assert B.typecode() == 'd'       
218        assert A == B   
219        assert A is not B   #Not the same object       
220
221
222
223    def test_file_function_time1(self):
224        """Test that File function interpolates correctly
225        between given times. No x,y dependency here.
226        """
227
228        #Write file
229        import os, time
230        from config import time_format
231        from math import sin, pi
232
233        #Typical ASCII file
234        finaltime = 1200
235        filename = 'test_file_function'
236        fid = open(filename + '.txt', 'w')
237        start = time.mktime(time.strptime('2000', '%Y'))
238        dt = 60  #One minute intervals
239        t = 0.0
240        while t <= finaltime:
241            t_string = time.strftime(time_format, time.gmtime(t+start))
242            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
243            t += dt
244
245        fid.close()
246
247        #Convert ASCII file to NetCDF (Which is what we really like!)
248        timefile2netcdf(filename)
249
250
251        #Create file function from time series
252        F = file_function(filename + '.tms',
253                          quantities = ['Attribute0',
254                                        'Attribute1',
255                                        'Attribute2'])
256       
257        #Now try interpolation
258        for i in range(20):
259            t = i*10
260            q = F(t)
261
262            #Exact linear intpolation
263            assert allclose(q[0], 2*t)
264            if i%6 == 0:
265                assert allclose(q[1], t**2)
266                assert allclose(q[2], sin(t*pi/600))
267
268        #Check non-exact
269
270        t = 90 #Halfway between 60 and 120
271        q = F(t)
272        assert allclose( (120**2 + 60**2)/2, q[1] )
273        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
274
275
276        t = 100 #Two thirds of the way between between 60 and 120
277        q = F(t)
278        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
279        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
280
281        os.remove(filename + '.txt')
282        os.remove(filename + '.tms')       
283
284
285       
286    def test_spatio_temporal_file_function(self):
287        """Test that spatio temporal file function performs the correct
288        interpolations in both time and space
289        NetCDF version (x,y,t dependency)       
290        """
291        import time
292
293        #Create sww file of simple propagation from left to right
294        #through rectangular domain
295        from shallow_water import Domain, Dirichlet_boundary
296        from mesh_factory import rectangular
297        from Numeric import take, concatenate, reshape
298
299        #Create basic mesh and shallow water domain
300        points, vertices, boundary = rectangular(3, 3)
301        domain1 = Domain(points, vertices, boundary)
302
303        from util import mean
304        domain1.reduction = mean
305        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
306                              # only one value.
307
308        domain1.default_order = 2
309        domain1.store = True
310        domain1.set_datadir('.')
311        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
312        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
313
314        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
315        domain1.set_quantity('elevation', 0)
316        domain1.set_quantity('friction', 0)
317        domain1.set_quantity('stage', 0)
318
319        # Boundary conditions
320        B0 = Dirichlet_boundary([0,0,0])
321        B6 = Dirichlet_boundary([0.6,0,0])
322        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
323        domain1.check_integrity()
324
325        finaltime = 8
326        #Evolution
327        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
328            pass
329            #domain1.write_time()
330
331
332        #Now read data from sww and check
333        from Scientific.IO.NetCDF import NetCDFFile
334        filename = domain1.get_name() + '.' + domain1.format
335        fid = NetCDFFile(filename)
336
337        x = fid.variables['x'][:]
338        y = fid.variables['y'][:]
339        stage = fid.variables['stage'][:]
340        xmomentum = fid.variables['xmomentum'][:]
341        ymomentum = fid.variables['ymomentum'][:]
342        time = fid.variables['time'][:]
343
344        #Take stage vertex values at last timestep on diagonal
345        #Diagonal is identified by vertices: 0, 5, 10, 15
346
347        timestep = len(time)-1 #Last timestep
348        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
349        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
350        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
351        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
352
353        #Reference interpolated values at midpoints on diagonal at
354        #this timestep are
355        r0 = (D[0] + D[1])/2
356        r1 = (D[1] + D[2])/2
357        r2 = (D[2] + D[3])/2
358
359        #And the midpoints are found now
360        Dx = take(reshape(x, (16,1)), [0,5,10,15])
361        Dy = take(reshape(y, (16,1)), [0,5,10,15])
362
363        diag = concatenate( (Dx, Dy), axis=1)
364        d_midpoints = (diag[1:] + diag[:-1])/2
365
366        #Let us see if the file function can find the correct
367        #values at the midpoints at the last timestep:
368        f = file_function(filename, domain1,
369                          interpolation_points = d_midpoints)
370
371        q = f(timestep/10., point_id=0); assert allclose(r0, q)
372        q = f(timestep/10., point_id=1); assert allclose(r1, q)
373        q = f(timestep/10., point_id=2); assert allclose(r2, q)
374
375
376        ##################
377        #Now do the same for the first timestep
378
379        timestep = 0 #First timestep
380        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
381        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
382        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
383        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
384
385        #Reference interpolated values at midpoints on diagonal at
386        #this timestep are
387        r0 = (D[0] + D[1])/2
388        r1 = (D[1] + D[2])/2
389        r2 = (D[2] + D[3])/2
390
391        #Let us see if the file function can find the correct
392        #values
393        q = f(0, point_id=0); assert allclose(r0, q)
394        q = f(0, point_id=1); assert allclose(r1, q)
395        q = f(0, point_id=2); assert allclose(r2, q)
396
397
398        ##################
399        #Now do it again for a timestep in the middle
400
401        timestep = 33
402        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
403        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
404        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
405        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
406
407        #Reference interpolated values at midpoints on diagonal at
408        #this timestep are
409        r0 = (D[0] + D[1])/2
410        r1 = (D[1] + D[2])/2
411        r2 = (D[2] + D[3])/2
412
413        q = f(timestep/10., point_id=0); assert allclose(r0, q)
414        q = f(timestep/10., point_id=1); assert allclose(r1, q)
415        q = f(timestep/10., point_id=2); assert allclose(r2, q)
416
417
418        ##################
419        #Now check temporal interpolation
420        #Halfway between timestep 15 and 16
421
422        timestep = 15
423        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
424        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
425        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
426        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
427
428        #Reference interpolated values at midpoints on diagonal at
429        #this timestep are
430        r0_0 = (D[0] + D[1])/2
431        r1_0 = (D[1] + D[2])/2
432        r2_0 = (D[2] + D[3])/2
433
434        #
435        timestep = 16
436        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
437        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
438        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
439        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
440
441        #Reference interpolated values at midpoints on diagonal at
442        #this timestep are
443        r0_1 = (D[0] + D[1])/2
444        r1_1 = (D[1] + D[2])/2
445        r2_1 = (D[2] + D[3])/2
446
447        # The reference values are
448        r0 = (r0_0 + r0_1)/2
449        r1 = (r1_0 + r1_1)/2
450        r2 = (r2_0 + r2_1)/2
451
452        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
453        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
454        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
455
456        ##################
457        #Finally check interpolation 2 thirds of the way
458        #between timestep 15 and 16
459
460        # The reference values are
461        r0 = (r0_0 + 2*r0_1)/3
462        r1 = (r1_0 + 2*r1_1)/3
463        r2 = (r2_0 + 2*r2_1)/3
464
465        #And the file function gives
466        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
467        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
468        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
469
470        fid.close()
471        import os
472        os.remove(filename)
473
474
475
476    def test_spatio_temporal_file_function_different_origin(self):
477        """Test that spatio temporal file function performs the correct
478        interpolations in both time and space where space is offset by
479        xllcorner and yllcorner
480        NetCDF version (x,y,t dependency)       
481        """
482        import time
483
484        #Create sww file of simple propagation from left to right
485        #through rectangular domain
486        from shallow_water import Domain, Dirichlet_boundary
487        from mesh_factory import rectangular
488        from Numeric import take, concatenate, reshape
489
490
491        from coordinate_transforms.geo_reference import Geo_reference
492        xllcorner = 2048
493        yllcorner = 11000
494        zone = 2
495
496        #Create basic mesh and shallow water domain
497        points, vertices, boundary = rectangular(3, 3)
498        domain1 = Domain(points, vertices, boundary,
499                         geo_reference = Geo_reference(xllcorner = xllcorner,
500                                                       yllcorner = yllcorner))
501       
502
503        from util import mean
504        domain1.reduction = mean
505        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
506                              # only one value.
507
508        domain1.default_order = 2
509        domain1.store = True
510        domain1.set_datadir('.')
511        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
512        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
513
514        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
515        domain1.set_quantity('elevation', 0)
516        domain1.set_quantity('friction', 0)
517        domain1.set_quantity('stage', 0)
518
519        # Boundary conditions
520        B0 = Dirichlet_boundary([0,0,0])
521        B6 = Dirichlet_boundary([0.6,0,0])
522        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
523        domain1.check_integrity()
524
525        finaltime = 8
526        #Evolution
527        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
528            pass
529            #domain1.write_time()
530
531
532        #Now read data from sww and check
533        from Scientific.IO.NetCDF import NetCDFFile
534        filename = domain1.get_name() + '.' + domain1.format
535        fid = NetCDFFile(filename)
536
537        x = fid.variables['x'][:]
538        y = fid.variables['y'][:]
539        stage = fid.variables['stage'][:]
540        xmomentum = fid.variables['xmomentum'][:]
541        ymomentum = fid.variables['ymomentum'][:]
542        time = fid.variables['time'][:]
543
544        #Take stage vertex values at last timestep on diagonal
545        #Diagonal is identified by vertices: 0, 5, 10, 15
546
547        timestep = len(time)-1 #Last timestep
548        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
549        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
550        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
551        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
552
553        #Reference interpolated values at midpoints on diagonal at
554        #this timestep are
555        r0 = (D[0] + D[1])/2
556        r1 = (D[1] + D[2])/2
557        r2 = (D[2] + D[3])/2
558
559        #And the midpoints are found now
560        Dx = take(reshape(x, (16,1)), [0,5,10,15])
561        Dy = take(reshape(y, (16,1)), [0,5,10,15])
562
563        diag = concatenate( (Dx, Dy), axis=1)
564        d_midpoints = (diag[1:] + diag[:-1])/2
565
566
567        #Adjust for georef - make interpolation points absolute
568        d_midpoints[:,0] += xllcorner
569        d_midpoints[:,1] += yllcorner               
570
571        #Let us see if the file function can find the correct
572        #values at the midpoints at the last timestep:
573        f = file_function(filename, domain1,
574                          interpolation_points = d_midpoints)
575
576        q = f(timestep/10., point_id=0); assert allclose(r0, q)
577        q = f(timestep/10., point_id=1); assert allclose(r1, q)
578        q = f(timestep/10., point_id=2); assert allclose(r2, q)
579
580
581        ##################
582        #Now do the same for the first timestep
583
584        timestep = 0 #First timestep
585        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
586        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
587        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
588        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
589
590        #Reference interpolated values at midpoints on diagonal at
591        #this timestep are
592        r0 = (D[0] + D[1])/2
593        r1 = (D[1] + D[2])/2
594        r2 = (D[2] + D[3])/2
595
596        #Let us see if the file function can find the correct
597        #values
598        q = f(0, point_id=0); assert allclose(r0, q)
599        q = f(0, point_id=1); assert allclose(r1, q)
600        q = f(0, point_id=2); assert allclose(r2, q)
601
602
603        ##################
604        #Now do it again for a timestep in the middle
605
606        timestep = 33
607        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
608        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
609        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
610        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
611
612        #Reference interpolated values at midpoints on diagonal at
613        #this timestep are
614        r0 = (D[0] + D[1])/2
615        r1 = (D[1] + D[2])/2
616        r2 = (D[2] + D[3])/2
617
618        q = f(timestep/10., point_id=0); assert allclose(r0, q)
619        q = f(timestep/10., point_id=1); assert allclose(r1, q)
620        q = f(timestep/10., point_id=2); assert allclose(r2, q)
621
622
623        ##################
624        #Now check temporal interpolation
625        #Halfway between timestep 15 and 16
626
627        timestep = 15
628        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
629        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
630        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
631        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
632
633        #Reference interpolated values at midpoints on diagonal at
634        #this timestep are
635        r0_0 = (D[0] + D[1])/2
636        r1_0 = (D[1] + D[2])/2
637        r2_0 = (D[2] + D[3])/2
638
639        #
640        timestep = 16
641        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
642        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
643        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
644        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
645
646        #Reference interpolated values at midpoints on diagonal at
647        #this timestep are
648        r0_1 = (D[0] + D[1])/2
649        r1_1 = (D[1] + D[2])/2
650        r2_1 = (D[2] + D[3])/2
651
652        # The reference values are
653        r0 = (r0_0 + r0_1)/2
654        r1 = (r1_0 + r1_1)/2
655        r2 = (r2_0 + r2_1)/2
656
657        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
658        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
659        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
660
661        ##################
662        #Finally check interpolation 2 thirds of the way
663        #between timestep 15 and 16
664
665        # The reference values are
666        r0 = (r0_0 + 2*r0_1)/3
667        r1 = (r1_0 + 2*r1_1)/3
668        r2 = (r2_0 + 2*r2_1)/3
669
670        #And the file function gives
671        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
672        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
673        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
674
675        fid.close()
676        import os
677        os.remove(filename)
678
679       
680
681
682    def test_spatio_temporal_file_function_time(self):
683        """Test that File function interpolates correctly
684        between given times.
685        NetCDF version (x,y,t dependency)
686        """
687
688        #Create NetCDF (sww) file to be read
689        # x: 0, 5, 10, 15
690        # y: -20, -10, 0, 10
691        # t: 0, 60, 120, ...., 1200
692        #
693        # test quantities (arbitrary but non-trivial expressions):
694        #
695        #   stage     = 3*x - y**2 + 2*t
696        #   xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2
697        #   ymomentum = x**2 + y**2 * sin(t*pi/600)
698
699        #NOTE: Nice test that may render some of the others redundant.
700
701        import os, time
702        from config import time_format
703        from Numeric import sin, pi, exp
704        from mesh_factory import rectangular
705        from shallow_water import Domain
706        import data_manager
707
708        finaltime = 1200
709        filename = 'test_file_function'
710
711        #Create a domain to hold test grid
712        #(0:15, -20:10)
713        points, vertices, boundary =\
714                rectangular(4, 4, 15, 30, origin = (0, -20))
715
716
717        #print 'Number of elements', len(vertices)
718        domain = Domain(points, vertices, boundary)
719        domain.smooth = False
720        domain.default_order = 2
721        domain.set_datadir('.')
722        domain.set_name(filename)
723        domain.store = True
724        domain.format = 'sww'   #Native netcdf visualisation format
725
726        #print points
727        start = time.mktime(time.strptime('2000', '%Y'))
728        domain.starttime = start
729
730
731        #Store structure
732        domain.initialise_storage()
733
734        #Compute artificial time steps and store
735        dt = 60  #One minute intervals
736        t = 0.0
737        while t <= finaltime:
738            #Compute quantities
739            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
740            domain.set_quantity('stage', f1)
741
742            f2 = lambda x,y: x+y+t**2
743            domain.set_quantity('xmomentum', f2)
744
745            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
746            domain.set_quantity('ymomentum', f3)
747
748            #Store and advance time
749            domain.time = t
750            domain.store_timestep(domain.conserved_quantities)
751            t += dt
752
753
754        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
755
756
757
758        #Deliberately set domain.starttime to too early
759        domain.starttime = start - 1
760
761        #Create file function
762        F = file_function(filename + '.sww', domain,
763                          quantities = domain.conserved_quantities,
764                          interpolation_points = interpolation_points)
765
766        #Check that FF updates fixes domain starttime
767        assert allclose(domain.starttime, start)
768
769        #Check that domain.starttime isn't updated if later
770        domain.starttime = start + 1
771        F = file_function(filename + '.sww', domain,
772                          quantities = domain.conserved_quantities,
773                          interpolation_points = interpolation_points)
774        assert allclose(domain.starttime, start+1)
775        domain.starttime = start
776
777
778        #Check linear interpolation in time
779        F = file_function(filename + '.sww', domain,
780                          quantities = domain.conserved_quantities,
781                          interpolation_points = interpolation_points)               
782        for id in range(len(interpolation_points)):
783            x = interpolation_points[id][0]
784            y = interpolation_points[id][1]
785
786            for i in range(20):
787                t = i*10
788                k = i%6
789
790                if k == 0:
791                    q0 = F(t, point_id=id)
792                    q1 = F(t+60, point_id=id)
793
794
795                q = F(t, point_id=id)
796                #print i, k, t, q
797                #print ' ', q0
798                #print ' ', q1
799                #print 's', (k*q1 + (6-k)*q0)/6
800                #print
801                assert allclose(q, (k*q1 + (6-k)*q0)/6)
802
803
804        #Another check of linear interpolation in time
805        for id in range(len(interpolation_points)):
806            q60 = F(60, point_id=id)
807            q120 = F(120, point_id=id)
808
809            t = 90 #Halfway between 60 and 120
810            q = F(t, point_id=id)
811            assert allclose( (q120+q60)/2, q )
812
813            t = 100 #Two thirds of the way between between 60 and 120
814            q = F(t, point_id=id)
815            assert allclose(q60/3 + 2*q120/3, q)
816
817
818
819        #Check that domain.starttime isn't updated if later than file starttime but earlier
820        #than file end time
821        delta = 23
822        domain.starttime = start + delta
823        F = file_function(filename + '.sww', domain,
824                          quantities = domain.conserved_quantities,
825                          interpolation_points = interpolation_points)
826        assert allclose(domain.starttime, start+delta)
827
828
829
830
831        #Now try interpolation with delta offset
832        for id in range(len(interpolation_points)):           
833            x = interpolation_points[id][0]
834            y = interpolation_points[id][1]
835
836            for i in range(20):
837                t = i*10
838                k = i%6
839
840                if k == 0:
841                    q0 = F(t-delta, point_id=id)
842                    q1 = F(t+60-delta, point_id=id)
843
844                q = F(t-delta, point_id=id)
845                assert allclose(q, (k*q1 + (6-k)*q0)/6)
846
847
848        os.remove(filename + '.sww')
849
850
851
852    def test_file_function_time_with_domain(self):
853        """Test that File function interpolates correctly
854        between given times. No x,y dependency here.
855        Use domain with starttime
856        """
857
858        #Write file
859        import os, time, calendar
860        from config import time_format
861        from math import sin, pi
862        from domain import Domain
863
864        finaltime = 1200
865        filename = 'test_file_function'
866        fid = open(filename + '.txt', 'w')
867        start = time.mktime(time.strptime('2000', '%Y'))
868        dt = 60  #One minute intervals
869        t = 0.0
870        while t <= finaltime:
871            t_string = time.strftime(time_format, time.gmtime(t+start))
872            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
873            t += dt
874
875        fid.close()
876
877
878        #Convert ASCII file to NetCDF (Which is what we really like!)
879        timefile2netcdf(filename)
880
881
882
883        a = [0.0, 0.0]
884        b = [4.0, 0.0]
885        c = [0.0, 3.0]
886
887        points = [a, b, c]
888        vertices = [[0,1,2]]
889        domain = Domain(points, vertices)
890
891        #Check that domain.starttime is updated if non-existing
892        F = file_function(filename + '.tms', domain)
893       
894        assert allclose(domain.starttime, start)
895
896        #Check that domain.starttime is updated if too early
897        domain.starttime = start - 1
898        F = file_function(filename + '.tms', domain)
899        assert allclose(domain.starttime, start)
900
901        #Check that domain.starttime isn't updated if later
902        domain.starttime = start + 1
903        F = file_function(filename + '.tms', domain)       
904        assert allclose(domain.starttime, start+1)
905
906        domain.starttime = start
907        F = file_function(filename + '.tms', domain,
908                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
909       
910
911        #print F.T
912        #print F.precomputed_values
913        #print 'F(60)', F(60)
914       
915        #Now try interpolation
916        for i in range(20):
917            t = i*10
918            q = F(t)
919
920            #Exact linear intpolation
921            assert allclose(q[0], 2*t)
922            if i%6 == 0:
923                assert allclose(q[1], t**2)
924                assert allclose(q[2], sin(t*pi/600))
925
926        #Check non-exact
927
928        t = 90 #Halfway between 60 and 120
929        q = F(t)
930        assert allclose( (120**2 + 60**2)/2, q[1] )
931        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
932
933
934        t = 100 #Two thirds of the way between between 60 and 120
935        q = F(t)
936        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
937        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
938
939        os.remove(filename + '.tms')
940        os.remove(filename + '.txt')       
941
942    def test_file_function_time_with_domain_different_start(self):
943        """Test that File function interpolates correctly
944        between given times. No x,y dependency here.
945        Use domain with a starttime later than that of file
946
947        ASCII version
948        """
949
950        #Write file
951        import os, time, calendar
952        from config import time_format
953        from math import sin, pi
954        from domain import Domain
955
956        finaltime = 1200
957        filename = 'test_file_function'
958        fid = open(filename + '.txt', 'w')
959        start = time.mktime(time.strptime('2000', '%Y'))
960        dt = 60  #One minute intervals
961        t = 0.0
962        while t <= finaltime:
963            t_string = time.strftime(time_format, time.gmtime(t+start))
964            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
965            t += dt
966
967        fid.close()
968
969        #Convert ASCII file to NetCDF (Which is what we really like!)
970        timefile2netcdf(filename)       
971
972        a = [0.0, 0.0]
973        b = [4.0, 0.0]
974        c = [0.0, 3.0]
975
976        points = [a, b, c]
977        vertices = [[0,1,2]]
978        domain = Domain(points, vertices)
979
980        #Check that domain.starttime isn't updated if later than file starttime but earlier
981        #than file end time
982        delta = 23
983        domain.starttime = start + delta
984        F = file_function(filename + '.tms', domain,
985                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])       
986        assert allclose(domain.starttime, start+delta)
987
988
989
990
991        #Now try interpolation with delta offset
992        for i in range(20):
993            t = i*10
994            q = F(t-delta)
995
996            #Exact linear intpolation
997            assert allclose(q[0], 2*t)
998            if i%6 == 0:
999                assert allclose(q[1], t**2)
1000                assert allclose(q[2], sin(t*pi/600))
1001
1002        #Check non-exact
1003
1004        t = 90 #Halfway between 60 and 120
1005        q = F(t-delta)
1006        assert allclose( (120**2 + 60**2)/2, q[1] )
1007        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
1008
1009
1010        t = 100 #Two thirds of the way between between 60 and 120
1011        q = F(t-delta)
1012        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
1013        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
1014
1015
1016        os.remove(filename + '.tms')
1017        os.remove(filename + '.txt')               
1018
1019
1020
1021
1022    #Polygon stuff
1023    def test_polygon_function_constants(self):
1024        p1 = [[0,0], [10,0], [10,10], [0,10]]
1025        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1026
1027        f = Polygon_function( [(p1, 1.0)] )
1028        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
1029        assert allclose(z, [1,1,0,0])
1030
1031
1032        f = Polygon_function( [(p2, 2.0)] )
1033        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
1034        assert allclose(z, [2,0,0,2])
1035
1036
1037        #Combined
1038        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
1039        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1040        assert allclose(z, [2,1,0,2])
1041
1042
1043    def test_polygon_function_callable(self):
1044        """Check that values passed into Polygon_function can be callable
1045        themselves.
1046        """
1047        p1 = [[0,0], [10,0], [10,10], [0,10]]
1048        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1049
1050        f = Polygon_function( [(p1, test_function)] )
1051        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
1052        assert allclose(z, [10,14,0,0])
1053
1054        #Combined
1055        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
1056        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1057        assert allclose(z, [2,14,0,2])
1058
1059
1060        #Combined w default
1061        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
1062        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1063        assert allclose(z, [2,14,3.14,2])
1064
1065
1066        #Combined w default func
1067        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
1068                              default = test_function)
1069        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1070        assert allclose(z, [2,14,35,2])
1071
1072
1073    def test_point_on_line(self):
1074
1075        #Endpoints first
1076        assert point_on_line( 0, 0, 0,0, 1,0 )
1077        assert point_on_line( 1, 0, 0,0, 1,0 )
1078
1079        #Then points on line
1080        assert point_on_line( 0.5, 0, 0,0, 1,0 )
1081        assert point_on_line( 0, 0.5, 0,1, 0,0 )
1082        assert point_on_line( 1, 0.5, 1,1, 1,0 )
1083        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
1084
1085        #Then points not on line
1086        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
1087        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
1088
1089        #From real example that failed
1090        assert not point_on_line( 40,50, 40,20, 40,40 )
1091
1092
1093        #From real example that failed
1094        assert not point_on_line( 40,19, 40,20, 40,40 )
1095
1096
1097
1098
1099    def test_inside_polygon_main(self):
1100
1101
1102        #Simplest case: Polygon is the unit square
1103        polygon = [[0,0], [1,0], [1,1], [0,1]]
1104
1105        assert inside_polygon( (0.5, 0.5), polygon )
1106        assert not inside_polygon( (0.5, 1.5), polygon )
1107        assert not inside_polygon( (0.5, -0.5), polygon )
1108        assert not inside_polygon( (-0.5, 0.5), polygon )
1109        assert not inside_polygon( (1.5, 0.5), polygon )
1110
1111        #Try point on borders
1112        assert inside_polygon( (1., 0.5), polygon, closed=True)
1113        assert inside_polygon( (0.5, 1), polygon, closed=True)
1114        assert inside_polygon( (0., 0.5), polygon, closed=True)
1115        assert inside_polygon( (0.5, 0.), polygon, closed=True)
1116
1117        assert not inside_polygon( (0.5, 1), polygon, closed=False)
1118        assert not inside_polygon( (0., 0.5), polygon, closed=False)
1119        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
1120        assert not inside_polygon( (1., 0.5), polygon, closed=False)
1121
1122
1123
1124        #From real example (that failed)
1125        polygon = [[20,20], [40,20], [40,40], [20,40]]
1126        points = [ [40, 50] ]
1127        res = inside_polygon(points, polygon)
1128        assert len(res) == 0
1129
1130        polygon = [[20,20], [40,20], [40,40], [20,40]]
1131        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
1132        res = inside_polygon(points, polygon)
1133        assert len(res) == 2
1134        assert allclose(res, [0,1])
1135
1136
1137
1138        #More convoluted and non convex polygon
1139        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1140        assert inside_polygon( (0.5, 0.5), polygon )
1141        assert inside_polygon( (1, -0.5), polygon )
1142        assert inside_polygon( (1.5, 0), polygon )
1143
1144        assert not inside_polygon( (0.5, 1.5), polygon )
1145        assert not inside_polygon( (0.5, -0.5), polygon )
1146
1147
1148        #Very convoluted polygon
1149        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1150        assert inside_polygon( (5, 5), polygon )
1151        assert inside_polygon( (17, 7), polygon )
1152        assert inside_polygon( (27, 2), polygon )
1153        assert inside_polygon( (35, -5), polygon )
1154        assert not inside_polygon( (15, 7), polygon )
1155        assert not inside_polygon( (24, 3), polygon )
1156        assert not inside_polygon( (25, -10), polygon )
1157
1158
1159
1160        #Another combination (that failed)
1161        polygon = [[0,0], [10,0], [10,10], [0,10]]
1162        assert inside_polygon( (5, 5), polygon )
1163        assert inside_polygon( (7, 7), polygon )
1164        assert not inside_polygon( (-17, 7), polygon )
1165        assert not inside_polygon( (7, 17), polygon )
1166        assert not inside_polygon( (17, 7), polygon )
1167        assert not inside_polygon( (27, 8), polygon )
1168        assert not inside_polygon( (35, -5), polygon )
1169
1170
1171
1172
1173        #Multiple polygons
1174
1175        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
1176                   [10,10], [11,10], [11,11], [10,11], [10,10]]
1177        assert inside_polygon( (0.5, 0.5), polygon )
1178        assert inside_polygon( (10.5, 10.5), polygon )
1179
1180        #FIXME: Fails if point is 5.5, 5.5
1181        assert not inside_polygon( (0, 5.5), polygon )
1182
1183        #Polygon with a hole
1184        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
1185                   [0,0], [1,0], [1,1], [0,1], [0,0]]
1186
1187        assert inside_polygon( (0, -0.5), polygon )
1188        assert not inside_polygon( (0.5, 0.5), polygon )
1189
1190    def test_inside_polygon_vector_version(self):
1191        #Now try the vector formulation returning indices
1192        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1193        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1194        res = inside_polygon( points, polygon, verbose=False )
1195
1196        assert allclose( res, [0,1,2] )
1197
1198    def test_outside_polygon(self):
1199        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1200
1201        assert not outside_polygon( [0.5, 0.5], U )
1202        #evaluate to False as the point 0.5, 0.5 is inside the unit square
1203       
1204        assert outside_polygon( [1.5, 0.5], U )
1205        #evaluate to True as the point 1.5, 0.5 is outside the unit square
1206       
1207        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1208        assert allclose( indices, [1] )
1209       
1210        #One more test of vector formulation returning indices
1211        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1212        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1213        res = outside_polygon( points, polygon )
1214
1215        assert allclose( res, [3, 4] )
1216
1217
1218
1219        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1220        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1221        res = outside_polygon( points, polygon )
1222
1223        assert allclose( res, [0, 4, 5] )       
1224     
1225    def test_outside_polygon2(self):
1226        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1227   
1228        assert not outside_polygon( [0.5, 1.0], U, closed = True )
1229        #evaluate to False as the point 0.5, 1.0 is inside the unit square
1230       
1231        assert outside_polygon( [0.5, 1.0], U, closed = False )
1232        #evaluate to True as the point 0.5, 1.0 is outside the unit square
1233
1234    def test_separate_points_by_polygon(self):
1235        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1236
1237        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1238        assert allclose( indices, [0,2,1] )
1239        assert count == 2
1240       
1241        #One more test of vector formulation returning indices
1242        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1243        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1244        res, count = separate_points_by_polygon( points, polygon )
1245
1246        assert allclose( res, [0,1,2,4,3] )
1247        assert count == 3
1248
1249
1250        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1251        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1252        res, count = separate_points_by_polygon( points, polygon )
1253
1254        assert allclose( res, [1,2,3,5,4,0] )       
1255        assert count == 3
1256       
1257
1258    def test_populate_polygon(self):
1259
1260        polygon = [[0,0], [1,0], [1,1], [0,1]]
1261        points = populate_polygon(polygon, 5)
1262
1263        assert len(points) == 5
1264        for point in points:
1265            assert inside_polygon(point, polygon)
1266
1267
1268        #Very convoluted polygon
1269        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1270
1271        points = populate_polygon(polygon, 5)
1272
1273        assert len(points) == 5
1274        for point in points:
1275            assert inside_polygon(point, polygon)
1276
1277
1278    def test_populate_polygon_with_exclude(self):
1279       
1280
1281        polygon = [[0,0], [1,0], [1,1], [0,1]]
1282        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
1283        points = populate_polygon(polygon, 5, exclude = [ex_poly])
1284
1285        assert len(points) == 5
1286        for point in points:
1287            assert inside_polygon(point, polygon)
1288            assert not inside_polygon(point, ex_poly)           
1289
1290
1291        #overlap
1292        polygon = [[0,0], [1,0], [1,1], [0,1]]
1293        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
1294        points = populate_polygon(polygon, 5, exclude = [ex_poly])
1295
1296        assert len(points) == 5
1297        for point in points:
1298            assert inside_polygon(point, polygon)
1299            assert not inside_polygon(point, ex_poly)                       
1300       
1301        #Multiple
1302        polygon = [[0,0], [1,0], [1,1], [0,1]]
1303        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
1304        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
1305       
1306        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
1307
1308        assert len(points) == 20
1309        for point in points:
1310            assert inside_polygon(point, polygon)
1311            assert not inside_polygon(point, ex_poly1)
1312            assert not inside_polygon(point, ex_poly2)                               
1313       
1314
1315        #Very convoluted polygon
1316        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1317        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
1318        points = populate_polygon(polygon, 20, exclude = [ex_poly])
1319       
1320        assert len(points) == 20
1321        for point in points:
1322            assert inside_polygon(point, polygon)
1323            assert not inside_polygon(point, ex_poly), '%s' %str(point)                       
1324
1325
1326
1327#-------------------------------------------------------------
1328if __name__ == "__main__":
1329    suite = unittest.makeSuite(Test_Util,'test')
1330    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
1331    runner = unittest.TextTestRunner()
1332    runner.run(suite)
1333
1334
1335
1336
Note: See TracBrowser for help on using the repository browser.