source: inundation/pyvolution/test_util.py @ 1884

Last change on this file since 1884 was 1884, checked in by ole, 18 years ago

Made interpolation_points in file_function absolute UTM coordinates and wrote test
File size: 43.2 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    def test_xya_ascii(self):
1022        import time, os
1023        FN = 'xyatest' + str(time.time()) + '.xya'
1024        fid = open(FN, 'w')
1025        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
1026        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
1027        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
1028        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
1029        fid.close()
1030
1031        points, attributes = read_xya(FN, format = 'asc')
1032
1033        assert allclose(points, [ [0,1], [1,0], [1,1] ])
1034        assert allclose(attributes['a1'], [10,30,40.2])
1035        assert allclose(attributes['a2'], [20,20,40.3])
1036        assert allclose(attributes['a3'], [30,10,40.4])
1037
1038        os.remove(FN)
1039
1040    def test_xya_ascii_w_names(self):
1041        import time, os
1042        FN = 'xyatest' + str(time.time()) + '.xya'
1043        fid = open(FN, 'w')
1044        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
1045        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
1046        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
1047        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
1048        fid.close()
1049
1050        points, attributes = read_xya(FN, format = 'asc')
1051
1052        assert allclose(points, [ [0,1], [1,0], [1,1] ])
1053
1054        assert allclose(attributes['a1'], [10,30,40.2])
1055        assert allclose(attributes['a2'], [20,20,40.3])
1056        assert allclose(attributes['a3'], [30,10,40.4])
1057
1058
1059        os.remove(FN)
1060
1061
1062
1063
1064    #Polygon stuff
1065    def test_polygon_function_constants(self):
1066        p1 = [[0,0], [10,0], [10,10], [0,10]]
1067        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1068
1069        f = Polygon_function( [(p1, 1.0)] )
1070        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
1071        assert allclose(z, [1,1,0,0])
1072
1073
1074        f = Polygon_function( [(p2, 2.0)] )
1075        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
1076        assert allclose(z, [2,0,0,2])
1077
1078
1079        #Combined
1080        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
1081        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1082        assert allclose(z, [2,1,0,2])
1083
1084
1085    def test_polygon_function_callable(self):
1086        """Check that values passed into Polygon_function can be callable
1087        themselves.
1088        """
1089        p1 = [[0,0], [10,0], [10,10], [0,10]]
1090        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1091
1092        f = Polygon_function( [(p1, test_function)] )
1093        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
1094        assert allclose(z, [10,14,0,0])
1095
1096        #Combined
1097        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
1098        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1099        assert allclose(z, [2,14,0,2])
1100
1101
1102        #Combined w default
1103        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
1104        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1105        assert allclose(z, [2,14,3.14,2])
1106
1107
1108        #Combined w default func
1109        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
1110                              default = test_function)
1111        z = f([5, 5, 27, 35], [5, 9, 8, -5])
1112        assert allclose(z, [2,14,35,2])
1113
1114
1115    def test_point_on_line(self):
1116
1117        #Endpoints first
1118        assert point_on_line( 0, 0, 0,0, 1,0 )
1119        assert point_on_line( 1, 0, 0,0, 1,0 )
1120
1121        #Then points on line
1122        assert point_on_line( 0.5, 0, 0,0, 1,0 )
1123        assert point_on_line( 0, 0.5, 0,1, 0,0 )
1124        assert point_on_line( 1, 0.5, 1,1, 1,0 )
1125        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
1126
1127        #Then points not on line
1128        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
1129        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
1130
1131        #From real example that failed
1132        assert not point_on_line( 40,50, 40,20, 40,40 )
1133
1134
1135        #From real example that failed
1136        assert not point_on_line( 40,19, 40,20, 40,40 )
1137
1138
1139
1140
1141    def test_inside_polygon_main(self):
1142
1143
1144        #Simplest case: Polygon is the unit square
1145        polygon = [[0,0], [1,0], [1,1], [0,1]]
1146
1147        assert inside_polygon( (0.5, 0.5), polygon )
1148        assert not inside_polygon( (0.5, 1.5), polygon )
1149        assert not inside_polygon( (0.5, -0.5), polygon )
1150        assert not inside_polygon( (-0.5, 0.5), polygon )
1151        assert not inside_polygon( (1.5, 0.5), polygon )
1152
1153        #Try point on borders
1154        assert inside_polygon( (1., 0.5), polygon, closed=True)
1155        assert inside_polygon( (0.5, 1), polygon, closed=True)
1156        assert inside_polygon( (0., 0.5), polygon, closed=True)
1157        assert inside_polygon( (0.5, 0.), polygon, closed=True)
1158
1159        assert not inside_polygon( (0.5, 1), polygon, closed=False)
1160        assert not inside_polygon( (0., 0.5), polygon, closed=False)
1161        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
1162        assert not inside_polygon( (1., 0.5), polygon, closed=False)
1163
1164
1165
1166        #From real example (that failed)
1167        polygon = [[20,20], [40,20], [40,40], [20,40]]
1168        points = [ [40, 50] ]
1169        res = inside_polygon(points, polygon)
1170        assert len(res) == 0
1171
1172        polygon = [[20,20], [40,20], [40,40], [20,40]]
1173        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
1174        res = inside_polygon(points, polygon)
1175        assert len(res) == 2
1176        assert allclose(res, [0,1])
1177
1178
1179
1180        #More convoluted and non convex polygon
1181        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1182        assert inside_polygon( (0.5, 0.5), polygon )
1183        assert inside_polygon( (1, -0.5), polygon )
1184        assert inside_polygon( (1.5, 0), polygon )
1185
1186        assert not inside_polygon( (0.5, 1.5), polygon )
1187        assert not inside_polygon( (0.5, -0.5), polygon )
1188
1189
1190        #Very convoluted polygon
1191        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1192        assert inside_polygon( (5, 5), polygon )
1193        assert inside_polygon( (17, 7), polygon )
1194        assert inside_polygon( (27, 2), polygon )
1195        assert inside_polygon( (35, -5), polygon )
1196        assert not inside_polygon( (15, 7), polygon )
1197        assert not inside_polygon( (24, 3), polygon )
1198        assert not inside_polygon( (25, -10), polygon )
1199
1200
1201
1202        #Another combination (that failed)
1203        polygon = [[0,0], [10,0], [10,10], [0,10]]
1204        assert inside_polygon( (5, 5), polygon )
1205        assert inside_polygon( (7, 7), polygon )
1206        assert not inside_polygon( (-17, 7), polygon )
1207        assert not inside_polygon( (7, 17), polygon )
1208        assert not inside_polygon( (17, 7), polygon )
1209        assert not inside_polygon( (27, 8), polygon )
1210        assert not inside_polygon( (35, -5), polygon )
1211
1212
1213
1214
1215        #Multiple polygons
1216
1217        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
1218                   [10,10], [11,10], [11,11], [10,11], [10,10]]
1219        assert inside_polygon( (0.5, 0.5), polygon )
1220        assert inside_polygon( (10.5, 10.5), polygon )
1221
1222        #FIXME: Fails if point is 5.5, 5.5
1223        assert not inside_polygon( (0, 5.5), polygon )
1224
1225        #Polygon with a hole
1226        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
1227                   [0,0], [1,0], [1,1], [0,1], [0,0]]
1228
1229        assert inside_polygon( (0, -0.5), polygon )
1230        assert not inside_polygon( (0.5, 0.5), polygon )
1231
1232    def test_inside_polygon_vector_version(self):
1233        #Now try the vector formulation returning indices
1234        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1235        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1236        res = inside_polygon( points, polygon, verbose=False )
1237
1238        assert allclose( res, [0,1,2] )
1239
1240    def test_outside_polygon(self):
1241        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1242
1243        assert not outside_polygon( [0.5, 0.5], U )
1244        #evaluate to False as the point 0.5, 0.5 is inside the unit square
1245       
1246        assert outside_polygon( [1.5, 0.5], U )
1247        #evaluate to True as the point 1.5, 0.5 is outside the unit square
1248       
1249        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1250        assert allclose( indices, [1] )
1251       
1252        #One more test of vector formulation returning indices
1253        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1254        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1255        res = outside_polygon( points, polygon )
1256
1257        assert allclose( res, [3, 4] )
1258
1259
1260
1261        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1262        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1263        res = outside_polygon( points, polygon )
1264
1265        assert allclose( res, [0, 4, 5] )       
1266     
1267    def test_outside_polygon2(self):
1268        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1269   
1270        assert not outside_polygon( [0.5, 1.0], U, closed = True )
1271        #evaluate to False as the point 0.5, 1.0 is inside the unit square
1272       
1273        assert outside_polygon( [0.5, 1.0], U, closed = False )
1274        #evaluate to True as the point 0.5, 1.0 is outside the unit square
1275
1276    def test_separate_points_by_polygon(self):
1277        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1278
1279        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1280        assert allclose( indices, [0,2,1] )
1281        assert count == 2
1282       
1283        #One more test of vector formulation returning indices
1284        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1285        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1286        res, count = separate_points_by_polygon( points, polygon )
1287
1288        assert allclose( res, [0,1,2,4,3] )
1289        assert count == 3
1290
1291
1292        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1293        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1294        res, count = separate_points_by_polygon( points, polygon )
1295
1296        assert allclose( res, [1,2,3,5,4,0] )       
1297        assert count == 3
1298       
1299
1300    def test_populate_polygon(self):
1301
1302        polygon = [[0,0], [1,0], [1,1], [0,1]]
1303        points = populate_polygon(polygon, 5)
1304
1305        assert len(points) == 5
1306        for point in points:
1307            assert inside_polygon(point, polygon)
1308
1309
1310        #Very convoluted polygon
1311        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1312
1313        points = populate_polygon(polygon, 5)
1314
1315        assert len(points) == 5
1316        for point in points:
1317            assert inside_polygon(point, polygon)
1318
1319
1320    def test_populate_polygon_with_exclude(self):
1321       
1322
1323        polygon = [[0,0], [1,0], [1,1], [0,1]]
1324        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
1325        points = populate_polygon(polygon, 5, exclude = [ex_poly])
1326
1327        assert len(points) == 5
1328        for point in points:
1329            assert inside_polygon(point, polygon)
1330            assert not inside_polygon(point, ex_poly)           
1331
1332
1333        #overlap
1334        polygon = [[0,0], [1,0], [1,1], [0,1]]
1335        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
1336        points = populate_polygon(polygon, 5, exclude = [ex_poly])
1337
1338        assert len(points) == 5
1339        for point in points:
1340            assert inside_polygon(point, polygon)
1341            assert not inside_polygon(point, ex_poly)                       
1342       
1343        #Multiple
1344        polygon = [[0,0], [1,0], [1,1], [0,1]]
1345        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
1346        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
1347       
1348        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
1349
1350        assert len(points) == 20
1351        for point in points:
1352            assert inside_polygon(point, polygon)
1353            assert not inside_polygon(point, ex_poly1)
1354            assert not inside_polygon(point, ex_poly2)                               
1355       
1356
1357        #Very convoluted polygon
1358        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1359        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
1360        points = populate_polygon(polygon, 20, exclude = [ex_poly])
1361       
1362        assert len(points) == 20
1363        for point in points:
1364            assert inside_polygon(point, polygon)
1365            assert not inside_polygon(point, ex_poly), '%s' %str(point)                       
1366
1367
1368
1369#-------------------------------------------------------------
1370if __name__ == "__main__":
1371    suite = unittest.makeSuite(Test_Util,'test')
1372    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
1373    runner = unittest.TextTestRunner()
1374    runner.run(suite)
1375
1376
1377
1378
Note: See TracBrowser for help on using the repository browser.