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

Last change on this file since 4659 was 4659, checked in by duncan, 17 years ago

remove the .xya file format.

File size: 31.9 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(points, elements, txt_file,
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(vertices, triangles,points,atts,
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(vertices, triangles,fileName_pts,
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(self):
364
365        a = [-1.0, 0.0]
366        b = [3.0, 4.0]
367        c = [4.0,1.0]
368        d = [-3.0, 2.0] #3
369        e = [-1.0,-2.0]
370        f = [1.0, -2.0] #5
371
372        vertices = [a, b, c, d,e,f]
373        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
374
375
376        fileName = tempfile.mktemp(".ddd")
377        file = open(fileName,"w")
378        file.write(" x,y, elevation \n\
379-2.0, 2.0, 0.\n\
380-1.0, 1.0, 0.\n\
3810.0, 2.0 , 2.\n\
3821.0, 1.0 , 2.\n\
383 2.0,  1.0 ,3. \n\
384 0.0,  0.0 , 0.\n\
385 1.0,  0.0 , 1.\n\
386 0.0,  -1.0, -1.\n\
387 -0.2, -0.5, -0.7\n\
388 -0.9, -1.5, -2.4\n\
389 0.5,  -1.9, -1.4\n\
390 3.0,  1.0 , 4.\n")
391        file.close()
392       
393        f = fit_to_mesh(vertices, triangles,fileName,
394                                alpha=0.0, max_read_lines=2)
395                        #use_cache=True, verbose=True)
396        answer = linear_function(vertices)
397        #print "f\n",f
398        #print "answer\n",answer
399        assert allclose(f, answer)
400   
401        os.remove(fileName)
402       
403    def test_fit_to_mesh_2_atts(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        # the 2nd att name is wacky so it's the first key off a hash table
419        file.write(" x,y, elevation, afriqction \n\
420-2.0, 2.0, 0., 0.\n\
421-1.0, 1.0, 0., 0.\n\
4220.0, 2.0 , 2., 20.\n\
4231.0, 1.0 , 2., 20.\n\
424 2.0,  1.0 ,3., 30. \n\
425 0.0,  0.0 , 0., 0.\n\
426 1.0,  0.0 , 1., 10.\n\
427 0.0,  -1.0, -1., -10.\n\
428 -0.2, -0.5, -0.7, -7.\n\
429 -0.9, -1.5, -2.4, -24. \n\
430 0.5,  -1.9, -1.4, -14. \n\
431 3.0,  1.0 , 4., 40. \n")
432        file.close()
433       
434        f = fit_to_mesh(vertices, triangles,fileName,
435                        alpha=0.0, 
436                        attribute_name='elevation', max_read_lines=2)
437        answer = linear_function(vertices)
438        #print "f\n",f
439        #print "answer\n",answer
440        assert allclose(f, answer)
441        os.remove(fileName)
442       
443    def test_fit_and_interpolation(self):
444
445        a = [0.0, 0.0]
446        b = [0.0, 2.0]
447        c = [2.0, 0.0]
448        d = [0.0, 4.0]
449        e = [2.0, 2.0]
450        f = [4.0, 0.0]
451
452        points = [a, b, c, d, e, f]
453        #bac, bce, ecf, dbe, daf, dae
454        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
455
456        #Get (enough) datapoints
457        data_points = [[ 0.66666667, 0.66666667],
458                       [ 1.33333333, 1.33333333],
459                       [ 2.66666667, 0.66666667],
460                       [ 0.66666667, 2.66666667],
461                       [ 0.0, 1.0],
462                       [ 0.0, 3.0],
463                       [ 1.0, 0.0],
464                       [ 1.0, 1.0],
465                       [ 1.0, 2.0],
466                       [ 1.0, 3.0],
467                       [ 2.0, 1.0],
468                       [ 3.0, 0.0],
469                       [ 3.0, 1.0]]
470
471        z = linear_function(data_points)
472        interp = Fit(points, triangles,
473                                alpha=0.0)
474
475        answer = linear_function(points)
476
477        f = interp.fit(data_points, z)
478
479        #print "f",f
480        #print "answer",answer
481        assert allclose(f, answer)
482
483       
484    def test_smoothing_and_interpolation(self):
485
486        a = [0.0, 0.0]
487        b = [0.0, 2.0]
488        c = [2.0, 0.0]
489        d = [0.0, 4.0]
490        e = [2.0, 2.0]
491        f = [4.0, 0.0]
492
493        points = [a, b, c, d, e, f]
494        #bac, bce, ecf, dbe, daf, dae
495        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
496
497        #Get (too few!) datapoints
498        data_points = [[ 0.66666667, 0.66666667],
499                       [ 1.33333333, 1.33333333],
500                       [ 2.66666667, 0.66666667],
501                       [ 0.66666667, 2.66666667]]
502
503        z = linear_function(data_points)
504        answer = linear_function(points)
505
506        #Make interpolator with too few data points and no smoothing
507        interp = Fit(points, triangles, alpha=0.0)
508        #Must raise an exception
509        try:
510            f = interp.fit(data_points,z)
511        except ToFewPointsError:
512            pass
513
514        #Now try with smoothing parameter
515        interp = Fit(points, triangles, alpha=1.0e-13)
516
517        f = interp.fit(data_points,z)
518        #f will be different from answer due to smoothing
519        assert allclose(f, answer,atol=5)
520
521
522    #Tests of smoothing matrix
523    def test_smoothing_matrix_one_triangle(self):
524        from Numeric import dot
525        a = [0.0, 0.0]
526        b = [0.0, 2.0]
527        c = [2.0,0.0]
528        points = [a, b, c]
529
530        vertices = [ [1,0,2] ]   #bac
531
532        interp = Fit(points, vertices)
533
534        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
535                                   [-0.5, 0.5, 0],
536                                   [-0.5, 0, 0.5]])
537
538        #Define f(x,y) = x
539        f = array([0,0,2]) #Value at global vertex 2
540
541        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
542        #           int 1 dx dy = area = 2
543        assert dot(dot(f, interp.get_D()), f) == 2
544
545        #Define f(x,y) = y
546        f = array([0,2,0])  #Value at global vertex 1
547
548        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
549        #           int 1 dx dy = area = 2
550        assert dot(dot(f, interp.get_D()), f) == 2
551
552        #Define f(x,y) = x+y
553        f = array([0,2,2])  #Values at global vertex 1 and 2
554
555        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
556        #           int 2 dx dy = 2*area = 4
557        assert dot(dot(f, interp.get_D()), f) == 4
558
559
560    def test_smoothing_matrix_more_triangles(self):
561        from Numeric import dot
562
563        a = [0.0, 0.0]
564        b = [0.0, 2.0]
565        c = [2.0,0.0]
566        d = [0.0, 4.0]
567        e = [2.0, 2.0]
568        f = [4.0,0.0]
569
570        points = [a, b, c, d, e, f]
571        #bac, bce, ecf, dbe, daf, dae
572        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
573
574        interp = Fit(points, vertices)
575
576
577        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
578        #                           [-0.5, 0.5, 0],
579        #                           [-0.5, 0, 0.5]])
580
581        #Define f(x,y) = x
582        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
583
584        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
585        #           int 1 dx dy = total area = 8
586        assert dot(dot(f, interp.get_D()), f) == 8
587
588        #Define f(x,y) = y
589        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
590
591        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
592        #           int 1 dx dy = area = 8
593        assert dot(dot(f, interp.get_D()), f) == 8
594
595        #Define f(x,y) = x+y
596        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
597
598        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
599        #           int 2 dx dy = 2*area = 16
600        assert dot(dot(f, interp.get_D()), f) == 16
601
602
603
604    def test_fit_and_interpolation_with_different_origins(self):
605        """Fit a surface to one set of points. Then interpolate that surface
606        using another set of points.
607        This test tests situtaion where points and mesh belong to a different
608        coordinate system as defined by origin.
609        """
610
611        #Setup mesh used to represent fitted function
612        a = [0.0, 0.0]
613        b = [0.0, 2.0]
614        c = [2.0, 0.0]
615        d = [0.0, 4.0]
616        e = [2.0, 2.0]
617        f = [4.0, 0.0]
618
619        points = [a, b, c, d, e, f]
620        #bac, bce, ecf, dbe, daf, dae
621        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
622
623        #Datapoints to fit from
624        data_points1 = [[ 0.66666667, 0.66666667],
625                        [ 1.33333333, 1.33333333],
626                        [ 2.66666667, 0.66666667],
627                        [ 0.66666667, 2.66666667],
628                        [ 0.0, 1.0],
629                        [ 0.0, 3.0],
630                        [ 1.0, 0.0],
631                        [ 1.0, 1.0],
632                        [ 1.0, 2.0],
633                        [ 1.0, 3.0],
634                        [ 2.0, 1.0],
635                        [ 3.0, 0.0],
636                        [ 3.0, 1.0]]
637
638
639        #First check that things are OK when using same origin
640        mesh_origin = (56, 290000, 618000) #zone, easting, northing
641        data_origin = (56, 290000, 618000) #zone, easting, northing
642
643
644        #Fit surface to mesh
645        interp = Fit(points, triangles, 
646                               alpha=0.0,
647                               mesh_origin = mesh_origin)
648
649        data_geo_spatial = Geospatial_data(data_points1,
650                         geo_reference = Geo_reference(56, 290000, 618000))
651        z = linear_function(data_points1) #Example z-values
652        f = interp.fit(data_geo_spatial, z)   #Fitted values at vertices
653
654        #Shift datapoints according to new origins
655        for k in range(len(data_points1)):
656            data_points1[k][0] += mesh_origin[1] - data_origin[1]
657            data_points1[k][1] += mesh_origin[2] - data_origin[2]
658
659
660
661        #Fit surface to mesh
662        interp = Fit(points, triangles, alpha=0.0)
663        #Fitted values at vertices (using same z as before)
664        f1 = interp.fit(data_points1,z) 
665
666        assert allclose(f,f1), 'Fit should have been unaltered'
667
668
669    def test_smooth_attributes_to_mesh_function(self):
670        """ Testing 2 attributes smoothed to the mesh
671        """
672
673        a = [0.0, 0.0]
674        b = [0.0, 5.0]
675        c = [5.0, 0.0]
676        points = [a, b, c]
677        triangles = [ [1,0,2] ]   #bac
678
679        d1 = [1.0, 1.0]
680        d2 = [1.0, 3.0]
681        d3 = [3.0, 1.0]
682        z1 = [2, 4]
683        z2 = [4, 8]
684        z3 = [4, 8]
685        data_coords = [d1, d2, d3]
686        z = [z1, z2, z3]
687
688        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
689        answer = [[0, 0], [5., 10.], [5., 10.]]
690
691        assert allclose(f, answer)
692
693    def test_fit_to_mesh_w_georef(self):
694        """Simple check that georef works at the fit_to_mesh level
695        """
696       
697        from anuga.coordinate_transforms.geo_reference import Geo_reference
698
699        #Mesh
700        vertex_coordinates = [[0.76, 0.76],
701                              [0.76, 5.76],
702                              [5.76, 0.76]]
703        triangles = [[0,2,1]]
704
705        mesh_geo = Geo_reference(56,-0.76,-0.76)
706        #print "mesh_geo.get_absolute(vertex_coordinates)", \
707         #     mesh_geo.get_absolute(vertex_coordinates)
708
709        #Data                       
710        data_points = [[ 201.0, 401.0],
711                       [ 201.0, 403.0],
712                       [ 203.0, 401.0]]
713
714        z = [2, 4, 4]
715       
716        data_geo = Geo_reference(56,-200,-400)
717
718        #print "data_geo.get_absolute(data_points)", \
719        #      data_geo.get_absolute(data_points)
720       
721        #Fit
722        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
723                         data_origin = data_geo.get_origin(),
724                         mesh_origin = mesh_geo.get_origin(),
725                         alpha = 0)
726        assert allclose( zz, [0,5,5] )
727
728
729    def test_fit_to_mesh_file2domain(self):
730        from load_mesh.loadASCII import import_mesh_file, \
731             export_mesh_file
732        import tempfile
733        import os
734
735        # create a .tsh file, no user outline
736        mesh_dic = {}
737        mesh_dic['vertices'] = [[0.0, 0.0],
738                                          [0.0, 5.0],
739                                          [5.0, 0.0]]
740        mesh_dic['triangles'] =  [[0, 2, 1]]
741        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
742        mesh_dic['triangle_tags'] = ['']
743        mesh_dic['vertex_attributes'] = [[], [], []]
744        mesh_dic['vertiex_attribute_titles'] = []
745        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
746        mesh_dic['segment_tags'] = ['external',
747                                                  'external',
748                                                  'external']
749        mesh_file = tempfile.mktemp(".tsh")
750        export_mesh_file(mesh_file,mesh_dic)
751
752        # create a points .csv file
753        point_file = tempfile.mktemp(".csv")
754        fd = open(point_file,'w')
755        fd.write("x,y, elevation, stage \n\
756        1.0, 1.0,2.,4 \n\
757        1.0, 3.0,4,8 \n\
758        3.0,1.0,4.,8 \n")
759        fd.close()
760
761        mesh_output_file = tempfile.mktemp(".tsh") 
762        fit_to_mesh_file(mesh_file,
763                         point_file,
764                         mesh_output_file,
765                         alpha = 0.0)
766        # load in the .tsh file we just wrote
767        mesh_dic = import_mesh_file(mesh_output_file)
768        #print "mesh_dic",mesh_dic
769        ans =[[0.0, 0.0],
770              [5.0, 10.0],
771              [5.0,10.0]]
772        assert allclose(mesh_dic['vertex_attributes'],ans)
773
774        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
775                        ['elevation','stage'],
776                        'test_fit_to_mesh_file failed')
777        domain = Domain(mesh_output_file, use_cache=True, verbose=False)
778       
779        answer = [0., 5., 5.]
780        assert allclose(domain.quantities['elevation'].vertex_values,
781                        answer)
782        #clean up
783        os.remove(mesh_file)
784        os.remove(point_file)
785        os.remove(mesh_output_file)
786
787    def test_fit_to_mesh_file3(self):
788        from load_mesh.loadASCII import import_mesh_file, \
789             export_mesh_file
790        import tempfile
791        import os
792
793        # create a .tsh file, no user outline
794        mesh_dic = {}
795        mesh_dic['vertices'] = [[0.76, 0.76],
796                                          [0.76, 5.76],
797                                          [5.76, 0.76]]
798        mesh_dic['triangles'] =  [[0, 2, 1]]
799        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
800        mesh_dic['triangle_tags'] = ['']
801        mesh_dic['vertex_attributes'] = [[], [], []]
802        mesh_dic['vertiex_attribute_titles'] = []
803        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
804        mesh_dic['segment_tags'] = ['external',
805                                                  'external',
806                                                  'external']
807        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
808        mesh_file = tempfile.mktemp(".tsh")
809        export_mesh_file(mesh_file,mesh_dic)
810
811        # create a points .csv file
812        point_file = tempfile.mktemp(".csv")
813        fd = open(point_file,'w')
814        fd.write("x,y, elevation, stage \n\
815        1.0, 1.0,2.,4 \n\
816        1.0, 3.0,4,8 \n\
817        3.0,1.0,4.,8 \n")
818        fd.close()
819       
820        mesh_output_file = tempfile.mktemp(".tsh")
821        fit_to_mesh_file(mesh_file,
822                         point_file,
823                         mesh_output_file,
824                         alpha = 0.0)
825        # load in the .tsh file we just wrote
826        mesh_dic = import_mesh_file(mesh_output_file)
827        #print "mesh_dic",mesh_dic
828        ans =[[0.0, 0.0],
829              [5.0, 10.0],
830              [5.0,10.0]]
831        assert allclose(mesh_dic['vertex_attributes'],ans)
832
833        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
834                        ['elevation','stage'],
835                        'test_fit_to_mesh_file failed')
836
837        #clean up
838        os.remove(mesh_file)
839        os.remove(point_file)
840        os.remove(mesh_output_file)
841
842    def test_fit_to_mesh_fileII(self):
843        from load_mesh.loadASCII import import_mesh_file, \
844             export_mesh_file
845        import tempfile
846        import os
847
848        # create a .tsh file, no user outline
849        mesh_dic = {}
850        mesh_dic['vertices'] = [[0.0, 0.0],
851                                [0.0, 5.0],
852                                [5.0, 0.0]]
853        mesh_dic['triangles'] =  [[0, 2, 1]]
854        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
855        mesh_dic['triangle_tags'] = ['']
856        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
857        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
858        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
859        mesh_dic['segment_tags'] = ['external',
860                                                  'external',
861                                                  'external']
862        mesh_file = tempfile.mktemp(".tsh")
863        export_mesh_file(mesh_file,mesh_dic)
864
865        # create a points .csv file
866        point_file = tempfile.mktemp(".csv")
867        fd = open(point_file,'w')
868        fd.write("x,y,elevation, stage \n\
869        1.0, 1.0,2.,4 \n\
870        1.0, 3.0,4,8 \n\
871        3.0,1.0,4.,8 \n")
872        fd.close()
873
874        mesh_output_file = "new_triangle.tsh"
875        fit_to_mesh_file(mesh_file,
876                         point_file,
877                         mesh_output_file,
878                         alpha = 0.0)
879        # load in the .tsh file we just wrote
880        mesh_dic = import_mesh_file(mesh_output_file)
881
882        assert allclose(mesh_dic['vertex_attributes'],
883                        [[1.0, 2.0,0.0, 0.0],
884                         [1.0, 2.0,5.0, 10.0],
885                         [1.0, 2.0,5.0,10.0]])
886
887        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
888                        ['density', 'temp','elevation','stage'],
889                        'test_fit_to_mesh_file failed')
890
891        #clean up
892        os.remove(mesh_file)
893        os.remove(mesh_output_file)
894        os.remove(point_file)
895
896    def test_fit_to_mesh_file_errors(self):
897        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
898        import tempfile
899        import os
900
901        # create a .tsh file, no user outline
902        mesh_dic = {}
903        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
904        mesh_dic['triangles'] =  [[0, 2, 1]]
905        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
906        mesh_dic['triangle_tags'] = ['']
907        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
908        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
909        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
910        mesh_dic['segment_tags'] = ['external', 'external','external']
911        mesh_file = tempfile.mktemp(".tsh")
912        export_mesh_file(mesh_file,mesh_dic)
913
914        # create a bad points .csv file
915        point_file = tempfile.mktemp(".csv")
916        fd = open(point_file,'w')
917        fd.write("x,y,elevation stage \n\
918        1.0, 1.0,2.,4 \n\
919        1.0, 3.0,4,8 \n\
920        3.0,1.0,4.,8 \n")
921        fd.close()
922
923        mesh_output_file = "new_triangle.tsh"
924        try:
925            fit_to_mesh_file(mesh_file, point_file,
926                             mesh_output_file, display_errors = False)
927        except SyntaxError:
928            pass
929        else:
930            #self.failUnless(0 ==1,  'Bad file did not raise error!')
931            raise 'Bad file did not raise error!'
932           
933        #clean up
934        os.remove(mesh_file)
935        os.remove(point_file)
936
937    def test_fit_to_mesh_file_errorsII(self):
938        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
939        import tempfile
940        import os
941
942        # create a .tsh file, no user outline
943        mesh_file = tempfile.mktemp(".tsh")
944        fd = open(mesh_file,'w')
945        fd.write("unit testing a bad .tsh file \n")
946        fd.close()
947
948        # create a points .csv file
949        point_file = tempfile.mktemp(".csv")
950        fd = open(point_file,'w')
951        fd.write("x,y,elevation, stage \n\
952        1.0, 1.0,2.,4 \n\
953        1.0, 3.0,4,8 \n\
954        3.0,1.0,4.,8 \n")
955        fd.close()
956
957        mesh_output_file = "new_triangle.tsh"
958        try:
959            fit_to_mesh_file(mesh_file, point_file,
960                             mesh_output_file, display_errors = False)
961        except IOError:
962            pass
963        else:
964            raise 'Bad file did not raise error!'
965           
966        #clean up
967        os.remove(mesh_file)
968        os.remove(point_file)
969
970    def test_fit_to_mesh_file_errorsIII(self):
971        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
972        import tempfile
973        import os
974
975        # create a .tsh file, no user outline
976        mesh_dic = {}
977        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
978        mesh_dic['triangles'] =  [[0, 2, 1]]
979        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
980        mesh_dic['triangle_tags'] = ['']
981        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
982        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
983        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
984        mesh_dic['segment_tags'] = ['external', 'external','external']
985        mesh_file = tempfile.mktemp(".tsh")
986        export_mesh_file(mesh_file,mesh_dic)
987
988
989        # create a points .csv file
990        point_file = tempfile.mktemp(".csv")
991        fd = open(point_file,'w')
992        fd.write("x,y,elevation, stage \n\
993        1.0, 1.0,2.,4 \n\
994        1.0, 3.0,4,8 \n\
995        3.0,1.0,4.,8 \n")
996        fd.close()
997
998        #This a deliberately illegal filename to invoke the error.
999        mesh_output_file = ".../\z\z:ya.tsh"       
1000
1001        try:
1002            fit_to_mesh_file(mesh_file, point_file,
1003                             mesh_output_file, display_errors = False)
1004        except IOError:
1005            pass
1006        else:
1007            raise 'Bad file did not raise error!'
1008       
1009        #clean up
1010        os.remove(mesh_file)
1011        os.remove(point_file)
1012 
1013    def Not_yet_test_smooth_att_to_mesh_with_excess_verts(self):
1014
1015        a = [0.0, 0.0]
1016        b = [0.0, 5.0]
1017        c = [5.0, 0.0]
1018        d = [1.0, 1.0]
1019        e = [18.0, 1000.0]
1020       
1021        points = [a, b, c, d, e]
1022        triangles = [ [1,0,2] ]   #bac
1023
1024        d1 = [1.0, 1.0]
1025        d2 = [1.0, 2.0]
1026        d3 = [3.0,1.0]
1027        d4 = [2.0,3.0]
1028        d5 = [2.0,2.0]
1029        d6 = [1.0,3.0]
1030        data_coords = [d1, d2, d3, d4, d5, d6]
1031        z = linear_function(data_coords)
1032        #print "z",z
1033
1034        interp = Fit(points, triangles, alpha=0.0)
1035       
1036        try:
1037            f = interp.fit(data_coords, z)
1038        except VertsWithNoTrianglesError:
1039            pass
1040        else:
1041            raise 'Verts with no triangles did not raise error!'
1042       
1043        #f = interp.fit(data_coords, z)
1044        #answer = linear_function(points)
1045
1046        #  Removing the bad verts that we don't care about
1047        #f = f[0:3]
1048        #answer = answer[0:3]
1049        #assert allclose(f, answer)
1050
1051#-------------------------------------------------------------
1052if __name__ == "__main__":
1053    suite = unittest.makeSuite(Test_Fit,'test')
1054    #suite = unittest.makeSuite(Test_Fit,'cache_test_fit_to_mesh_pts')
1055    #suite = unittest.makeSuite(Test_Fit,'')
1056    runner = unittest.TextTestRunner(verbosity=1)
1057    runner.run(suite)
1058
1059
1060
1061
1062
Note: See TracBrowser for help on using the repository browser.