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

Last change on this file since 6536 was 6174, checked in by rwilson, 15 years ago

Changed .array(A) to .array(A, num.Int) where appropriate. Helps when going to numpy.

File size: 33.5 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 Numeric 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], num.Int) #Value at global vertex 2      #array default#
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], num.Int)  #Value at global vertex 1      #array default#
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], num.Int)  #Values at global vertex 1 and 2      #array default#
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], num.Int) #f evaluated at points a-f      #array default#
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], num.Int) #f evaluated at points a-f      #array default#
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], num.Int)  #f evaluated at points a-f     #array default#
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    #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_function')
1102    #suite = unittest.makeSuite(Test_Fit,'')
1103    runner = unittest.TextTestRunner(verbosity=1)
1104    runner.run(suite)
1105
1106
1107
1108
1109
Note: See TracBrowser for help on using the repository browser.