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

Last change on this file since 1122 was 1097, checked in by duncan, 20 years ago

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