source: inundation/pyvolution/test_data_manager.py @ 2296

Last change on this file since 2296 was 2296, checked in by duncan, 18 years ago

comments, trying to get it to delete files it produces. Not fully done.

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