source: inundation/pyvolution/test_data_manager.py @ 2267

Last change on this file since 2267 was 2263, checked in by ole, 19 years ago

Work towards a test_all.py at the inundation directory level. Most tests pass except a few in least_squares, data_manager and cg_solve...
File size: 117.2 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import copy
6from Numeric import zeros, array, allclose, Float
7from util import mean
8import tempfile
9import os
10from Scientific.IO.NetCDF import NetCDFFile
11
12from data_manager import *
13from shallow_water import *
14from config import epsilon
15import data_manager
16
17from coordinate_transforms.geo_reference import Geo_reference
18
19class Test_Data_Manager(unittest.TestCase):
20    def setUp(self):
21        import time
22        from mesh_factory import rectangular
23
24        #Create basic mesh
25        points, vertices, boundary = rectangular(2, 2)
26
27        #Create shallow water domain
28        domain = Domain(points, vertices, boundary)
29        domain.default_order=2
30
31
32        #Set some field values
33        domain.set_quantity('elevation', lambda x,y: -x)
34        domain.set_quantity('friction', 0.03)
35
36
37        ######################
38        # Boundary conditions
39        B = Transmissive_boundary(domain)
40        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
41
42
43        ######################
44        #Initial condition - with jumps
45
46
47        bed = domain.quantities['elevation'].vertex_values
48        stage = zeros(bed.shape, Float)
49
50        h = 0.3
51        for i in range(stage.shape[0]):
52            if i % 2 == 0:
53                stage[i,:] = bed[i,:] + h
54            else:
55                stage[i,:] = bed[i,:]
56
57        domain.set_quantity('stage', stage)
58        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
59
60        domain.distribute_to_vertices_and_edges()
61
62
63        self.domain = domain
64
65        C = domain.get_vertex_coordinates()
66        self.X = C[:,0:6:2].copy()
67        self.Y = C[:,1:6:2].copy()
68
69        self.F = bed
70
71
72
73
74        #Write A testfile
75
76        self.test_MOST_file = 'xxxsmall_ha.nc'
77
78        longitudes = [150.66667, 150.83334, 151., 151.16667]
79        latitudes = [-34.5, -34.33333, -34.16667, -34]
80
81        long_name = 'LON'
82        lat_name = 'LAT'
83
84        nx = 4
85        ny = 4
86        six = 6
87
88
89        fid = NetCDFFile(self.test_MOST_file, 'w')
90       
91        fid.createDimension(long_name,nx)
92        fid.createVariable(long_name,'d',(long_name,))
93        fid.variables[long_name].point_spacing='uneven'
94        fid.variables[long_name].units='degrees_east'
95        fid.variables[long_name].assignValue(longitudes)
96
97        fid.createDimension(lat_name,ny)
98        fid.createVariable(lat_name,'d',(lat_name,))
99        fid.variables[lat_name].point_spacing='uneven'
100        fid.variables[lat_name].units='degrees_north'
101        fid.variables[lat_name].assignValue(latitudes)
102
103        fid.createDimension('TIME',six)
104        fid.createVariable('TIME','d',('TIME',))
105        fid.variables['TIME'].point_spacing='uneven'
106        fid.variables['TIME'].units='seconds'
107        fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
108
109        fid.createVariable('HA','d',('TIME', lat_name, long_name))
110        #fid.createVariable('HA','d',(lat_name, long_name))
111        fid.variables['HA'].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
112                                          [-0.1214216, 0, 0, 0],
113                                          [0, 0, 0, 0],
114                                          [0, 0, 0, 0]],
115                                         [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
116                                          [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
117                                          [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
118                                          [0, 2.291054e-005, 2.291054e-005, 0]],
119                                         [[0.3400644, 0.0001374632, -23.31503, -6.50198],
120                                          [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
121                                          [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
122                                          [0, 0.0001374632, 0.0001374632, 0]],
123                                         [[0.3400644, 0.0002520159, -23.29672, -6.50198],
124                                          [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
125                                          [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
126                                          [0, 0.0002520159, 0.0002520159, 0]],
127                                         [[0.3400644, 0.0003665686, -23.27842, -6.50198],
128                                          [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
129                                          [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
130                                          [0, 0.0003665686, 0.0003665686, 0]],
131                                         [[0.3400644, 0.0004811212, -23.26012, -6.50198],
132                                          [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
133                                          [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
134                                          [0, 0.0004811212, 0.0004811212, 0]]])
135           
136       
137        fid.close()
138       
139
140
141    def tearDown(self):
142        import os
143        #os.remove(self.test_MOST_file)
144       
145
146    def test_sww_constant(self):
147        """Test that constant sww information can be written correctly
148        (non smooth)
149        """
150
151        import time, os
152        from Numeric import array, zeros, allclose, Float, concatenate
153        from Scientific.IO.NetCDF import NetCDFFile
154
155        self.domain.filename = 'datatest' + str(id(self))
156        self.domain.format = 'sww'
157        self.domain.smooth = False
158
159        sww = get_dataobject(self.domain)
160        sww.store_connectivity()
161
162        #Check contents
163        #Get NetCDF
164        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
165
166        # Get the variables
167        x = fid.variables['x']
168        y = fid.variables['y']
169        z = fid.variables['elevation']
170
171        volumes = fid.variables['volumes']
172
173
174        assert allclose (x[:], self.X.flat)
175        assert allclose (y[:], self.Y.flat)
176        assert allclose (z[:], self.F.flat)
177
178        V = volumes
179
180        P = len(self.domain)
181        for k in range(P):
182            assert V[k, 0] == 3*k
183            assert V[k, 1] == 3*k+1
184            assert V[k, 2] == 3*k+2
185
186
187        fid.close()
188
189        #Cleanup
190        os.remove(sww.filename)
191
192
193    def test_sww_constant_smooth(self):
194        """Test that constant sww information can be written correctly
195        (non smooth)
196        """
197
198        import time, os
199        from Numeric import array, zeros, allclose, Float, concatenate
200        from Scientific.IO.NetCDF import NetCDFFile
201
202        self.domain.filename = 'datatest' + str(id(self))
203        self.domain.format = 'sww'
204        self.domain.smooth = True
205
206        sww = get_dataobject(self.domain)
207        sww.store_connectivity()
208
209        #Check contents
210        #Get NetCDF
211        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
212
213        # Get the variables
214        x = fid.variables['x']
215        y = fid.variables['y']
216        z = fid.variables['elevation']
217
218        volumes = fid.variables['volumes']
219
220        X = x[:]
221        Y = y[:]
222
223        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
224        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
225        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
226
227        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
228
229        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
230
231        Z = z[:]
232        assert Z[4] == -0.5
233
234        V = volumes
235        assert V[2,0] == 4
236        assert V[2,1] == 5
237        assert V[2,2] == 1
238
239        assert V[4,0] == 6
240        assert V[4,1] == 7
241        assert V[4,2] == 3
242
243
244        fid.close()
245
246        #Cleanup
247        os.remove(sww.filename)
248
249
250
251    def test_sww_variable(self):
252        """Test that sww information can be written correctly
253        """
254
255        import time, os
256        from Numeric import array, zeros, allclose, Float, concatenate
257        from Scientific.IO.NetCDF import NetCDFFile
258
259        self.domain.filename = 'datatest' + str(id(self))
260        self.domain.format = 'sww'
261        self.domain.smooth = True
262        self.domain.reduction = mean
263
264        sww = get_dataobject(self.domain)
265        sww.store_connectivity()
266        sww.store_timestep('stage')
267
268        #Check contents
269        #Get NetCDF
270        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
271
272
273        # Get the variables
274        x = fid.variables['x']
275        y = fid.variables['y']
276        z = fid.variables['elevation']
277        time = fid.variables['time']
278        stage = fid.variables['stage']
279
280
281        Q = self.domain.quantities['stage']
282        Q0 = Q.vertex_values[:,0]
283        Q1 = Q.vertex_values[:,1]
284        Q2 = Q.vertex_values[:,2]
285
286        A = stage[0,:]
287        #print A[0], (Q2[0,0] + Q1[1,0])/2
288        assert allclose(A[0], (Q2[0] + Q1[1])/2)
289        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
290        assert allclose(A[2], Q0[3])
291        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
292
293        #Center point
294        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
295                                 Q0[5] + Q2[6] + Q1[7])/6)
296
297
298
299        fid.close()
300
301        #Cleanup
302        os.remove(sww.filename)
303
304
305    def test_sww_variable2(self):
306        """Test that sww information can be written correctly
307        multiple timesteps. Use average as reduction operator
308        """
309
310        import time, os
311        from Numeric import array, zeros, allclose, Float, concatenate
312        from Scientific.IO.NetCDF import NetCDFFile
313
314        self.domain.filename = 'datatest' + str(id(self))
315        self.domain.format = 'sww'
316        self.domain.smooth = True
317
318        self.domain.reduction = mean
319
320        sww = get_dataobject(self.domain)
321        sww.store_connectivity()
322        sww.store_timestep('stage')
323        self.domain.evolve_to_end(finaltime = 0.01)
324        sww.store_timestep('stage')
325
326
327        #Check contents
328        #Get NetCDF
329        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
330
331        # Get the variables
332        x = fid.variables['x']
333        y = fid.variables['y']
334        z = fid.variables['elevation']
335        time = fid.variables['time']
336        stage = fid.variables['stage']
337
338        #Check values
339        Q = self.domain.quantities['stage']
340        Q0 = Q.vertex_values[:,0]
341        Q1 = Q.vertex_values[:,1]
342        Q2 = Q.vertex_values[:,2]
343
344        A = stage[1,:]
345        assert allclose(A[0], (Q2[0] + Q1[1])/2)
346        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
347        assert allclose(A[2], Q0[3])
348        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
349
350        #Center point
351        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
352                                 Q0[5] + Q2[6] + Q1[7])/6)
353
354
355        fid.close()
356
357        #Cleanup
358        os.remove(sww.filename)
359
360    def test_sww_variable3(self):
361        """Test that sww information can be written correctly
362        multiple timesteps using a different reduction operator (min)
363        """
364
365        import time, os
366        from Numeric import array, zeros, allclose, Float, concatenate
367        from Scientific.IO.NetCDF import NetCDFFile
368
369        self.domain.filename = 'datatest' + str(id(self))
370        self.domain.format = 'sww'
371        self.domain.smooth = True
372        self.domain.reduction = min
373
374        sww = get_dataobject(self.domain)
375        sww.store_connectivity()
376        sww.store_timestep('stage')
377
378        self.domain.evolve_to_end(finaltime = 0.01)
379        sww.store_timestep('stage')
380
381
382        #Check contents
383        #Get NetCDF
384        fid = NetCDFFile(sww.filename, 'r')
385
386
387        # Get the variables
388        x = fid.variables['x']
389        y = fid.variables['y']
390        z = fid.variables['elevation']
391        time = fid.variables['time']
392        stage = fid.variables['stage']
393
394        #Check values
395        Q = self.domain.quantities['stage']
396        Q0 = Q.vertex_values[:,0]
397        Q1 = Q.vertex_values[:,1]
398        Q2 = Q.vertex_values[:,2]
399
400        A = stage[1,:]
401        assert allclose(A[0], min(Q2[0], Q1[1]))
402        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
403        assert allclose(A[2], Q0[3])
404        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
405
406        #Center point
407        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
408                                  Q0[5], Q2[6], Q1[7]))
409
410
411        fid.close()
412
413        #Cleanup
414        os.remove(sww.filename)
415
416
417    def test_sync(self):
418        """Test info stored at each timestep is as expected (incl initial condition)
419        """
420
421        import time, os, config
422        from Numeric import array, zeros, allclose, Float, concatenate
423        from Scientific.IO.NetCDF import NetCDFFile
424
425        self.domain.filename = 'synctest'
426        self.domain.format = 'sww'
427        self.domain.smooth = False
428        self.domain.store = True
429        self.domain.beta_h = 0
430
431        #Evolution
432        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
433            stage = self.domain.quantities['stage'].vertex_values
434
435            #Get NetCDF
436            fid = NetCDFFile(self.domain.writer.filename, 'r')
437            stage_file = fid.variables['stage']
438
439            if t == 0.0:
440                assert allclose(stage, self.initial_stage)
441                assert allclose(stage_file[:], stage.flat)
442            else:
443                assert not allclose(stage, self.initial_stage)
444                assert not allclose(stage_file[:], stage.flat)
445
446            fid.close()
447
448        os.remove(self.domain.writer.filename)
449
450
451
452    def test_sww_DSG(self):
453        """Not a test, rather a look at the sww format
454        """
455
456        import time, os
457        from Numeric import array, zeros, allclose, Float, concatenate
458        from Scientific.IO.NetCDF import NetCDFFile
459
460        self.domain.filename = 'datatest' + str(id(self))
461        self.domain.format = 'sww'
462        self.domain.smooth = True
463        self.domain.reduction = mean
464
465        sww = get_dataobject(self.domain)
466        sww.store_connectivity()
467        sww.store_timestep('stage')
468
469        #Check contents
470        #Get NetCDF
471        fid = NetCDFFile(sww.filename, 'r')
472
473        # Get the variables
474        x = fid.variables['x']
475        y = fid.variables['y']
476        z = fid.variables['elevation']
477
478        volumes = fid.variables['volumes']
479        time = fid.variables['time']
480
481        # 2D
482        stage = fid.variables['stage']
483
484        X = x[:]
485        Y = y[:]
486        Z = z[:]
487        V = volumes[:]
488        T = time[:]
489        S = stage[:,:]
490
491#         print "****************************"
492#         print "X ",X
493#         print "****************************"
494#         print "Y ",Y
495#         print "****************************"
496#         print "Z ",Z
497#         print "****************************"
498#         print "V ",V
499#         print "****************************"
500#         print "Time ",T
501#         print "****************************"
502#         print "Stage ",S
503#         print "****************************"
504
505
506        fid.close()
507
508        #Cleanup
509        os.remove(sww.filename)
510
511
512    #def test_write_pts(self):
513    #    #Obsolete
514    #   
515    #    #Get (enough) datapoints
516    #
517    #    from Numeric import array
518    #    points = array([[ 0.66666667, 0.66666667],
519    #                    [ 1.33333333, 1.33333333],
520    #                    [ 2.66666667, 0.66666667],
521    #                    [ 0.66666667, 2.66666667],
522    #                    [ 0.0, 1.0],
523    #                    [ 0.0, 3.0],
524    #                    [ 1.0, 0.0],
525    #                    [ 1.0, 1.0],
526    #                    [ 1.0, 2.0],
527    #                    [ 1.0, 3.0],
528    #                    [ 2.0, 1.0],
529    #                    [ 3.0, 0.0],
530    #                    [ 3.0, 1.0]])
531    #
532    #    z = points[:,0] + 2*points[:,1]
533    #
534    #    ptsfile = 'testptsfile.pts'
535    #    write_ptsfile(ptsfile, points, z,
536    #                  attribute_name = 'linear_combination')
537    #
538    #    #Check contents
539    #    #Get NetCDF
540    #    from Scientific.IO.NetCDF import NetCDFFile       
541    #    fid = NetCDFFile(ptsfile, 'r')
542    #
543    #    # Get the variables
544    #    #print fid.variables.keys()
545    #    points1 = fid.variables['points']
546    #    z1 = fid.variables['linear_combination']
547    #
548    #    #Check values#
549    #
550    #    #print points[:]
551    #    #print ref_points
552    #    assert allclose(points, points1) 
553    #
554    #    #print attributes[:]
555    #    #print ref_elevation
556    #    assert allclose(z, z1)
557    #
558    #    #Cleanup
559    #    fid.close()
560    #
561    #    import os
562    #    os.remove(ptsfile)
563       
564
565
566
567    def test_dem2pts(self):
568        """Test conversion from dem in ascii format to native NetCDF xya format
569        """
570
571        import time, os
572        from Numeric import array, zeros, allclose, Float, concatenate
573        from Scientific.IO.NetCDF import NetCDFFile
574
575        #Write test asc file
576        root = 'demtest'
577
578        filename = root+'.asc'
579        fid = open(filename, 'w')
580        fid.write("""ncols         5
581nrows         6
582xllcorner     2000.5
583yllcorner     3000.5
584cellsize      25
585NODATA_value  -9999
586""")
587        #Create linear function
588
589        ref_points = []
590        ref_elevation = []
591        for i in range(6):
592            y = (6-i)*25.0
593            for j in range(5):
594                x = j*25.0
595                z = x+2*y
596
597                ref_points.append( [x,y] )
598                ref_elevation.append(z)
599                fid.write('%f ' %z)
600            fid.write('\n')
601
602        fid.close()
603
604        #Write prj file with metadata
605        metafilename = root+'.prj'
606        fid = open(metafilename, 'w')
607
608
609        fid.write("""Projection UTM
610Zone 56
611Datum WGS84
612Zunits NO
613Units METERS
614Spheroid WGS84
615Xshift 0.0000000000
616Yshift 10000000.0000000000
617Parameters
618""")
619        fid.close()
620
621        #Convert to NetCDF pts
622        convert_dem_from_ascii2netcdf(root)
623        dem2pts(root)
624
625        #Check contents
626        #Get NetCDF
627        fid = NetCDFFile(root+'.pts', 'r')
628
629        # Get the variables
630        #print fid.variables.keys()
631        points = fid.variables['points']
632        elevation = fid.variables['elevation']
633
634        #Check values
635
636        #print points[:]
637        #print ref_points
638        assert allclose(points, ref_points)
639
640        #print attributes[:]
641        #print ref_elevation
642        assert allclose(elevation, ref_elevation)
643
644        #Cleanup
645        fid.close()
646
647
648        os.remove(root + '.pts')
649        os.remove(root + '.dem')
650        os.remove(root + '.asc')
651        os.remove(root + '.prj')
652
653
654
655    def test_dem2pts_bounding_box(self):
656        """Test conversion from dem in ascii format to native NetCDF xya format
657        """
658
659        import time, os
660        from Numeric import array, zeros, allclose, Float, concatenate
661        from Scientific.IO.NetCDF import NetCDFFile
662
663        #Write test asc file
664        root = 'demtest'
665
666        filename = root+'.asc'
667        fid = open(filename, 'w')
668        fid.write("""ncols         5
669nrows         6
670xllcorner     2000.5
671yllcorner     3000.5
672cellsize      25
673NODATA_value  -9999
674""")
675        #Create linear function
676
677        ref_points = []
678        ref_elevation = []
679        for i in range(6):
680            y = (6-i)*25.0
681            for j in range(5):
682                x = j*25.0
683                z = x+2*y
684
685                ref_points.append( [x,y] )
686                ref_elevation.append(z)
687                fid.write('%f ' %z)
688            fid.write('\n')
689
690        fid.close()
691
692        #Write prj file with metadata
693        metafilename = root+'.prj'
694        fid = open(metafilename, 'w')
695
696
697        fid.write("""Projection UTM
698Zone 56
699Datum WGS84
700Zunits NO
701Units METERS
702Spheroid WGS84
703Xshift 0.0000000000
704Yshift 10000000.0000000000
705Parameters
706""")
707        fid.close()
708
709        #Convert to NetCDF pts
710        convert_dem_from_ascii2netcdf(root)
711        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
712                northing_min=3035.0, northing_max=3125.5)
713
714        #Check contents
715        #Get NetCDF
716        fid = NetCDFFile(root+'.pts', 'r')
717
718        # Get the variables
719        #print fid.variables.keys()
720        points = fid.variables['points']
721        elevation = fid.variables['elevation']
722
723        #Check values
724        assert fid.xllcorner[0] == 2010.0
725        assert fid.yllcorner[0] == 3035.0
726
727        #create new reference points
728        ref_points = []
729        ref_elevation = []
730        for i in range(4):
731            y = (4-i)*25.0 + 25.0
732            y_new = y + 3000.5 - 3035.0
733            for j in range(4):
734                x = j*25.0 + 25.0
735                x_new = x + 2000.5 - 2010.0
736                z = x+2*y
737
738                ref_points.append( [x_new,y_new] )
739                ref_elevation.append(z)
740
741        #print points[:]
742        #print ref_points
743        assert allclose(points, ref_points)
744
745        #print attributes[:]
746        #print ref_elevation
747        assert allclose(elevation, ref_elevation)
748
749        #Cleanup
750        fid.close()
751
752
753        os.remove(root + '.pts')
754        os.remove(root + '.dem')
755        os.remove(root + '.asc')
756        os.remove(root + '.prj')
757
758
759
760    def test_hecras_cross_sections2pts(self):
761        """Test conversion from HECRAS cross sections in ascii format
762        to native NetCDF pts format
763        """
764
765        import time, os
766        from Numeric import array, zeros, allclose, Float, concatenate
767        from Scientific.IO.NetCDF import NetCDFFile
768
769        #Write test asc file
770        root = 'hecrastest'
771
772        filename = root+'.sdf'
773        fid = open(filename, 'w')
774        fid.write("""
775# RAS export file created on Mon 15Aug2005 11:42
776# by HEC-RAS Version 3.1.1
777
778BEGIN HEADER:
779  UNITS: METRIC
780  DTM TYPE: TIN
781  DTM: v:\1\cit\perth_topo\river_tin
782  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
783  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
784  MAP PROJECTION: UTM
785  PROJECTION ZONE: 50
786  DATUM: AGD66
787  VERTICAL DATUM:
788  NUMBER OF REACHES:  19
789  NUMBER OF CROSS-SECTIONS:  2
790END HEADER:
791
792
793BEGIN CROSS-SECTIONS:
794
795  CROSS-SECTION:
796    STREAM ID:Southern-Wungong
797    REACH ID:Southern-Wungong
798    STATION:21410   
799    CUT LINE:
800      407546.08 , 6437277.542
801      407329.32 , 6437489.482
802      407283.11 , 6437541.232
803    SURFACE LINE:
804     407546.08,   6437277.54,   52.14
805     407538.88,   6437284.58,   51.07
806     407531.68,   6437291.62,   50.56
807     407524.48,   6437298.66,   49.58
808     407517.28,   6437305.70,   49.09
809     407510.08,   6437312.74,   48.76
810  END:
811
812  CROSS-SECTION:
813    STREAM ID:Swan River     
814    REACH ID:Swan Mouth     
815    STATION:840.*   
816    CUT LINE:
817      381178.0855 , 6452559.0685
818      380485.4755 , 6453169.272
819    SURFACE LINE:
820     381178.09,   6452559.07,   4.17
821     381169.49,   6452566.64,   4.26
822     381157.78,   6452576.96,   4.34
823     381155.97,   6452578.56,   4.35
824     381143.72,   6452589.35,   4.43
825     381136.69,   6452595.54,   4.58
826     381114.74,   6452614.88,   4.41
827     381075.53,   6452649.43,   4.17
828     381071.47,   6452653.00,   3.99
829     381063.46,   6452660.06,   3.67
830     381054.41,   6452668.03,   3.67
831  END: 
832END CROSS-SECTIONS:
833""")
834
835        fid.close()
836
837
838        #Convert to NetCDF pts
839        hecras_cross_sections2pts(root)
840
841        #Check contents
842        #Get NetCDF
843        fid = NetCDFFile(root+'.pts', 'r')
844
845        # Get the variables
846        #print fid.variables.keys()
847        points = fid.variables['points']
848        elevation = fid.variables['elevation']
849
850        #Check values
851        ref_points = [[407546.08, 6437277.54],
852                      [407538.88, 6437284.58],
853                      [407531.68, 6437291.62],
854                      [407524.48, 6437298.66],
855                      [407517.28, 6437305.70],
856                      [407510.08, 6437312.74]]
857
858        ref_points += [[381178.09, 6452559.07], 
859                       [381169.49, 6452566.64], 
860                       [381157.78, 6452576.96], 
861                       [381155.97, 6452578.56], 
862                       [381143.72, 6452589.35], 
863                       [381136.69, 6452595.54], 
864                       [381114.74, 6452614.88], 
865                       [381075.53, 6452649.43], 
866                       [381071.47, 6452653.00], 
867                       [381063.46, 6452660.06], 
868                       [381054.41, 6452668.03]] 
869       
870                     
871        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
872        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
873
874        #print points[:]
875        #print ref_points
876        assert allclose(points, ref_points)
877
878        #print attributes[:]
879        #print ref_elevation
880        assert allclose(elevation, ref_elevation)
881
882        #Cleanup
883        fid.close()
884
885
886        os.remove(root + '.sdf')
887        os.remove(root + '.pts')
888
889
890
891    def test_sww2dem_asc_elevation(self):
892        """Test that sww information can be converted correctly to asc/prj
893        format readable by e.g. ArcView
894        """
895
896        import time, os
897        from Numeric import array, zeros, allclose, Float, concatenate
898        from Scientific.IO.NetCDF import NetCDFFile
899
900        #Setup
901        self.domain.filename = 'datatest'
902
903        prjfile = self.domain.filename + '_elevation.prj'
904        ascfile = self.domain.filename + '_elevation.asc'
905        swwfile = self.domain.filename + '.sww'
906
907        self.domain.set_datadir('.')
908        self.domain.format = 'sww'
909        self.domain.smooth = True
910        self.domain.set_quantity('elevation', lambda x,y: -x-y)
911
912        self.domain.geo_reference = Geo_reference(56,308500,6189000)
913
914        sww = get_dataobject(self.domain)
915        sww.store_connectivity()
916        sww.store_timestep('stage')
917
918        self.domain.evolve_to_end(finaltime = 0.01)
919        sww.store_timestep('stage')
920
921        cellsize = 0.25
922        #Check contents
923        #Get NetCDF
924
925        fid = NetCDFFile(sww.filename, 'r')
926
927        # Get the variables
928        x = fid.variables['x'][:]
929        y = fid.variables['y'][:]
930        z = fid.variables['elevation'][:]
931        time = fid.variables['time'][:]
932        stage = fid.variables['stage'][:]
933
934
935        #Export to ascii/prj files
936        sww2dem(self.domain.filename,
937                quantity = 'elevation',
938                cellsize = cellsize,
939                verbose = False,
940                format = 'asc')
941               
942        #Check prj (meta data)
943        prjid = open(prjfile)
944        lines = prjid.readlines()
945        prjid.close()
946
947        L = lines[0].strip().split()
948        assert L[0].strip().lower() == 'projection'
949        assert L[1].strip().lower() == 'utm'
950
951        L = lines[1].strip().split()
952        assert L[0].strip().lower() == 'zone'
953        assert L[1].strip().lower() == '56'
954
955        L = lines[2].strip().split()
956        assert L[0].strip().lower() == 'datum'
957        assert L[1].strip().lower() == 'wgs84'
958
959        L = lines[3].strip().split()
960        assert L[0].strip().lower() == 'zunits'
961        assert L[1].strip().lower() == 'no'
962
963        L = lines[4].strip().split()
964        assert L[0].strip().lower() == 'units'
965        assert L[1].strip().lower() == 'meters'
966
967        L = lines[5].strip().split()
968        assert L[0].strip().lower() == 'spheroid'
969        assert L[1].strip().lower() == 'wgs84'
970
971        L = lines[6].strip().split()
972        assert L[0].strip().lower() == 'xshift'
973        assert L[1].strip().lower() == '500000'
974
975        L = lines[7].strip().split()
976        assert L[0].strip().lower() == 'yshift'
977        assert L[1].strip().lower() == '10000000'
978
979        L = lines[8].strip().split()
980        assert L[0].strip().lower() == 'parameters'
981
982
983        #Check asc file
984        ascid = open(ascfile)
985        lines = ascid.readlines()
986        ascid.close()
987
988        L = lines[0].strip().split()
989        assert L[0].strip().lower() == 'ncols'
990        assert L[1].strip().lower() == '5'
991
992        L = lines[1].strip().split()
993        assert L[0].strip().lower() == 'nrows'
994        assert L[1].strip().lower() == '5'
995
996        L = lines[2].strip().split()
997        assert L[0].strip().lower() == 'xllcorner'
998        assert allclose(float(L[1].strip().lower()), 308500)
999
1000        L = lines[3].strip().split()
1001        assert L[0].strip().lower() == 'yllcorner'
1002        assert allclose(float(L[1].strip().lower()), 6189000)
1003
1004        L = lines[4].strip().split()
1005        assert L[0].strip().lower() == 'cellsize'
1006        assert allclose(float(L[1].strip().lower()), cellsize)
1007
1008        L = lines[5].strip().split()
1009        assert L[0].strip() == 'NODATA_value'
1010        assert L[1].strip().lower() == '-9999'
1011
1012        #Check grid values
1013        for j in range(5):
1014            L = lines[6+j].strip().split()
1015            y = (4-j) * cellsize
1016            for i in range(5):
1017                assert allclose(float(L[i]), -i*cellsize - y)
1018
1019
1020        fid.close()
1021
1022        #Cleanup
1023        os.remove(prjfile)
1024        os.remove(ascfile)
1025        os.remove(swwfile)
1026
1027
1028
1029    def test_sww2dem_larger(self):
1030        """Test that sww information can be converted correctly to asc/prj
1031        format readable by e.g. ArcView. Here:
1032
1033        ncols         11
1034        nrows         11
1035        xllcorner     308500
1036        yllcorner     6189000
1037        cellsize      10.000000
1038        NODATA_value  -9999
1039        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1040         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1041         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1042         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
1043         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
1044         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
1045         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
1046         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
1047         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
1048         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
1049           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
1050       
1051        """
1052
1053        import time, os
1054        from Numeric import array, zeros, allclose, Float, concatenate
1055        from Scientific.IO.NetCDF import NetCDFFile
1056
1057        #Setup
1058
1059        from mesh_factory import rectangular
1060
1061        #Create basic mesh (100m x 100m)
1062        points, vertices, boundary = rectangular(2, 2, 100, 100)
1063
1064        #Create shallow water domain
1065        domain = Domain(points, vertices, boundary)
1066        domain.default_order = 2
1067
1068        domain.filename = 'datatest'
1069
1070        prjfile = domain.filename + '_elevation.prj'
1071        ascfile = domain.filename + '_elevation.asc'
1072        swwfile = domain.filename + '.sww'
1073
1074        domain.set_datadir('.')
1075        domain.format = 'sww'
1076        domain.smooth = True
1077        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1078
1079        #
1080        domain.set_quantity('elevation', lambda x,y: -x-y)
1081        domain.set_quantity('stage', 0)       
1082
1083        B = Transmissive_boundary(domain)
1084        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1085
1086
1087        #
1088        sww = get_dataobject(domain)
1089        sww.store_connectivity()
1090        sww.store_timestep('stage')
1091
1092        domain.evolve_to_end(finaltime = 0.01)
1093        sww.store_timestep('stage')
1094
1095        cellsize = 10  #10m grid
1096
1097       
1098        #Check contents
1099        #Get NetCDF
1100
1101        fid = NetCDFFile(sww.filename, 'r')
1102
1103        # Get the variables
1104        x = fid.variables['x'][:]
1105        y = fid.variables['y'][:]
1106        z = fid.variables['elevation'][:]
1107        time = fid.variables['time'][:]
1108        stage = fid.variables['stage'][:]
1109
1110
1111        #Export to ascii/prj files
1112        sww2dem(domain.filename,
1113                quantity = 'elevation',
1114                cellsize = cellsize,
1115                verbose = False,
1116                format = 'asc')
1117       
1118               
1119        #Check prj (meta data)
1120        prjid = open(prjfile)
1121        lines = prjid.readlines()
1122        prjid.close()
1123
1124        L = lines[0].strip().split()
1125        assert L[0].strip().lower() == 'projection'
1126        assert L[1].strip().lower() == 'utm'
1127
1128        L = lines[1].strip().split()
1129        assert L[0].strip().lower() == 'zone'
1130        assert L[1].strip().lower() == '56'
1131
1132        L = lines[2].strip().split()
1133        assert L[0].strip().lower() == 'datum'
1134        assert L[1].strip().lower() == 'wgs84'
1135
1136        L = lines[3].strip().split()
1137        assert L[0].strip().lower() == 'zunits'
1138        assert L[1].strip().lower() == 'no'
1139
1140        L = lines[4].strip().split()
1141        assert L[0].strip().lower() == 'units'
1142        assert L[1].strip().lower() == 'meters'
1143
1144        L = lines[5].strip().split()
1145        assert L[0].strip().lower() == 'spheroid'
1146        assert L[1].strip().lower() == 'wgs84'
1147
1148        L = lines[6].strip().split()
1149        assert L[0].strip().lower() == 'xshift'
1150        assert L[1].strip().lower() == '500000'
1151
1152        L = lines[7].strip().split()
1153        assert L[0].strip().lower() == 'yshift'
1154        assert L[1].strip().lower() == '10000000'
1155
1156        L = lines[8].strip().split()
1157        assert L[0].strip().lower() == 'parameters'
1158
1159
1160        #Check asc file
1161        ascid = open(ascfile)
1162        lines = ascid.readlines()
1163        ascid.close()
1164
1165        L = lines[0].strip().split()
1166        assert L[0].strip().lower() == 'ncols'
1167        assert L[1].strip().lower() == '11'
1168
1169        L = lines[1].strip().split()
1170        assert L[0].strip().lower() == 'nrows'
1171        assert L[1].strip().lower() == '11'
1172
1173        L = lines[2].strip().split()
1174        assert L[0].strip().lower() == 'xllcorner'
1175        assert allclose(float(L[1].strip().lower()), 308500)
1176
1177        L = lines[3].strip().split()
1178        assert L[0].strip().lower() == 'yllcorner'
1179        assert allclose(float(L[1].strip().lower()), 6189000)
1180
1181        L = lines[4].strip().split()
1182        assert L[0].strip().lower() == 'cellsize'
1183        assert allclose(float(L[1].strip().lower()), cellsize)
1184
1185        L = lines[5].strip().split()
1186        assert L[0].strip() == 'NODATA_value'
1187        assert L[1].strip().lower() == '-9999'
1188
1189        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1190        for i, line in enumerate(lines[6:]):
1191            for j, value in enumerate( line.split() ):
1192                #assert allclose(float(value), -(10-i+j)*cellsize)
1193                assert float(value) == -(10-i+j)*cellsize
1194       
1195
1196        fid.close()
1197
1198        #Cleanup
1199        os.remove(prjfile)
1200        os.remove(ascfile)
1201        os.remove(swwfile)
1202
1203
1204
1205    def test_sww2dem_boundingbox(self):
1206        """Test that sww information can be converted correctly to asc/prj
1207        format readable by e.g. ArcView.
1208        This will test that mesh can be restricted by bounding box
1209
1210        Original extent is 100m x 100m:
1211
1212        Eastings:   308500 -  308600
1213        Northings: 6189000 - 6189100
1214
1215        Bounding box changes this to the 50m x 50m square defined by
1216
1217        Eastings:   308530 -  308570
1218        Northings: 6189050 - 6189100
1219
1220        The cropped values should be
1221
1222         -130 -140 -150 -160 -170 
1223         -120 -130 -140 -150 -160 
1224         -110 -120 -130 -140 -150 
1225         -100 -110 -120 -130 -140 
1226          -90 -100 -110 -120 -130 
1227          -80  -90 -100 -110 -120
1228
1229        and the new lower reference point should be   
1230        Eastings:   308530
1231        Northings: 6189050
1232       
1233        Original dataset is the same as in test_sww2dem_larger()
1234       
1235        """
1236
1237        import time, os
1238        from Numeric import array, zeros, allclose, Float, concatenate
1239        from Scientific.IO.NetCDF import NetCDFFile
1240
1241        #Setup
1242
1243        from mesh_factory import rectangular
1244
1245        #Create basic mesh (100m x 100m)
1246        points, vertices, boundary = rectangular(2, 2, 100, 100)
1247
1248        #Create shallow water domain
1249        domain = Domain(points, vertices, boundary)
1250        domain.default_order = 2
1251
1252        domain.filename = 'datatest'
1253
1254        prjfile = domain.filename + '_elevation.prj'
1255        ascfile = domain.filename + '_elevation.asc'
1256        swwfile = domain.filename + '.sww'
1257
1258        domain.set_datadir('.')
1259        domain.format = 'sww'
1260        domain.smooth = True
1261        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1262
1263        #
1264        domain.set_quantity('elevation', lambda x,y: -x-y)
1265        domain.set_quantity('stage', 0)       
1266
1267        B = Transmissive_boundary(domain)
1268        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1269
1270
1271        #
1272        sww = get_dataobject(domain)
1273        sww.store_connectivity()
1274        sww.store_timestep('stage')
1275
1276        domain.evolve_to_end(finaltime = 0.01)
1277        sww.store_timestep('stage')
1278
1279        cellsize = 10  #10m grid
1280
1281       
1282        #Check contents
1283        #Get NetCDF
1284
1285        fid = NetCDFFile(sww.filename, 'r')
1286
1287        # Get the variables
1288        x = fid.variables['x'][:]
1289        y = fid.variables['y'][:]
1290        z = fid.variables['elevation'][:]
1291        time = fid.variables['time'][:]
1292        stage = fid.variables['stage'][:]
1293
1294
1295        #Export to ascii/prj files
1296        sww2dem(domain.filename,
1297                quantity = 'elevation',
1298                cellsize = cellsize,
1299                easting_min = 308530,
1300                easting_max = 308570,
1301                northing_min = 6189050, 
1302                northing_max = 6189100,                               
1303                verbose = False,
1304                format = 'asc')
1305       
1306        fid.close()
1307
1308       
1309        #Check prj (meta data)
1310        prjid = open(prjfile)
1311        lines = prjid.readlines()
1312        prjid.close()
1313
1314        L = lines[0].strip().split()
1315        assert L[0].strip().lower() == 'projection'
1316        assert L[1].strip().lower() == 'utm'
1317
1318        L = lines[1].strip().split()
1319        assert L[0].strip().lower() == 'zone'
1320        assert L[1].strip().lower() == '56'
1321
1322        L = lines[2].strip().split()
1323        assert L[0].strip().lower() == 'datum'
1324        assert L[1].strip().lower() == 'wgs84'
1325
1326        L = lines[3].strip().split()
1327        assert L[0].strip().lower() == 'zunits'
1328        assert L[1].strip().lower() == 'no'
1329
1330        L = lines[4].strip().split()
1331        assert L[0].strip().lower() == 'units'
1332        assert L[1].strip().lower() == 'meters'
1333
1334        L = lines[5].strip().split()
1335        assert L[0].strip().lower() == 'spheroid'
1336        assert L[1].strip().lower() == 'wgs84'
1337
1338        L = lines[6].strip().split()
1339        assert L[0].strip().lower() == 'xshift'
1340        assert L[1].strip().lower() == '500000'
1341
1342        L = lines[7].strip().split()
1343        assert L[0].strip().lower() == 'yshift'
1344        assert L[1].strip().lower() == '10000000'
1345
1346        L = lines[8].strip().split()
1347        assert L[0].strip().lower() == 'parameters'
1348
1349
1350        #Check asc file
1351        ascid = open(ascfile)
1352        lines = ascid.readlines()
1353        ascid.close()
1354
1355        L = lines[0].strip().split()
1356        assert L[0].strip().lower() == 'ncols'
1357        assert L[1].strip().lower() == '5'
1358
1359        L = lines[1].strip().split()
1360        assert L[0].strip().lower() == 'nrows'
1361        assert L[1].strip().lower() == '6'
1362
1363        L = lines[2].strip().split()
1364        assert L[0].strip().lower() == 'xllcorner'
1365        assert allclose(float(L[1].strip().lower()), 308530)
1366
1367        L = lines[3].strip().split()
1368        assert L[0].strip().lower() == 'yllcorner'
1369        assert allclose(float(L[1].strip().lower()), 6189050)
1370
1371        L = lines[4].strip().split()
1372        assert L[0].strip().lower() == 'cellsize'
1373        assert allclose(float(L[1].strip().lower()), cellsize)
1374
1375        L = lines[5].strip().split()
1376        assert L[0].strip() == 'NODATA_value'
1377        assert L[1].strip().lower() == '-9999'
1378
1379        #Check grid values
1380        for i, line in enumerate(lines[6:]):
1381            for j, value in enumerate( line.split() ):
1382                #assert float(value) == -(10-i+j)*cellsize               
1383                assert float(value) == -(10-i+j+3)*cellsize
1384       
1385
1386
1387        #Cleanup
1388        os.remove(prjfile)
1389        os.remove(ascfile)
1390        os.remove(swwfile)
1391
1392
1393
1394    def test_sww2dem_asc_stage_reduction(self):
1395        """Test that sww information can be converted correctly to asc/prj
1396        format readable by e.g. ArcView
1397
1398        This tests the reduction of quantity stage using min
1399        """
1400
1401        import time, os
1402        from Numeric import array, zeros, allclose, Float, concatenate
1403        from Scientific.IO.NetCDF import NetCDFFile
1404
1405        #Setup
1406        self.domain.filename = 'datatest'
1407
1408        prjfile = self.domain.filename + '_stage.prj'
1409        ascfile = self.domain.filename + '_stage.asc'
1410        swwfile = self.domain.filename + '.sww'
1411
1412        self.domain.set_datadir('.')
1413        self.domain.format = 'sww'
1414        self.domain.smooth = True
1415        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1416
1417        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1418
1419
1420        sww = get_dataobject(self.domain)
1421        sww.store_connectivity()
1422        sww.store_timestep('stage')
1423
1424        self.domain.evolve_to_end(finaltime = 0.01)
1425        sww.store_timestep('stage')
1426
1427        cellsize = 0.25
1428        #Check contents
1429        #Get NetCDF
1430
1431        fid = NetCDFFile(sww.filename, 'r')
1432
1433        # Get the variables
1434        x = fid.variables['x'][:]
1435        y = fid.variables['y'][:]
1436        z = fid.variables['elevation'][:]
1437        time = fid.variables['time'][:]
1438        stage = fid.variables['stage'][:]
1439
1440
1441        #Export to ascii/prj files
1442        sww2dem(self.domain.filename,
1443                quantity = 'stage',
1444                cellsize = cellsize,
1445                reduction = min,
1446                format = 'asc')
1447
1448
1449        #Check asc file
1450        ascid = open(ascfile)
1451        lines = ascid.readlines()
1452        ascid.close()
1453
1454        L = lines[0].strip().split()
1455        assert L[0].strip().lower() == 'ncols'
1456        assert L[1].strip().lower() == '5'
1457
1458        L = lines[1].strip().split()
1459        assert L[0].strip().lower() == 'nrows'
1460        assert L[1].strip().lower() == '5'
1461
1462        L = lines[2].strip().split()
1463        assert L[0].strip().lower() == 'xllcorner'
1464        assert allclose(float(L[1].strip().lower()), 308500)
1465
1466        L = lines[3].strip().split()
1467        assert L[0].strip().lower() == 'yllcorner'
1468        assert allclose(float(L[1].strip().lower()), 6189000)
1469
1470        L = lines[4].strip().split()
1471        assert L[0].strip().lower() == 'cellsize'
1472        assert allclose(float(L[1].strip().lower()), cellsize)
1473
1474        L = lines[5].strip().split()
1475        assert L[0].strip() == 'NODATA_value'
1476        assert L[1].strip().lower() == '-9999'
1477
1478
1479        #Check grid values (where applicable)
1480        for j in range(5):
1481            if j%2 == 0:
1482                L = lines[6+j].strip().split()
1483                jj = 4-j
1484                for i in range(5):
1485                    if i%2 == 0:
1486                        index = jj/2 + i/2*3
1487                        val0 = stage[0,index]
1488                        val1 = stage[1,index]
1489
1490                        #print i, j, index, ':', L[i], val0, val1
1491                        assert allclose(float(L[i]), min(val0, val1))
1492
1493
1494        fid.close()
1495
1496        #Cleanup
1497        os.remove(prjfile)
1498        os.remove(ascfile)
1499        #os.remove(swwfile)
1500
1501
1502
1503    def test_sww2dem_asc_derived_quantity(self):
1504        """Test that sww information can be converted correctly to asc/prj
1505        format readable by e.g. ArcView
1506
1507        This tests the use of derived quantities
1508        """
1509
1510        import time, os
1511        from Numeric import array, zeros, allclose, Float, concatenate
1512        from Scientific.IO.NetCDF import NetCDFFile
1513
1514        #Setup
1515        self.domain.filename = 'datatest'
1516
1517        prjfile = self.domain.filename + '_depth.prj'
1518        ascfile = self.domain.filename + '_depth.asc'
1519        swwfile = self.domain.filename + '.sww'
1520
1521        self.domain.set_datadir('.')
1522        self.domain.format = 'sww'
1523        self.domain.smooth = True
1524        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1525        self.domain.set_quantity('stage', 0.0)
1526
1527        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1528
1529
1530        sww = get_dataobject(self.domain)
1531        sww.store_connectivity()
1532        sww.store_timestep('stage')
1533
1534        self.domain.evolve_to_end(finaltime = 0.01)
1535        sww.store_timestep('stage')
1536
1537        cellsize = 0.25
1538        #Check contents
1539        #Get NetCDF
1540
1541        fid = NetCDFFile(sww.filename, 'r')
1542
1543        # Get the variables
1544        x = fid.variables['x'][:]
1545        y = fid.variables['y'][:]
1546        z = fid.variables['elevation'][:]
1547        time = fid.variables['time'][:]
1548        stage = fid.variables['stage'][:]
1549
1550
1551        #Export to ascii/prj files
1552        sww2dem(self.domain.filename,
1553                basename_out = 'datatest_depth',
1554                quantity = 'stage - elevation',
1555                cellsize = cellsize,
1556                reduction = min,
1557                format = 'asc',
1558                verbose = False)
1559
1560
1561        #Check asc file
1562        ascid = open(ascfile)
1563        lines = ascid.readlines()
1564        ascid.close()
1565
1566        L = lines[0].strip().split()
1567        assert L[0].strip().lower() == 'ncols'
1568        assert L[1].strip().lower() == '5'
1569
1570        L = lines[1].strip().split()
1571        assert L[0].strip().lower() == 'nrows'
1572        assert L[1].strip().lower() == '5'
1573
1574        L = lines[2].strip().split()
1575        assert L[0].strip().lower() == 'xllcorner'
1576        assert allclose(float(L[1].strip().lower()), 308500)
1577
1578        L = lines[3].strip().split()
1579        assert L[0].strip().lower() == 'yllcorner'
1580        assert allclose(float(L[1].strip().lower()), 6189000)
1581
1582        L = lines[4].strip().split()
1583        assert L[0].strip().lower() == 'cellsize'
1584        assert allclose(float(L[1].strip().lower()), cellsize)
1585
1586        L = lines[5].strip().split()
1587        assert L[0].strip() == 'NODATA_value'
1588        assert L[1].strip().lower() == '-9999'
1589
1590
1591        #Check grid values (where applicable)
1592        for j in range(5):
1593            if j%2 == 0:
1594                L = lines[6+j].strip().split()
1595                jj = 4-j
1596                for i in range(5):
1597                    if i%2 == 0:
1598                        index = jj/2 + i/2*3
1599                        val0 = stage[0,index] - z[index]
1600                        val1 = stage[1,index] - z[index]
1601
1602                        #print i, j, index, ':', L[i], val0, val1
1603                        assert allclose(float(L[i]), min(val0, val1))
1604
1605
1606        fid.close()
1607
1608        #Cleanup
1609        os.remove(prjfile)
1610        os.remove(ascfile)
1611        #os.remove(swwfile)
1612
1613
1614
1615
1616
1617    def test_sww2dem_asc_missing_points(self):
1618        """Test that sww information can be converted correctly to asc/prj
1619        format readable by e.g. ArcView
1620
1621        This test includes the writing of missing values
1622        """
1623
1624        import time, os
1625        from Numeric import array, zeros, allclose, Float, concatenate
1626        from Scientific.IO.NetCDF import NetCDFFile
1627
1628        #Setup mesh not coinciding with rectangle.
1629        #This will cause missing values to occur in gridded data
1630
1631
1632        points = [                        [1.0, 1.0],
1633                              [0.5, 0.5], [1.0, 0.5],
1634                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1635
1636        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1637
1638        #Create shallow water domain
1639        domain = Domain(points, vertices)
1640        domain.default_order=2
1641
1642
1643        #Set some field values
1644        domain.set_quantity('elevation', lambda x,y: -x-y)
1645        domain.set_quantity('friction', 0.03)
1646
1647
1648        ######################
1649        # Boundary conditions
1650        B = Transmissive_boundary(domain)
1651        domain.set_boundary( {'exterior': B} )
1652
1653
1654        ######################
1655        #Initial condition - with jumps
1656
1657        bed = domain.quantities['elevation'].vertex_values
1658        stage = zeros(bed.shape, Float)
1659
1660        h = 0.3
1661        for i in range(stage.shape[0]):
1662            if i % 2 == 0:
1663                stage[i,:] = bed[i,:] + h
1664            else:
1665                stage[i,:] = bed[i,:]
1666
1667        domain.set_quantity('stage', stage)
1668        domain.distribute_to_vertices_and_edges()
1669
1670        domain.filename = 'datatest'
1671
1672        prjfile = domain.filename + '_elevation.prj'
1673        ascfile = domain.filename + '_elevation.asc'
1674        swwfile = domain.filename + '.sww'
1675
1676        domain.set_datadir('.')
1677        domain.format = 'sww'
1678        domain.smooth = True
1679
1680        domain.geo_reference = Geo_reference(56,308500,6189000)
1681
1682        sww = get_dataobject(domain)
1683        sww.store_connectivity()
1684        sww.store_timestep('stage')
1685
1686        cellsize = 0.25
1687        #Check contents
1688        #Get NetCDF
1689
1690        fid = NetCDFFile(swwfile, 'r')
1691
1692        # Get the variables
1693        x = fid.variables['x'][:]
1694        y = fid.variables['y'][:]
1695        z = fid.variables['elevation'][:]
1696        time = fid.variables['time'][:]
1697
1698        try:
1699            geo_reference = Geo_reference(NetCDFObject=fid)
1700        except AttributeError, e:
1701            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1702
1703        #Export to ascii/prj files
1704        sww2dem(domain.filename,
1705                quantity = 'elevation',
1706                cellsize = cellsize,
1707                verbose = False,
1708                format = 'asc')
1709
1710
1711        #Check asc file
1712        ascid = open(ascfile)
1713        lines = ascid.readlines()
1714        ascid.close()
1715
1716        L = lines[0].strip().split()
1717        assert L[0].strip().lower() == 'ncols'
1718        assert L[1].strip().lower() == '5'
1719
1720        L = lines[1].strip().split()
1721        assert L[0].strip().lower() == 'nrows'
1722        assert L[1].strip().lower() == '5'
1723
1724        L = lines[2].strip().split()
1725        assert L[0].strip().lower() == 'xllcorner'
1726        assert allclose(float(L[1].strip().lower()), 308500)
1727
1728        L = lines[3].strip().split()
1729        assert L[0].strip().lower() == 'yllcorner'
1730        assert allclose(float(L[1].strip().lower()), 6189000)
1731
1732        L = lines[4].strip().split()
1733        assert L[0].strip().lower() == 'cellsize'
1734        assert allclose(float(L[1].strip().lower()), cellsize)
1735
1736        L = lines[5].strip().split()
1737        assert L[0].strip() == 'NODATA_value'
1738        assert L[1].strip().lower() == '-9999'
1739
1740        #Check grid values
1741        for j in range(5):
1742            L = lines[6+j].strip().split()
1743            assert len(L) == 5
1744            y = (4-j) * cellsize
1745
1746            for i in range(5):
1747                #print i
1748                if i+j >= 4:
1749                    assert allclose(float(L[i]), -i*cellsize - y)
1750                else:
1751                    #Missing values
1752                    assert allclose(float(L[i]), -9999)
1753
1754
1755
1756        fid.close()
1757
1758        #Cleanup
1759        os.remove(prjfile)
1760        os.remove(ascfile)
1761        os.remove(swwfile)
1762
1763    def test_sww2ers_simple(self):
1764        """Test that sww information can be converted correctly to asc/prj
1765        format readable by e.g. ArcView
1766        """
1767
1768        import time, os
1769        from Numeric import array, zeros, allclose, Float, concatenate
1770        from Scientific.IO.NetCDF import NetCDFFile
1771
1772
1773        NODATA_value = 1758323
1774
1775        #Setup
1776        self.domain.filename = 'datatest'
1777
1778        headerfile = self.domain.filename + '.ers'
1779        swwfile = self.domain.filename + '.sww'
1780
1781        self.domain.set_datadir('.')
1782        self.domain.format = 'sww'
1783        self.domain.smooth = True
1784        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1785
1786        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1787
1788        sww = get_dataobject(self.domain)
1789        sww.store_connectivity()
1790        sww.store_timestep('stage')
1791
1792        self.domain.evolve_to_end(finaltime = 0.01)
1793        sww.store_timestep('stage')
1794
1795        cellsize = 0.25
1796        #Check contents
1797        #Get NetCDF
1798
1799        fid = NetCDFFile(sww.filename, 'r')
1800
1801        # Get the variables
1802        x = fid.variables['x'][:]
1803        y = fid.variables['y'][:]
1804        z = fid.variables['elevation'][:]
1805        time = fid.variables['time'][:]
1806        stage = fid.variables['stage'][:]
1807
1808
1809        #Export to ers files
1810        #sww2ers(self.domain.filename,
1811        #        quantity = 'elevation',
1812        #        cellsize = cellsize,
1813        #        verbose = False)
1814               
1815        sww2dem(self.domain.filename,
1816                quantity = 'elevation',
1817                cellsize = cellsize,
1818                NODATA_value = NODATA_value,
1819                verbose = False,
1820                format = 'ers')
1821
1822        #Check header data
1823        from ermapper_grids import read_ermapper_header, read_ermapper_data
1824       
1825        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
1826        #print header
1827        assert header['projection'].lower() == '"utm-56"'
1828        assert header['datum'].lower() == '"wgs84"'
1829        assert header['units'].lower() == '"meters"'   
1830        assert header['value'].lower() == '"elevation"'         
1831        assert header['xdimension'] == '0.25'
1832        assert header['ydimension'] == '0.25'   
1833        assert float(header['eastings']) == 308500.0   #xllcorner
1834        assert float(header['northings']) == 6189000.0 #yllcorner       
1835        assert int(header['nroflines']) == 5
1836        assert int(header['nrofcellsperline']) == 5     
1837        assert int(header['nullcellvalue']) == NODATA_value
1838        #FIXME - there is more in the header                   
1839
1840               
1841        #Check grid data               
1842        grid = read_ermapper_data(self.domain.filename + '_elevation') 
1843       
1844        #FIXME (Ole): Why is this the desired reference grid for -x-y?
1845        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
1846                    -1,    -1.25, -1.5,  -1.75, -2.0,
1847                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
1848                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1849                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
1850                                         
1851
1852        #print grid
1853        assert allclose(grid, ref_grid)
1854
1855        fid.close()
1856       
1857        #Cleanup
1858        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1859        #Done (Ole) - it was because sww2ers didn't close it's sww file
1860        os.remove(sww.filename)
1861
1862
1863
1864    def test_ferret2sww(self):
1865        """Test that georeferencing etc works when converting from
1866        ferret format (lat/lon) to sww format (UTM)
1867        """
1868        from Scientific.IO.NetCDF import NetCDFFile
1869        import os, sys
1870
1871        #The test file has
1872        # LON = 150.66667, 150.83334, 151, 151.16667
1873        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1874        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1875        #
1876        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1877        # Fourth value (index==3) is -6.50198 cm
1878
1879       
1880
1881        #Read
1882        from coordinate_transforms.redfearn import redfearn
1883        fid = NetCDFFile(self.test_MOST_file)
1884        first_value = fid.variables['HA'][:][0,0,0]
1885        fourth_value = fid.variables['HA'][:][0,0,3]
1886
1887
1888        #Call conversion (with zero origin)
1889        ferret2sww('small', verbose=False,
1890                   origin = (56, 0, 0))
1891
1892
1893        #Work out the UTM coordinates for first point
1894        zone, e, n = redfearn(-34.5, 150.66667)
1895        #print zone, e, n
1896
1897        #Read output file 'small.sww'
1898        fid = NetCDFFile('small.sww')
1899
1900        x = fid.variables['x'][:]
1901        y = fid.variables['y'][:]
1902
1903        #Check that first coordinate is correctly represented
1904        assert allclose(x[0], e)
1905        assert allclose(y[0], n)
1906
1907        #Check first value
1908        stage = fid.variables['stage'][:]
1909        xmomentum = fid.variables['xmomentum'][:]
1910        ymomentum = fid.variables['ymomentum'][:]
1911
1912        #print ymomentum
1913
1914        assert allclose(stage[0,0], first_value/100)  #Meters
1915
1916        #Check fourth value
1917        assert allclose(stage[0,3], fourth_value/100)  #Meters
1918
1919        fid.close()
1920
1921        #Cleanup
1922        import os
1923        os.remove('small.sww')
1924
1925
1926
1927    def test_ferret2sww_2(self):
1928        """Test that georeferencing etc works when converting from
1929        ferret format (lat/lon) to sww format (UTM)
1930        """
1931        from Scientific.IO.NetCDF import NetCDFFile
1932
1933        #The test file has
1934        # LON = 150.66667, 150.83334, 151, 151.16667
1935        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1936        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1937        #
1938        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1939        # Fourth value (index==3) is -6.50198 cm
1940
1941
1942        from coordinate_transforms.redfearn import redfearn
1943
1944        #fid = NetCDFFile('small_ha.nc')
1945        fid = NetCDFFile(self.test_MOST_file)
1946
1947        #Pick a coordinate and a value
1948
1949        time_index = 1
1950        lat_index = 0
1951        lon_index = 2
1952
1953        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
1954        test_time = fid.variables['TIME'][:][time_index]
1955        test_lat = fid.variables['LAT'][:][lat_index]
1956        test_lon = fid.variables['LON'][:][lon_index]
1957
1958        linear_point_index = lat_index*4 + lon_index
1959        fid.close()
1960
1961        #Call conversion (with zero origin)
1962        ferret2sww('small', verbose=False,
1963                   origin = (56, 0, 0))
1964
1965
1966        #Work out the UTM coordinates for test point
1967        zone, e, n = redfearn(test_lat, test_lon)
1968
1969        #Read output file 'small.sww'
1970        fid = NetCDFFile('small.sww')
1971
1972        x = fid.variables['x'][:]
1973        y = fid.variables['y'][:]
1974
1975        #Check that test coordinate is correctly represented
1976        assert allclose(x[linear_point_index], e)
1977        assert allclose(y[linear_point_index], n)
1978
1979        #Check test value
1980        stage = fid.variables['stage'][:]
1981
1982        assert allclose(stage[time_index, linear_point_index], test_value/100)
1983
1984        fid.close()
1985
1986        #Cleanup
1987        import os
1988        os.remove('small.sww')
1989
1990
1991
1992    def test_ferret2sww3(self):
1993        """Elevation included
1994        """
1995        from Scientific.IO.NetCDF import NetCDFFile
1996
1997        #The test file has
1998        # LON = 150.66667, 150.83334, 151, 151.16667
1999        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2000        # ELEVATION = [-1 -2 -3 -4
2001        #             -5 -6 -7 -8
2002        #              ...
2003        #              ...      -16]
2004        # where the top left corner is -1m,
2005        # and the ll corner is -13.0m
2006        #
2007        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2008        # Fourth value (index==3) is -6.50198 cm
2009
2010        from coordinate_transforms.redfearn import redfearn
2011        import os
2012        fid1 = NetCDFFile('test_ha.nc','w')
2013        fid2 = NetCDFFile('test_ua.nc','w')
2014        fid3 = NetCDFFile('test_va.nc','w')
2015        fid4 = NetCDFFile('test_e.nc','w')
2016
2017        h1_list = [150.66667,150.83334,151.]
2018        h2_list = [-34.5,-34.33333]
2019
2020        long_name = 'LON'
2021        lat_name = 'LAT'
2022
2023        nx = 3
2024        ny = 2
2025
2026        for fid in [fid1,fid2,fid3]:
2027            fid.createDimension(long_name,nx)
2028            fid.createVariable(long_name,'d',(long_name,))
2029            fid.variables[long_name].point_spacing='uneven'
2030            fid.variables[long_name].units='degrees_east'
2031            fid.variables[long_name].assignValue(h1_list)
2032
2033            fid.createDimension(lat_name,ny)
2034            fid.createVariable(lat_name,'d',(lat_name,))
2035            fid.variables[lat_name].point_spacing='uneven'
2036            fid.variables[lat_name].units='degrees_north'
2037            fid.variables[lat_name].assignValue(h2_list)
2038
2039            fid.createDimension('TIME',2)
2040            fid.createVariable('TIME','d',('TIME',))
2041            fid.variables['TIME'].point_spacing='uneven'
2042            fid.variables['TIME'].units='seconds'
2043            fid.variables['TIME'].assignValue([0.,1.])
2044            if fid == fid3: break
2045
2046
2047        for fid in [fid4]:
2048            fid.createDimension(long_name,nx)
2049            fid.createVariable(long_name,'d',(long_name,))
2050            fid.variables[long_name].point_spacing='uneven'
2051            fid.variables[long_name].units='degrees_east'
2052            fid.variables[long_name].assignValue(h1_list)
2053
2054            fid.createDimension(lat_name,ny)
2055            fid.createVariable(lat_name,'d',(lat_name,))
2056            fid.variables[lat_name].point_spacing='uneven'
2057            fid.variables[lat_name].units='degrees_north'
2058            fid.variables[lat_name].assignValue(h2_list)
2059
2060        name = {}
2061        name[fid1]='HA'
2062        name[fid2]='UA'
2063        name[fid3]='VA'
2064        name[fid4]='ELEVATION'
2065
2066        units = {}
2067        units[fid1]='cm'
2068        units[fid2]='cm/s'
2069        units[fid3]='cm/s'
2070        units[fid4]='m'
2071
2072        values = {}
2073        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2074        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2075        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2076        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2077
2078        for fid in [fid1,fid2,fid3]:
2079          fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
2080          fid.variables[name[fid]].point_spacing='uneven'
2081          fid.variables[name[fid]].units=units[fid]
2082          fid.variables[name[fid]].assignValue(values[fid])
2083          fid.variables[name[fid]].missing_value = -99999999.
2084          if fid == fid3: break
2085
2086        for fid in [fid4]:
2087            fid.createVariable(name[fid],'d',(lat_name,long_name))
2088            fid.variables[name[fid]].point_spacing='uneven'
2089            fid.variables[name[fid]].units=units[fid]
2090            fid.variables[name[fid]].assignValue(values[fid])
2091            fid.variables[name[fid]].missing_value = -99999999.
2092
2093
2094        fid1.sync(); fid1.close()
2095        fid2.sync(); fid2.close()
2096        fid3.sync(); fid3.close()
2097        fid4.sync(); fid4.close()
2098
2099        fid1 = NetCDFFile('test_ha.nc','r')
2100        fid2 = NetCDFFile('test_e.nc','r')
2101        fid3 = NetCDFFile('test_va.nc','r')
2102
2103
2104        first_amp = fid1.variables['HA'][:][0,0,0]
2105        third_amp = fid1.variables['HA'][:][0,0,2]
2106        first_elevation = fid2.variables['ELEVATION'][0,0]
2107        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2108        first_speed = fid3.variables['VA'][0,0,0]
2109        third_speed = fid3.variables['VA'][:][0,0,2]
2110
2111        fid1.close()
2112        fid2.close()
2113        fid3.close()
2114
2115        #Call conversion (with zero origin)
2116        ferret2sww('test', verbose=False,
2117                   origin = (56, 0, 0))
2118
2119        os.remove('test_va.nc')
2120        os.remove('test_ua.nc')
2121        os.remove('test_ha.nc')
2122        os.remove('test_e.nc')
2123
2124        #Read output file 'test.sww'
2125        fid = NetCDFFile('test.sww')
2126
2127
2128        #Check first value
2129        elevation = fid.variables['elevation'][:]
2130        stage = fid.variables['stage'][:]
2131        xmomentum = fid.variables['xmomentum'][:]
2132        ymomentum = fid.variables['ymomentum'][:]
2133
2134        #print ymomentum
2135        first_height = first_amp/100 - first_elevation
2136        third_height = third_amp/100 - third_elevation
2137        first_momentum=first_speed*first_height/100
2138        third_momentum=third_speed*third_height/100
2139
2140        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2141        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2142
2143        fid.close()
2144
2145        #Cleanup
2146        os.remove('test.sww')
2147
2148
2149
2150
2151    def test_ferret2sww_nz_origin(self):
2152        from Scientific.IO.NetCDF import NetCDFFile
2153        from coordinate_transforms.redfearn import redfearn
2154
2155        #Call conversion (with nonzero origin)
2156        ferret2sww('small', verbose=False,
2157                   origin = (56, 100000, 200000))
2158
2159
2160        #Work out the UTM coordinates for first point
2161        zone, e, n = redfearn(-34.5, 150.66667)
2162
2163        #Read output file 'small.sww'
2164        fid = NetCDFFile('small.sww', 'r')
2165
2166        x = fid.variables['x'][:]
2167        y = fid.variables['y'][:]
2168
2169        #Check that first coordinate is correctly represented
2170        assert allclose(x[0], e-100000)
2171        assert allclose(y[0], n-200000)
2172
2173        fid.close()
2174
2175        #Cleanup
2176        os.remove('small.sww')
2177
2178
2179
2180    def test_sww_extent(self):
2181        """Not a test, rather a look at the sww format
2182        """
2183
2184        import time, os
2185        from Numeric import array, zeros, allclose, Float, concatenate
2186        from Scientific.IO.NetCDF import NetCDFFile
2187
2188        self.domain.filename = 'datatest' + str(id(self))
2189        self.domain.format = 'sww'
2190        self.domain.smooth = True
2191        self.domain.reduction = mean
2192        self.domain.set_datadir('.')
2193
2194
2195        sww = get_dataobject(self.domain)
2196        sww.store_connectivity()
2197        sww.store_timestep('stage')
2198        self.domain.time = 2.
2199
2200        #Modify stage at second timestep
2201        stage = self.domain.quantities['stage'].vertex_values
2202        self.domain.set_quantity('stage', stage/2)
2203
2204        sww.store_timestep('stage')
2205
2206        file_and_extension_name = self.domain.filename + ".sww"
2207        #print "file_and_extension_name",file_and_extension_name
2208        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2209               extent_sww(file_and_extension_name )
2210
2211        assert allclose(xmin, 0.0)
2212        assert allclose(xmax, 1.0)
2213        assert allclose(ymin, 0.0)
2214        assert allclose(ymax, 1.0)
2215        assert allclose(stagemin, -0.85)
2216        assert allclose(stagemax, 0.15)
2217
2218
2219        #Cleanup
2220        os.remove(sww.filename)
2221
2222
2223
2224    def test_sww2domain(self):
2225        ################################################
2226        #Create a test domain, and evolve and save it.
2227        ################################################
2228        from mesh_factory import rectangular
2229        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2230             Constant_height, Time_boundary, Transmissive_boundary
2231        from Numeric import array
2232
2233        #Create basic mesh
2234
2235        yiel=0.01
2236        points, vertices, boundary = rectangular(10,10)
2237
2238        #Create shallow water domain
2239        domain = Domain(points, vertices, boundary)
2240        domain.geo_reference = Geo_reference(56,11,11)
2241        domain.smooth = False
2242        domain.visualise = False
2243        domain.store = True
2244        domain.filename = 'bedslope'
2245        domain.default_order=2
2246        #Bed-slope and friction
2247        domain.set_quantity('elevation', lambda x,y: -x/3)
2248        domain.set_quantity('friction', 0.1)
2249        # Boundary conditions
2250        from math import sin, pi
2251        Br = Reflective_boundary(domain)
2252        Bt = Transmissive_boundary(domain)
2253        Bd = Dirichlet_boundary([0.2,0.,0.])
2254        Bw = Time_boundary(domain=domain,
2255                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2256
2257        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2258        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2259
2260        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2261        #Initial condition
2262        h = 0.05
2263        elevation = domain.quantities['elevation'].vertex_values
2264        domain.set_quantity('stage', elevation + h)
2265        #elevation = domain.get_quantity('elevation')
2266        #domain.set_quantity('stage', elevation + h)
2267
2268        domain.check_integrity()
2269        #Evolution
2270        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2271        #    domain.write_time()
2272            pass
2273
2274
2275        ##########################################
2276        #Import the example's file as a new domain
2277        ##########################################
2278        from data_manager import sww2domain
2279        from Numeric import allclose
2280        import os
2281
2282        filename = domain.datadir+os.sep+domain.filename+'.sww'
2283        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2284        #points, vertices, boundary = rectangular(15,15)
2285        #domain2.boundary = boundary
2286        ###################
2287        ##NOW TEST IT!!!
2288        ###################
2289
2290        bits = ['vertex_coordinates']
2291        for quantity in ['elevation']+domain.quantities_to_be_stored:
2292            bits.append('quantities["%s"].get_integral()'%quantity)
2293            bits.append('get_quantity("%s")'%quantity)
2294
2295        for bit in bits:
2296            #print 'testing that domain.'+bit+' has been restored'
2297            #print bit
2298        #print 'done'
2299            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2300
2301        ######################################
2302        #Now evolve them both, just to be sure
2303        ######################################x
2304        visualise = False
2305        #visualise = True
2306        domain.visualise = visualise
2307        domain.time = 0.
2308        from time import sleep
2309
2310        final = .1
2311        domain.set_quantity('friction', 0.1)
2312        domain.store = False
2313        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2314
2315
2316        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2317            if visualise: sleep(1.)
2318            #domain.write_time()
2319            pass
2320
2321        final = final - (domain2.starttime-domain.starttime)
2322        #BUT since domain1 gets time hacked back to 0:
2323        final = final + (domain2.starttime-domain.starttime)
2324
2325        domain2.smooth = False
2326        domain2.visualise = visualise
2327        domain2.store = False
2328        domain2.default_order=2
2329        domain2.set_quantity('friction', 0.1)
2330        #Bed-slope and friction
2331        # Boundary conditions
2332        Bd2=Dirichlet_boundary([0.2,0.,0.])
2333        domain2.boundary = domain.boundary
2334        #print 'domain2.boundary'
2335        #print domain2.boundary
2336        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2337        #domain2.set_boundary({'exterior': Bd})
2338
2339        domain2.check_integrity()
2340
2341        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2342            if visualise: sleep(1.)
2343            #domain2.write_time()
2344            pass
2345
2346        ###################
2347        ##NOW TEST IT!!!
2348        ##################
2349
2350        bits = [ 'vertex_coordinates']
2351
2352        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2353            bits.append('quantities["%s"].get_integral()'%quantity)
2354            bits.append('get_quantity("%s")'%quantity)
2355
2356        for bit in bits:
2357            #print bit
2358            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2359
2360
2361    def test_sww2domain2(self):
2362        ##################################################################
2363        #Same as previous test, but this checks how NaNs are handled.
2364        ##################################################################
2365
2366
2367        from mesh_factory import rectangular
2368        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2369             Constant_height, Time_boundary, Transmissive_boundary
2370        from Numeric import array
2371
2372        #Create basic mesh
2373        points, vertices, boundary = rectangular(2,2)
2374
2375        #Create shallow water domain
2376        domain = Domain(points, vertices, boundary)
2377        domain.smooth = False
2378        domain.visualise = False
2379        domain.store = True
2380        domain.filename = 'bedslope'
2381        domain.default_order=2
2382        domain.quantities_to_be_stored=['stage']
2383
2384        domain.set_quantity('elevation', lambda x,y: -x/3)
2385        domain.set_quantity('friction', 0.1)
2386
2387        from math import sin, pi
2388        Br = Reflective_boundary(domain)
2389        Bt = Transmissive_boundary(domain)
2390        Bd = Dirichlet_boundary([0.2,0.,0.])
2391        Bw = Time_boundary(domain=domain,
2392                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2393
2394        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2395
2396        h = 0.05
2397        elevation = domain.quantities['elevation'].vertex_values
2398        domain.set_quantity('stage', elevation + h)
2399
2400        domain.check_integrity()
2401
2402        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2403            pass
2404            #domain.write_time()
2405
2406
2407
2408        ##################################
2409        #Import the file as a new domain
2410        ##################################
2411        from data_manager import sww2domain
2412        from Numeric import allclose
2413        import os
2414
2415        filename = domain.datadir+os.sep+domain.filename+'.sww'
2416
2417        #Fail because NaNs are present
2418        try:
2419            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2420            assert True == False
2421        except:
2422            #Now import it, filling NaNs to be 0
2423            filler = 0
2424            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2425        bits = [ 'geo_reference.get_xllcorner()',
2426                'geo_reference.get_yllcorner()',
2427                'vertex_coordinates']
2428
2429        for quantity in ['elevation']+domain.quantities_to_be_stored:
2430            bits.append('quantities["%s"].get_integral()'%quantity)
2431            bits.append('get_quantity("%s")'%quantity)
2432
2433        for bit in bits:
2434        #    print 'testing that domain.'+bit+' has been restored'
2435            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2436
2437        assert max(max(domain2.get_quantity('xmomentum')))==filler
2438        assert min(min(domain2.get_quantity('xmomentum')))==filler
2439        assert max(max(domain2.get_quantity('ymomentum')))==filler
2440        assert min(min(domain2.get_quantity('ymomentum')))==filler
2441
2442        #print 'passed'
2443
2444        #cleanup
2445        #import os
2446        #os.remove(domain.datadir+'/'+domain.filename+'.sww')
2447
2448
2449    #def test_weed(self):
2450        from data_manager import weed
2451
2452        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2453        volumes1 = [[0,1,2],[3,4,5]]
2454        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2455        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2456
2457        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2458
2459        assert len(points2)==len(coordinates2)
2460        for i in range(len(coordinates2)):
2461            coordinate = tuple(coordinates2[i])
2462            assert points2.has_key(coordinate)
2463            points2[coordinate]=i
2464
2465        for triangle in volumes1:
2466            for coordinate in triangle:
2467                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2468                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2469
2470
2471    #FIXME This fails - smooth makes the comparism too hard for allclose
2472    def ztest_sww2domain3(self):
2473        ################################################
2474        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2475        ################################################
2476        from mesh_factory import rectangular
2477        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2478             Constant_height, Time_boundary, Transmissive_boundary
2479        from Numeric import array
2480        #Create basic mesh
2481
2482        yiel=0.01
2483        points, vertices, boundary = rectangular(10,10)
2484
2485        #Create shallow water domain
2486        domain = Domain(points, vertices, boundary)
2487        domain.geo_reference = Geo_reference(56,11,11)
2488        domain.smooth = True
2489        domain.visualise = False
2490        domain.store = True
2491        domain.filename = 'bedslope'
2492        domain.default_order=2
2493        #Bed-slope and friction
2494        domain.set_quantity('elevation', lambda x,y: -x/3)
2495        domain.set_quantity('friction', 0.1)
2496        # Boundary conditions
2497        from math import sin, pi
2498        Br = Reflective_boundary(domain)
2499        Bt = Transmissive_boundary(domain)
2500        Bd = Dirichlet_boundary([0.2,0.,0.])
2501        Bw = Time_boundary(domain=domain,
2502                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2503
2504        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2505
2506        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2507        #Initial condition
2508        h = 0.05
2509        elevation = domain.quantities['elevation'].vertex_values
2510        domain.set_quantity('stage', elevation + h)
2511        #elevation = domain.get_quantity('elevation')
2512        #domain.set_quantity('stage', elevation + h)
2513
2514        domain.check_integrity()
2515        #Evolution
2516        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2517        #    domain.write_time()
2518            pass
2519
2520
2521        ##########################################
2522        #Import the example's file as a new domain
2523        ##########################################
2524        from data_manager import sww2domain
2525        from Numeric import allclose
2526        import os
2527
2528        filename = domain.datadir+os.sep+domain.filename+'.sww'
2529        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2530        #points, vertices, boundary = rectangular(15,15)
2531        #domain2.boundary = boundary
2532        ###################
2533        ##NOW TEST IT!!!
2534        ###################
2535
2536        #FIXME smooth domain so that they can be compared
2537
2538
2539        bits = []#'vertex_coordinates']
2540        for quantity in ['elevation']+domain.quantities_to_be_stored:
2541            bits.append('quantities["%s"].get_integral()'%quantity)
2542            #bits.append('get_quantity("%s")'%quantity)
2543
2544        for bit in bits:
2545            #print 'testing that domain.'+bit+' has been restored'
2546            #print bit
2547            #print 'done'
2548            #print ('domain.'+bit), eval('domain.'+bit)
2549            #print ('domain2.'+bit), eval('domain2.'+bit)
2550            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2551            pass
2552
2553        ######################################
2554        #Now evolve them both, just to be sure
2555        ######################################x
2556        visualise = False
2557        visualise = True
2558        domain.visualise = visualise
2559        domain.time = 0.
2560        from time import sleep
2561
2562        final = .5
2563        domain.set_quantity('friction', 0.1)
2564        domain.store = False
2565        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2566
2567        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2568            if visualise: sleep(.03)
2569            #domain.write_time()
2570            pass
2571
2572        domain2.smooth = True
2573        domain2.visualise = visualise
2574        domain2.store = False
2575        domain2.default_order=2
2576        domain2.set_quantity('friction', 0.1)
2577        #Bed-slope and friction
2578        # Boundary conditions
2579        Bd2=Dirichlet_boundary([0.2,0.,0.])
2580        Br2 = Reflective_boundary(domain2)
2581        domain2.boundary = domain.boundary
2582        #print 'domain2.boundary'
2583        #print domain2.boundary
2584        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2585        #domain2.boundary = domain.boundary
2586        #domain2.set_boundary({'exterior': Bd})
2587
2588        domain2.check_integrity()
2589
2590        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2591            if visualise: sleep(.03)
2592            #domain2.write_time()
2593            pass
2594
2595        ###################
2596        ##NOW TEST IT!!!
2597        ##################
2598
2599        bits = [ 'vertex_coordinates']
2600
2601        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2602            #bits.append('quantities["%s"].get_integral()'%quantity)
2603            bits.append('get_quantity("%s")'%quantity)
2604
2605        for bit in bits:
2606            print bit
2607            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2608
2609
2610    def test_decimate_dem(self):
2611        """Test decimation of dem file
2612        """
2613
2614        import os
2615        from Numeric import ones, allclose, Float, arange
2616        from Scientific.IO.NetCDF import NetCDFFile
2617
2618        #Write test dem file
2619        root = 'decdemtest'
2620
2621        filename = root + '.dem'
2622        fid = NetCDFFile(filename, 'w')
2623
2624        fid.institution = 'Geoscience Australia'
2625        fid.description = 'NetCDF DEM format for compact and portable ' +\
2626                          'storage of spatial point data'
2627
2628        nrows = 15
2629        ncols = 18
2630
2631        fid.ncols = ncols
2632        fid.nrows = nrows
2633        fid.xllcorner = 2000.5
2634        fid.yllcorner = 3000.5
2635        fid.cellsize = 25
2636        fid.NODATA_value = -9999
2637
2638        fid.zone = 56
2639        fid.false_easting = 0.0
2640        fid.false_northing = 0.0
2641        fid.projection = 'UTM'
2642        fid.datum = 'WGS84'
2643        fid.units = 'METERS'
2644
2645        fid.createDimension('number_of_points', nrows*ncols)
2646
2647        fid.createVariable('elevation', Float, ('number_of_points',))
2648
2649        elevation = fid.variables['elevation']
2650
2651        elevation[:] = (arange(nrows*ncols))
2652
2653        fid.close()
2654
2655        #generate the elevation values expected in the decimated file
2656        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2657                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2658                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2659                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2660                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2661                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2662                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2663                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2664                         (144+145+146+162+163+164+180+181+182) / 9.0,
2665                         (148+149+150+166+167+168+184+185+186) / 9.0,
2666                         (152+153+154+170+171+172+188+189+190) / 9.0,
2667                         (156+157+158+174+175+176+192+193+194) / 9.0,
2668                         (216+217+218+234+235+236+252+253+254) / 9.0,
2669                         (220+221+222+238+239+240+256+257+258) / 9.0,
2670                         (224+225+226+242+243+244+260+261+262) / 9.0,
2671                         (228+229+230+246+247+248+264+265+266) / 9.0]
2672
2673        #generate a stencil for computing the decimated values
2674        stencil = ones((3,3), Float) / 9.0
2675
2676        decimate_dem(root, stencil=stencil, cellsize_new=100)
2677
2678        #Open decimated NetCDF file
2679        fid = NetCDFFile(root + '_100.dem', 'r')
2680
2681        # Get decimated elevation
2682        elevation = fid.variables['elevation']
2683
2684        #Check values
2685        assert allclose(elevation, ref_elevation)
2686
2687        #Cleanup
2688        fid.close()
2689
2690        os.remove(root + '.dem')
2691        os.remove(root + '_100.dem')
2692
2693    def test_decimate_dem_NODATA(self):
2694        """Test decimation of dem file that includes NODATA values
2695        """
2696
2697        import os
2698        from Numeric import ones, allclose, Float, arange, reshape
2699        from Scientific.IO.NetCDF import NetCDFFile
2700
2701        #Write test dem file
2702        root = 'decdemtest'
2703
2704        filename = root + '.dem'
2705        fid = NetCDFFile(filename, 'w')
2706
2707        fid.institution = 'Geoscience Australia'
2708        fid.description = 'NetCDF DEM format for compact and portable ' +\
2709                          'storage of spatial point data'
2710
2711        nrows = 15
2712        ncols = 18
2713        NODATA_value = -9999
2714
2715        fid.ncols = ncols
2716        fid.nrows = nrows
2717        fid.xllcorner = 2000.5
2718        fid.yllcorner = 3000.5
2719        fid.cellsize = 25
2720        fid.NODATA_value = NODATA_value
2721
2722        fid.zone = 56
2723        fid.false_easting = 0.0
2724        fid.false_northing = 0.0
2725        fid.projection = 'UTM'
2726        fid.datum = 'WGS84'
2727        fid.units = 'METERS'
2728
2729        fid.createDimension('number_of_points', nrows*ncols)
2730
2731        fid.createVariable('elevation', Float, ('number_of_points',))
2732
2733        elevation = fid.variables['elevation']
2734
2735        #generate initial elevation values
2736        elevation_tmp = (arange(nrows*ncols))
2737        #add some NODATA values
2738        elevation_tmp[0]   = NODATA_value
2739        elevation_tmp[95]  = NODATA_value
2740        elevation_tmp[188] = NODATA_value
2741        elevation_tmp[189] = NODATA_value
2742        elevation_tmp[190] = NODATA_value
2743        elevation_tmp[209] = NODATA_value
2744        elevation_tmp[252] = NODATA_value
2745
2746        elevation[:] = elevation_tmp
2747
2748        fid.close()
2749
2750        #generate the elevation values expected in the decimated file
2751        ref_elevation = [NODATA_value,
2752                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2753                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2754                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2755                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2756                         NODATA_value,
2757                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2758                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2759                         (144+145+146+162+163+164+180+181+182) / 9.0,
2760                         (148+149+150+166+167+168+184+185+186) / 9.0,
2761                         NODATA_value,
2762                         (156+157+158+174+175+176+192+193+194) / 9.0,
2763                         NODATA_value,
2764                         (220+221+222+238+239+240+256+257+258) / 9.0,
2765                         (224+225+226+242+243+244+260+261+262) / 9.0,
2766                         (228+229+230+246+247+248+264+265+266) / 9.0]
2767
2768        #generate a stencil for computing the decimated values
2769        stencil = ones((3,3), Float) / 9.0
2770
2771        decimate_dem(root, stencil=stencil, cellsize_new=100)
2772
2773        #Open decimated NetCDF file
2774        fid = NetCDFFile(root + '_100.dem', 'r')
2775
2776        # Get decimated elevation
2777        elevation = fid.variables['elevation']
2778
2779        #Check values
2780        assert allclose(elevation, ref_elevation)
2781
2782        #Cleanup
2783        fid.close()
2784
2785        os.remove(root + '.dem')
2786        os.remove(root + '_100.dem')
2787
2788    def xxxtestz_sww2ers_real(self):
2789        """Test that sww information can be converted correctly to asc/prj
2790        format readable by e.g. ArcView
2791        """
2792
2793        import time, os
2794        from Numeric import array, zeros, allclose, Float, concatenate
2795        from Scientific.IO.NetCDF import NetCDFFile
2796
2797        # the memory optimised least squares
2798        #  cellsize = 20,   # this one seems to hang
2799        #  cellsize = 200000, # Ran 1 test in 269.703s
2800                                #Ran 1 test in 267.344s
2801        #  cellsize = 20000,  # Ran 1 test in 460.922s
2802        #  cellsize = 2000   #Ran 1 test in 5340.250s
2803        #  cellsize = 200   #this one seems to hang, building matirx A
2804
2805        # not optimised
2806        # seems to hang
2807        #  cellsize = 2000   # Ran 1 test in 5334.563s
2808        #Export to ascii/prj files
2809        sww2dem('karratha_100m',
2810                quantity = 'depth',
2811                cellsize = 200000,
2812                verbose = True)
2813
2814    def test_read_asc(self):
2815        """Test conversion from dem in ascii format to native NetCDF xya format
2816        """
2817
2818        import time, os
2819        from Numeric import array, zeros, allclose, Float, concatenate
2820        from Scientific.IO.NetCDF import NetCDFFile
2821
2822        import data_manager
2823        #Write test asc file
2824        filename = tempfile.mktemp(".000")
2825        fid = open(filename, 'w')
2826        fid.write("""ncols         7
2827nrows         4
2828xllcorner     2000.5
2829yllcorner     3000.5
2830cellsize      25
2831NODATA_value  -9999
2832    97.921    99.285   125.588   180.830   258.645   342.872   415.836
2833   473.157   514.391   553.893   607.120   678.125   777.283   883.038
2834   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
2835   502.645   516.230   504.739   450.604   388.500   338.097   514.980
2836""")
2837        fid.close()
2838        bath_metadata, grid = \
2839                   data_manager._read_asc(filename, verbose=False)
2840        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
2841        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
2842        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
2843        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
2844        self.failUnless(grid[0][0]  == 97.921,  'Failed')
2845        self.failUnless(grid[3][6]  == 514.980,  'Failed')
2846
2847        os.remove(filename)
2848       
2849    def test_asc_csiro2sww(self):
2850        import tempfile
2851
2852        bath_dir = tempfile.mkdtemp()
2853        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
2854        #bath_dir = 'bath_data_manager_test'
2855        #print "os.getcwd( )",os.getcwd( )
2856        elevation_dir =  tempfile.mkdtemp()
2857        #elevation_dir = 'elev_expanded'
2858        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
2859        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
2860       
2861        fid = open(bath_dir_filename, 'w')
2862        fid.write(""" ncols             3
2863 nrows             2
2864 xllcorner    148.00000
2865 yllcorner    -38.00000
2866 cellsize       0.25
2867 nodata_value   -9999.0
2868    9000.000    -1000.000    3000.0
2869   -1000.000    9000.000  -1000.000
2870""")
2871        fid.close()
2872       
2873        fid = open(elevation_dir_filename1, 'w')
2874        fid.write(""" ncols             3
2875 nrows             2
2876 xllcorner    148.00000
2877 yllcorner    -38.00000
2878 cellsize       0.25
2879 nodata_value   -9999.0
2880    9000.000    0.000    3000.0
2881     0.000     9000.000     0.000
2882""")
2883        fid.close()
2884
2885        fid = open(elevation_dir_filename2, 'w')
2886        fid.write(""" ncols             3
2887 nrows             2
2888 xllcorner    148.00000
2889 yllcorner    -38.00000
2890 cellsize       0.25
2891 nodata_value   -9999.0
2892    9000.000    4000.000    4000.0
2893    4000.000    9000.000    4000.000
2894""")
2895        fid.close()
2896         
2897        ucur_dir =  tempfile.mkdtemp()
2898        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
2899        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
2900       
2901        fid = open(ucur_dir_filename1, 'w')
2902        fid.write(""" ncols             3
2903 nrows             2
2904 xllcorner    148.00000
2905 yllcorner    -38.00000
2906 cellsize       0.25
2907 nodata_value   -9999.0
2908    90.000    60.000    30.0
2909    10.000    10.000    10.000
2910""")
2911        fid.close()
2912        fid = open(ucur_dir_filename2, 'w')
2913        fid.write(""" ncols             3
2914 nrows             2
2915 xllcorner    148.00000
2916 yllcorner    -38.00000
2917 cellsize       0.25
2918 nodata_value   -9999.0
2919    90.000    60.000    30.0
2920    10.000    10.000    10.000
2921""")
2922        fid.close()
2923       
2924        vcur_dir =  tempfile.mkdtemp()
2925        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
2926        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
2927       
2928        fid = open(vcur_dir_filename1, 'w')
2929        fid.write(""" ncols             3
2930 nrows             2
2931 xllcorner    148.00000
2932 yllcorner    -38.00000
2933 cellsize       0.25
2934 nodata_value   -9999.0
2935    90.000    60.000    30.0
2936    10.000    10.000    10.000
2937""")
2938        fid.close()
2939        fid = open(vcur_dir_filename2, 'w')
2940        fid.write(""" ncols             3
2941 nrows             2
2942 xllcorner    148.00000
2943 yllcorner    -38.00000
2944 cellsize       0.25
2945 nodata_value   -9999.0
2946    90.000    60.000    30.0
2947    10.000    10.000    10.000
2948""")
2949        fid.close()
2950       
2951        sww_file = 'a_test.sww'
2952        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
2953
2954        # check the sww file
2955       
2956        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
2957        x = fid.variables['x'][:]
2958        y = fid.variables['y'][:]
2959        z = fid.variables['z'][:]
2960        stage = fid.variables['stage'][:]
2961        xmomentum = fid.variables['xmomentum'][:]
2962        geo_ref = Geo_reference(NetCDFObject=fid)
2963        #print "geo_ref",geo_ref
2964        x_ref = geo_ref.get_xllcorner()
2965        y_ref = geo_ref.get_yllcorner()
2966        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
2967        assert allclose(x_ref, 587798.418) # (-38, 148)
2968        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
2969       
2970        #Zone:   55   
2971        #Easting:  588095.674  Northing: 5821451.722
2972        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
2973        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
2974
2975        #Zone:   55   
2976        #Easting:  632145.632  Northing: 5820863.269
2977        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
2978        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
2979
2980        #Zone:   55   
2981        #Easting:  609748.788  Northing: 5793447.860
2982        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
2983        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
2984
2985        assert allclose(z[0],9000.0 )
2986        assert allclose(stage[0][1],0.0 )
2987
2988        #(4000+1000)*60
2989        assert allclose(xmomentum[1][1],300000.0 )
2990       
2991       
2992        fid.close()
2993   
2994        #tidy up
2995        os.remove(bath_dir_filename)
2996        os.rmdir(bath_dir)
2997
2998        os.remove(elevation_dir_filename1)
2999        os.remove(elevation_dir_filename2)
3000        os.rmdir(elevation_dir)
3001       
3002        os.remove(ucur_dir_filename1)
3003        os.remove(ucur_dir_filename2)
3004        os.rmdir(ucur_dir)
3005       
3006        os.remove(vcur_dir_filename1)
3007        os.remove(vcur_dir_filename2)
3008        os.rmdir(vcur_dir)
3009
3010
3011        # remove sww file
3012        os.remove(sww_file)
3013       
3014    def test_asc_csiro2sww2(self):
3015        import tempfile
3016
3017        bath_dir = tempfile.mkdtemp()
3018        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3019        #bath_dir = 'bath_data_manager_test'
3020        #print "os.getcwd( )",os.getcwd( )
3021        elevation_dir =  tempfile.mkdtemp()
3022        #elevation_dir = 'elev_expanded'
3023        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3024        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3025       
3026        fid = open(bath_dir_filename, 'w')
3027        fid.write(""" ncols             3
3028 nrows             2
3029 xllcorner    148.00000
3030 yllcorner    -38.00000
3031 cellsize       0.25
3032 nodata_value   -9999.0
3033    9000.000    -1000.000    3000.0
3034   -1000.000    9000.000  -1000.000
3035""")
3036        fid.close()
3037       
3038        fid = open(elevation_dir_filename1, 'w')
3039        fid.write(""" ncols             3
3040 nrows             2
3041 xllcorner    148.00000
3042 yllcorner    -38.00000
3043 cellsize       0.25
3044 nodata_value   -9999.0
3045    9000.000    0.000    3000.0
3046     0.000     -9999.000     -9999.000
3047""")
3048        fid.close()
3049
3050        fid = open(elevation_dir_filename2, 'w')
3051        fid.write(""" ncols             3
3052 nrows             2
3053 xllcorner    148.00000
3054 yllcorner    -38.00000
3055 cellsize       0.25
3056 nodata_value   -9999.0
3057    9000.000    4000.000    4000.0
3058    4000.000    9000.000    4000.000
3059""")
3060        fid.close()
3061         
3062        ucur_dir =  tempfile.mkdtemp()
3063        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3064        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3065       
3066        fid = open(ucur_dir_filename1, 'w')
3067        fid.write(""" ncols             3
3068 nrows             2
3069 xllcorner    148.00000
3070 yllcorner    -38.00000
3071 cellsize       0.25
3072 nodata_value   -9999.0
3073    90.000    60.000    30.0
3074    10.000    10.000    10.000
3075""")
3076        fid.close()
3077        fid = open(ucur_dir_filename2, 'w')
3078        fid.write(""" ncols             3
3079 nrows             2
3080 xllcorner    148.00000
3081 yllcorner    -38.00000
3082 cellsize       0.25
3083 nodata_value   -9999.0
3084    90.000    60.000    30.0
3085    10.000    10.000    10.000
3086""")
3087        fid.close()
3088       
3089        vcur_dir =  tempfile.mkdtemp()
3090        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3091        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3092       
3093        fid = open(vcur_dir_filename1, 'w')
3094        fid.write(""" ncols             3
3095 nrows             2
3096 xllcorner    148.00000
3097 yllcorner    -38.00000
3098 cellsize       0.25
3099 nodata_value   -9999.0
3100    90.000    60.000    30.0
3101    10.000    10.000    10.000
3102""")
3103        fid.close()
3104        fid = open(vcur_dir_filename2, 'w')
3105        fid.write(""" ncols             3
3106 nrows             2
3107 xllcorner    148.00000
3108 yllcorner    -38.00000
3109 cellsize       0.25
3110 nodata_value   -9999.0
3111    90.000    60.000    30.0
3112    10.000    10.000    10.000
3113""")
3114        fid.close()
3115       
3116        try:
3117            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3118                          vcur_dir, sww_file)
3119        except:
3120            #tidy up
3121            os.remove(bath_dir_filename)
3122            os.rmdir(bath_dir)
3123
3124            os.remove(elevation_dir_filename1)
3125            os.remove(elevation_dir_filename2)
3126            os.rmdir(elevation_dir)
3127       
3128            os.remove(ucur_dir_filename1)
3129            os.remove(ucur_dir_filename2)
3130            os.rmdir(ucur_dir)
3131       
3132            os.remove(vcur_dir_filename1)
3133            os.remove(vcur_dir_filename2)
3134            os.rmdir(vcur_dir)
3135        else:
3136            #tidy up
3137            os.remove(bath_dir_filename)
3138            os.rmdir(bath_dir)
3139
3140            os.remove(elevation_dir_filename1)
3141            os.remove(elevation_dir_filename2)
3142            os.rmdir(elevation_dir)
3143            raise 'Should raise exception'
3144       
3145            os.remove(ucur_dir_filename1)
3146            os.remove(ucur_dir_filename2)
3147            os.rmdir(ucur_dir)
3148       
3149            os.remove(vcur_dir_filename1)
3150            os.remove(vcur_dir_filename2)
3151            os.rmdir(vcur_dir)
3152
3153       
3154     
3155    def test_asc_csiro2sww3(self):
3156        import tempfile
3157
3158        bath_dir = tempfile.mkdtemp()
3159        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3160        #bath_dir = 'bath_data_manager_test'
3161        #print "os.getcwd( )",os.getcwd( )
3162        elevation_dir =  tempfile.mkdtemp()
3163        #elevation_dir = 'elev_expanded'
3164        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3165        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3166       
3167        fid = open(bath_dir_filename, 'w')
3168        fid.write(""" ncols             3
3169 nrows             2
3170 xllcorner    148.00000
3171 yllcorner    -38.00000
3172 cellsize       0.25
3173 nodata_value   -9999.0
3174    9000.000    -1000.000    3000.0
3175   -1000.000    9000.000  -1000.000
3176""")
3177        fid.close()
3178       
3179        fid = open(elevation_dir_filename1, 'w')
3180        fid.write(""" ncols             3
3181 nrows             2
3182 xllcorner    148.00000
3183 yllcorner    -38.00000
3184 cellsize       0.25
3185 nodata_value   -9999.0
3186    9000.000    0.000    3000.0
3187     0.000     -9999.000     -9999.000
3188""")
3189        fid.close()
3190
3191        fid = open(elevation_dir_filename2, 'w')
3192        fid.write(""" ncols             3
3193 nrows             2
3194 xllcorner    148.00000
3195 yllcorner    -38.00000
3196 cellsize       0.25
3197 nodata_value   -9999.0
3198    9000.000    4000.000    4000.0
3199    4000.000    9000.000    4000.000
3200""")
3201        fid.close()
3202         
3203        ucur_dir =  tempfile.mkdtemp()
3204        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3205        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3206       
3207        fid = open(ucur_dir_filename1, 'w')
3208        fid.write(""" ncols             3
3209 nrows             2
3210 xllcorner    148.00000
3211 yllcorner    -38.00000
3212 cellsize       0.25
3213 nodata_value   -9999.0
3214    90.000    60.000    30.0
3215    10.000    10.000    10.000
3216""")
3217        fid.close()
3218        fid = open(ucur_dir_filename2, 'w')
3219        fid.write(""" ncols             3
3220 nrows             2
3221 xllcorner    148.00000
3222 yllcorner    -38.00000
3223 cellsize       0.25
3224 nodata_value   -9999.0
3225    90.000    60.000    30.0
3226    10.000    10.000    10.000
3227""")
3228        fid.close()
3229       
3230        vcur_dir =  tempfile.mkdtemp()
3231        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3232        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3233       
3234        fid = open(vcur_dir_filename1, 'w')
3235        fid.write(""" ncols             3
3236 nrows             2
3237 xllcorner    148.00000
3238 yllcorner    -38.00000
3239 cellsize       0.25
3240 nodata_value   -9999.0
3241    90.000    60.000    30.0
3242    10.000    10.000    10.000
3243""")
3244        fid.close()
3245        fid = open(vcur_dir_filename2, 'w')
3246        fid.write(""" ncols             3
3247 nrows             2
3248 xllcorner    148.00000
3249 yllcorner    -38.00000
3250 cellsize       0.25
3251 nodata_value   -9999.0
3252    90.000    60.000    30.0
3253    10.000    10.000    10.000
3254""")
3255        fid.close()
3256       
3257        sww_file = 'a_test.sww'
3258        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3259                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3260                      mean_stage = 100)
3261
3262        # check the sww file
3263       
3264        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3265        x = fid.variables['x'][:]
3266        y = fid.variables['y'][:]
3267        z = fid.variables['z'][:]
3268        stage = fid.variables['stage'][:]
3269        xmomentum = fid.variables['xmomentum'][:]
3270        geo_ref = Geo_reference(NetCDFObject=fid)
3271        #print "geo_ref",geo_ref
3272        x_ref = geo_ref.get_xllcorner()
3273        y_ref = geo_ref.get_yllcorner()
3274        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3275        assert allclose(x_ref, 587798.418) # (-38, 148)
3276        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3277       
3278        #Zone:   55   
3279        #Easting:  588095.674  Northing: 5821451.722
3280        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3281        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3282
3283        #Zone:   55   
3284        #Easting:  632145.632  Northing: 5820863.269
3285        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3286        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3287
3288        #Zone:   55   
3289        #Easting:  609748.788  Northing: 5793447.860
3290        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3291        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3292
3293        assert allclose(z[0],9000.0 )
3294        assert allclose(stage[0][4],100.0 )
3295        assert allclose(stage[0][5],100.0 )
3296
3297        #(100.0 - 9000)*10
3298        assert allclose(xmomentum[0][4], -89000.0 )
3299       
3300        #(100.0 - -1000.000)*10
3301        assert allclose(xmomentum[0][5], 11000.0 )
3302       
3303        fid.close()
3304   
3305        #tidy up
3306        os.remove(bath_dir_filename)
3307        os.rmdir(bath_dir)
3308
3309        os.remove(elevation_dir_filename1)
3310        os.remove(elevation_dir_filename2)
3311        os.rmdir(elevation_dir)
3312       
3313        os.remove(ucur_dir_filename1)
3314        os.remove(ucur_dir_filename2)
3315        os.rmdir(ucur_dir)
3316       
3317        os.remove(vcur_dir_filename1)
3318        os.remove(vcur_dir_filename2)
3319        os.rmdir(vcur_dir)
3320
3321        # remove sww file
3322        os.remove(sww_file)
3323   
3324     
3325    def test_asc_csiro2sww4(self):
3326        """
3327        Test specifying the extent
3328        """
3329       
3330        import tempfile
3331
3332        bath_dir = tempfile.mkdtemp()
3333        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3334        #bath_dir = 'bath_data_manager_test'
3335        #print "os.getcwd( )",os.getcwd( )
3336        elevation_dir =  tempfile.mkdtemp()
3337        #elevation_dir = 'elev_expanded'
3338        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3339        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3340       
3341        fid = open(bath_dir_filename, 'w')
3342        fid.write(""" ncols             4
3343 nrows             4
3344 xllcorner    148.00000
3345 yllcorner    -38.00000
3346 cellsize       0.25
3347 nodata_value   -9999.0
3348   -9000.000    -1000.000   -3000.0 -2000.000
3349   -1000.000    9000.000  -1000.000 -3000.000
3350   -4000.000    6000.000   2000.000 -5000.000
3351   -9000.000    -1000.000   -3000.0 -2000.000
3352""")
3353        fid.close()
3354       
3355        fid = open(elevation_dir_filename1, 'w')
3356        fid.write(""" ncols             4
3357 nrows             4
3358 xllcorner    148.00000
3359 yllcorner    -38.00000
3360 cellsize       0.25
3361 nodata_value   -9999.0
3362   -900.000    -100.000   -300.0 -200.000
3363   -100.000    900.000  -100.000 -300.000
3364   -400.000    600.000   200.000 -500.000
3365   -900.000    -100.000   -300.0 -200.000
3366""")
3367        fid.close()
3368
3369        fid = open(elevation_dir_filename2, 'w')
3370        fid.write(""" ncols             4
3371 nrows             4
3372 xllcorner    148.00000
3373 yllcorner    -38.00000
3374 cellsize       0.25
3375 nodata_value   -9999.0
3376   -990.000    -110.000   -330.0 -220.000
3377   -110.000    990.000  -110.000 -330.000
3378   -440.000    660.000   220.000 -550.000
3379   -990.000    -110.000   -330.0 -220.000
3380""")
3381        fid.close()
3382         
3383        ucur_dir =  tempfile.mkdtemp()
3384        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3385        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3386       
3387        fid = open(ucur_dir_filename1, 'w')
3388        fid.write(""" ncols             4
3389 nrows             4
3390 xllcorner    148.00000
3391 yllcorner    -38.00000
3392 cellsize       0.25
3393 nodata_value   -9999.0
3394   -90.000    -10.000   -30.0 -20.000
3395   -10.000    90.000  -10.000 -30.000
3396   -40.000    60.000   20.000 -50.000
3397   -90.000    -10.000   -30.0 -20.000
3398""")
3399        fid.close()
3400        fid = open(ucur_dir_filename2, 'w')
3401        fid.write(""" ncols             4
3402 nrows             4
3403 xllcorner    148.00000
3404 yllcorner    -38.00000
3405 cellsize       0.25
3406 nodata_value   -9999.0
3407   -90.000    -10.000   -30.0 -20.000
3408   -10.000    99.000  -11.000 -30.000
3409   -40.000    66.000   22.000 -50.000
3410   -90.000    -10.000   -30.0 -20.000
3411""")
3412        fid.close()
3413       
3414        vcur_dir =  tempfile.mkdtemp()
3415        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3416        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3417       
3418        fid = open(vcur_dir_filename1, 'w')
3419        fid.write(""" ncols             4
3420 nrows             4
3421 xllcorner    148.00000
3422 yllcorner    -38.00000
3423 cellsize       0.25
3424 nodata_value   -9999.0
3425   -90.000    -10.000   -30.0 -20.000
3426   -10.000    80.000  -20.000 -30.000
3427   -40.000    50.000   10.000 -50.000
3428   -90.000    -10.000   -30.0 -20.000
3429""")
3430        fid.close()
3431        fid = open(vcur_dir_filename2, 'w')
3432        fid.write(""" ncols             4
3433 nrows             4
3434 xllcorner    148.00000
3435 yllcorner    -38.00000
3436 cellsize       0.25
3437 nodata_value   -9999.0
3438   -90.000    -10.000   -30.0 -20.000
3439   -10.000    88.000  -22.000 -30.000
3440   -40.000    55.000   11.000 -50.000
3441   -90.000    -10.000   -30.0 -20.000
3442""")
3443        fid.close()
3444       
3445        sww_file = tempfile.mktemp(".sww")
3446        #sww_file = 'a_test.sww'
3447        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3448                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3449                      mean_stage = 100,
3450                       minlat = -37.6, maxlat = -37.6,
3451                  minlon = 148.3, maxlon = 148.3
3452                      #,verbose = True
3453                      )
3454
3455        # check the sww file
3456       
3457        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3458        x = fid.variables['x'][:]
3459        y = fid.variables['y'][:]
3460        z = fid.variables['z'][:]
3461        stage = fid.variables['stage'][:]
3462        xmomentum = fid.variables['xmomentum'][:]
3463        ymomentum = fid.variables['ymomentum'][:]
3464        geo_ref = Geo_reference(NetCDFObject=fid)
3465        #print "geo_ref",geo_ref
3466        x_ref = geo_ref.get_xllcorner()
3467        y_ref = geo_ref.get_yllcorner()
3468        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3469
3470        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3471        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3472        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3473
3474        #Easting:  632145.632  Northing: 5820863.269
3475        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3476       
3477        #print "x",x
3478        #print "y",y
3479        self.failUnless(len(x) == 4,'failed') # 2*2
3480        self.failUnless(len(x) == 4,'failed') # 2*2
3481       
3482        #Zone:   55
3483        #Easting:  632145.632  Northing: 5820863.269
3484        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3485        # magic number - y is close enough for me.
3486        assert allclose(x[3], 632145.63 - x_ref)
3487        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3488
3489        assert allclose(z[0],9000.0 ) #z is elevation info
3490        #print "z",z
3491        # 2 time steps, 4 points
3492        self.failUnless(xmomentum.shape == (2,4), 'failed') 
3493        self.failUnless(ymomentum.shape == (2,4), 'failed') 
3494       
3495        #(100.0 - -1000.000)*10
3496        #assert allclose(xmomentum[0][5], 11000.0 )
3497       
3498        fid.close()
3499       
3500        # is the sww file readable?
3501        #Lets see if we can convert it to a dem!
3502        #print "sww_file",sww_file
3503        #dem_file = tempfile.mktemp(".dem")
3504        domain = sww2domain(sww_file) ###, dem_file)
3505        domain.check_integrity()
3506
3507        #tidy up
3508        os.remove(bath_dir_filename)
3509        os.rmdir(bath_dir)
3510
3511        os.remove(elevation_dir_filename1)
3512        os.remove(elevation_dir_filename2)
3513        os.rmdir(elevation_dir)
3514       
3515        os.remove(ucur_dir_filename1)
3516        os.remove(ucur_dir_filename2)
3517        os.rmdir(ucur_dir)
3518       
3519        os.remove(vcur_dir_filename1)
3520        os.remove(vcur_dir_filename2)
3521        os.rmdir(vcur_dir)
3522
3523
3524     
3525       
3526        # remove sww file
3527        os.remove(sww_file)
3528       
3529        # remove dem file
3530        #os.remove(dem_file)
3531
3532    def test_get_min_max_indexes(self):
3533        latitudes = [3,2,1,0]
3534        longitudes = [0,10,20,30]
3535
3536        # k - lat
3537        # l - lon
3538        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3539            latitudes,longitudes,
3540            -10,4,-10,31)
3541
3542        #print "kmin",kmin;print "kmax",kmax
3543        #print "lmin",lmin;print "lmax",lmax
3544        latitudes_new = latitudes[kmin:kmax]
3545        longitudes_news = longitudes[lmin:lmax]
3546        #print "latitudes_new", latitudes_new
3547        #print "longitudes_news",longitudes_news
3548        self.failUnless(latitudes == latitudes_new and \
3549                        longitudes == longitudes_news,
3550                         'failed')
3551       
3552        ## 2nd test
3553        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3554            latitudes,longitudes,
3555            0.5,2.5,5,25)
3556        #print "kmin",kmin;print "kmax",kmax
3557        #print "lmin",lmin;print "lmax",lmax
3558        latitudes_new = latitudes[kmin:kmax]
3559        longitudes_news = longitudes[lmin:lmax]
3560        #print "latitudes_new", latitudes_new
3561        #print "longitudes_news",longitudes_news
3562       
3563        self.failUnless(latitudes == latitudes_new and \
3564                        longitudes == longitudes_news,
3565                         'failed')
3566       
3567        ## 3rd test
3568        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3569                                                                    longitudes,
3570                                                      1.1,1.9,12,17)
3571        #print "kmin",kmin;print "kmax",kmax
3572        #print "lmin",lmin;print "lmax",lmax
3573        latitudes_new = latitudes[kmin:kmax]
3574        longitudes_news = longitudes[lmin:lmax]
3575        #print "latitudes_new", latitudes_new
3576        #print "longitudes_news",longitudes_news
3577       
3578        self.failUnless(latitudes_new == [2, 1] and \
3579                        longitudes_news == [10, 20],
3580                         'failed')
3581
3582       
3583        ## 4th test
3584        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3585            latitudes,longitudes,
3586                                                      -0.1,1.9,-2,17)
3587        #print "kmin",kmin;print "kmax",kmax
3588        #print "lmin",lmin;print "lmax",lmax
3589        latitudes_new = latitudes[kmin:kmax]
3590        longitudes_news = longitudes[lmin:lmax]
3591        #print "latitudes_new", latitudes_new
3592        #print "longitudes_news",longitudes_news
3593       
3594        self.failUnless(latitudes_new == [2, 1, 0] and \
3595                        longitudes_news == [0, 10, 20],
3596                         'failed')
3597        ## 5th test
3598        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3599            latitudes,longitudes,
3600            0.1,1.9,2,17)
3601        #print "kmin",kmin;print "kmax",kmax
3602        #print "lmin",lmin;print "lmax",lmax
3603        latitudes_new = latitudes[kmin:kmax]
3604        longitudes_news = longitudes[lmin:lmax]
3605        #print "latitudes_new", latitudes_new
3606        #print "longitudes_news",longitudes_news
3607       
3608        self.failUnless(latitudes_new == [2, 1, 0] and \
3609                        longitudes_news == [0, 10, 20],
3610                         'failed')
3611       
3612        ## 6th test
3613       
3614        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3615            latitudes,longitudes,
3616            1.5,4,18,32)
3617        #print "kmin",kmin;print "kmax",kmax
3618        #print "lmin",lmin;print "lmax",lmax
3619        latitudes_new = latitudes[kmin:kmax]
3620        longitudes_news = longitudes[lmin:lmax]
3621        #print "latitudes_new", latitudes_new
3622        #print "longitudes_news",longitudes_news
3623       
3624        self.failUnless(latitudes_new == [3, 2, 1] and \
3625                        longitudes_news == [10, 20, 30],
3626                         'failed')
3627       
3628       
3629        ## 7th test
3630        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3631        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3632            latitudes,longitudes,
3633            1.5,1.5,15,15)
3634        #print "kmin",kmin;print "kmax",kmax
3635        #print "lmin",lmin;print "lmax",lmax
3636        latitudes_new = latitudes[kmin:kmax]
3637        longitudes_news = longitudes[lmin:lmax]
3638        m2d = m2d[kmin:kmax,lmin:lmax]
3639        #print "m2d", m2d
3640        #print "latitudes_new", latitudes_new
3641        #print "longitudes_news",longitudes_news
3642       
3643        self.failUnless(latitudes_new == [2, 1] and \
3644                        longitudes_news == [10, 20],
3645                         'failed')
3646       
3647        self.failUnless(m2d == [[5,6],[9,10]],
3648                         'failed')
3649       
3650    def test_get_min_max_indexes2(self):
3651        latitudes = [-30,-35,-40,-45]
3652        longitudes = [148,149,150,151]
3653
3654        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3655       
3656        # k - lat
3657        # l - lon
3658        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3659            latitudes,longitudes,
3660            -37,-27,147,149.5)
3661
3662        #print "kmin",kmin;print "kmax",kmax
3663        #print "lmin",lmin;print "lmax",lmax
3664        #print "m2d", m2d
3665        #print "latitudes", latitudes
3666        #print "longitudes",longitudes
3667        #print "latitudes[kmax]", latitudes[kmax]
3668        latitudes_new = latitudes[kmin:kmax]
3669        longitudes_new = longitudes[lmin:lmax]
3670        m2d = m2d[kmin:kmax,lmin:lmax]
3671        #print "m2d", m2d
3672        #print "latitudes_new", latitudes_new
3673        #print "longitudes_new",longitudes_new
3674       
3675        self.failUnless(latitudes_new == [-30, -35, -40] and \
3676                        longitudes_new == [148, 149,150],
3677                         'failed')
3678        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
3679                         'failed')
3680       
3681    def test_get_min_max_indexes3(self):
3682        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3683        longitudes = [148,149,150,151]
3684
3685        # k - lat
3686        # l - lon
3687        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3688            latitudes,longitudes,
3689            -43,-37,148.5,149.5)
3690
3691       
3692        #print "kmin",kmin;print "kmax",kmax
3693        #print "lmin",lmin;print "lmax",lmax
3694        #print "latitudes", latitudes
3695        #print "longitudes",longitudes
3696        latitudes_new = latitudes[kmin:kmax]
3697        longitudes_news = longitudes[lmin:lmax]
3698        #print "latitudes_new", latitudes_new
3699        #print "longitudes_news",longitudes_news
3700       
3701        self.failUnless(latitudes_new == [-35, -40, -45] and \
3702                        longitudes_news == [148, 149,150],
3703                         'failed')
3704       
3705    def test_get_min_max_indexes4(self):
3706        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3707        longitudes = [148,149,150,151]
3708
3709        # k - lat
3710        # l - lon
3711        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3712            latitudes,longitudes)
3713       
3714        #print "kmin",kmin;print "kmax",kmax
3715        #print "lmin",lmin;print "lmax",lmax
3716        #print "latitudes", latitudes
3717        #print "longitudes",longitudes
3718        latitudes_new = latitudes[kmin:kmax]
3719        longitudes_news = longitudes[lmin:lmax]
3720        #print "latitudes_new", latitudes_new
3721        #print "longitudes_news",longitudes_news
3722       
3723        self.failUnless(latitudes_new == latitudes  and \
3724                        longitudes_news == longitudes,
3725                         'failed')
3726       
3727    def test_tsh2sww(self):
3728        import os
3729        import tempfile
3730       
3731        tsh_file = tempfile.mktemp(".tsh")
3732        file = open(tsh_file,"w")
3733        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
37340 0.0 0.0 0.0 0.0 0.01 \n \
37351 1.0 0.0 10.0 10.0 0.02  \n \
37362 0.0 1.0 0.0 10.0 0.03  \n \
37373 0.5 0.25 8.0 12.0 0.04  \n \
3738# Vert att title  \n \
3739elevation  \n \
3740stage  \n \
3741friction  \n \
37422 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
37430 0 3 2 -1  -1  1 dsg\n\
37441 0 1 3 -1  0 -1   ole nielsen\n\
37454 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
37460 1 0 2 \n\
37471 0 2 3 \n\
37482 2 3 \n\
37493 3 1 1 \n\
37503 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
37510 216.0 -86.0 \n \
37521 160.0 -167.0 \n \
37532 114.0 -91.0 \n \
37543 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
37550 0 1 0 \n \
37561 1 2 0 \n \
37572 2 0 0 \n \
37580 # <x> <y> ...Mesh Holes... \n \
37590 # <x> <y> <attribute>...Mesh Regions... \n \
37600 # <x> <y> <attribute>...Mesh Regions, area... \n\
3761#Geo reference \n \
376256 \n \
3763140 \n \
3764120 \n")
3765        file.close()
3766
3767        #sww_file = tempfile.mktemp(".sww")
3768        #print "sww_file",sww_file
3769        #print "sww_file",tsh_file
3770        tsh2sww(tsh_file)
3771
3772        os.remove(tsh_file)
3773        os.remove(tsh_file[:-4] + '.sww')
3774#-------------------------------------------------------------
3775if __name__ == "__main__":
3776    #suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
3777    suite = unittest.makeSuite(Test_Data_Manager,'test')
3778    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
3779    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')   
3780    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
3781    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
3782    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
3783    runner = unittest.TextTestRunner()
3784    runner.run(suite)
3785 
Note: See TracBrowser for help on using the repository browser.