source: anuga_core/source/anuga/fit_interpolate/test_fit.py @ 7276

Last change on this file since 7276 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

File size: 33.2 KB
Line 
1#!/usr/bin/env python
2
3"""
4Note, fitting/blocking is also tested in
5test_quantity.test_set_values_from_UTM_pts.
6"""
7
8#TEST
9import sys
10import unittest
11from math import sqrt
12import tempfile
13import os
14
15from fit import *
16from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
17from anuga.utilities.sparse import Sparse, Sparse_CSR
18from anuga.coordinate_transforms.geo_reference import Geo_reference
19from anuga.utilities.numerical_tools import ensure_numeric
20from anuga.geospatial_data.geospatial_data import Geospatial_data
21from anuga.shallow_water import Domain
22
23import numpy as num
24
25
26def distance(x, y):
27    return sqrt(num.sum((num.array(x)-num.array(y))**2))
28
29def linear_function(point):
30    point = num.array(point)
31    return point[:,0]+point[:,1]
32
33
34class Test_Fit(unittest.TestCase):
35
36    def setUp(self):
37        pass
38
39    def tearDown(self):
40        pass
41
42
43    def test_smooth_attributes_to_mesh(self):
44        a = [0.0, 0.0]
45        b = [0.0, 5.0]
46        c = [5.0, 0.0]
47        points = [a, b, c]
48        triangles = [ [1,0,2] ]   #bac
49
50        d1 = [1.0, 1.0]
51        d2 = [1.0, 3.0]
52        d3 = [3.0,1.0]
53        z1 = 2
54        z2 = 4
55        z3 = 4
56        data_coords = [d1, d2, d3]
57
58        z = [z1, z2, z3]
59        fit = Fit(points, triangles, alpha=0)
60        #print "interp.get_A()", interp.get_A()
61        fit._build_matrix_AtA_Atz(ensure_numeric(data_coords),
62                                  ensure_numeric(z))
63        #print "Atz - from fit", fit.Atz
64        #print "AtA - from fit", fit.AtA.todense()
65        #print "z",z
66
67        assert num.allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7)
68
69        f = fit.fit()
70       
71        answer = [0, 5., 5.]
72
73        #print "f\n",f
74        #print "answer\n",answer
75
76        assert num.allclose(f, answer, atol=1e-7)
77
78    def test_smooth_att_to_meshII(self):
79
80        a = [0.0, 0.0]
81        b = [0.0, 5.0]
82        c = [5.0, 0.0]
83        points = [a, b, c]
84        triangles = [ [1,0,2] ]   #bac
85
86        d1 = [1.0, 1.0]
87        d2 = [1.0, 2.0]
88        d3 = [3.0,1.0]
89        data_coords = [d1, d2, d3]
90        z = linear_function(data_coords)
91        #print "z",z
92
93        interp = Fit(points, triangles, alpha=0.0)
94        f = interp.fit(data_coords, z)
95        answer = linear_function(points)
96        #print "f\n",f
97        #print "answer\n",answer
98
99        assert num.allclose(f, answer)
100
101    def test_smooth_attributes_to_meshIII(self):
102
103        a = [-1.0, 0.0]
104        b = [3.0, 4.0]
105        c = [4.0,1.0]
106        d = [-3.0, 2.0] #3
107        e = [-1.0,-2.0]
108        f = [1.0, -2.0] #5
109
110        vertices = [a, b, c, d,e,f]
111        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
112
113        point_coords = [[-2.0, 2.0],
114                        [-1.0, 1.0],
115                        [0.0,2.0],
116                        [1.0, 1.0],
117                        [2.0, 1.0],
118                        [0.0,0.0],
119                        [1.0, 0.0],
120                        [0.0, -1.0],
121                        [-0.2,-0.5],
122                        [-0.9, -1.5],
123                        [0.5, -1.9],
124                        [3.0,1.0]]
125
126        z = linear_function(point_coords)
127        interp = Fit(vertices, triangles,
128                                alpha=0.0)
129
130        #print 'z',z
131        f = interp.fit(point_coords,z)
132        answer = linear_function(vertices)
133        #print "f\n",f
134        #print "answer\n",answer
135        assert num.allclose(f, answer)
136
137   
138    def test_smooth_attributes_to_meshIV(self):
139        # Testing 2 attributes smoothed to the mesh
140       
141
142        a = [0.0, 0.0]
143        b = [0.0, 5.0]
144        c = [5.0, 0.0]
145        points = [a, b, c]
146        triangles = [ [1,0,2] ]   #bac
147
148        d1 = [1.0, 1.0]
149        d2 = [1.0, 3.0]
150        d3 = [3.0, 1.0]
151        z1 = [2, 4]
152        z2 = [4, 8]
153        z3 = [4, 8]
154        data_coords = [d1, d2, d3]
155
156        z = [z1, z2, z3]
157        fit = Fit(points, triangles, alpha=0)
158       
159        f =  fit.fit(data_coords,z)
160        answer = [[0,0], [5., 10.], [5., 10.]]
161        assert num.allclose(f, answer)
162
163    def test_smooth_attributes_to_mesh_build_fit_subset(self):
164
165        a = [-1.0, 0.0]
166        b = [3.0, 4.0]
167        c = [4.0,1.0]
168        d = [-3.0, 2.0] #3
169        e = [-1.0,-2.0]
170        f = [1.0, -2.0] #5
171
172        vertices = [a, b, c, d,e,f]
173        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
174
175        interp = Fit(vertices, triangles,
176                                alpha=0.0)
177
178        point_coords = [[-2.0, 2.0],
179                        [-1.0, 1.0],
180                        [0.0,2.0],
181                        [1.0, 1.0],
182                       ]
183
184        z = linear_function(point_coords)
185
186        f = interp.build_fit_subset(point_coords,z)
187       
188        point_coords = [
189                        [2.0, 1.0],
190                        [0.0,0.0],
191                        [1.0, 0.0],
192                        [0.0, -1.0],
193                        [-0.2,-0.5],
194                        [-0.9, -1.5],
195                        [0.5, -1.9],
196                        [3.0,1.0]]
197       
198        z = linear_function(point_coords)
199
200        f = interp.build_fit_subset(point_coords,z)
201       
202        #print 'z',z
203        f = interp.fit()
204        answer = linear_function(vertices)
205        #print "f\n",f
206        #print "answer\n",answer
207        assert num.allclose(f, answer)
208
209        # test fit 2 mesh as well.
210    def test_fit_file_blocking(self):
211
212        a = [-1.0, 0.0]
213        b = [3.0, 4.0]
214        c = [4.0,1.0]
215        d = [-3.0, 2.0] #3
216        e = [-1.0,-2.0]
217        f = [1.0, -2.0] #5
218
219        vertices = [a, b, c, d,e,f]
220        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
221
222        interp = Fit(vertices, triangles,
223                     alpha=0.0)
224
225
226        fileName = tempfile.mktemp(".ddd")
227        file = open(fileName,"w")
228        file.write(" x,y, elevation \n\
229-2.0, 2.0, 0.\n\
230-1.0, 1.0, 0.\n\
2310.0, 2.0 , 2.\n\
2321.0, 1.0 , 2.\n\
233 2.0,  1.0 ,3. \n\
234 0.0,  0.0 , 0.\n\
235 1.0,  0.0 , 1.\n\
236 0.0,  -1.0, -1.\n\
237 -0.2, -0.5, -0.7\n\
238 -0.9, -1.5, -2.4\n\
239 0.5,  -1.9, -1.4\n\
240 3.0,  1.0 , 4.\n")
241        file.close()
242       
243        f = interp.fit(fileName, max_read_lines=2)
244        answer = linear_function(vertices)
245        #print "f\n",f
246        #print "answer\n",answer
247        assert num.allclose(f, answer)
248        os.remove(fileName)
249
250    def test_fit_to_mesh_UTM_file(self):
251        #Get (enough) datapoints
252        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
253                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
254
255        data_geo_spatial = Geospatial_data(data_points,
256                                           points_are_lats_longs=True)
257        points_UTM = data_geo_spatial.get_data_points(absolute=True)
258        attributes = linear_function(points_UTM)
259        att = 'elevation'
260       
261        #Create .txt file
262        txt_file = tempfile.mktemp(".txt")
263        file = open(txt_file,"w")
264        file.write(" x,y," + att + " \n")
265        for data_point, attribute in map(None, points_UTM, attributes):
266            row = str(data_point[0]) + ',' + str(data_point[1]) \
267                  + ',' + str(attribute)
268            #print "row", row
269            file.write(row + "\n")
270        file.close()
271
272        # setting up the mesh
273        a = [240000, 7620000]
274        b = [240000, 7680000]
275        c = [300000, 7620000]
276        points = [a, b, c]
277        elements = [[0,2,1]]
278        f = fit_to_mesh(txt_file, points, elements, 
279                        alpha=0.0, max_read_lines=2)
280        answer = linear_function(points)
281        #print "f",f
282        #print "answer",answer
283        assert num.allclose(f, answer)
284
285        # Delete file!
286        os.remove(txt_file)
287       
288    def cache_test_fit_to_mesh_pts(self):
289        a = [-1.0, 0.0]
290        b = [3.0, 4.0]
291        c = [4.0,1.0]
292        d = [-3.0, 2.0] #3
293        e = [-1.0,-2.0]
294        f = [1.0, -2.0] #5
295
296        vertices = [a, b, c, d,e,f]
297        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
298
299
300        fileName = tempfile.mktemp(".txt")
301        file = open(fileName,"w")
302        file.write(" x, y, elevation \n\
303-2.0, 2.0, 0.\n\
304-1.0, 1.0, 0.\n\
3050.0, 2.0 , 2.\n\
3061.0, 1.0 , 2.\n\
307 2.0,  1.0 ,3. \n\
308 0.0,  0.0 , 0.\n\
309 1.0,  0.0 , 1.\n\
310 0.0,  -1.0, -1.\n\
311 -0.2, -0.5, -0.7\n\
312 -0.9, -1.5, -2.4\n\
313 0.5,  -1.9, -1.4\n\
314 3.0,  1.0 , 4.\n")
315        file.close()
316        geo = Geospatial_data(fileName)
317        fileName_pts = tempfile.mktemp(".pts")
318        points = geo.get_data_points(absolute=True)
319        atts = geo.get_attributes()
320        f = fit_to_mesh(points,atts, vertices, triangles,
321                                alpha=0.0, max_read_lines=2,
322                        use_cache=True, verbose=True)
323        answer = linear_function(vertices)
324        #print "f\n",f
325        #print "answer\n",answer
326        assert num.allclose(f, answer)
327        os.remove(fileName)
328       
329    def test_fit_to_mesh_pts(self):
330        a = [-1.0, 0.0]
331        b = [3.0, 4.0]
332        c = [4.0,1.0]
333        d = [-3.0, 2.0] #3
334        e = [-1.0,-2.0]
335        f = [1.0, -2.0] #5
336
337        vertices = [a, b, c, d,e,f]
338        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
339
340
341        fileName = tempfile.mktemp(".txt")
342        file = open(fileName,"w")
343        file.write(" x, y, elevation \n\
344-2.0, 2.0, 0.\n\
345-1.0, 1.0, 0.\n\
3460.0, 2.0 , 2.\n\
3471.0, 1.0 , 2.\n\
348 2.0,  1.0 ,3. \n\
349 0.0,  0.0 , 0.\n\
350 1.0,  0.0 , 1.\n\
351 0.0,  -1.0, -1.\n\
352 -0.2, -0.5, -0.7\n\
353 -0.9, -1.5, -2.4\n\
354 0.5,  -1.9, -1.4\n\
355 3.0,  1.0 , 4.\n")
356        file.close()
357        geo = Geospatial_data(fileName)
358        fileName_pts = tempfile.mktemp(".pts")
359        geo.export_points_file(fileName_pts)
360        f = fit_to_mesh(fileName_pts, vertices, triangles,
361                                alpha=0.0, max_read_lines=2)
362        answer = linear_function(vertices)
363        #print "f\n",f
364        #print "answer\n",answer
365        assert num.allclose(f, answer)
366        os.remove(fileName)
367        os.remove(fileName_pts)
368       
369    def test_fit_to_mesh_pts_passing_mesh_in(self):
370        a = [-1.0, 0.0]
371        b = [3.0, 4.0]
372        c = [4.0,1.0]
373        d = [-3.0, 2.0] #3
374        e = [-1.0,-2.0]
375        f = [1.0, -2.0] #5
376
377        vertices = [a, b, c, d,e,f]
378        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
379
380
381        fileName = tempfile.mktemp(".txt")
382        file = open(fileName,"w")
383        file.write(" x, y, elevation \n\
384-2.0, 2.0, 0.\n\
385-1.0, 1.0, 0.\n\
3860.0, 2.0 , 2.\n\
3871.0, 1.0 , 2.\n\
388 2.0,  1.0 ,3. \n\
389 0.0,  0.0 , 0.\n\
390 1.0,  0.0 , 1.\n\
391 0.0,  -1.0, -1.\n\
392 -0.2, -0.5, -0.7\n\
393 -0.9, -1.5, -2.4\n\
394 0.5,  -1.9, -1.4\n\
395 3.0,  1.0 , 4.\n")
396        file.close()
397        geo = Geospatial_data(fileName)
398        fileName_pts = tempfile.mktemp(".pts")
399        geo.export_points_file(fileName_pts)
400        f = fit_to_mesh(fileName_pts, vertices, triangles,
401                                alpha=0.0, max_read_lines=2)
402        answer = linear_function(vertices)
403        #print "f\n",f
404        #print "answer\n",answer
405        assert num.allclose(f, answer)
406        os.remove(fileName)
407        os.remove(fileName_pts)
408       
409    def test_fit_to_mesh(self):
410
411        a = [-1.0, 0.0]
412        b = [3.0, 4.0]
413        c = [4.0,1.0]
414        d = [-3.0, 2.0] #3
415        e = [-1.0,-2.0]
416        f = [1.0, -2.0] #5
417
418        vertices = [a, b, c, d,e,f]
419        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
420
421
422        fileName = tempfile.mktemp(".ddd")
423        file = open(fileName,"w")
424        file.write(" x,y, elevation \n\
425-2.0, 2.0, 0.\n\
426-1.0, 1.0, 0.\n\
4270.0, 2.0 , 2.\n\
4281.0, 1.0 , 2.\n\
429 2.0,  1.0 ,3. \n\
430 0.0,  0.0 , 0.\n\
431 1.0,  0.0 , 1.\n\
432 0.0,  -1.0, -1.\n\
433 -0.2, -0.5, -0.7\n\
434 -0.9, -1.5, -2.4\n\
435 0.5,  -1.9, -1.4\n\
436 3.0,  1.0 , 4.\n")
437        file.close()
438       
439        f = fit_to_mesh(fileName, vertices, triangles,
440                                alpha=0.0, max_read_lines=2)
441                        #use_cache=True, verbose=True)
442        answer = linear_function(vertices)
443        #print "f\n",f
444        #print "answer\n",answer
445        assert num.allclose(f, answer)
446   
447        os.remove(fileName)
448       
449    def test_fit_to_mesh_2_atts(self):
450
451        a = [-1.0, 0.0]
452        b = [3.0, 4.0]
453        c = [4.0,1.0]
454        d = [-3.0, 2.0] #3
455        e = [-1.0,-2.0]
456        f = [1.0, -2.0] #5
457
458        vertices = [a, b, c, d,e,f]
459        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
460
461
462        fileName = tempfile.mktemp(".ddd")
463        file = open(fileName,"w")
464        # the 2nd att name is wacky so it's the first key off a hash table
465        file.write(" x,y, elevation, afriqction \n\
466-2.0, 2.0, 0., 0.\n\
467-1.0, 1.0, 0., 0.\n\
4680.0, 2.0 , 2., 20.\n\
4691.0, 1.0 , 2., 20.\n\
470 2.0,  1.0 ,3., 30. \n\
471 0.0,  0.0 , 0., 0.\n\
472 1.0,  0.0 , 1., 10.\n\
473 0.0,  -1.0, -1., -10.\n\
474 -0.2, -0.5, -0.7, -7.\n\
475 -0.9, -1.5, -2.4, -24. \n\
476 0.5,  -1.9, -1.4, -14. \n\
477 3.0,  1.0 , 4., 40. \n")
478        file.close()
479       
480        f = fit_to_mesh(fileName, vertices, triangles,
481                        alpha=0.0, 
482                        attribute_name='elevation', max_read_lines=2)
483        answer = linear_function(vertices)
484        #print "f\n",f
485        #print "answer\n",answer
486        assert num.allclose(f, answer)
487        os.remove(fileName)
488       
489    def test_fit_and_interpolation(self):
490
491        a = [0.0, 0.0]
492        b = [0.0, 2.0]
493        c = [2.0, 0.0]
494        d = [0.0, 4.0]
495        e = [2.0, 2.0]
496        f = [4.0, 0.0]
497
498        points = [a, b, c, d, e, f]
499        #bac, bce, ecf, dbe, daf, dae
500        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
501
502        #Get (enough) datapoints
503        data_points = [[ 0.66666667, 0.66666667],
504                       [ 1.33333333, 1.33333333],
505                       [ 2.66666667, 0.66666667],
506                       [ 0.66666667, 2.66666667],
507                       [ 0.0, 1.0],
508                       [ 0.0, 3.0],
509                       [ 1.0, 0.0],
510                       [ 1.0, 1.0],
511                       [ 1.0, 2.0],
512                       [ 1.0, 3.0],
513                       [ 2.0, 1.0],
514                       [ 3.0, 0.0],
515                       [ 3.0, 1.0]]
516
517        z = linear_function(data_points)
518        interp = Fit(points, triangles,
519                                alpha=0.0)
520
521        answer = linear_function(points)
522
523        f = interp.fit(data_points, z)
524
525        #print "f",f
526        #print "answer",answer
527        assert num.allclose(f, answer)
528
529       
530    def test_smoothing_and_interpolation(self):
531
532        a = [0.0, 0.0]
533        b = [0.0, 2.0]
534        c = [2.0, 0.0]
535        d = [0.0, 4.0]
536        e = [2.0, 2.0]
537        f = [4.0, 0.0]
538
539        points = [a, b, c, d, e, f]
540        #bac, bce, ecf, dbe, daf, dae
541        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
542
543        #Get (too few!) datapoints
544        data_points = [[ 0.66666667, 0.66666667],
545                       [ 1.33333333, 1.33333333],
546                       [ 2.66666667, 0.66666667],
547                       [ 0.66666667, 2.66666667]]
548
549        z = linear_function(data_points)
550        answer = linear_function(points)
551
552        #Make interpolator with too few data points and no smoothing
553        interp = Fit(points, triangles, alpha=0.0)
554        #Must raise an exception
555        try:
556            f = interp.fit(data_points,z)
557        except TooFewPointsError:
558            pass
559
560        #Now try with smoothing parameter
561        interp = Fit(points, triangles, alpha=1.0e-13)
562
563        f = interp.fit(data_points,z)
564        #f will be different from answer due to smoothing
565        assert num.allclose(f, answer,atol=5)
566
567
568    #Tests of smoothing matrix
569    def test_smoothing_matrix_one_triangle(self):
570        a = [0.0, 0.0]
571        b = [0.0, 2.0]
572        c = [2.0,0.0]
573        points = [a, b, c]
574
575        vertices = [ [1,0,2] ]   #bac
576
577        interp = Fit(points, vertices)
578
579        assert num.allclose(interp.get_D(), [[1, -0.5, -0.5],
580                                             [-0.5, 0.5, 0],
581                                             [-0.5, 0, 0.5]])
582
583        #Define f(x,y) = x
584        f = num.array([0,0,2]) #Value at global vertex 2
585
586        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
587        #           int 1 dx dy = area = 2
588        assert num.dot(num.dot(f, interp.get_D()), f) == 2
589
590        #Define f(x,y) = y
591        f = num.array([0,2,0])  #Value at global vertex 1
592
593        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
594        #           int 1 dx dy = area = 2
595        assert num.dot(num.dot(f, interp.get_D()), f) == 2
596
597        #Define f(x,y) = x+y
598        f = num.array([0,2,2])  #Values at global vertex 1 and 2
599
600        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
601        #           int 2 dx dy = 2*area = 4
602        assert num.dot(num.dot(f, interp.get_D()), f) == 4
603
604
605    def test_smoothing_matrix_more_triangles(self):
606        a = [0.0, 0.0]
607        b = [0.0, 2.0]
608        c = [2.0,0.0]
609        d = [0.0, 4.0]
610        e = [2.0, 2.0]
611        f = [4.0,0.0]
612
613        points = [a, b, c, d, e, f]
614        #bac, bce, ecf, dbe, daf, dae
615        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
616
617        interp = Fit(points, vertices)
618
619
620        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
621        #                           [-0.5, 0.5, 0],
622        #                           [-0.5, 0, 0.5]])
623
624        #Define f(x,y) = x
625        f = num.array([0,0,2,0,2,4]) #f evaluated at points a-f
626
627        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
628        #           int 1 dx dy = total area = 8
629        assert num.dot(num.dot(f, interp.get_D()), f) == 8
630
631        #Define f(x,y) = y
632        f = num.array([0,2,0,4,2,0]) #f evaluated at points a-f
633
634        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
635        #           int 1 dx dy = area = 8
636        assert num.dot(num.dot(f, interp.get_D()), f) == 8
637
638        #Define f(x,y) = x+y
639        f = num.array([0,2,2,4,4,4])  #f evaluated at points a-f
640
641        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
642        #           int 2 dx dy = 2*area = 16
643        assert num.dot(num.dot(f, interp.get_D()), f) == 16
644
645
646
647    def test_fit_and_interpolation_with_different_origins(self):
648        """Fit a surface to one set of points. Then interpolate that surface
649        using another set of points.
650        This test tests situtaion where points and mesh belong to a different
651        coordinate system as defined by origin.
652        """
653
654        #Setup mesh used to represent fitted function
655        a = [0.0, 0.0]
656        b = [0.0, 2.0]
657        c = [2.0, 0.0]
658        d = [0.0, 4.0]
659        e = [2.0, 2.0]
660        f = [4.0, 0.0]
661
662        points = [a, b, c, d, e, f]
663        #bac, bce, ecf, dbe, daf, dae
664        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
665
666        #Datapoints to fit from
667        data_points1 = [[ 0.66666667, 0.66666667],
668                        [ 1.33333333, 1.33333333],
669                        [ 2.66666667, 0.66666667],
670                        [ 0.66666667, 2.66666667],
671                        [ 0.0, 1.0],
672                        [ 0.0, 3.0],
673                        [ 1.0, 0.0],
674                        [ 1.0, 1.0],
675                        [ 1.0, 2.0],
676                        [ 1.0, 3.0],
677                        [ 2.0, 1.0],
678                        [ 3.0, 0.0],
679                        [ 3.0, 1.0]]
680
681
682        #First check that things are OK when using same origin
683        mesh_origin = (56, 290000, 618000) #zone, easting, northing
684        data_origin = (56, 290000, 618000) #zone, easting, northing
685
686
687        #Fit surface to mesh
688        interp = Fit(points, triangles, 
689                               alpha=0.0,
690                               mesh_origin = mesh_origin)
691
692        data_geo_spatial = Geospatial_data(data_points1,
693                         geo_reference = Geo_reference(56, 290000, 618000))
694        z = linear_function(data_points1) #Example z-values
695        f = interp.fit(data_geo_spatial, z)   #Fitted values at vertices
696
697        #Shift datapoints according to new origins
698        for k in range(len(data_points1)):
699            data_points1[k][0] += mesh_origin[1] - data_origin[1]
700            data_points1[k][1] += mesh_origin[2] - data_origin[2]
701
702
703
704        #Fit surface to mesh
705        interp = Fit(points, triangles, alpha=0.0)
706        #Fitted values at vertices (using same z as before)
707        f1 = interp.fit(data_points1,z) 
708
709        assert num.allclose(f,f1), 'Fit should have been unaltered'
710
711
712    def test_smooth_attributes_to_mesh_function(self):
713        #Testing 2 attributes smoothed to the mesh
714       
715
716        a = [0.0, 0.0]
717        b = [0.0, 5.0]
718        c = [5.0, 0.0]
719        points = [a, b, c]
720        triangles = [ [1,0,2] ]   #bac
721
722        d1 = [1.0, 1.0]
723        d2 = [1.0, 3.0]
724        d3 = [3.0, 1.0]
725        z1 = [2, 4]
726        z2 = [4, 8]
727        z3 = [4, 8]
728        data_coords = [d1, d2, d3]
729        z = [z1, z2, z3]
730
731        f = fit_to_mesh(data_coords, vertex_coordinates=points,
732                        triangles=triangles,
733                        point_attributes=z, alpha=0.0)
734        answer = [[0, 0], [5., 10.], [5., 10.]]
735
736        assert num.allclose(f, answer)
737
738    def test_fit_to_mesh_w_georef(self):
739        """Simple check that georef works at the fit_to_mesh level
740        """
741       
742        from anuga.coordinate_transforms.geo_reference import Geo_reference
743
744        #Mesh
745        vertex_coordinates = [[0.76, 0.76],
746                              [0.76, 5.76],
747                              [5.76, 0.76]]
748        triangles = [[0,2,1]]
749
750        mesh_geo = Geo_reference(56,-0.76,-0.76)
751        #print "mesh_geo.get_absolute(vertex_coordinates)", \
752         #     mesh_geo.get_absolute(vertex_coordinates)
753
754        #Data                       
755        data_points = [[ 201.0, 401.0],
756                       [ 201.0, 403.0],
757                       [ 203.0, 401.0]]
758
759        z = [2, 4, 4]
760       
761        data_geo = Geo_reference(56,-200,-400)
762
763        #print "data_geo.get_absolute(data_points)", \
764        #      data_geo.get_absolute(data_points)
765       
766        #Fit
767        zz = fit_to_mesh(data_points, vertex_coordinates=vertex_coordinates,
768                         triangles=triangles,
769                         point_attributes=z,
770                         data_origin = data_geo.get_origin(),
771                         mesh_origin = mesh_geo.get_origin(),
772                         alpha = 0)
773        assert num.allclose( zz, [0,5,5] )
774
775
776    def test_fit_to_mesh_file2domain(self):
777        from load_mesh.loadASCII import import_mesh_file, \
778             export_mesh_file
779        import tempfile
780        import os
781
782        # create a .tsh file, no user outline
783        mesh_dic = {}
784        mesh_dic['vertices'] = [[0.0, 0.0],
785                                          [0.0, 5.0],
786                                          [5.0, 0.0]]
787        mesh_dic['triangles'] =  [[0, 2, 1]]
788        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
789        mesh_dic['triangle_tags'] = ['']
790        mesh_dic['vertex_attributes'] = [[], [], []]
791        mesh_dic['vertiex_attribute_titles'] = []
792        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
793        mesh_dic['segment_tags'] = ['external',
794                                                  'external',
795                                                  'external']
796        mesh_file = tempfile.mktemp(".tsh")
797        export_mesh_file(mesh_file,mesh_dic)
798
799        # create a points .csv file
800        point_file = tempfile.mktemp(".csv")
801        fd = open(point_file,'w')
802        fd.write("x,y, elevation, stage \n\
803        1.0, 1.0,2.,4 \n\
804        1.0, 3.0,4,8 \n\
805        3.0,1.0,4.,8 \n")
806        fd.close()
807
808        mesh_output_file = tempfile.mktemp(".tsh") 
809        fit_to_mesh_file(mesh_file,
810                         point_file,
811                         mesh_output_file,
812                         alpha = 0.0)
813        # load in the .tsh file we just wrote
814        mesh_dic = import_mesh_file(mesh_output_file)
815        #print "mesh_dic",mesh_dic
816        ans =[[0.0, 0.0],
817              [5.0, 10.0],
818              [5.0,10.0]]
819        assert num.allclose(mesh_dic['vertex_attributes'],ans)
820
821        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
822                        ['elevation','stage'],
823                        'test_fit_to_mesh_file failed')
824        domain = Domain(mesh_output_file, use_cache=True, verbose=False)
825       
826        answer = [0., 5., 5.]
827        assert num.allclose(domain.quantities['elevation'].vertex_values,
828                            answer)
829        #clean up
830        os.remove(mesh_file)
831        os.remove(point_file)
832        os.remove(mesh_output_file)
833
834    def test_fit_to_mesh_file3(self):
835        from load_mesh.loadASCII import import_mesh_file, \
836             export_mesh_file
837        import tempfile
838        import os
839
840        # create a .tsh file, no user outline
841        mesh_dic = {}
842        mesh_dic['vertices'] = [[0.76, 0.76],
843                                          [0.76, 5.76],
844                                          [5.76, 0.76]]
845        mesh_dic['triangles'] =  [[0, 2, 1]]
846        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
847        mesh_dic['triangle_tags'] = ['']
848        mesh_dic['vertex_attributes'] = [[], [], []]
849        mesh_dic['vertiex_attribute_titles'] = []
850        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
851        mesh_dic['segment_tags'] = ['external',
852                                                  'external',
853                                                  'external']
854        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
855        mesh_file = tempfile.mktemp(".tsh")
856        export_mesh_file(mesh_file,mesh_dic)
857
858        # create a points .csv file
859        point_file = tempfile.mktemp(".csv")
860        fd = open(point_file,'w')
861        fd.write("x,y, elevation, stage \n\
862        1.0, 1.0,2.,4 \n\
863        1.0, 3.0,4,8 \n\
864        3.0,1.0,4.,8 \n")
865        fd.close()
866       
867        mesh_output_file = tempfile.mktemp(".tsh")
868        fit_to_mesh_file(mesh_file,
869                         point_file,
870                         mesh_output_file,
871                         alpha = 0.0)
872        # load in the .tsh file we just wrote
873        mesh_dic = import_mesh_file(mesh_output_file)
874        #print "mesh_dic",mesh_dic
875        ans =[[0.0, 0.0],
876              [5.0, 10.0],
877              [5.0,10.0]]
878        assert num.allclose(mesh_dic['vertex_attributes'],ans)
879
880        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
881                        ['elevation','stage'],
882                        'test_fit_to_mesh_file failed')
883
884        #clean up
885        os.remove(mesh_file)
886        os.remove(point_file)
887        os.remove(mesh_output_file)
888
889    def test_fit_to_mesh_fileII(self):
890        from load_mesh.loadASCII import import_mesh_file, \
891             export_mesh_file
892        import tempfile
893        import os
894
895        # create a .tsh file, no user outline
896        mesh_dic = {}
897        mesh_dic['vertices'] = [[0.0, 0.0],
898                                [0.0, 5.0],
899                                [5.0, 0.0]]
900        mesh_dic['triangles'] =  [[0, 2, 1]]
901        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
902        mesh_dic['triangle_tags'] = ['']
903        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
904        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
905        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
906        mesh_dic['segment_tags'] = ['external',
907                                                  'external',
908                                                  'external']
909        mesh_file = tempfile.mktemp(".tsh")
910        export_mesh_file(mesh_file,mesh_dic)
911
912        # create a points .csv file
913        point_file = tempfile.mktemp(".csv")
914        fd = open(point_file,'w')
915        fd.write("x,y,elevation, stage \n\
916        1.0, 1.0,2.,4 \n\
917        1.0, 3.0,4,8 \n\
918        3.0,1.0,4.,8 \n")
919        fd.close()
920
921        mesh_output_file = "new_triangle.tsh"
922        fit_to_mesh_file(mesh_file,
923                         point_file,
924                         mesh_output_file,
925                         alpha = 0.0)
926        # load in the .tsh file we just wrote
927        mesh_dic = import_mesh_file(mesh_output_file)
928
929        assert num.allclose(mesh_dic['vertex_attributes'],
930                            [[1.0, 2.0,0.0, 0.0],
931                             [1.0, 2.0,5.0, 10.0],
932                             [1.0, 2.0,5.0,10.0]])
933
934        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
935                        ['density', 'temp','elevation','stage'],
936                        'test_fit_to_mesh_file failed')
937
938        #clean up
939        os.remove(mesh_file)
940        os.remove(mesh_output_file)
941        os.remove(point_file)
942
943    def test_fit_to_mesh_file_errors(self):
944        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
945        import tempfile
946        import os
947
948        # create a .tsh file, no user outline
949        mesh_dic = {}
950        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
951        mesh_dic['triangles'] =  [[0, 2, 1]]
952        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
953        mesh_dic['triangle_tags'] = ['']
954        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
955        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
956        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
957        mesh_dic['segment_tags'] = ['external', 'external','external']
958        mesh_file = tempfile.mktemp(".tsh")
959        export_mesh_file(mesh_file,mesh_dic)
960
961        # create a bad points .csv file
962        point_file = tempfile.mktemp(".csv")
963        fd = open(point_file,'w')
964        fd.write("x,y,elevation stage \n\
965        1.0, 1.0,2.,4 \n\
966        1.0, 3.0,4,8 \n\
967        3.0,1.0,4.,8 \n")
968        fd.close()
969
970        mesh_output_file = "new_triangle.tsh"
971        try:
972            fit_to_mesh_file(mesh_file, point_file,
973                             mesh_output_file, display_errors = False)
974        except SyntaxError:
975            pass
976        else:
977            #self.failUnless(0 ==1,  'Bad file did not raise error!')
978            raise 'Bad file did not raise error!'
979           
980        #clean up
981        os.remove(mesh_file)
982        os.remove(point_file)
983
984    def test_fit_to_mesh_file_errorsII(self):
985        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
986        import tempfile
987        import os
988
989        # create a .tsh file, no user outline
990        mesh_file = tempfile.mktemp(".tsh")
991        fd = open(mesh_file,'w')
992        fd.write("unit testing a bad .tsh file \n")
993        fd.close()
994
995        # create a points .csv file
996        point_file = tempfile.mktemp(".csv")
997        fd = open(point_file,'w')
998        fd.write("x,y,elevation, stage \n\
999        1.0, 1.0,2.,4 \n\
1000        1.0, 3.0,4,8 \n\
1001        3.0,1.0,4.,8 \n")
1002        fd.close()
1003
1004        mesh_output_file = "new_triangle.tsh"
1005        try:
1006            fit_to_mesh_file(mesh_file, point_file,
1007                             mesh_output_file, display_errors = False)
1008        except IOError:
1009            pass
1010        else:
1011            raise 'Bad file did not raise error!'
1012           
1013        #clean up
1014        os.remove(mesh_file)
1015        os.remove(point_file)
1016
1017    def test_fit_to_mesh_file_errorsIII(self):
1018        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
1019        import tempfile
1020        import os
1021
1022        # create a .tsh file, no user outline
1023        mesh_dic = {}
1024        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
1025        mesh_dic['triangles'] =  [[0, 2, 1]]
1026        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
1027        mesh_dic['triangle_tags'] = ['']
1028        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
1029        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
1030        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
1031        mesh_dic['segment_tags'] = ['external', 'external','external']
1032        mesh_file = tempfile.mktemp(".tsh")
1033        export_mesh_file(mesh_file,mesh_dic)
1034
1035
1036        # create a points .csv file
1037        point_file = tempfile.mktemp(".csv")
1038        fd = open(point_file,'w')
1039        fd.write("x,y,elevation, stage \n\
1040        1.0, 1.0,2.,4 \n\
1041        1.0, 3.0,4,8 \n\
1042        3.0,1.0,4.,8 \n")
1043        fd.close()
1044
1045        #This a deliberately illegal filename to invoke the error.
1046        mesh_output_file = ".../\z\z:ya.tsh"       
1047
1048        try:
1049            fit_to_mesh_file(mesh_file, point_file,
1050                             mesh_output_file, display_errors = False)
1051        except IOError:
1052            pass
1053        else:
1054            raise 'Bad file did not raise error!'
1055       
1056        #clean up
1057        os.remove(mesh_file)
1058        os.remove(point_file)
1059 
1060    def Not_yet_test_smooth_att_to_mesh_with_excess_verts(self):
1061
1062        a = [0.0, 0.0]
1063        b = [0.0, 5.0]
1064        c = [5.0, 0.0]
1065        d = [1.0, 1.0]
1066        e = [18.0, 1000.0]
1067       
1068        points = [a, b, c, d, e]
1069        triangles = [ [1,0,2] ]   #bac
1070
1071        d1 = [1.0, 1.0]
1072        d2 = [1.0, 2.0]
1073        d3 = [3.0,1.0]
1074        d4 = [2.0,3.0]
1075        d5 = [2.0,2.0]
1076        d6 = [1.0,3.0]
1077        data_coords = [d1, d2, d3, d4, d5, d6]
1078        z = linear_function(data_coords)
1079        #print "z",z
1080
1081        interp = Fit(points, triangles, alpha=0.0)
1082       
1083        try:
1084            f = interp.fit(data_coords, z)
1085        except VertsWithNoTrianglesError:
1086            pass
1087        else:
1088            raise 'Verts with no triangles did not raise error!'
1089       
1090        #f = interp.fit(data_coords, z)
1091        #answer = linear_function(points)
1092
1093        #  Removing the bad verts that we don't care about
1094        #f = f[0:3]
1095        #answer = answer[0:3]
1096        #assert allclose(f, answer)
1097
1098#-------------------------------------------------------------
1099if __name__ == "__main__":
1100    suite = unittest.makeSuite(Test_Fit,'test')
1101    runner = unittest.TextTestRunner() #verbosity=1)
1102    runner.run(suite)
1103
Note: See TracBrowser for help on using the repository browser.