source: inundation/pyvolution/test_data_manager.py @ 2492

Last change on this file since 2492 was 2492, checked in by sexton, 18 years ago

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