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

Last change on this file since 5181 was 4859, checked in by duncan, 16 years ago

check the last triangle fit/interp speed up. Upgrading the code to time runs

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