source: inundation/ga/storm_surge/pyvolution/test_least_squares.py @ 1697

Last change on this file since 1697 was 1670, checked in by ole, 20 years ago

Refactored Interpolation_function and tested
File size: 46.4 KB
RevLine 
[1018]1#!/usr/bin/env python
2
3#TEST
4
5import unittest
6from math import sqrt
7
8
9from least_squares import *
10from Numeric import allclose, array, transpose
11
[1097]12from coordinate_transforms.geo_reference import Geo_reference
13
[1018]14def distance(x, y):
15    return sqrt( sum( (array(x)-array(y))**2 ))
16
17def linear_function(point):
18    point = array(point)
19    return point[:,0]+point[:,1]
20
21
22class Test_Least_Squares(unittest.TestCase):
23
24    def setUp(self):
25        pass
26
27    def tearDown(self):
28        pass
29
30    def test_datapoint_at_centroid(self):
31        a = [0.0, 0.0]
32        b = [0.0, 2.0]
33        c = [2.0,0.0]
34        points = [a, b, c]
35        vertices = [ [1,0,2] ]   #bac
36
37        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
38
39        interp = Interpolation(points, vertices, data)
40        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
41
42
43    def test_quad_tree(self):
44        p0 = [-10.0, -10.0]
45        p1 = [20.0, -10.0]
46        p2 = [-10.0, 20.0]
47        p3 = [10.0, 50.0]
48        p4 = [30.0, 30.0]
49        p5 = [50.0, 10.0]
50        p6 = [40.0, 60.0]
51        p7 = [60.0, 40.0]
52        p8 = [-66.0, 20.0]
53        p9 = [10.0, -66.0]
54
55        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
56        triangles = [ [0, 1, 2],
57                      [3, 2, 4],
58                      [4, 2, 1],
59                      [4, 1, 5],
60                      [3, 4, 6],
61                      [6, 4, 7],
62                      [7, 4, 5],
63                      [8, 0, 2],
64                      [0, 9, 1]]
65
66        data = [ [4,4] ]
67        interp = Interpolation(points, triangles, data, alpha = 0.0,
68                               max_points_per_cell = 4)
69        #print "PDSG - interp.get_A()", interp.get_A()
70        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
71                      0., 0. , 0., 0., 0., 0.]]
72        assert allclose(interp.get_A(), answer)
73        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
74        #print "PDSG - interp.get_A()", interp.get_A()
75        answer =  [ [ 0.0,  0.0,  0.0,  0.,
76                      0., 0. , 0., 0., 0., 0.]]
77        assert allclose(interp.get_A(), answer)
78
79
80        #point outside of quad tree root cell
81        interp.set_point_coordinates([[-70, -70]])
82        #print "PDSG - interp.get_A()", interp.get_A()
83        answer =  [ [ 0.0,  0.0,  0.0,  0.,
84                      0., 0. , 0., 0., 0., 0.]]
85        assert allclose(interp.get_A(), answer)
86
87    def test_expand_search(self):
88        p0 = [-10.0, -10.0]
89        p1 = [20.0, -10.0]
90        p2 = [-10.0, 20.0]
91        p3 = [10.0, 50.0]
92        p4 = [30.0, 30.0]
93        p5 = [50.0, 10.0]
94        p6 = [40.0, 60.0]
95        p7 = [60.0, 40.0]
96        p8 = [-66.0, 20.0]
97        p9 = [10.0, -66.0]
98
99        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
100        triangles = [ [0, 1, 2],
101                      [3, 2, 4],
102                      [4, 2, 1],
103                      [4, 1, 5],
104                      [3, 4, 6],
105                      [6, 4, 7],
106                      [7, 4, 5],
107                      [8, 0, 2],
108                      [0, 9, 1]]
109
110        data = [ [4,4],
111                 [-30,10],
112                 [-20,0],
113                 [-20,10],
114                 [0,30],
115                 [10,-40],
116                 [10,-30],
117                 [10,-20],
118                 [10,10],
119                 [10,20],
120                 [10,30],
121                 [10,40],
122                 [20,10],
123                 [25,45],
124                 [30,0],
125                 [30,10],
126                 [30,30],
127                 [30,40],
128                 [30,50],
129                 [40,10],
130                 [40,30],
131                 [40,40],
132                 [40,50],
133                 [50,20],
134                 [50,30],
135                 [50,40],
136                 [50,50],
137                 [30,0],
138                 [-20,-20]]
139        point_attributes = [ -400000,
140                     10,
141                     10,
142                     10,
143                    10,
144                     10,
145                     10,
146                     10,
147                     10,
148                     10,
149                     10,
150                     10,
151                     10,
152                     10,
153                    10,
154                     10,
155                     10,
156                     10,
157                     10,
158                     10,
159                     10,
160                     10,
161                     10,
162                     10,
163                     10,
164                     10,
165                     10,
166                   10,
167                   99]
168
169        interp = Interpolation(points, triangles, data,
170                               alpha=0.0, expand_search=False, #verbose = True, #False,
171                               max_points_per_cell = 4)
172        calc = interp.fit_points(point_attributes, )
173        #print "calc",calc
174
175        # the point at 4,4 is ignored.  An expanded search has to be done
176        # to fine which triangel it's in.
177        # An expanded search isn't done to find that the last point
178        # isn't in the mesh.  But this isn't tested.
179        answer= [ 10,
180                     10,
181                     10,
182                     10,
183                    10,
184                     10,
185                     10,
186                     10,
187                     10,
188                     10]
189        assert allclose(calc, answer)
190
191    def test_quad_treeII(self):
192        p0 = [-66.0, 14.0]
193        p1 = [14.0, -66.0]
194        p2 = [14.0, 14.0]
195        p3 = [60.0, 20.0]
196        p4 = [10.0, 60.0]
197        p5 = [60.0, 60.0]
198
199        points = [p0, p1, p2, p3, p4, p5]
200        triangles = [ [0, 1, 2],
201                      [3, 2, 1],
202                      [0, 2, 4],
203                      [4, 2, 5],
204                      [5, 2, 3]]
205
206        data = [ [-26.0,-26.0] ]
207        interp = Interpolation(points, triangles, data, alpha = 0.0,
208                 max_points_per_cell = 4)
209        #print "PDSG - interp.get_A()", interp.get_A()
210        answer =  [ [ 0.5,  0.5,  0.0,  0.,
211                      0., 0.]]
212        assert allclose(interp.get_A(), answer)
213        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
214        #print "PDSG -30,-30 - interp.get_A()", interp.get_A()
215        answer =  [ [ 0.0,  0.0,  0.0,  0.,
216                      0., 0.]]
217        assert allclose(interp.get_A(), answer)
218
219
220        #point outside of quad tree root cell
221        interp.set_point_coordinates([[-70, -70]])
222        #print "PDSG -70,-70 interp.get_A()", interp.get_A()
223        answer =  [ [ 0.0,  0.0,  0.0,  0.,
224                      0., 0. ]]
225        assert allclose(interp.get_A(), answer)
226
227
228    def test_datapoints_at_vertices(self):
229        """Test that data points coinciding with vertices yield a diagonal matrix
230        """
231
232        a = [0.0, 0.0]
233        b = [0.0, 2.0]
234        c = [2.0,0.0]
235        points = [a, b, c]
236        vertices = [ [1,0,2] ]   #bac
237
238        data = points #Use data at vertices
239
240        interp = Interpolation(points, vertices, data)
241        assert allclose(interp.get_A(), [[1., 0., 0.],
242                                   [0., 1., 0.],
243                                   [0., 0., 1.]])
244
245
246
247    def test_datapoints_on_edge_midpoints(self):
248        """Try datapoints midway on edges -
249        each point should affect two matrix entries equally
250        """
251
252        a = [0.0, 0.0]
253        b = [0.0, 2.0]
254        c = [2.0,0.0]
255        points = [a, b, c]
256        vertices = [ [1,0,2] ]   #bac
257
258        data = [ [0., 1.], [1., 0.], [1., 1.] ]
259
260        interp = Interpolation(points, vertices, data)
261
262        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
263                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
264                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
265
266
267    def test_datapoints_on_edges(self):
268        """Try datapoints on edges -
269        each point should affect two matrix entries in proportion
270        """
271
272        a = [0.0, 0.0]
273        b = [0.0, 2.0]
274        c = [2.0,0.0]
275        points = [a, b, c]
276        vertices = [ [1,0,2] ]   #bac
277
278        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
279
280        interp = Interpolation(points, vertices, data)
281
282        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
283                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
284                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
285
286    def test_arbitrary_datapoints(self):
287        """Try arbitrary datapoints
288        """
289
290        from Numeric import sum
291
292        a = [0.0, 0.0]
293        b = [0.0, 2.0]
294        c = [2.0,0.0]
295        points = [a, b, c]
296        vertices = [ [1,0,2] ]   #bac
297
298        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
299
300        interp = Interpolation(points, vertices, data)
301        #print "interp.get_A()", interp.get_A()
302        assert allclose(sum(interp.get_A(), axis=1), 1.0)
303
304    def test_arbitrary_datapoints_some_outside(self):
305        """Try arbitrary datapoints one outside the triangle.
306        That one should be ignored
307        """
308
309        from Numeric import sum
310
311        a = [0.0, 0.0]
312        b = [0.0, 2.0]
313        c = [2.0,0.0]
314        points = [a, b, c]
315        vertices = [ [1,0,2] ]   #bac
316
317        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
318
319
320        interp = Interpolation(points, vertices, data, precrop = True)
321        assert allclose(sum(interp.get_A(), axis=1), 1.0)
322
323        interp = Interpolation(points, vertices, data, precrop = False)
324        assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
325
326
327
328    # this causes a memory error in scipy.sparse
329    def test_more_triangles(self):
330
331        a = [-1.0, 0.0]
332        b = [3.0, 4.0]
333        c = [4.0,1.0]
334        d = [-3.0, 2.0] #3
335        e = [-1.0,-2.0]
336        f = [1.0, -2.0] #5
337
338        points = [a, b, c, d,e,f]
339        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
340
341        #Data points
342        data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
343        interp = Interpolation(points, triangles, data_points)
344
345        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
346                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
347                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
348                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
349                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
350                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
351
352
353        A = interp.get_A()
354        for i in range(A.shape[0]):
355            for j in range(A.shape[1]):
356                if not allclose(A[i,j], answer[i][j]):
357                    print i,j,':',A[i,j], answer[i][j]
358
359
360        assert allclose(interp.get_A(), answer)
361
362
363
364
365    def test_smooth_attributes_to_mesh(self):
366        a = [0.0, 0.0]
367        b = [0.0, 5.0]
368        c = [5.0, 0.0]
369        points = [a, b, c]
370        triangles = [ [1,0,2] ]   #bac
371
372        d1 = [1.0, 1.0]
373        d2 = [1.0, 3.0]
374        d3 = [3.0,1.0]
375        z1 = 2
376        z2 = 4
377        z3 = 4
378        data_coords = [d1, d2, d3]
379
380        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
381        z = [z1, z2, z3]
382        f = interp.fit(z)
383        answer = [0, 5., 5.]
384
385        #print "f\n",f
386        #print "answer\n",answer
387
388        assert allclose(f, answer, atol=1e-7)
389
390
391    def test_smooth_att_to_meshII(self):
392
393        a = [0.0, 0.0]
394        b = [0.0, 5.0]
395        c = [5.0, 0.0]
396        points = [a, b, c]
397        triangles = [ [1,0,2] ]   #bac
398
399        d1 = [1.0, 1.0]
400        d2 = [1.0, 2.0]
401        d3 = [3.0,1.0]
402        data_coords = [d1, d2, d3]
403        z = linear_function(data_coords)
404        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
405        f = interp.fit(z)
406        answer = linear_function(points)
407
408        assert allclose(f, answer)
409
410    def test_smooth_attributes_to_meshIII(self):
411
412        a = [-1.0, 0.0]
413        b = [3.0, 4.0]
414        c = [4.0,1.0]
415        d = [-3.0, 2.0] #3
416        e = [-1.0,-2.0]
417        f = [1.0, -2.0] #5
418
419        vertices = [a, b, c, d,e,f]
420        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
421
422        point_coords = [[-2.0, 2.0],
423                        [-1.0, 1.0],
424                        [0.0,2.0],
425                        [1.0, 1.0],
426                        [2.0, 1.0],
427                        [0.0,0.0],
428                        [1.0, 0.0],
429                        [0.0, -1.0],
430                        [-0.2,-0.5],
431                        [-0.9, -1.5],
432                        [0.5, -1.9],
433                        [3.0,1.0]]
434
435        z = linear_function(point_coords)
436        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
437
438        #print 'z',z
439        f = interp.fit(z)
440        answer = linear_function(vertices)
441        #print "f\n",f
442        #print "answer\n",answer
443        assert allclose(f, answer)
444
445
446    def test_smooth_attributes_to_meshIV(self):
447        """ Testing 2 attributes smoothed to the mesh
448        """
449
450        a = [0.0, 0.0]
451        b = [0.0, 5.0]
452        c = [5.0, 0.0]
453        points = [a, b, c]
454        triangles = [ [1,0,2] ]   #bac
455
456        d1 = [1.0, 1.0]
457        d2 = [1.0, 3.0]
458        d3 = [3.0, 1.0]
459        z1 = [2, 4]
460        z2 = [4, 8]
461        z3 = [4, 8]
462        data_coords = [d1, d2, d3]
463
464        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
465        z = [z1, z2, z3]
466        f =  interp.fit_points(z)
467        answer = [[0,0], [5., 10.], [5., 10.]]
468        assert allclose(f, answer)
469
470    def test_interpolate_attributes_to_points(self):
471        v0 = [0.0, 0.0]
472        v1 = [0.0, 5.0]
473        v2 = [5.0, 0.0]
474
475        vertices = [v0, v1, v2]
476        triangles = [ [1,0,2] ]   #bac
477
478        d0 = [1.0, 1.0]
479        d1 = [1.0, 2.0]
480        d2 = [3.0, 1.0]
481        point_coords = [ d0, d1, d2]
482
483        interp = Interpolation(vertices, triangles, point_coords)
484        f = linear_function(vertices)
485        z = interp.interpolate(f)
486        answer = linear_function(point_coords)
487
488
489        assert allclose(z, answer)
490
491
492    def test_interpolate_attributes_to_pointsII(self):
493        a = [-1.0, 0.0]
494        b = [3.0, 4.0]
495        c = [4.0, 1.0]
496        d = [-3.0, 2.0] #3
497        e = [-1.0, -2.0]
498        f = [1.0, -2.0] #5
499
500        vertices = [a, b, c, d,e,f]
501        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
502
503
504        point_coords = [[-2.0, 2.0],
505                        [-1.0, 1.0],
506                        [0.0, 2.0],
507                        [1.0, 1.0],
508                        [2.0, 1.0],
509                        [0.0, 0.0],
510                        [1.0, 0.0],
511                        [0.0, -1.0],
512                        [-0.2, -0.5],
513                        [-0.9, -1.5],
514                        [0.5, -1.9],
515                        [3.0, 1.0]]
516
517        interp = Interpolation(vertices, triangles, point_coords)
518        f = linear_function(vertices)
519        z = interp.interpolate(f)
520        answer = linear_function(point_coords)
521        #print "z",z
522        #print "answer",answer
523        assert allclose(z, answer)
524
525    def test_interpolate_attributes_to_pointsIII(self):
526        """Test linear interpolation of known values at vertices to
527        new points inside a triangle
528        """
529        a = [0.0, 0.0]
530        b = [0.0, 5.0]
531        c = [5.0, 0.0]
532        d = [5.0, 5.0]
533
534        vertices = [a, b, c, d]
[1158]535        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
[1018]536
537        #Points within triangle 1
538        d0 = [1.0, 1.0]
539        d1 = [1.0, 2.0]
540        d2 = [3.0, 1.0]
541
542        #Point within triangle 2
543        d3 = [4.0, 3.0]
544
545        #Points on common edge
546        d4 = [2.5, 2.5]
547        d5 = [4.0, 1.0]
548
549        #Point on common vertex
550        d6 = [0., 5.]
551
552
553        point_coords = [d0, d1, d2, d3, d4, d5, d6]
554
555        interp = Interpolation(vertices, triangles, point_coords)
556
557        #Known values at vertices
558        #Functions are x+y, x+2y, 2x+y, x-y-5
559        f = [ [0., 0., 0., -5.],        # (0,0)
560              [5., 10., 5., -10.],      # (0,5)
561              [5., 5., 10.0, 0.],       # (5,0)
562              [10., 15., 15., -5.]]     # (5,5)
563
564        z = interp.interpolate(f)
565        answer = [ [2., 3., 3., -5.],   # (1,1)
566                   [3., 5., 4., -6.],   # (1,2)
567                   [4., 5., 7., -3.],   # (3,1)
568                   [7., 10., 11., -4.], # (4,3)
569                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
570                   [5., 6., 9., -2.],   # (4,1)
571                   [5., 10., 5., -10.]] # (0,5)
572
573        #print "***********"
574        #print "z",z
575        #print "answer",answer
576        #print "***********"
577
578        assert allclose(z, answer)
579
580    def test_interpolate_attributes_to_pointsIV(self):
581        a = [-1.0, 0.0]
582        b = [3.0, 4.0]
583        c = [4.0, 1.0]
584        d = [-3.0, 2.0] #3
585        e = [-1.0, -2.0]
586        f = [1.0, -2.0] #5
587
588        vertices = [a, b, c, d,e,f]
589        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
590
591
592        point_coords = [[-2.0, 2.0],
593                        [-1.0, 1.0],
594                        [0.0, 2.0],
595                        [1.0, 1.0],
596                        [2.0, 1.0],
597                        [0.0, 0.0],
598                        [1.0, 0.0],
599                        [0.0, -1.0],
600                        [-0.2, -0.5],
601                        [-0.9, -1.5],
602                        [0.5, -1.9],
603                        [3.0, 1.0]]
604
605        interp = Interpolation(vertices, triangles, point_coords)
606        f = array([linear_function(vertices),2*linear_function(vertices) ])
607        f = transpose(f)
608        #print "f",f
609        z = interp.interpolate(f)
610        answer = [linear_function(point_coords),
611                  2*linear_function(point_coords) ]
612        answer = transpose(answer)
613        #print "z",z
614        #print "answer",answer
615        assert allclose(z, answer)
616
617    def test_smooth_attributes_to_mesh_function(self):
618        """ Testing 2 attributes smoothed to the mesh
619        """
620
621        a = [0.0, 0.0]
622        b = [0.0, 5.0]
623        c = [5.0, 0.0]
624        points = [a, b, c]
625        triangles = [ [1,0,2] ]   #bac
626
627        d1 = [1.0, 1.0]
628        d2 = [1.0, 3.0]
629        d3 = [3.0, 1.0]
630        z1 = [2, 4]
631        z2 = [4, 8]
632        z3 = [4, 8]
633        data_coords = [d1, d2, d3]
634        z = [z1, z2, z3]
635
636        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
637        answer = [[0, 0], [5., 10.], [5., 10.]]
638
639        assert allclose(f, answer)
640
641
642
643    def test_pts2rectangular(self):
644
645        import time, os
646        FN = 'xyatest' + str(time.time()) + '.xya'
647        fid = open(FN, 'w')
648        fid.write(' %s \n' %('elevation'))
649        fid.write('%f %f %f\n' %(1,1,2) )
650        fid.write('%f %f %f\n' %(1,3,4) )
651        fid.write('%f %f %f\n' %(3,1,4) )
652        fid.close()
653
654        points, triangles, boundary, attributes =\
655                pts2rectangular(FN, 4, 4, format = 'asc')
656
657
658        data_coords = [ [1,1], [1,3], [3,1] ]
659        z = [2, 4, 4]
660
661        ref = fit_to_mesh(points, triangles, data_coords, z)
662
663        #print attributes
664        #print ref
665        assert allclose(attributes, ref)
666
667        os.remove(FN)
668
669
670    #Tests of smoothing matrix
671    def test_smoothing_matrix_one_triangle(self):
672        from Numeric import dot
673        a = [0.0, 0.0]
674        b = [0.0, 2.0]
675        c = [2.0,0.0]
676        points = [a, b, c]
677
678        vertices = [ [1,0,2] ]   #bac
679
680        interp = Interpolation(points, vertices)
681
682        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
683                                   [-0.5, 0.5, 0],
684                                   [-0.5, 0, 0.5]])
685
686        #Define f(x,y) = x
687        f = array([0,0,2]) #Value at global vertex 2
688
689        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
690        #           int 1 dx dy = area = 2
691        assert dot(dot(f, interp.get_D()), f) == 2
692
693        #Define f(x,y) = y
694        f = array([0,2,0])  #Value at global vertex 1
695
696        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
697        #           int 1 dx dy = area = 2
698        assert dot(dot(f, interp.get_D()), f) == 2
699
700        #Define f(x,y) = x+y
701        f = array([0,2,2])  #Values at global vertex 1 and 2
702
703        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
704        #           int 2 dx dy = 2*area = 4
705        assert dot(dot(f, interp.get_D()), f) == 4
706
707
708
709    def test_smoothing_matrix_more_triangles(self):
710        from Numeric import dot
711
712        a = [0.0, 0.0]
713        b = [0.0, 2.0]
714        c = [2.0,0.0]
715        d = [0.0, 4.0]
716        e = [2.0, 2.0]
717        f = [4.0,0.0]
718
719        points = [a, b, c, d, e, f]
720        #bac, bce, ecf, dbe, daf, dae
721        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
722
723        interp = Interpolation(points, vertices)
724
725
726        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
727        #                           [-0.5, 0.5, 0],
728        #                           [-0.5, 0, 0.5]])
729
730        #Define f(x,y) = x
731        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
732
733        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
734        #           int 1 dx dy = total area = 8
735        assert dot(dot(f, interp.get_D()), f) == 8
736
737        #Define f(x,y) = y
738        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
739
740        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
741        #           int 1 dx dy = area = 8
742        assert dot(dot(f, interp.get_D()), f) == 8
743
744        #Define f(x,y) = x+y
745        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
746
747        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
748        #           int 2 dx dy = 2*area = 16
749        assert dot(dot(f, interp.get_D()), f) == 16
750
751
752    def test_fit_and_interpolation(self):
753        from mesh import Mesh
754
755        a = [0.0, 0.0]
756        b = [0.0, 2.0]
757        c = [2.0, 0.0]
758        d = [0.0, 4.0]
759        e = [2.0, 2.0]
760        f = [4.0, 0.0]
761
762        points = [a, b, c, d, e, f]
763        #bac, bce, ecf, dbe, daf, dae
764        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
765
766        #Get (enough) datapoints
767        data_points = [[ 0.66666667, 0.66666667],
768                       [ 1.33333333, 1.33333333],
769                       [ 2.66666667, 0.66666667],
770                       [ 0.66666667, 2.66666667],
771                       [ 0.0, 1.0],
772                       [ 0.0, 3.0],
773                       [ 1.0, 0.0],
774                       [ 1.0, 1.0],
775                       [ 1.0, 2.0],
776                       [ 1.0, 3.0],
777                       [ 2.0, 1.0],
778                       [ 3.0, 0.0],
779                       [ 3.0, 1.0]]
780
781        interp = Interpolation(points, triangles, data_points, alpha=0.0)
782
783        z = linear_function(data_points)
784        answer = linear_function(points)
785
786        f = interp.fit(z)
787
788        #print "f",f
789        #print "answer",answer
790        assert allclose(f, answer)
791
792        #Map back
793        z1 = interp.interpolate(f)
794        #print "z1\n", z1
795        #print "z\n",z
796        assert allclose(z, z1)
797
798
799    def test_smoothing_and_interpolation(self):
800
801        a = [0.0, 0.0]
802        b = [0.0, 2.0]
803        c = [2.0, 0.0]
804        d = [0.0, 4.0]
805        e = [2.0, 2.0]
806        f = [4.0, 0.0]
807
808        points = [a, b, c, d, e, f]
809        #bac, bce, ecf, dbe, daf, dae
810        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
811
812        #Get (too few!) datapoints
813        data_points = [[ 0.66666667, 0.66666667],
814                       [ 1.33333333, 1.33333333],
815                       [ 2.66666667, 0.66666667],
816                       [ 0.66666667, 2.66666667]]
817
818        z = linear_function(data_points)
819        answer = linear_function(points)
820
821        #Make interpolator with too few data points and no smoothing
822        interp = Interpolation(points, triangles, data_points, alpha=0.0)
823        #Must raise an exception
824        try:
825            f = interp.fit(z)
826        except:
827            pass
828
829        #Now try with smoothing parameter
830        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
831
832        f = interp.fit(z)
[1657]833        #f will be different from answer due to smoothing
[1018]834        assert allclose(f, answer,atol=5)
835
836        #Map back
837        z1 = interp.interpolate(f)
838        assert allclose(z, z1)
839
840
841
842    def test_fit_and_interpolation_with_new_points(self):
843        """Fit a surface to one set of points. Then interpolate that surface
844        using another set of points.
845        """
846        from mesh import Mesh
847
848
849        #Setup mesh used to represent fitted function
850        a = [0.0, 0.0]
851        b = [0.0, 2.0]
852        c = [2.0, 0.0]
853        d = [0.0, 4.0]
854        e = [2.0, 2.0]
855        f = [4.0, 0.0]
856
857        points = [a, b, c, d, e, f]
858        #bac, bce, ecf, dbe, daf, dae
859        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
860
861        #Datapoints to fit from
862        data_points1 = [[ 0.66666667, 0.66666667],
863                        [ 1.33333333, 1.33333333],
864                        [ 2.66666667, 0.66666667],
865                        [ 0.66666667, 2.66666667],
866                        [ 0.0, 1.0],
867                        [ 0.0, 3.0],
868                        [ 1.0, 0.0],
869                        [ 1.0, 1.0],
870                        [ 15, -17],   #Outside mesh
871                        [ 1.0, 2.0],
872                        [ 1.0, 3.0],
873                        [ 2.0, 1.0],
874                        [ 3.0, 0.0],
875                        [ 3.0, 1.0]]
876
877        #Fit surface to mesh
878        interp = Interpolation(points, triangles, data_points1, alpha=0.0,
879                               precrop = True)
880        z = linear_function(data_points1) #Example z-values
881        f = interp.fit(z)                 #Fitted values at vertices
882
883
884
885        #New datapoints where interpolated values are sought
886        data_points2 = [[ 0.0, 0.0],
887                        [ 0.5, 0.5],
888                        [ 0.7, 0.7],
889                        [-13, 65],  #Outside
890                        [ 1.0, 0.5],
891                        [ 2.0, 0.4],
892                        [ 2.8, 1.2]]
893
894
895
896        #Build new A matrix based on new points (without precrop)
897        interp.build_interpolation_matrix_A(data_points2, precrop = False)
898
899        #Interpolate using fitted surface
900        z1 = interp.interpolate(f)
901
902        #import Numeric
903        #data_points2 = Numeric.take(data_points2, interp.point_indices)
904
905        #Desired result (OK for points inside)
906
907        answer = linear_function(data_points2)
908        import Numeric
909        z1 = Numeric.take(z1, [0,1,2,4,5,6])
910        answer = Numeric.take(answer, [0,1,2,4,5,6])
911        assert allclose(z1, answer)
912
913        #Build new A matrix based on new points (with precrop)
914        interp.build_interpolation_matrix_A(data_points2, precrop = True)
915
916        #Interpolate using fitted surface
917        z1 = interp.interpolate(f)
918
919        import Numeric
920        data_points2 = Numeric.take(data_points2, interp.point_indices)
921
922        #Desired result
923        answer = linear_function(data_points2)
924        assert allclose(z1, answer)
925
926
[1657]927
928
929
930
[1632]931    def test_interpolation_from_discontinuous_vertex_values(self):
932        """test_interpolation_from_discontinuous_vertex_values.
933        This will test the format used internally in pyvolution and also
934        interpolation from sww files
935        """
936       
937        from mesh import Mesh
[1018]938
[1632]939
940        #Setup mesh used to represent discontinuous function
941        a = [0.0, 0.0]
942        b = [0.0, 2.0]
943        c = [2.0, 0.0]
944        d = [0.0, 4.0]
945        e = [2.0, 2.0]
946        f = [4.0, 0.0]
947
948        points = [b, a, c,
949                  b, c, e,
950                  e, c, f,
951                  d, b, e]
952       
953        #bac, bce, ecf, dbe
954        triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
955
956       
957        vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
958                         
959       
960
961        #New datapoints where interpolated values are sought
962        data_points = [[0.0, 0.0],  #T0
963                       [0.5, 0.5],  #T0
964                       [1.5, 1.5],  #T1
965                       [2.5, 0.5],  #T2
966                       [0.0, 3.0],  #T3
967                       [1.0, 2.0],  #In between T1 and T3 (T1 is used) FIXME?
968                       [2.0, 1.0],  #In between T1 and T2 (T1 is used) FIXME?
969                       [1.0, 1.0]]  #In between T1 and T0 (T0 is used) FIXME?
970       
971
972
973
974        #Build interpolation matrix
975        interp = Interpolation(points, triangles, data_points)
976                               #, alpha=0.0, precrop = True)
977
978        #print interp.A.todense()
979        #print vertex_values
980
981        #Interpolate using fitted surface
982        z = interp.interpolate(vertex_values)
983
984        #print z
985
986        assert allclose(z, [0,0,1,2,5,1,1,0])
987
988
989
990
[1670]991    def test_interpolation_function_time_only(self):
992        """Test spatio-temporal interpolation
993        Test that spatio temporal function performs the correct
994        interpolations in both time and space
995        """
[1632]996
997
[1670]998        #Three timesteps
999        time = [1.0, 5.0, 6.0]
1000       
[1632]1001
[1670]1002        #One quantity
1003        Q = zeros( (3,6), Float )
[1632]1004
[1670]1005        #Linear in time and space
1006        a = [0.0, 0.0]
1007        b = [0.0, 2.0]
1008        c = [2.0, 0.0]
1009        d = [0.0, 4.0]
1010        e = [2.0, 2.0]
1011        f = [4.0, 0.0]
1012
1013        points = [a, b, c, d, e, f]
1014       
1015        for i, t in enumerate(time):
1016            Q[i, :] = t*linear_function(points)
1017
1018           
1019        #Check basic interpolation of one quantity using averaging
1020        #(no interpolation points or spatial info)
1021        from util import mean       
1022        I = Interpolation_function(time, [mean(Q[0,:]),
1023                                          mean(Q[1,:]),
1024                                          mean(Q[2,:])])
1025
1026
1027
1028        #Check temporal interpolation
1029        for i in [0,1,2]:
1030            assert allclose(I(time[i]), mean(Q[i,:]))
1031
1032        #Midway   
1033        assert allclose(I( (time[0] + time[1])/2 ),
1034                        (I(time[0]) + I(time[1]))/2 )
1035
1036        assert allclose(I( (time[1] + time[2])/2 ),
1037                        (I(time[1]) + I(time[2]))/2 )
1038
1039        assert allclose(I( (time[0] + time[2])/2 ),
1040                        (I(time[0]) + I(time[2]))/2 )                 
1041
1042        #1/3
1043        assert allclose(I( (time[0] + time[2])/3 ),
1044                        (I(time[0]) + I(time[2]))/3 )                         
1045
1046
1047        #Out of bounds checks
1048        try:
1049            I(time[0]-1) 
1050        except:
1051            pass
1052        else:
1053            raise 'Should raise exception'
1054
1055        try:
1056            I(time[-1]+1) 
1057        except:
1058            pass
1059        else:
1060            raise 'Should raise exception'       
1061
1062
1063       
1064
1065    def test_interpolation_function(self):
1066        """Test spatio-temporal interpolation
1067        Test that spatio temporal function performs the correct
1068        interpolations in both time and space
1069        """
1070
1071
1072        #Three timesteps
1073        time = [1.0, 5.0, 6.0]
1074       
1075
1076        #Setup mesh used to represent fitted function
1077        a = [0.0, 0.0]
1078        b = [0.0, 2.0]
1079        c = [2.0, 0.0]
1080        d = [0.0, 4.0]
1081        e = [2.0, 2.0]
1082        f = [4.0, 0.0]
1083
1084        points = [a, b, c, d, e, f]
1085        #bac, bce, ecf, dbe
1086        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1087
1088
1089        #New datapoints where interpolated values are sought
1090        interpolation_points = [[ 0.0, 0.0],
1091                                [ 0.5, 0.5],
1092                                [ 0.7, 0.7],
1093                                [ 1.0, 0.5],
1094                                [ 2.0, 0.4],
1095                                [ 2.8, 1.2]]
1096
1097
1098        #One quantity
1099        Q = zeros( (3,6), Float )
1100
1101        #Linear in time and space
1102        for i, t in enumerate(time):
1103            Q[i, :] = t*linear_function(points)
1104
1105
1106        #Check interpolation of one quantity using interpolaton points)
1107        I = Interpolation_function(time, Q,
1108                                   vertex_coordinates = points,
1109                                   triangles = triangles, 
1110                                   interpolation_points = interpolation_points,
1111                                   verbose = False)
1112
1113
1114        answer = linear_function(interpolation_points)
1115
1116        t = time[0]
1117        for j in range(50): #t in [1, 6]
1118            for id in range(len(interpolation_points)):
1119                assert allclose(I(t, id), t*answer[id])
1120
1121            t += 0.1   
1122
1123
1124        try:   
1125            I(1)
1126        except:
1127            pass
1128        else:
1129            raise 'Should raise exception'
1130           
1131        #
1132        #interpolation_points = [[ 0.0, 0.0],
1133        #                        [ 0.5, 0.5],
1134        #                        [ 0.7, 0.7],
1135        #                        [-13, 65],  #Outside
1136        #                        [ 1.0, 0.5],
1137        #                        [ 2.0, 0.4],
1138        #                        [ 2.8, 1.2]]
1139        #
1140        #try:
1141        #    I = Interpolation_function(time, Q,
1142        #                               vertex_coordinates = points,
1143        #                               triangles = triangles,
1144        #                               interpolation_points = interpolation_points,
1145        #                               verbose = False)
1146        #except:
1147        #    pass
1148        #else:
1149        #    raise 'Should raise exception'
1150
1151
1152
1153
[1018]1154    def test_fit_and_interpolation_with_different_origins(self):
1155        """Fit a surface to one set of points. Then interpolate that surface
1156        using another set of points.
1157        This test tests situtaion where points and mesh belong to a different
1158        coordinate system as defined by origin.
1159        """
1160        from mesh import Mesh
1161
1162        #Setup mesh used to represent fitted function
1163        a = [0.0, 0.0]
1164        b = [0.0, 2.0]
1165        c = [2.0, 0.0]
1166        d = [0.0, 4.0]
1167        e = [2.0, 2.0]
1168        f = [4.0, 0.0]
1169
1170        points = [a, b, c, d, e, f]
1171        #bac, bce, ecf, dbe, daf, dae
1172        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1173
1174        #Datapoints to fit from
1175        data_points1 = [[ 0.66666667, 0.66666667],
1176                        [ 1.33333333, 1.33333333],
1177                        [ 2.66666667, 0.66666667],
1178                        [ 0.66666667, 2.66666667],
1179                        [ 0.0, 1.0],
1180                        [ 0.0, 3.0],
1181                        [ 1.0, 0.0],
1182                        [ 1.0, 1.0],
1183                        [ 1.0, 2.0],
1184                        [ 1.0, 3.0],
1185                        [ 2.0, 1.0],
1186                        [ 3.0, 0.0],
1187                        [ 3.0, 1.0]]
1188
1189
1190        #First check that things are OK when using same origin
1191        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1192        data_origin = (56, 290000, 618000) #zone, easting, northing
1193
1194
1195        #Fit surface to mesh
1196        interp = Interpolation(points, triangles, data_points1,
1197                               alpha=0.0,
1198                               data_origin = data_origin,
1199                               mesh_origin = mesh_origin)
1200
1201        z = linear_function(data_points1) #Example z-values
1202        f = interp.fit(z)                 #Fitted values at vertices
1203
1204
1205        #New datapoints where interpolated values are sought
1206        data_points2 = [[ 0.0, 0.0],
1207                        [ 0.5, 0.5],
1208                        [ 0.7, 0.7],
1209                        [ 1.0, 0.5],
1210                        [ 2.0, 0.4],
1211                        [ 2.8, 1.2]]
1212
1213
1214        #Build new A matrix based on new points
1215        interp.build_interpolation_matrix_A(data_points2)
1216
1217        #Interpolate using fitted surface
1218        z1 = interp.interpolate(f)
1219
1220        #Desired result
1221        answer = linear_function(data_points2)
1222        assert allclose(z1, answer)
1223
1224
1225        ##############################################
1226
[1096]1227        #Then check situation where points are relative to a different
[1018]1228        #origin (same zone, though, until we figure that out (FIXME))
1229
1230        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1231        data_origin = (56, 10000, 10000) #zone, easting, northing
1232
1233        #Shift datapoints according to new origin
1234
1235        for k in range(len(data_points1)):
1236            data_points1[k][0] += mesh_origin[1] - data_origin[1]
1237            data_points1[k][1] += mesh_origin[2] - data_origin[2]
1238
1239        for k in range(len(data_points2)):
1240            data_points2[k][0] += mesh_origin[1] - data_origin[1]
1241            data_points2[k][1] += mesh_origin[2] - data_origin[2]
1242
1243
1244
1245        #Fit surface to mesh
1246        interp = Interpolation(points, triangles, data_points1,
1247                               alpha=0.0,
1248                               data_origin = data_origin,
1249                               mesh_origin = mesh_origin)
1250
1251        f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
1252
1253        assert allclose(f,f1), 'Fit should have been unaltered'
1254
1255
1256        #Build new A matrix based on new points
1257        interp.build_interpolation_matrix_A(data_points2)
1258
1259        #Interpolate using fitted surface
1260        z1 = interp.interpolate(f)
1261        assert allclose(z1, answer)
1262
1263
1264        #########################################################
1265        #Finally try to relate data_points2 to new origin without
1266        #rebuilding matrix
1267
1268        data_origin = (56, 2000, 2000) #zone, easting, northing
1269        for k in range(len(data_points2)):
1270            data_points2[k][0] += 8000
1271            data_points2[k][1] += 8000
1272
1273        #Build new A matrix based on new points
1274        interp.build_interpolation_matrix_A(data_points2,
1275                                            data_origin = data_origin)
1276
1277        #Interpolate using fitted surface
1278        z1 = interp.interpolate(f)
1279        assert allclose(z1, answer)
1280
1281
1282    def test_fit_to_mesh_file(self):
[1423]1283        from load_mesh.loadASCII import import_mesh_file, \
[1018]1284             export_mesh_file
1285        import tempfile
1286        import os
1287
1288        # create a .tsh file, no user outline
1289        mesh_dic = {}
1290        mesh_dic['vertices'] = [[0.0, 0.0],
1291                                          [0.0, 5.0],
1292                                          [5.0, 0.0]]
1293        mesh_dic['triangles'] =  [[0, 2, 1]]
1294        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[1183]1295        mesh_dic['triangle_tags'] = ['']
[1018]1296        mesh_dic['vertex_attributes'] = [[], [], []]
1297        mesh_dic['vertiex_attribute_titles'] = []
1298        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1299        mesh_dic['segment_tags'] = ['external',
1300                                                  'external',
1301                                                  'external']
1302        mesh_file = tempfile.mktemp(".tsh")
1303        export_mesh_file(mesh_file,mesh_dic)
1304
1305        # create an .xya file
1306        point_file = tempfile.mktemp(".xya")
1307        fd = open(point_file,'w')
1308        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
1309        fd.close()
1310
[1469]1311        mesh_output_file = tempfile.mktemp(".tsh") 
[1018]1312        fit_to_mesh_file(mesh_file,
1313                         point_file,
1314                         mesh_output_file,
1315                         alpha = 0.0)
1316        # load in the .tsh file we just wrote
[1423]1317        mesh_dic = import_mesh_file(mesh_output_file)
[1018]1318        #print "mesh_dic",mesh_dic
1319        ans =[[0.0, 0.0],
1320              [5.0, 10.0],
1321              [5.0,10.0]]
1322        assert allclose(mesh_dic['vertex_attributes'],ans)
1323
1324        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1325                        ['elevation','stage'],
1326                        'test_fit_to_mesh_file failed')
1327
1328        #clean up
1329        os.remove(mesh_file)
1330        os.remove(point_file)
1331        os.remove(mesh_output_file)
1332
[1097]1333    def test_fit_to_mesh_file3(self):
[1423]1334        from load_mesh.loadASCII import import_mesh_file, \
[1097]1335             export_mesh_file
1336        import tempfile
1337        import os
1338
1339        # create a .tsh file, no user outline
1340        mesh_dic = {}
1341        mesh_dic['vertices'] = [[0.76, 0.76],
1342                                          [0.76, 5.76],
1343                                          [5.76, 0.76]]
1344        mesh_dic['triangles'] =  [[0, 2, 1]]
1345        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[1183]1346        mesh_dic['triangle_tags'] = ['']
[1097]1347        mesh_dic['vertex_attributes'] = [[], [], []]
1348        mesh_dic['vertiex_attribute_titles'] = []
1349        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1350        mesh_dic['segment_tags'] = ['external',
1351                                                  'external',
1352                                                  'external']
1353        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
1354        mesh_file = tempfile.mktemp(".tsh")
1355        export_mesh_file(mesh_file,mesh_dic)
1356
1357        #FIXME - make this test the georef in the points file as well.
1358        # create an .xya file
1359        point_file = tempfile.mktemp(".xya")
1360        fd = open(point_file,'w')
1361        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
1362        fd.close()
1363
[1469]1364        mesh_output_file = tempfile.mktemp(".tsh")
[1097]1365        fit_to_mesh_file(mesh_file,
1366                         point_file,
1367                         mesh_output_file,
1368                         alpha = 0.0)
1369        # load in the .tsh file we just wrote
[1423]1370        mesh_dic = import_mesh_file(mesh_output_file)
[1097]1371        #print "mesh_dic",mesh_dic
1372        ans =[[0.0, 0.0],
1373              [5.0, 10.0],
1374              [5.0,10.0]]
1375        assert allclose(mesh_dic['vertex_attributes'],ans)
1376
1377        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1378                        ['elevation','stage'],
1379                        'test_fit_to_mesh_file failed')
1380
1381        #clean up
1382        os.remove(mesh_file)
1383        os.remove(point_file)
1384        os.remove(mesh_output_file)
1385
[1018]1386    def test_fit_to_mesh_fileII(self):
[1423]1387        from load_mesh.loadASCII import import_mesh_file, \
[1018]1388             export_mesh_file
1389        import tempfile
1390        import os
1391
1392        # create a .tsh file, no user outline
1393        mesh_dic = {}
1394        mesh_dic['vertices'] = [[0.0, 0.0],
[1632]1395                                [0.0, 5.0],
1396                                [5.0, 0.0]]
[1018]1397        mesh_dic['triangles'] =  [[0, 2, 1]]
1398        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[1183]1399        mesh_dic['triangle_tags'] = ['']
[1018]1400        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1401        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1402        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1403        mesh_dic['segment_tags'] = ['external',
1404                                                  'external',
1405                                                  'external']
1406        mesh_file = tempfile.mktemp(".tsh")
1407        export_mesh_file(mesh_file,mesh_dic)
1408
1409        # create an .xya file
1410        point_file = tempfile.mktemp(".xya")
1411        fd = open(point_file,'w')
1412        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
1413        fd.close()
1414
1415        mesh_output_file = "new_triangle.tsh"
1416        fit_to_mesh_file(mesh_file,
1417                         point_file,
1418                         mesh_output_file,
1419                         alpha = 0.0)
1420        # load in the .tsh file we just wrote
[1423]1421        mesh_dic = import_mesh_file(mesh_output_file)
[1018]1422
1423        assert allclose(mesh_dic['vertex_attributes'],
1424                        [[1.0, 2.0,0.0, 0.0],
1425                         [1.0, 2.0,5.0, 10.0],
1426                         [1.0, 2.0,5.0,10.0]])
1427
1428        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1429                        ['density', 'temp','elevation','stage'],
1430                        'test_fit_to_mesh_file failed')
1431
1432        #clean up
1433        os.remove(mesh_file)
1434        os.remove(mesh_output_file)
1435        os.remove(point_file)
1436
1437    def test_fit_to_msh_netcdf_fileII(self):
[1632]1438        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
[1018]1439        import tempfile
1440        import os
1441
1442        # create a .tsh file, no user outline
1443        mesh_dic = {}
1444        mesh_dic['vertices'] = [[0.0, 0.0],
[1632]1445                                [0.0, 5.0],
1446                                [5.0, 0.0]]
[1018]1447        mesh_dic['triangles'] =  [[0, 2, 1]]
1448        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[1183]1449        mesh_dic['triangle_tags'] = ['']
[1018]1450        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1451        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1452        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1453        mesh_dic['segment_tags'] = ['external',
1454                                                  'external',
1455                                                  'external']
1456        mesh_file = tempfile.mktemp(".msh")
1457        export_mesh_file(mesh_file,mesh_dic)
1458
1459        # create an .xya file
1460        point_file = tempfile.mktemp(".xya")
1461        fd = open(point_file,'w')
1462        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
1463        fd.close()
1464
1465        mesh_output_file = "new_triangle.msh"
1466        fit_to_mesh_file(mesh_file,
1467                         point_file,
1468                         mesh_output_file,
1469                         alpha = 0.0)
1470        # load in the .tsh file we just wrote
[1423]1471        mesh_dic = import_mesh_file(mesh_output_file)
[1018]1472
1473        assert allclose(mesh_dic['vertex_attributes'],
1474                        [[1.0, 2.0,0.0, 0.0],
1475                         [1.0, 2.0,5.0, 10.0],
1476                         [1.0, 2.0,5.0,10.0]])
1477
1478        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
1479                        ['density', 'temp','elevation','stage'],
1480                        'test_fit_to_mesh_file failed')
1481
1482        #clean up
1483        os.remove(mesh_file)
1484        os.remove(mesh_output_file)
1485        os.remove(point_file)
1486
1487#-------------------------------------------------------------
1488if __name__ == "__main__":
1489    suite = unittest.makeSuite(Test_Least_Squares,'test')
1490
1491    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
1492    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
1493    runner = unittest.TextTestRunner(verbosity=1)
1494    runner.run(suite)
1495
1496
1497
1498
1499
Note: See TracBrowser for help on using the repository browser.