source: anuga_core/source/anuga/shallow_water/test_data_manager.py @ 3750

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

more error handling

File size: 161.7 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import copy
6from Numeric import zeros, array, allclose, Float
7from anuga.utilities.numerical_tools import mean
8import tempfile
9import os
10from Scientific.IO.NetCDF import NetCDFFile
11from struct import pack
12
13from anuga.shallow_water import *
14from anuga.shallow_water.data_manager import *
15from anuga.config import epsilon
16from anuga.utilities.anuga_exceptions import ANUGAError
17from anuga.utilities.numerical_tools import ensure_numeric
18
19# This is needed to run the tests of local functions
20import data_manager 
21from anuga.coordinate_transforms.redfearn import redfearn
22from anuga.coordinate_transforms.geo_reference import Geo_reference
23
24class Test_Data_Manager(unittest.TestCase):
25    def setUp(self):
26        import time
27        from mesh_factory import rectangular
28
29        #Create basic mesh
30        points, vertices, boundary = rectangular(2, 2)
31
32        #Create shallow water domain
33        domain = Domain(points, vertices, boundary)
34        domain.default_order=2
35
36
37        #Set some field values
38        domain.set_quantity('elevation', lambda x,y: -x)
39        domain.set_quantity('friction', 0.03)
40
41
42        ######################
43        # Boundary conditions
44        B = Transmissive_boundary(domain)
45        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
46
47
48        ######################
49        #Initial condition - with jumps
50
51
52        bed = domain.quantities['elevation'].vertex_values
53        stage = zeros(bed.shape, Float)
54
55        h = 0.3
56        for i in range(stage.shape[0]):
57            if i % 2 == 0:
58                stage[i,:] = bed[i,:] + h
59            else:
60                stage[i,:] = bed[i,:]
61
62        domain.set_quantity('stage', stage)
63        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
64
65        domain.distribute_to_vertices_and_edges()
66
67
68        self.domain = domain
69
70        C = domain.get_vertex_coordinates()
71        self.X = C[:,0:6:2].copy()
72        self.Y = C[:,1:6:2].copy()
73
74        self.F = bed
75
76        #Write A testfile (not realistic. Values aren't realistic)
77        self.test_MOST_file = 'most_small'
78
79        longitudes = [150.66667, 150.83334, 151., 151.16667]
80        latitudes = [-34.5, -34.33333, -34.16667, -34]
81
82        long_name = 'LON'
83        lat_name = 'LAT'
84
85        nx = 4
86        ny = 4
87        six = 6
88
89
90        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
91            fid = NetCDFFile(self.test_MOST_file + ext, 'w')
92
93            fid.createDimension(long_name,nx)
94            fid.createVariable(long_name,'d',(long_name,))
95            fid.variables[long_name].point_spacing='uneven'
96            fid.variables[long_name].units='degrees_east'
97            fid.variables[long_name].assignValue(longitudes)
98
99            fid.createDimension(lat_name,ny)
100            fid.createVariable(lat_name,'d',(lat_name,))
101            fid.variables[lat_name].point_spacing='uneven'
102            fid.variables[lat_name].units='degrees_north'
103            fid.variables[lat_name].assignValue(latitudes)
104
105            fid.createDimension('TIME',six)
106            fid.createVariable('TIME','d',('TIME',))
107            fid.variables['TIME'].point_spacing='uneven'
108            fid.variables['TIME'].units='seconds'
109            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
110
111
112            name = ext[1:3].upper()
113            if name == 'E.': name = 'ELEVATION'
114            fid.createVariable(name,'d',('TIME', lat_name, long_name))
115            fid.variables[name].units='CENTIMETERS'
116            fid.variables[name].missing_value=-1.e+034
117
118            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
119                                              [-0.1214216, 0, 0, 0],
120                                              [0, 0, 0, 0],
121                                              [0, 0, 0, 0]],
122                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
123                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
124                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
125                                              [0, 2.291054e-005, 2.291054e-005, 0]],
126                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
127                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
128                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
129                                              [0, 0.0001374632, 0.0001374632, 0]],
130                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
131                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
132                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
133                                              [0, 0.0002520159, 0.0002520159, 0]],
134                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
135                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
136                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
137                                              [0, 0.0003665686, 0.0003665686, 0]],
138                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
139                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
140                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
141                                              [0, 0.0004811212, 0.0004811212, 0]]])
142
143
144            fid.close()
145
146
147
148
149    def tearDown(self):
150        import os
151        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
152            #print 'Trying to remove', self.test_MOST_file + ext
153            os.remove(self.test_MOST_file + ext)
154
155    def test_sww_constant(self):
156        """Test that constant sww information can be written correctly
157        (non smooth)
158        """
159
160        import time, os
161        from Numeric import array, zeros, allclose, Float, concatenate
162        from Scientific.IO.NetCDF import NetCDFFile
163
164        self.domain.filename = 'datatest' + str(id(self))
165        self.domain.format = 'sww'
166        self.domain.smooth = False
167
168        sww = get_dataobject(self.domain)
169        sww.store_connectivity()
170
171        #Check contents
172        #Get NetCDF
173        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
174
175        # Get the variables
176        x = fid.variables['x']
177        y = fid.variables['y']
178        z = fid.variables['elevation']
179
180        volumes = fid.variables['volumes']
181
182
183        assert allclose (x[:], self.X.flat)
184        assert allclose (y[:], self.Y.flat)
185        assert allclose (z[:], self.F.flat)
186
187        V = volumes
188
189        P = len(self.domain)
190        for k in range(P):
191            assert V[k, 0] == 3*k
192            assert V[k, 1] == 3*k+1
193            assert V[k, 2] == 3*k+2
194
195
196        fid.close()
197
198        #Cleanup
199        os.remove(sww.filename)
200
201
202    def test_sww_constant_smooth(self):
203        """Test that constant sww information can be written correctly
204        (non smooth)
205        """
206
207        import time, os
208        from Numeric import array, zeros, allclose, Float, concatenate
209        from Scientific.IO.NetCDF import NetCDFFile
210
211        self.domain.filename = 'datatest' + str(id(self))
212        self.domain.format = 'sww'
213        self.domain.smooth = True
214
215        sww = get_dataobject(self.domain)
216        sww.store_connectivity()
217
218        #Check contents
219        #Get NetCDF
220        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
221
222        # Get the variables
223        x = fid.variables['x']
224        y = fid.variables['y']
225        z = fid.variables['elevation']
226
227        volumes = fid.variables['volumes']
228
229        X = x[:]
230        Y = y[:]
231
232        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
233        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
234        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
235
236        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
237
238        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
239
240        Z = z[:]
241        assert Z[4] == -0.5
242
243        V = volumes
244        assert V[2,0] == 4
245        assert V[2,1] == 5
246        assert V[2,2] == 1
247
248        assert V[4,0] == 6
249        assert V[4,1] == 7
250        assert V[4,2] == 3
251
252
253        fid.close()
254
255        #Cleanup
256        os.remove(sww.filename)
257
258
259
260    def test_sww_variable(self):
261        """Test that sww information can be written correctly
262        """
263
264        import time, os
265        from Numeric import array, zeros, allclose, Float, concatenate
266        from Scientific.IO.NetCDF import NetCDFFile
267
268        self.domain.filename = 'datatest' + str(id(self))
269        self.domain.format = 'sww'
270        self.domain.smooth = True
271        self.domain.reduction = mean
272
273        sww = get_dataobject(self.domain)
274        sww.store_connectivity()
275        sww.store_timestep('stage')
276
277        #Check contents
278        #Get NetCDF
279        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
280
281
282        # Get the variables
283        x = fid.variables['x']
284        y = fid.variables['y']
285        z = fid.variables['elevation']
286        time = fid.variables['time']
287        stage = fid.variables['stage']
288
289
290        Q = self.domain.quantities['stage']
291        Q0 = Q.vertex_values[:,0]
292        Q1 = Q.vertex_values[:,1]
293        Q2 = Q.vertex_values[:,2]
294
295        A = stage[0,:]
296        #print A[0], (Q2[0,0] + Q1[1,0])/2
297        assert allclose(A[0], (Q2[0] + Q1[1])/2)
298        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
299        assert allclose(A[2], Q0[3])
300        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
301
302        #Center point
303        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
304                                 Q0[5] + Q2[6] + Q1[7])/6)
305
306
307
308        fid.close()
309
310        #Cleanup
311        os.remove(sww.filename)
312
313
314    def test_sww_variable2(self):
315        """Test that sww information can be written correctly
316        multiple timesteps. Use average as reduction operator
317        """
318
319        import time, os
320        from Numeric import array, zeros, allclose, Float, concatenate
321        from Scientific.IO.NetCDF import NetCDFFile
322
323        self.domain.filename = 'datatest' + str(id(self))
324        self.domain.format = 'sww'
325        self.domain.smooth = True
326
327        self.domain.reduction = mean
328
329        sww = get_dataobject(self.domain)
330        sww.store_connectivity()
331        sww.store_timestep('stage')
332        self.domain.evolve_to_end(finaltime = 0.01)
333        sww.store_timestep('stage')
334
335
336        #Check contents
337        #Get NetCDF
338        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
339
340        # Get the variables
341        x = fid.variables['x']
342        y = fid.variables['y']
343        z = fid.variables['elevation']
344        time = fid.variables['time']
345        stage = fid.variables['stage']
346
347        #Check values
348        Q = self.domain.quantities['stage']
349        Q0 = Q.vertex_values[:,0]
350        Q1 = Q.vertex_values[:,1]
351        Q2 = Q.vertex_values[:,2]
352
353        A = stage[1,:]
354        assert allclose(A[0], (Q2[0] + Q1[1])/2)
355        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
356        assert allclose(A[2], Q0[3])
357        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
358
359        #Center point
360        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
361                                 Q0[5] + Q2[6] + Q1[7])/6)
362
363
364        fid.close()
365
366        #Cleanup
367        os.remove(sww.filename)
368
369    def test_sww_variable3(self):
370        """Test that sww information can be written correctly
371        multiple timesteps using a different reduction operator (min)
372        """
373
374        import time, os
375        from Numeric import array, zeros, allclose, Float, concatenate
376        from Scientific.IO.NetCDF import NetCDFFile
377
378        self.domain.filename = 'datatest' + str(id(self))
379        self.domain.format = 'sww'
380        self.domain.smooth = True
381        self.domain.reduction = min
382
383        sww = get_dataobject(self.domain)
384        sww.store_connectivity()
385        sww.store_timestep('stage')
386
387        self.domain.evolve_to_end(finaltime = 0.01)
388        sww.store_timestep('stage')
389
390
391        #Check contents
392        #Get NetCDF
393        fid = NetCDFFile(sww.filename, 'r')
394
395
396        # Get the variables
397        x = fid.variables['x']
398        y = fid.variables['y']
399        z = fid.variables['elevation']
400        time = fid.variables['time']
401        stage = fid.variables['stage']
402
403        #Check values
404        Q = self.domain.quantities['stage']
405        Q0 = Q.vertex_values[:,0]
406        Q1 = Q.vertex_values[:,1]
407        Q2 = Q.vertex_values[:,2]
408
409        A = stage[1,:]
410        assert allclose(A[0], min(Q2[0], Q1[1]))
411        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
412        assert allclose(A[2], Q0[3])
413        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
414
415        #Center point
416        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
417                                  Q0[5], Q2[6], Q1[7]))
418
419
420        fid.close()
421
422        #Cleanup
423        os.remove(sww.filename)
424
425
426    def test_sync(self):
427        """Test info stored at each timestep is as expected (incl initial condition)
428        """
429
430        import time, os, config
431        from Numeric import array, zeros, allclose, Float, concatenate
432        from Scientific.IO.NetCDF import NetCDFFile
433
434        self.domain.filename = 'synctest'
435        self.domain.format = 'sww'
436        self.domain.smooth = False
437        self.domain.store = True
438        self.domain.beta_h = 0
439
440        #Evolution
441        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
442            stage = self.domain.quantities['stage'].vertex_values
443
444            #Get NetCDF
445            fid = NetCDFFile(self.domain.writer.filename, 'r')
446            stage_file = fid.variables['stage']
447
448            if t == 0.0:
449                assert allclose(stage, self.initial_stage)
450                assert allclose(stage_file[:], stage.flat)
451            else:
452                assert not allclose(stage, self.initial_stage)
453                assert not allclose(stage_file[:], stage.flat)
454
455            fid.close()
456
457        os.remove(self.domain.writer.filename)
458
459
460    def test_sww_minimum_storable_height(self):
461        """Test that sww information can be written correctly
462        multiple timesteps using a different reduction operator (min)
463        """
464
465        import time, os
466        from Numeric import array, zeros, allclose, Float, concatenate
467        from Scientific.IO.NetCDF import NetCDFFile
468
469        self.domain.filename = 'datatest' + str(id(self))
470        self.domain.format = 'sww'
471        self.domain.smooth = True
472        self.domain.reduction = min
473        self.domain.minimum_storable_height = 100
474
475        sww = get_dataobject(self.domain)
476        sww.store_connectivity()
477        sww.store_timestep('stage')
478
479        self.domain.evolve_to_end(finaltime = 0.01)
480        sww.store_timestep('stage')
481
482
483        #Check contents
484        #Get NetCDF
485        fid = NetCDFFile(sww.filename, 'r')
486
487
488        # Get the variables
489        x = fid.variables['x']
490        y = fid.variables['y']
491        z = fid.variables['elevation']
492        time = fid.variables['time']
493        stage = fid.variables['stage']
494
495        #Check values
496        Q = self.domain.quantities['stage']
497        Q0 = Q.vertex_values[:,0]
498        Q1 = Q.vertex_values[:,1]
499        Q2 = Q.vertex_values[:,2]
500
501        A = stage[1,:]
502        assert allclose(stage[1,:], z[:])
503        fid.close()
504
505        #Cleanup
506        os.remove(sww.filename)
507
508
509    def Not_a_test_sww_DSG(self):
510        """Not a test, rather a look at the sww format
511        """
512
513        import time, os
514        from Numeric import array, zeros, allclose, Float, concatenate
515        from Scientific.IO.NetCDF import NetCDFFile
516
517        self.domain.filename = 'datatest' + str(id(self))
518        self.domain.format = 'sww'
519        self.domain.smooth = True
520        self.domain.reduction = mean
521
522        sww = get_dataobject(self.domain)
523        sww.store_connectivity()
524        sww.store_timestep('stage')
525
526        #Check contents
527        #Get NetCDF
528        fid = NetCDFFile(sww.filename, 'r')
529
530        # Get the variables
531        x = fid.variables['x']
532        y = fid.variables['y']
533        z = fid.variables['elevation']
534
535        volumes = fid.variables['volumes']
536        time = fid.variables['time']
537
538        # 2D
539        stage = fid.variables['stage']
540
541        X = x[:]
542        Y = y[:]
543        Z = z[:]
544        V = volumes[:]
545        T = time[:]
546        S = stage[:,:]
547
548#         print "****************************"
549#         print "X ",X
550#         print "****************************"
551#         print "Y ",Y
552#         print "****************************"
553#         print "Z ",Z
554#         print "****************************"
555#         print "V ",V
556#         print "****************************"
557#         print "Time ",T
558#         print "****************************"
559#         print "Stage ",S
560#         print "****************************"
561
562
563        fid.close()
564
565        #Cleanup
566        os.remove(sww.filename)
567
568
569    #def test_write_pts(self):
570    #    #Obsolete
571    #
572    #    #Get (enough) datapoints
573    #
574    #    from Numeric import array
575    #    points = array([[ 0.66666667, 0.66666667],
576    #                    [ 1.33333333, 1.33333333],
577    #                    [ 2.66666667, 0.66666667],
578    #                    [ 0.66666667, 2.66666667],
579    #                    [ 0.0, 1.0],
580    #                    [ 0.0, 3.0],
581    #                    [ 1.0, 0.0],
582    #                    [ 1.0, 1.0],
583    #                    [ 1.0, 2.0],
584    #                    [ 1.0, 3.0],
585    #                    [ 2.0, 1.0],
586    #                    [ 3.0, 0.0],
587    #                    [ 3.0, 1.0]])
588    #
589    #    z = points[:,0] + 2*points[:,1]
590    #
591    #    ptsfile = 'testptsfile.pts'
592    #    write_ptsfile(ptsfile, points, z,
593    #                  attribute_name = 'linear_combination')
594    #
595    #    #Check contents
596    #    #Get NetCDF
597    #    from Scientific.IO.NetCDF import NetCDFFile
598    #    fid = NetCDFFile(ptsfile, 'r')
599    #
600    #    # Get the variables
601    #    #print fid.variables.keys()
602    #    points1 = fid.variables['points']
603    #    z1 = fid.variables['linear_combination']
604    #
605    #    #Check values#
606    #
607    #    #print points[:]
608    #    #print ref_points
609    #    assert allclose(points, points1)
610    #
611    #    #print attributes[:]
612    #    #print ref_elevation
613    #    assert allclose(z, z1)
614    #
615    #    #Cleanup
616    #    fid.close()
617    #
618    #    import os
619    #    os.remove(ptsfile)
620
621
622    def test_dem2pts_bounding_box_v2(self):
623        """Test conversion from dem in ascii format to native NetCDF xya format
624        """
625
626        import time, os
627        from Numeric import array, zeros, allclose, Float, concatenate, ones
628        from Scientific.IO.NetCDF import NetCDFFile
629
630        #Write test asc file
631        root = 'demtest'
632
633        filename = root+'.asc'
634        fid = open(filename, 'w')
635        fid.write("""ncols         10
636nrows         10
637xllcorner     2000
638yllcorner     3000
639cellsize      1
640NODATA_value  -9999
641""")
642        #Create linear function
643        ref_points = []
644        ref_elevation = []
645        x0 = 2000
646        y = 3010
647        yvec = range(10)
648        xvec = range(10)
649        z = -1
650        for i in range(10):
651            y = y - 1
652            for j in range(10):
653                x = x0 + xvec[j]
654                z += 1
655                ref_points.append ([x,y])
656                ref_elevation.append(z)
657                fid.write('%f ' %z)
658            fid.write('\n')
659
660        fid.close()
661
662        #print 'sending pts', ref_points
663        #print 'sending elev', ref_elevation
664
665        #Write prj file with metadata
666        metafilename = root+'.prj'
667        fid = open(metafilename, 'w')
668
669
670        fid.write("""Projection UTM
671Zone 56
672Datum WGS84
673Zunits NO
674Units METERS
675Spheroid WGS84
676Xshift 0.0000000000
677Yshift 10000000.0000000000
678Parameters
679""")
680        fid.close()
681
682        #Convert to NetCDF pts
683        convert_dem_from_ascii2netcdf(root)
684        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
685                northing_min=3003.0, northing_max=3006.0)
686
687        #Check contents
688        #Get NetCDF
689        fid = NetCDFFile(root+'.pts', 'r')
690
691        # Get the variables
692        #print fid.variables.keys()
693        points = fid.variables['points']
694        elevation = fid.variables['elevation']
695
696        #Check values
697        assert fid.xllcorner[0] == 2002.0
698        assert fid.yllcorner[0] == 3003.0
699
700        #create new reference points
701        newz = []
702        newz[0:5] = ref_elevation[32:38]
703        newz[6:11] = ref_elevation[42:48]
704        newz[12:17] = ref_elevation[52:58]
705        newz[18:23] = ref_elevation[62:68]
706        ref_elevation = []
707        ref_elevation = newz
708        ref_points = []
709        x0 = 2002
710        y = 3007
711        yvec = range(4)
712        xvec = range(6)
713        for i in range(4):
714            y = y - 1
715            ynew = y - 3003.0
716            for j in range(6):
717                x = x0 + xvec[j]
718                xnew = x - 2002.0
719                ref_points.append ([xnew,ynew]) #Relative point values
720
721        assert allclose(points, ref_points)
722
723        assert allclose(elevation, ref_elevation)
724
725        #Cleanup
726        fid.close()
727
728
729        os.remove(root + '.pts')
730        os.remove(root + '.dem')
731        os.remove(root + '.asc')
732        os.remove(root + '.prj')
733
734
735    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
736        """Test conversion from dem in ascii format to native NetCDF xya format
737        """
738
739        import time, os
740        from Numeric import array, zeros, allclose, Float, concatenate, ones
741        from Scientific.IO.NetCDF import NetCDFFile
742
743        #Write test asc file
744        root = 'demtest'
745
746        filename = root+'.asc'
747        fid = open(filename, 'w')
748        fid.write("""ncols         10
749nrows         10
750xllcorner     2000
751yllcorner     3000
752cellsize      1
753NODATA_value  -9999
754""")
755        #Create linear function
756        ref_points = []
757        ref_elevation = []
758        x0 = 2000
759        y = 3010
760        yvec = range(10)
761        xvec = range(10)
762        #z = range(100)
763        z = zeros(100)
764        NODATA_value = -9999
765        count = -1
766        for i in range(10):
767            y = y - 1
768            for j in range(10):
769                x = x0 + xvec[j]
770                ref_points.append ([x,y])
771                count += 1
772                z[count] = (4*i - 3*j)%13
773                if j == 4: z[count] = NODATA_value #column inside clipping region
774                if j == 8: z[count] = NODATA_value #column outside clipping region
775                if i == 9: z[count] = NODATA_value #row outside clipping region
776                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
777                ref_elevation.append( z[count] )
778                fid.write('%f ' %z[count])
779            fid.write('\n')
780
781        fid.close()
782
783        #print 'sending elev', ref_elevation
784
785        #Write prj file with metadata
786        metafilename = root+'.prj'
787        fid = open(metafilename, 'w')
788
789
790        fid.write("""Projection UTM
791Zone 56
792Datum WGS84
793Zunits NO
794Units METERS
795Spheroid WGS84
796Xshift 0.0000000000
797Yshift 10000000.0000000000
798Parameters
799""")
800        fid.close()
801
802        #Convert to NetCDF pts
803        convert_dem_from_ascii2netcdf(root)
804        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
805                northing_min=3003.0, northing_max=3006.0)
806
807        #Check contents
808        #Get NetCDF
809        fid = NetCDFFile(root+'.pts', 'r')
810
811        # Get the variables
812        #print fid.variables.keys()
813        points = fid.variables['points']
814        elevation = fid.variables['elevation']
815
816        #Check values
817        assert fid.xllcorner[0] == 2002.0
818        assert fid.yllcorner[0] == 3003.0
819
820        #create new reference points
821        newz = zeros(19)
822        newz[0:2] = ref_elevation[32:34]
823        newz[2:5] = ref_elevation[35:38]
824        newz[5:7] = ref_elevation[42:44]
825        newz[7] = ref_elevation[45]
826        newz[8] = ref_elevation[47]
827        newz[9:11] = ref_elevation[52:54]
828        newz[11:14] = ref_elevation[55:58]
829        newz[14:16] = ref_elevation[62:64]
830        newz[16:19] = ref_elevation[65:68]
831
832
833        ref_elevation = newz
834        ref_points = []
835        new_ref_points = []
836        x0 = 2002
837        y = 3007
838        yvec = range(4)
839        xvec = range(6)
840        for i in range(4):
841            y = y - 1
842            ynew = y - 3003.0
843            for j in range(6):
844                x = x0 + xvec[j]
845                xnew = x - 2002.0
846                if j <> 2 and (i<>1 or j<>4):
847                    ref_points.append([x,y])
848                    new_ref_points.append ([xnew,ynew])
849
850
851        assert allclose(points, new_ref_points)
852        assert allclose(elevation, ref_elevation)
853
854        #Cleanup
855        fid.close()
856
857
858        os.remove(root + '.pts')
859        os.remove(root + '.dem')
860        os.remove(root + '.asc')
861        os.remove(root + '.prj')
862
863
864    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
865        """Test conversion from dem in ascii format to native NetCDF xya format
866        Check missing values on clipping boundary
867        """
868
869        import time, os
870        from Numeric import array, zeros, allclose, Float, concatenate, ones
871        from Scientific.IO.NetCDF import NetCDFFile
872
873        #Write test asc file
874        root = 'demtest'
875
876        filename = root+'.asc'
877        fid = open(filename, 'w')
878        fid.write("""ncols         10
879nrows         10
880xllcorner     2000
881yllcorner     3000
882cellsize      1
883NODATA_value  -9999
884""")
885        #Create linear function
886        ref_points = []
887        ref_elevation = []
888        x0 = 2000
889        y = 3010
890        yvec = range(10)
891        xvec = range(10)
892        #z = range(100)
893        z = zeros(100)
894        NODATA_value = -9999
895        count = -1
896        for i in range(10):
897            y = y - 1
898            for j in range(10):
899                x = x0 + xvec[j]
900                ref_points.append ([x,y])
901                count += 1
902                z[count] = (4*i - 3*j)%13
903                if j == 4: z[count] = NODATA_value #column inside clipping region
904                if j == 8: z[count] = NODATA_value #column outside clipping region
905                if i == 6: z[count] = NODATA_value #row on clipping boundary
906                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
907                ref_elevation.append( z[count] )
908                fid.write('%f ' %z[count])
909            fid.write('\n')
910
911        fid.close()
912
913        #print 'sending elev', ref_elevation
914
915        #Write prj file with metadata
916        metafilename = root+'.prj'
917        fid = open(metafilename, 'w')
918
919
920        fid.write("""Projection UTM
921Zone 56
922Datum WGS84
923Zunits NO
924Units METERS
925Spheroid WGS84
926Xshift 0.0000000000
927Yshift 10000000.0000000000
928Parameters
929""")
930        fid.close()
931
932        #Convert to NetCDF pts
933        convert_dem_from_ascii2netcdf(root)
934        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
935                northing_min=3003.0, northing_max=3006.0)
936
937        #Check contents
938        #Get NetCDF
939        fid = NetCDFFile(root+'.pts', 'r')
940
941        # Get the variables
942        #print fid.variables.keys()
943        points = fid.variables['points']
944        elevation = fid.variables['elevation']
945
946        #Check values
947        assert fid.xllcorner[0] == 2002.0
948        assert fid.yllcorner[0] == 3003.0
949
950        #create new reference points
951        newz = zeros(14)
952        newz[0:2] = ref_elevation[32:34]
953        newz[2:5] = ref_elevation[35:38]
954        newz[5:7] = ref_elevation[42:44]
955        newz[7] = ref_elevation[45]
956        newz[8] = ref_elevation[47]
957        newz[9:11] = ref_elevation[52:54]
958        newz[11:14] = ref_elevation[55:58]
959
960
961
962        ref_elevation = newz
963        ref_points = []
964        new_ref_points = []
965        x0 = 2002
966        y = 3007
967        yvec = range(4)
968        xvec = range(6)
969        for i in range(4):
970            y = y - 1
971            ynew = y - 3003.0
972            for j in range(6):
973                x = x0 + xvec[j]
974                xnew = x - 2002.0
975                if j <> 2 and (i<>1 or j<>4) and i<>3:
976                    ref_points.append([x,y])
977                    new_ref_points.append ([xnew,ynew])
978
979
980        #print points[:],points[:].shape
981        #print new_ref_points, len(new_ref_points)
982
983        assert allclose(elevation, ref_elevation)
984        assert allclose(points, new_ref_points)
985
986
987        #Cleanup
988        fid.close()
989
990
991        os.remove(root + '.pts')
992        os.remove(root + '.dem')
993        os.remove(root + '.asc')
994        os.remove(root + '.prj')
995
996
997    def test_hecras_cross_sections2pts(self):
998        """Test conversion from HECRAS cross sections in ascii format
999        to native NetCDF pts format
1000        """
1001
1002        import time, os
1003        from Numeric import array, zeros, allclose, Float, concatenate
1004        from Scientific.IO.NetCDF import NetCDFFile
1005
1006        #Write test asc file
1007        root = 'hecrastest'
1008
1009        filename = root+'.sdf'
1010        fid = open(filename, 'w')
1011        fid.write("""
1012# RAS export file created on Mon 15Aug2005 11:42
1013# by HEC-RAS Version 3.1.1
1014
1015BEGIN HEADER:
1016  UNITS: METRIC
1017  DTM TYPE: TIN
1018  DTM: v:\1\cit\perth_topo\river_tin
1019  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
1020  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
1021  MAP PROJECTION: UTM
1022  PROJECTION ZONE: 50
1023  DATUM: AGD66
1024  VERTICAL DATUM:
1025  NUMBER OF REACHES:  19
1026  NUMBER OF CROSS-SECTIONS:  2
1027END HEADER:
1028
1029
1030BEGIN CROSS-SECTIONS:
1031
1032  CROSS-SECTION:
1033    STREAM ID:Southern-Wungong
1034    REACH ID:Southern-Wungong
1035    STATION:21410
1036    CUT LINE:
1037      407546.08 , 6437277.542
1038      407329.32 , 6437489.482
1039      407283.11 , 6437541.232
1040    SURFACE LINE:
1041     407546.08,   6437277.54,   52.14
1042     407538.88,   6437284.58,   51.07
1043     407531.68,   6437291.62,   50.56
1044     407524.48,   6437298.66,   49.58
1045     407517.28,   6437305.70,   49.09
1046     407510.08,   6437312.74,   48.76
1047  END:
1048
1049  CROSS-SECTION:
1050    STREAM ID:Swan River
1051    REACH ID:Swan Mouth
1052    STATION:840.*
1053    CUT LINE:
1054      381178.0855 , 6452559.0685
1055      380485.4755 , 6453169.272
1056    SURFACE LINE:
1057     381178.09,   6452559.07,   4.17
1058     381169.49,   6452566.64,   4.26
1059     381157.78,   6452576.96,   4.34
1060     381155.97,   6452578.56,   4.35
1061     381143.72,   6452589.35,   4.43
1062     381136.69,   6452595.54,   4.58
1063     381114.74,   6452614.88,   4.41
1064     381075.53,   6452649.43,   4.17
1065     381071.47,   6452653.00,   3.99
1066     381063.46,   6452660.06,   3.67
1067     381054.41,   6452668.03,   3.67
1068  END:
1069END CROSS-SECTIONS:
1070""")
1071
1072        fid.close()
1073
1074
1075        #Convert to NetCDF pts
1076        hecras_cross_sections2pts(root)
1077
1078        #Check contents
1079        #Get NetCDF
1080        fid = NetCDFFile(root+'.pts', 'r')
1081
1082        # Get the variables
1083        #print fid.variables.keys()
1084        points = fid.variables['points']
1085        elevation = fid.variables['elevation']
1086
1087        #Check values
1088        ref_points = [[407546.08, 6437277.54],
1089                      [407538.88, 6437284.58],
1090                      [407531.68, 6437291.62],
1091                      [407524.48, 6437298.66],
1092                      [407517.28, 6437305.70],
1093                      [407510.08, 6437312.74]]
1094
1095        ref_points += [[381178.09, 6452559.07],
1096                       [381169.49, 6452566.64],
1097                       [381157.78, 6452576.96],
1098                       [381155.97, 6452578.56],
1099                       [381143.72, 6452589.35],
1100                       [381136.69, 6452595.54],
1101                       [381114.74, 6452614.88],
1102                       [381075.53, 6452649.43],
1103                       [381071.47, 6452653.00],
1104                       [381063.46, 6452660.06],
1105                       [381054.41, 6452668.03]]
1106
1107
1108        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
1109        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
1110
1111        #print points[:]
1112        #print ref_points
1113        assert allclose(points, ref_points)
1114
1115        #print attributes[:]
1116        #print ref_elevation
1117        assert allclose(elevation, ref_elevation)
1118
1119        #Cleanup
1120        fid.close()
1121
1122
1123        os.remove(root + '.sdf')
1124        os.remove(root + '.pts')
1125
1126
1127
1128    def test_sww2dem_asc_elevation(self):
1129        """Test that sww information can be converted correctly to asc/prj
1130        format readable by e.g. ArcView
1131        """
1132
1133        import time, os
1134        from Numeric import array, zeros, allclose, Float, concatenate
1135        from Scientific.IO.NetCDF import NetCDFFile
1136
1137        #Setup
1138        self.domain.filename = 'datatest'
1139
1140        prjfile = self.domain.filename + '_elevation.prj'
1141        ascfile = self.domain.filename + '_elevation.asc'
1142        swwfile = self.domain.filename + '.sww'
1143
1144        self.domain.set_datadir('.')
1145        self.domain.format = 'sww'
1146        self.domain.smooth = True
1147        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1148
1149        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1150
1151        sww = get_dataobject(self.domain)
1152        sww.store_connectivity()
1153        sww.store_timestep('stage')
1154
1155        self.domain.evolve_to_end(finaltime = 0.01)
1156        sww.store_timestep('stage')
1157
1158        cellsize = 0.25
1159        #Check contents
1160        #Get NetCDF
1161
1162        fid = NetCDFFile(sww.filename, 'r')
1163
1164        # Get the variables
1165        x = fid.variables['x'][:]
1166        y = fid.variables['y'][:]
1167        z = fid.variables['elevation'][:]
1168        time = fid.variables['time'][:]
1169        stage = fid.variables['stage'][:]
1170
1171
1172        #Export to ascii/prj files
1173        sww2dem(self.domain.filename,
1174                quantity = 'elevation',
1175                cellsize = cellsize,
1176                verbose = False,
1177                format = 'asc')
1178
1179        #Check prj (meta data)
1180        prjid = open(prjfile)
1181        lines = prjid.readlines()
1182        prjid.close()
1183
1184        L = lines[0].strip().split()
1185        assert L[0].strip().lower() == 'projection'
1186        assert L[1].strip().lower() == 'utm'
1187
1188        L = lines[1].strip().split()
1189        assert L[0].strip().lower() == 'zone'
1190        assert L[1].strip().lower() == '56'
1191
1192        L = lines[2].strip().split()
1193        assert L[0].strip().lower() == 'datum'
1194        assert L[1].strip().lower() == 'wgs84'
1195
1196        L = lines[3].strip().split()
1197        assert L[0].strip().lower() == 'zunits'
1198        assert L[1].strip().lower() == 'no'
1199
1200        L = lines[4].strip().split()
1201        assert L[0].strip().lower() == 'units'
1202        assert L[1].strip().lower() == 'meters'
1203
1204        L = lines[5].strip().split()
1205        assert L[0].strip().lower() == 'spheroid'
1206        assert L[1].strip().lower() == 'wgs84'
1207
1208        L = lines[6].strip().split()
1209        assert L[0].strip().lower() == 'xshift'
1210        assert L[1].strip().lower() == '500000'
1211
1212        L = lines[7].strip().split()
1213        assert L[0].strip().lower() == 'yshift'
1214        assert L[1].strip().lower() == '10000000'
1215
1216        L = lines[8].strip().split()
1217        assert L[0].strip().lower() == 'parameters'
1218
1219
1220        #Check asc file
1221        ascid = open(ascfile)
1222        lines = ascid.readlines()
1223        ascid.close()
1224
1225        L = lines[0].strip().split()
1226        assert L[0].strip().lower() == 'ncols'
1227        assert L[1].strip().lower() == '5'
1228
1229        L = lines[1].strip().split()
1230        assert L[0].strip().lower() == 'nrows'
1231        assert L[1].strip().lower() == '5'
1232
1233        L = lines[2].strip().split()
1234        assert L[0].strip().lower() == 'xllcorner'
1235        assert allclose(float(L[1].strip().lower()), 308500)
1236
1237        L = lines[3].strip().split()
1238        assert L[0].strip().lower() == 'yllcorner'
1239        assert allclose(float(L[1].strip().lower()), 6189000)
1240
1241        L = lines[4].strip().split()
1242        assert L[0].strip().lower() == 'cellsize'
1243        assert allclose(float(L[1].strip().lower()), cellsize)
1244
1245        L = lines[5].strip().split()
1246        assert L[0].strip() == 'NODATA_value'
1247        assert L[1].strip().lower() == '-9999'
1248
1249        #Check grid values
1250        for j in range(5):
1251            L = lines[6+j].strip().split()
1252            y = (4-j) * cellsize
1253            for i in range(5):
1254                assert allclose(float(L[i]), -i*cellsize - y)
1255
1256
1257        fid.close()
1258
1259        #Cleanup
1260        os.remove(prjfile)
1261        os.remove(ascfile)
1262        os.remove(swwfile)
1263
1264
1265
1266    def test_sww2dem_larger(self):
1267        """Test that sww information can be converted correctly to asc/prj
1268        format readable by e.g. ArcView. Here:
1269
1270        ncols         11
1271        nrows         11
1272        xllcorner     308500
1273        yllcorner     6189000
1274        cellsize      10.000000
1275        NODATA_value  -9999
1276        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1277         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1278         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1279         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
1280         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
1281         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
1282         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
1283         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
1284         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
1285         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
1286           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
1287
1288        """
1289
1290        import time, os
1291        from Numeric import array, zeros, allclose, Float, concatenate
1292        from Scientific.IO.NetCDF import NetCDFFile
1293
1294        #Setup
1295
1296        from mesh_factory import rectangular
1297
1298        #Create basic mesh (100m x 100m)
1299        points, vertices, boundary = rectangular(2, 2, 100, 100)
1300
1301        #Create shallow water domain
1302        domain = Domain(points, vertices, boundary)
1303        domain.default_order = 2
1304
1305        domain.filename = 'datatest'
1306
1307        prjfile = domain.filename + '_elevation.prj'
1308        ascfile = domain.filename + '_elevation.asc'
1309        swwfile = domain.filename + '.sww'
1310
1311        domain.set_datadir('.')
1312        domain.format = 'sww'
1313        domain.smooth = True
1314        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1315
1316        #
1317        domain.set_quantity('elevation', lambda x,y: -x-y)
1318        domain.set_quantity('stage', 0)
1319
1320        B = Transmissive_boundary(domain)
1321        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1322
1323
1324        #
1325        sww = get_dataobject(domain)
1326        sww.store_connectivity()
1327        sww.store_timestep('stage')
1328
1329        domain.evolve_to_end(finaltime = 0.01)
1330        sww.store_timestep('stage')
1331
1332        cellsize = 10  #10m grid
1333
1334
1335        #Check contents
1336        #Get NetCDF
1337
1338        fid = NetCDFFile(sww.filename, 'r')
1339
1340        # Get the variables
1341        x = fid.variables['x'][:]
1342        y = fid.variables['y'][:]
1343        z = fid.variables['elevation'][:]
1344        time = fid.variables['time'][:]
1345        stage = fid.variables['stage'][:]
1346
1347
1348        #Export to ascii/prj files
1349        sww2dem(domain.filename,
1350                quantity = 'elevation',
1351                cellsize = cellsize,
1352                verbose = False,
1353                format = 'asc')
1354
1355
1356        #Check prj (meta data)
1357        prjid = open(prjfile)
1358        lines = prjid.readlines()
1359        prjid.close()
1360
1361        L = lines[0].strip().split()
1362        assert L[0].strip().lower() == 'projection'
1363        assert L[1].strip().lower() == 'utm'
1364
1365        L = lines[1].strip().split()
1366        assert L[0].strip().lower() == 'zone'
1367        assert L[1].strip().lower() == '56'
1368
1369        L = lines[2].strip().split()
1370        assert L[0].strip().lower() == 'datum'
1371        assert L[1].strip().lower() == 'wgs84'
1372
1373        L = lines[3].strip().split()
1374        assert L[0].strip().lower() == 'zunits'
1375        assert L[1].strip().lower() == 'no'
1376
1377        L = lines[4].strip().split()
1378        assert L[0].strip().lower() == 'units'
1379        assert L[1].strip().lower() == 'meters'
1380
1381        L = lines[5].strip().split()
1382        assert L[0].strip().lower() == 'spheroid'
1383        assert L[1].strip().lower() == 'wgs84'
1384
1385        L = lines[6].strip().split()
1386        assert L[0].strip().lower() == 'xshift'
1387        assert L[1].strip().lower() == '500000'
1388
1389        L = lines[7].strip().split()
1390        assert L[0].strip().lower() == 'yshift'
1391        assert L[1].strip().lower() == '10000000'
1392
1393        L = lines[8].strip().split()
1394        assert L[0].strip().lower() == 'parameters'
1395
1396
1397        #Check asc file
1398        ascid = open(ascfile)
1399        lines = ascid.readlines()
1400        ascid.close()
1401
1402        L = lines[0].strip().split()
1403        assert L[0].strip().lower() == 'ncols'
1404        assert L[1].strip().lower() == '11'
1405
1406        L = lines[1].strip().split()
1407        assert L[0].strip().lower() == 'nrows'
1408        assert L[1].strip().lower() == '11'
1409
1410        L = lines[2].strip().split()
1411        assert L[0].strip().lower() == 'xllcorner'
1412        assert allclose(float(L[1].strip().lower()), 308500)
1413
1414        L = lines[3].strip().split()
1415        assert L[0].strip().lower() == 'yllcorner'
1416        assert allclose(float(L[1].strip().lower()), 6189000)
1417
1418        L = lines[4].strip().split()
1419        assert L[0].strip().lower() == 'cellsize'
1420        assert allclose(float(L[1].strip().lower()), cellsize)
1421
1422        L = lines[5].strip().split()
1423        assert L[0].strip() == 'NODATA_value'
1424        assert L[1].strip().lower() == '-9999'
1425
1426        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1427        for i, line in enumerate(lines[6:]):
1428            for j, value in enumerate( line.split() ):
1429                #assert allclose(float(value), -(10-i+j)*cellsize)
1430                assert float(value) == -(10-i+j)*cellsize
1431
1432
1433        fid.close()
1434
1435        #Cleanup
1436        os.remove(prjfile)
1437        os.remove(ascfile)
1438        os.remove(swwfile)
1439
1440
1441
1442    def test_sww2dem_boundingbox(self):
1443        """Test that sww information can be converted correctly to asc/prj
1444        format readable by e.g. ArcView.
1445        This will test that mesh can be restricted by bounding box
1446
1447        Original extent is 100m x 100m:
1448
1449        Eastings:   308500 -  308600
1450        Northings: 6189000 - 6189100
1451
1452        Bounding box changes this to the 50m x 50m square defined by
1453
1454        Eastings:   308530 -  308570
1455        Northings: 6189050 - 6189100
1456
1457        The cropped values should be
1458
1459         -130 -140 -150 -160 -170
1460         -120 -130 -140 -150 -160
1461         -110 -120 -130 -140 -150
1462         -100 -110 -120 -130 -140
1463          -90 -100 -110 -120 -130
1464          -80  -90 -100 -110 -120
1465
1466        and the new lower reference point should be
1467        Eastings:   308530
1468        Northings: 6189050
1469
1470        Original dataset is the same as in test_sww2dem_larger()
1471
1472        """
1473
1474        import time, os
1475        from Numeric import array, zeros, allclose, Float, concatenate
1476        from Scientific.IO.NetCDF import NetCDFFile
1477
1478        #Setup
1479
1480        from mesh_factory import rectangular
1481
1482        #Create basic mesh (100m x 100m)
1483        points, vertices, boundary = rectangular(2, 2, 100, 100)
1484
1485        #Create shallow water domain
1486        domain = Domain(points, vertices, boundary)
1487        domain.default_order = 2
1488
1489        domain.filename = 'datatest'
1490
1491        prjfile = domain.filename + '_elevation.prj'
1492        ascfile = domain.filename + '_elevation.asc'
1493        swwfile = domain.filename + '.sww'
1494
1495        domain.set_datadir('.')
1496        domain.format = 'sww'
1497        domain.smooth = True
1498        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1499
1500        #
1501        domain.set_quantity('elevation', lambda x,y: -x-y)
1502        domain.set_quantity('stage', 0)
1503
1504        B = Transmissive_boundary(domain)
1505        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1506
1507
1508        #
1509        sww = get_dataobject(domain)
1510        sww.store_connectivity()
1511        sww.store_timestep('stage')
1512
1513        domain.evolve_to_end(finaltime = 0.01)
1514        sww.store_timestep('stage')
1515
1516        cellsize = 10  #10m grid
1517
1518
1519        #Check contents
1520        #Get NetCDF
1521
1522        fid = NetCDFFile(sww.filename, 'r')
1523
1524        # Get the variables
1525        x = fid.variables['x'][:]
1526        y = fid.variables['y'][:]
1527        z = fid.variables['elevation'][:]
1528        time = fid.variables['time'][:]
1529        stage = fid.variables['stage'][:]
1530
1531
1532        #Export to ascii/prj files
1533        sww2dem(domain.filename,
1534                quantity = 'elevation',
1535                cellsize = cellsize,
1536                easting_min = 308530,
1537                easting_max = 308570,
1538                northing_min = 6189050,
1539                northing_max = 6189100,
1540                verbose = False,
1541                format = 'asc')
1542
1543        fid.close()
1544
1545
1546        #Check prj (meta data)
1547        prjid = open(prjfile)
1548        lines = prjid.readlines()
1549        prjid.close()
1550
1551        L = lines[0].strip().split()
1552        assert L[0].strip().lower() == 'projection'
1553        assert L[1].strip().lower() == 'utm'
1554
1555        L = lines[1].strip().split()
1556        assert L[0].strip().lower() == 'zone'
1557        assert L[1].strip().lower() == '56'
1558
1559        L = lines[2].strip().split()
1560        assert L[0].strip().lower() == 'datum'
1561        assert L[1].strip().lower() == 'wgs84'
1562
1563        L = lines[3].strip().split()
1564        assert L[0].strip().lower() == 'zunits'
1565        assert L[1].strip().lower() == 'no'
1566
1567        L = lines[4].strip().split()
1568        assert L[0].strip().lower() == 'units'
1569        assert L[1].strip().lower() == 'meters'
1570
1571        L = lines[5].strip().split()
1572        assert L[0].strip().lower() == 'spheroid'
1573        assert L[1].strip().lower() == 'wgs84'
1574
1575        L = lines[6].strip().split()
1576        assert L[0].strip().lower() == 'xshift'
1577        assert L[1].strip().lower() == '500000'
1578
1579        L = lines[7].strip().split()
1580        assert L[0].strip().lower() == 'yshift'
1581        assert L[1].strip().lower() == '10000000'
1582
1583        L = lines[8].strip().split()
1584        assert L[0].strip().lower() == 'parameters'
1585
1586
1587        #Check asc file
1588        ascid = open(ascfile)
1589        lines = ascid.readlines()
1590        ascid.close()
1591
1592        L = lines[0].strip().split()
1593        assert L[0].strip().lower() == 'ncols'
1594        assert L[1].strip().lower() == '5'
1595
1596        L = lines[1].strip().split()
1597        assert L[0].strip().lower() == 'nrows'
1598        assert L[1].strip().lower() == '6'
1599
1600        L = lines[2].strip().split()
1601        assert L[0].strip().lower() == 'xllcorner'
1602        assert allclose(float(L[1].strip().lower()), 308530)
1603
1604        L = lines[3].strip().split()
1605        assert L[0].strip().lower() == 'yllcorner'
1606        assert allclose(float(L[1].strip().lower()), 6189050)
1607
1608        L = lines[4].strip().split()
1609        assert L[0].strip().lower() == 'cellsize'
1610        assert allclose(float(L[1].strip().lower()), cellsize)
1611
1612        L = lines[5].strip().split()
1613        assert L[0].strip() == 'NODATA_value'
1614        assert L[1].strip().lower() == '-9999'
1615
1616        #Check grid values
1617        for i, line in enumerate(lines[6:]):
1618            for j, value in enumerate( line.split() ):
1619                #assert float(value) == -(10-i+j)*cellsize
1620                assert float(value) == -(10-i+j+3)*cellsize
1621
1622
1623
1624        #Cleanup
1625        os.remove(prjfile)
1626        os.remove(ascfile)
1627        os.remove(swwfile)
1628
1629
1630
1631    def test_sww2dem_asc_stage_reduction(self):
1632        """Test that sww information can be converted correctly to asc/prj
1633        format readable by e.g. ArcView
1634
1635        This tests the reduction of quantity stage using min
1636        """
1637
1638        import time, os
1639        from Numeric import array, zeros, allclose, Float, concatenate
1640        from Scientific.IO.NetCDF import NetCDFFile
1641
1642        #Setup
1643        self.domain.filename = 'datatest'
1644
1645        prjfile = self.domain.filename + '_stage.prj'
1646        ascfile = self.domain.filename + '_stage.asc'
1647        swwfile = self.domain.filename + '.sww'
1648
1649        self.domain.set_datadir('.')
1650        self.domain.format = 'sww'
1651        self.domain.smooth = True
1652        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1653
1654        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1655
1656
1657        sww = get_dataobject(self.domain)
1658        sww.store_connectivity()
1659        sww.store_timestep('stage')
1660
1661        self.domain.evolve_to_end(finaltime = 0.01)
1662        sww.store_timestep('stage')
1663
1664        cellsize = 0.25
1665        #Check contents
1666        #Get NetCDF
1667
1668        fid = NetCDFFile(sww.filename, 'r')
1669
1670        # Get the variables
1671        x = fid.variables['x'][:]
1672        y = fid.variables['y'][:]
1673        z = fid.variables['elevation'][:]
1674        time = fid.variables['time'][:]
1675        stage = fid.variables['stage'][:]
1676
1677
1678        #Export to ascii/prj files
1679        sww2dem(self.domain.filename,
1680                quantity = 'stage',
1681                cellsize = cellsize,
1682                reduction = min,
1683                format = 'asc')
1684
1685
1686        #Check asc file
1687        ascid = open(ascfile)
1688        lines = ascid.readlines()
1689        ascid.close()
1690
1691        L = lines[0].strip().split()
1692        assert L[0].strip().lower() == 'ncols'
1693        assert L[1].strip().lower() == '5'
1694
1695        L = lines[1].strip().split()
1696        assert L[0].strip().lower() == 'nrows'
1697        assert L[1].strip().lower() == '5'
1698
1699        L = lines[2].strip().split()
1700        assert L[0].strip().lower() == 'xllcorner'
1701        assert allclose(float(L[1].strip().lower()), 308500)
1702
1703        L = lines[3].strip().split()
1704        assert L[0].strip().lower() == 'yllcorner'
1705        assert allclose(float(L[1].strip().lower()), 6189000)
1706
1707        L = lines[4].strip().split()
1708        assert L[0].strip().lower() == 'cellsize'
1709        assert allclose(float(L[1].strip().lower()), cellsize)
1710
1711        L = lines[5].strip().split()
1712        assert L[0].strip() == 'NODATA_value'
1713        assert L[1].strip().lower() == '-9999'
1714
1715
1716        #Check grid values (where applicable)
1717        for j in range(5):
1718            if j%2 == 0:
1719                L = lines[6+j].strip().split()
1720                jj = 4-j
1721                for i in range(5):
1722                    if i%2 == 0:
1723                        index = jj/2 + i/2*3
1724                        val0 = stage[0,index]
1725                        val1 = stage[1,index]
1726
1727                        #print i, j, index, ':', L[i], val0, val1
1728                        assert allclose(float(L[i]), min(val0, val1))
1729
1730
1731        fid.close()
1732
1733        #Cleanup
1734        os.remove(prjfile)
1735        os.remove(ascfile)
1736        #os.remove(swwfile)
1737
1738
1739
1740    def test_sww2dem_asc_derived_quantity(self):
1741        """Test that sww information can be converted correctly to asc/prj
1742        format readable by e.g. ArcView
1743
1744        This tests the use of derived quantities
1745        """
1746
1747        import time, os
1748        from Numeric import array, zeros, allclose, Float, concatenate
1749        from Scientific.IO.NetCDF import NetCDFFile
1750
1751        #Setup
1752        self.domain.filename = 'datatest'
1753
1754        prjfile = self.domain.filename + '_depth.prj'
1755        ascfile = self.domain.filename + '_depth.asc'
1756        swwfile = self.domain.filename + '.sww'
1757
1758        self.domain.set_datadir('.')
1759        self.domain.format = 'sww'
1760        self.domain.smooth = True
1761        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1762        self.domain.set_quantity('stage', 0.0)
1763
1764        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1765
1766
1767        sww = get_dataobject(self.domain)
1768        sww.store_connectivity()
1769        sww.store_timestep('stage')
1770
1771        self.domain.evolve_to_end(finaltime = 0.01)
1772        sww.store_timestep('stage')
1773
1774        cellsize = 0.25
1775        #Check contents
1776        #Get NetCDF
1777
1778        fid = NetCDFFile(sww.filename, 'r')
1779
1780        # Get the variables
1781        x = fid.variables['x'][:]
1782        y = fid.variables['y'][:]
1783        z = fid.variables['elevation'][:]
1784        time = fid.variables['time'][:]
1785        stage = fid.variables['stage'][:]
1786
1787
1788        #Export to ascii/prj files
1789        sww2dem(self.domain.filename,
1790                basename_out = 'datatest_depth',
1791                quantity = 'stage - elevation',
1792                cellsize = cellsize,
1793                reduction = min,
1794                format = 'asc',
1795                verbose = False)
1796
1797
1798        #Check asc file
1799        ascid = open(ascfile)
1800        lines = ascid.readlines()
1801        ascid.close()
1802
1803        L = lines[0].strip().split()
1804        assert L[0].strip().lower() == 'ncols'
1805        assert L[1].strip().lower() == '5'
1806
1807        L = lines[1].strip().split()
1808        assert L[0].strip().lower() == 'nrows'
1809        assert L[1].strip().lower() == '5'
1810
1811        L = lines[2].strip().split()
1812        assert L[0].strip().lower() == 'xllcorner'
1813        assert allclose(float(L[1].strip().lower()), 308500)
1814
1815        L = lines[3].strip().split()
1816        assert L[0].strip().lower() == 'yllcorner'
1817        assert allclose(float(L[1].strip().lower()), 6189000)
1818
1819        L = lines[4].strip().split()
1820        assert L[0].strip().lower() == 'cellsize'
1821        assert allclose(float(L[1].strip().lower()), cellsize)
1822
1823        L = lines[5].strip().split()
1824        assert L[0].strip() == 'NODATA_value'
1825        assert L[1].strip().lower() == '-9999'
1826
1827
1828        #Check grid values (where applicable)
1829        for j in range(5):
1830            if j%2 == 0:
1831                L = lines[6+j].strip().split()
1832                jj = 4-j
1833                for i in range(5):
1834                    if i%2 == 0:
1835                        index = jj/2 + i/2*3
1836                        val0 = stage[0,index] - z[index]
1837                        val1 = stage[1,index] - z[index]
1838
1839                        #print i, j, index, ':', L[i], val0, val1
1840                        assert allclose(float(L[i]), min(val0, val1))
1841
1842
1843        fid.close()
1844
1845        #Cleanup
1846        os.remove(prjfile)
1847        os.remove(ascfile)
1848        #os.remove(swwfile)
1849
1850
1851
1852
1853
1854    def test_sww2dem_asc_missing_points(self):
1855        """Test that sww information can be converted correctly to asc/prj
1856        format readable by e.g. ArcView
1857
1858        This test includes the writing of missing values
1859        """
1860
1861        import time, os
1862        from Numeric import array, zeros, allclose, Float, concatenate
1863        from Scientific.IO.NetCDF import NetCDFFile
1864
1865        #Setup mesh not coinciding with rectangle.
1866        #This will cause missing values to occur in gridded data
1867
1868
1869        points = [                        [1.0, 1.0],
1870                              [0.5, 0.5], [1.0, 0.5],
1871                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1872
1873        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1874
1875        #Create shallow water domain
1876        domain = Domain(points, vertices)
1877        domain.default_order=2
1878
1879
1880        #Set some field values
1881        domain.set_quantity('elevation', lambda x,y: -x-y)
1882        domain.set_quantity('friction', 0.03)
1883
1884
1885        ######################
1886        # Boundary conditions
1887        B = Transmissive_boundary(domain)
1888        domain.set_boundary( {'exterior': B} )
1889
1890
1891        ######################
1892        #Initial condition - with jumps
1893
1894        bed = domain.quantities['elevation'].vertex_values
1895        stage = zeros(bed.shape, Float)
1896
1897        h = 0.3
1898        for i in range(stage.shape[0]):
1899            if i % 2 == 0:
1900                stage[i,:] = bed[i,:] + h
1901            else:
1902                stage[i,:] = bed[i,:]
1903
1904        domain.set_quantity('stage', stage)
1905        domain.distribute_to_vertices_and_edges()
1906
1907        domain.filename = 'datatest'
1908
1909        prjfile = domain.filename + '_elevation.prj'
1910        ascfile = domain.filename + '_elevation.asc'
1911        swwfile = domain.filename + '.sww'
1912
1913        domain.set_datadir('.')
1914        domain.format = 'sww'
1915        domain.smooth = True
1916
1917        domain.geo_reference = Geo_reference(56,308500,6189000)
1918
1919        sww = get_dataobject(domain)
1920        sww.store_connectivity()
1921        sww.store_timestep('stage')
1922
1923        cellsize = 0.25
1924        #Check contents
1925        #Get NetCDF
1926
1927        fid = NetCDFFile(swwfile, 'r')
1928
1929        # Get the variables
1930        x = fid.variables['x'][:]
1931        y = fid.variables['y'][:]
1932        z = fid.variables['elevation'][:]
1933        time = fid.variables['time'][:]
1934
1935        try:
1936            geo_reference = Geo_reference(NetCDFObject=fid)
1937        except AttributeError, e:
1938            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1939
1940        #Export to ascii/prj files
1941        sww2dem(domain.filename,
1942                quantity = 'elevation',
1943                cellsize = cellsize,
1944                verbose = False,
1945                format = 'asc')
1946
1947
1948        #Check asc file
1949        ascid = open(ascfile)
1950        lines = ascid.readlines()
1951        ascid.close()
1952
1953        L = lines[0].strip().split()
1954        assert L[0].strip().lower() == 'ncols'
1955        assert L[1].strip().lower() == '5'
1956
1957        L = lines[1].strip().split()
1958        assert L[0].strip().lower() == 'nrows'
1959        assert L[1].strip().lower() == '5'
1960
1961        L = lines[2].strip().split()
1962        assert L[0].strip().lower() == 'xllcorner'
1963        assert allclose(float(L[1].strip().lower()), 308500)
1964
1965        L = lines[3].strip().split()
1966        assert L[0].strip().lower() == 'yllcorner'
1967        assert allclose(float(L[1].strip().lower()), 6189000)
1968
1969        L = lines[4].strip().split()
1970        assert L[0].strip().lower() == 'cellsize'
1971        assert allclose(float(L[1].strip().lower()), cellsize)
1972
1973        L = lines[5].strip().split()
1974        assert L[0].strip() == 'NODATA_value'
1975        assert L[1].strip().lower() == '-9999'
1976
1977        #Check grid values
1978        for j in range(5):
1979            L = lines[6+j].strip().split()
1980            assert len(L) == 5
1981            y = (4-j) * cellsize
1982
1983            for i in range(5):
1984                #print i
1985                if i+j >= 4:
1986                    assert allclose(float(L[i]), -i*cellsize - y)
1987                else:
1988                    #Missing values
1989                    assert allclose(float(L[i]), -9999)
1990
1991
1992
1993        fid.close()
1994
1995        #Cleanup
1996        os.remove(prjfile)
1997        os.remove(ascfile)
1998        os.remove(swwfile)
1999
2000    def test_sww2ers_simple(self):
2001        """Test that sww information can be converted correctly to asc/prj
2002        format readable by e.g. ArcView
2003        """
2004
2005        import time, os
2006        from Numeric import array, zeros, allclose, Float, concatenate
2007        from Scientific.IO.NetCDF import NetCDFFile
2008
2009
2010        NODATA_value = 1758323
2011
2012        #Setup
2013        self.domain.filename = 'datatest'
2014
2015        headerfile = self.domain.filename + '.ers'
2016        swwfile = self.domain.filename + '.sww'
2017
2018        self.domain.set_datadir('.')
2019        self.domain.format = 'sww'
2020        self.domain.smooth = True
2021        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2022
2023        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2024
2025        sww = get_dataobject(self.domain)
2026        sww.store_connectivity()
2027        sww.store_timestep('stage')
2028
2029        self.domain.evolve_to_end(finaltime = 0.01)
2030        sww.store_timestep('stage')
2031
2032        cellsize = 0.25
2033        #Check contents
2034        #Get NetCDF
2035
2036        fid = NetCDFFile(sww.filename, 'r')
2037
2038        # Get the variables
2039        x = fid.variables['x'][:]
2040        y = fid.variables['y'][:]
2041        z = fid.variables['elevation'][:]
2042        time = fid.variables['time'][:]
2043        stage = fid.variables['stage'][:]
2044
2045
2046        #Export to ers files
2047        sww2dem(self.domain.filename,
2048                quantity = 'elevation',
2049                cellsize = cellsize,
2050                NODATA_value = NODATA_value,
2051                verbose = False,
2052                format = 'ers')
2053
2054        #Check header data
2055        from ermapper_grids import read_ermapper_header, read_ermapper_data
2056
2057        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
2058        #print header
2059        assert header['projection'].lower() == '"utm-56"'
2060        assert header['datum'].lower() == '"wgs84"'
2061        assert header['units'].lower() == '"meters"'
2062        assert header['value'].lower() == '"elevation"'
2063        assert header['xdimension'] == '0.25'
2064        assert header['ydimension'] == '0.25'
2065        assert float(header['eastings']) == 308500.0   #xllcorner
2066        assert float(header['northings']) == 6189000.0 #yllcorner
2067        assert int(header['nroflines']) == 5
2068        assert int(header['nrofcellsperline']) == 5
2069        assert int(header['nullcellvalue']) == NODATA_value
2070        #FIXME - there is more in the header
2071
2072
2073        #Check grid data
2074        grid = read_ermapper_data(self.domain.filename + '_elevation')
2075
2076        #FIXME (Ole): Why is this the desired reference grid for -x-y?
2077        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
2078                    -1,    -1.25, -1.5,  -1.75, -2.0,
2079                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
2080                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
2081                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
2082
2083
2084        #print grid
2085        assert allclose(grid, ref_grid)
2086
2087        fid.close()
2088
2089        #Cleanup
2090        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
2091        #Done (Ole) - it was because sww2ers didn't close it's sww file
2092        os.remove(sww.filename)
2093        os.remove(self.domain.filename + '_elevation')
2094        os.remove(self.domain.filename + '_elevation.ers')
2095
2096
2097
2098    def test_sww2pts_centroids(self):
2099        """Test that sww information can be converted correctly to pts data at specified coordinates
2100        - in this case, the centroids.
2101        """
2102
2103        import time, os
2104        from Numeric import array, zeros, allclose, Float, concatenate, NewAxis
2105        from Scientific.IO.NetCDF import NetCDFFile
2106        from anuga.geospatial_data.geospatial_data import Geospatial_data
2107
2108        # Used for points that lie outside mesh
2109        NODATA_value = 1758323
2110
2111        # Setup
2112        self.domain.filename = 'datatest'
2113
2114        ptsfile = self.domain.filename + '_elevation.pts'
2115        swwfile = self.domain.filename + '.sww'
2116
2117        self.domain.set_datadir('.')
2118        self.domain.format = 'sww'
2119        self.smooth = True #self.set_store_vertices_uniquely(False)
2120        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2121
2122        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2123
2124        sww = get_dataobject(self.domain)
2125        sww.store_connectivity()
2126        sww.store_timestep('stage')
2127
2128        self.domain.evolve_to_end(finaltime = 0.01)
2129        sww.store_timestep('stage')
2130
2131        # Check contents in NetCDF
2132        fid = NetCDFFile(sww.filename, 'r')
2133
2134        # Get the variables
2135        x = fid.variables['x'][:]
2136        y = fid.variables['y'][:]
2137        elevation = fid.variables['elevation'][:]
2138        time = fid.variables['time'][:]
2139        stage = fid.variables['stage'][:]
2140
2141        volumes = fid.variables['volumes'][:]
2142
2143
2144        # Invoke interpolation for vertex points       
2145        points = concatenate( (x[:,NewAxis],y[:,NewAxis]), axis=1 )
2146        sww2pts(self.domain.filename,
2147                quantity = 'elevation',
2148                data_points = points,
2149                NODATA_value = NODATA_value,
2150                verbose = False)
2151        ref_point_values = elevation
2152        point_values = Geospatial_data(ptsfile).get_attributes()
2153        #print 'P', point_values
2154        #print 'Ref', ref_point_values       
2155        assert allclose(point_values, ref_point_values)       
2156
2157
2158
2159        # Invoke interpolation for centroids
2160        points = self.domain.get_centroid_coordinates()
2161        #print points
2162        sww2pts(self.domain.filename,
2163                quantity = 'elevation',
2164                data_points = points,
2165                NODATA_value = NODATA_value,
2166                verbose = False)
2167        ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids
2168
2169       
2170        point_values = Geospatial_data(ptsfile).get_attributes()
2171        #print 'P', point_values
2172        #print 'Ref', ref_point_values       
2173        assert allclose(point_values, ref_point_values)       
2174
2175
2176
2177        fid.close()
2178
2179        #Cleanup
2180        os.remove(sww.filename)
2181        os.remove(ptsfile)
2182
2183
2184
2185
2186    def test_ferret2sww1(self):
2187        """Test that georeferencing etc works when converting from
2188        ferret format (lat/lon) to sww format (UTM)
2189        """
2190        from Scientific.IO.NetCDF import NetCDFFile
2191        import os, sys
2192
2193        #The test file has
2194        # LON = 150.66667, 150.83334, 151, 151.16667
2195        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2196        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2197        #
2198        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2199        # Fourth value (index==3) is -6.50198 cm
2200
2201
2202
2203        #Read
2204        from anuga.coordinate_transforms.redfearn import redfearn
2205        #fid = NetCDFFile(self.test_MOST_file)
2206        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2207        first_value = fid.variables['HA'][:][0,0,0]
2208        fourth_value = fid.variables['HA'][:][0,0,3]
2209        fid.close()
2210
2211
2212        #Call conversion (with zero origin)
2213        #ferret2sww('small', verbose=False,
2214        #           origin = (56, 0, 0))
2215        ferret2sww(self.test_MOST_file, verbose=False,
2216                   origin = (56, 0, 0))
2217
2218        #Work out the UTM coordinates for first point
2219        zone, e, n = redfearn(-34.5, 150.66667)
2220        #print zone, e, n
2221
2222        #Read output file 'small.sww'
2223        #fid = NetCDFFile('small.sww')
2224        fid = NetCDFFile(self.test_MOST_file + '.sww')
2225
2226        x = fid.variables['x'][:]
2227        y = fid.variables['y'][:]
2228
2229        #Check that first coordinate is correctly represented
2230        assert allclose(x[0], e)
2231        assert allclose(y[0], n)
2232
2233        #Check first value
2234        stage = fid.variables['stage'][:]
2235        xmomentum = fid.variables['xmomentum'][:]
2236        ymomentum = fid.variables['ymomentum'][:]
2237
2238        #print ymomentum
2239
2240        assert allclose(stage[0,0], first_value/100)  #Meters
2241
2242        #Check fourth value
2243        assert allclose(stage[0,3], fourth_value/100)  #Meters
2244
2245        fid.close()
2246
2247        #Cleanup
2248        import os
2249        os.remove(self.test_MOST_file + '.sww')
2250
2251
2252    def test_ferret2sww_2(self):
2253        """Test that georeferencing etc works when converting from
2254        ferret format (lat/lon) to sww format (UTM)
2255        """
2256        from Scientific.IO.NetCDF import NetCDFFile
2257
2258        #The test file has
2259        # LON = 150.66667, 150.83334, 151, 151.16667
2260        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2261        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2262        #
2263        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2264        # Fourth value (index==3) is -6.50198 cm
2265
2266
2267        from anuga.coordinate_transforms.redfearn import redfearn
2268
2269        #fid = NetCDFFile('small_ha.nc')
2270        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2271
2272        #Pick a coordinate and a value
2273
2274        time_index = 1
2275        lat_index = 0
2276        lon_index = 2
2277
2278        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
2279        test_time = fid.variables['TIME'][:][time_index]
2280        test_lat = fid.variables['LAT'][:][lat_index]
2281        test_lon = fid.variables['LON'][:][lon_index]
2282
2283        linear_point_index = lat_index*4 + lon_index
2284        fid.close()
2285
2286        #Call conversion (with zero origin)
2287        ferret2sww(self.test_MOST_file, verbose=False,
2288                   origin = (56, 0, 0))
2289
2290
2291        #Work out the UTM coordinates for test point
2292        zone, e, n = redfearn(test_lat, test_lon)
2293
2294        #Read output file 'small.sww'
2295        fid = NetCDFFile(self.test_MOST_file + '.sww')
2296
2297        x = fid.variables['x'][:]
2298        y = fid.variables['y'][:]
2299
2300        #Check that test coordinate is correctly represented
2301        assert allclose(x[linear_point_index], e)
2302        assert allclose(y[linear_point_index], n)
2303
2304        #Check test value
2305        stage = fid.variables['stage'][:]
2306
2307        assert allclose(stage[time_index, linear_point_index], test_value/100)
2308
2309        fid.close()
2310
2311        #Cleanup
2312        import os
2313        os.remove(self.test_MOST_file + '.sww')
2314
2315
2316
2317    def test_ferret2sww3(self):
2318        """Elevation included
2319        """
2320        from Scientific.IO.NetCDF import NetCDFFile
2321
2322        #The test file has
2323        # LON = 150.66667, 150.83334, 151, 151.16667
2324        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2325        # ELEVATION = [-1 -2 -3 -4
2326        #             -5 -6 -7 -8
2327        #              ...
2328        #              ...      -16]
2329        # where the top left corner is -1m,
2330        # and the ll corner is -13.0m
2331        #
2332        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2333        # Fourth value (index==3) is -6.50198 cm
2334
2335        from anuga.coordinate_transforms.redfearn import redfearn
2336        import os
2337        fid1 = NetCDFFile('test_ha.nc','w')
2338        fid2 = NetCDFFile('test_ua.nc','w')
2339        fid3 = NetCDFFile('test_va.nc','w')
2340        fid4 = NetCDFFile('test_e.nc','w')
2341
2342        h1_list = [150.66667,150.83334,151.]
2343        h2_list = [-34.5,-34.33333]
2344
2345        long_name = 'LON'
2346        lat_name = 'LAT'
2347        time_name = 'TIME'
2348
2349        nx = 3
2350        ny = 2
2351
2352        for fid in [fid1,fid2,fid3]:
2353            fid.createDimension(long_name,nx)
2354            fid.createVariable(long_name,'d',(long_name,))
2355            fid.variables[long_name].point_spacing='uneven'
2356            fid.variables[long_name].units='degrees_east'
2357            fid.variables[long_name].assignValue(h1_list)
2358
2359            fid.createDimension(lat_name,ny)
2360            fid.createVariable(lat_name,'d',(lat_name,))
2361            fid.variables[lat_name].point_spacing='uneven'
2362            fid.variables[lat_name].units='degrees_north'
2363            fid.variables[lat_name].assignValue(h2_list)
2364
2365            fid.createDimension(time_name,2)
2366            fid.createVariable(time_name,'d',(time_name,))
2367            fid.variables[time_name].point_spacing='uneven'
2368            fid.variables[time_name].units='seconds'
2369            fid.variables[time_name].assignValue([0.,1.])
2370            if fid == fid3: break
2371
2372
2373        for fid in [fid4]:
2374            fid.createDimension(long_name,nx)
2375            fid.createVariable(long_name,'d',(long_name,))
2376            fid.variables[long_name].point_spacing='uneven'
2377            fid.variables[long_name].units='degrees_east'
2378            fid.variables[long_name].assignValue(h1_list)
2379
2380            fid.createDimension(lat_name,ny)
2381            fid.createVariable(lat_name,'d',(lat_name,))
2382            fid.variables[lat_name].point_spacing='uneven'
2383            fid.variables[lat_name].units='degrees_north'
2384            fid.variables[lat_name].assignValue(h2_list)
2385
2386        name = {}
2387        name[fid1]='HA'
2388        name[fid2]='UA'
2389        name[fid3]='VA'
2390        name[fid4]='ELEVATION'
2391
2392        units = {}
2393        units[fid1]='cm'
2394        units[fid2]='cm/s'
2395        units[fid3]='cm/s'
2396        units[fid4]='m'
2397
2398        values = {}
2399        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2400        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2401        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2402        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2403
2404        for fid in [fid1,fid2,fid3]:
2405          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
2406          fid.variables[name[fid]].point_spacing='uneven'
2407          fid.variables[name[fid]].units=units[fid]
2408          fid.variables[name[fid]].assignValue(values[fid])
2409          fid.variables[name[fid]].missing_value = -99999999.
2410          if fid == fid3: break
2411
2412        for fid in [fid4]:
2413            fid.createVariable(name[fid],'d',(lat_name,long_name))
2414            fid.variables[name[fid]].point_spacing='uneven'
2415            fid.variables[name[fid]].units=units[fid]
2416            fid.variables[name[fid]].assignValue(values[fid])
2417            fid.variables[name[fid]].missing_value = -99999999.
2418
2419
2420        fid1.sync(); fid1.close()
2421        fid2.sync(); fid2.close()
2422        fid3.sync(); fid3.close()
2423        fid4.sync(); fid4.close()
2424
2425        fid1 = NetCDFFile('test_ha.nc','r')
2426        fid2 = NetCDFFile('test_e.nc','r')
2427        fid3 = NetCDFFile('test_va.nc','r')
2428
2429
2430        first_amp = fid1.variables['HA'][:][0,0,0]
2431        third_amp = fid1.variables['HA'][:][0,0,2]
2432        first_elevation = fid2.variables['ELEVATION'][0,0]
2433        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2434        first_speed = fid3.variables['VA'][0,0,0]
2435        third_speed = fid3.variables['VA'][:][0,0,2]
2436
2437        fid1.close()
2438        fid2.close()
2439        fid3.close()
2440
2441        #Call conversion (with zero origin)
2442        ferret2sww('test', verbose=False,
2443                   origin = (56, 0, 0), inverted_bathymetry=False)
2444
2445        os.remove('test_va.nc')
2446        os.remove('test_ua.nc')
2447        os.remove('test_ha.nc')
2448        os.remove('test_e.nc')
2449
2450        #Read output file 'test.sww'
2451        fid = NetCDFFile('test.sww')
2452
2453
2454        #Check first value
2455        elevation = fid.variables['elevation'][:]
2456        stage = fid.variables['stage'][:]
2457        xmomentum = fid.variables['xmomentum'][:]
2458        ymomentum = fid.variables['ymomentum'][:]
2459
2460        #print ymomentum
2461        first_height = first_amp/100 - first_elevation
2462        third_height = third_amp/100 - third_elevation
2463        first_momentum=first_speed*first_height/100
2464        third_momentum=third_speed*third_height/100
2465
2466        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2467        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2468
2469        fid.close()
2470
2471        #Cleanup
2472        os.remove('test.sww')
2473
2474
2475
2476    def test_ferret2sww4(self):
2477        """Like previous but with augmented variable names as
2478        in files produced by ferret as opposed to MOST
2479        """
2480        from Scientific.IO.NetCDF import NetCDFFile
2481
2482        #The test file has
2483        # LON = 150.66667, 150.83334, 151, 151.16667
2484        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2485        # ELEVATION = [-1 -2 -3 -4
2486        #             -5 -6 -7 -8
2487        #              ...
2488        #              ...      -16]
2489        # where the top left corner is -1m,
2490        # and the ll corner is -13.0m
2491        #
2492        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2493        # Fourth value (index==3) is -6.50198 cm
2494
2495        from anuga.coordinate_transforms.redfearn import redfearn
2496        import os
2497        fid1 = NetCDFFile('test_ha.nc','w')
2498        fid2 = NetCDFFile('test_ua.nc','w')
2499        fid3 = NetCDFFile('test_va.nc','w')
2500        fid4 = NetCDFFile('test_e.nc','w')
2501
2502        h1_list = [150.66667,150.83334,151.]
2503        h2_list = [-34.5,-34.33333]
2504
2505#        long_name = 'LON961_1261'
2506#        lat_name = 'LAT481_841'
2507#        time_name = 'TIME1'
2508
2509        long_name = 'LON'
2510        lat_name = 'LAT'
2511        time_name = 'TIME'
2512
2513        nx = 3
2514        ny = 2
2515
2516        for fid in [fid1,fid2,fid3]:
2517            fid.createDimension(long_name,nx)
2518            fid.createVariable(long_name,'d',(long_name,))
2519            fid.variables[long_name].point_spacing='uneven'
2520            fid.variables[long_name].units='degrees_east'
2521            fid.variables[long_name].assignValue(h1_list)
2522
2523            fid.createDimension(lat_name,ny)
2524            fid.createVariable(lat_name,'d',(lat_name,))
2525            fid.variables[lat_name].point_spacing='uneven'
2526            fid.variables[lat_name].units='degrees_north'
2527            fid.variables[lat_name].assignValue(h2_list)
2528
2529            fid.createDimension(time_name,2)
2530            fid.createVariable(time_name,'d',(time_name,))
2531            fid.variables[time_name].point_spacing='uneven'
2532            fid.variables[time_name].units='seconds'
2533            fid.variables[time_name].assignValue([0.,1.])
2534            if fid == fid3: break
2535
2536
2537        for fid in [fid4]:
2538            fid.createDimension(long_name,nx)
2539            fid.createVariable(long_name,'d',(long_name,))
2540            fid.variables[long_name].point_spacing='uneven'
2541            fid.variables[long_name].units='degrees_east'
2542            fid.variables[long_name].assignValue(h1_list)
2543
2544            fid.createDimension(lat_name,ny)
2545            fid.createVariable(lat_name,'d',(lat_name,))
2546            fid.variables[lat_name].point_spacing='uneven'
2547            fid.variables[lat_name].units='degrees_north'
2548            fid.variables[lat_name].assignValue(h2_list)
2549
2550        name = {}
2551        name[fid1]='HA'
2552        name[fid2]='UA'
2553        name[fid3]='VA'
2554        name[fid4]='ELEVATION'
2555
2556        units = {}
2557        units[fid1]='cm'
2558        units[fid2]='cm/s'
2559        units[fid3]='cm/s'
2560        units[fid4]='m'
2561
2562        values = {}
2563        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2564        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2565        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2566        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2567
2568        for fid in [fid1,fid2,fid3]:
2569          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
2570          fid.variables[name[fid]].point_spacing='uneven'
2571          fid.variables[name[fid]].units=units[fid]
2572          fid.variables[name[fid]].assignValue(values[fid])
2573          fid.variables[name[fid]].missing_value = -99999999.
2574          if fid == fid3: break
2575
2576        for fid in [fid4]:
2577            fid.createVariable(name[fid],'d',(lat_name,long_name))
2578            fid.variables[name[fid]].point_spacing='uneven'
2579            fid.variables[name[fid]].units=units[fid]
2580            fid.variables[name[fid]].assignValue(values[fid])
2581            fid.variables[name[fid]].missing_value = -99999999.
2582
2583
2584        fid1.sync(); fid1.close()
2585        fid2.sync(); fid2.close()
2586        fid3.sync(); fid3.close()
2587        fid4.sync(); fid4.close()
2588
2589        fid1 = NetCDFFile('test_ha.nc','r')
2590        fid2 = NetCDFFile('test_e.nc','r')
2591        fid3 = NetCDFFile('test_va.nc','r')
2592
2593
2594        first_amp = fid1.variables['HA'][:][0,0,0]
2595        third_amp = fid1.variables['HA'][:][0,0,2]
2596        first_elevation = fid2.variables['ELEVATION'][0,0]
2597        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2598        first_speed = fid3.variables['VA'][0,0,0]
2599        third_speed = fid3.variables['VA'][:][0,0,2]
2600
2601        fid1.close()
2602        fid2.close()
2603        fid3.close()
2604
2605        #Call conversion (with zero origin)
2606        ferret2sww('test', verbose=False, origin = (56, 0, 0)
2607                   , inverted_bathymetry=False)
2608
2609        os.remove('test_va.nc')
2610        os.remove('test_ua.nc')
2611        os.remove('test_ha.nc')
2612        os.remove('test_e.nc')
2613
2614        #Read output file 'test.sww'
2615        fid = NetCDFFile('test.sww')
2616
2617
2618        #Check first value
2619        elevation = fid.variables['elevation'][:]
2620        stage = fid.variables['stage'][:]
2621        xmomentum = fid.variables['xmomentum'][:]
2622        ymomentum = fid.variables['ymomentum'][:]
2623
2624        #print ymomentum
2625        first_height = first_amp/100 - first_elevation
2626        third_height = third_amp/100 - third_elevation
2627        first_momentum=first_speed*first_height/100
2628        third_momentum=third_speed*third_height/100
2629
2630        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2631        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2632
2633        fid.close()
2634
2635        #Cleanup
2636        os.remove('test.sww')
2637
2638
2639
2640
2641    def test_ferret2sww_nz_origin(self):
2642        from Scientific.IO.NetCDF import NetCDFFile
2643        from anuga.coordinate_transforms.redfearn import redfearn
2644
2645        #Call conversion (with nonzero origin)
2646        ferret2sww(self.test_MOST_file, verbose=False,
2647                   origin = (56, 100000, 200000))
2648
2649
2650        #Work out the UTM coordinates for first point
2651        zone, e, n = redfearn(-34.5, 150.66667)
2652
2653        #Read output file 'small.sww'
2654        #fid = NetCDFFile('small.sww', 'r')
2655        fid = NetCDFFile(self.test_MOST_file + '.sww')
2656
2657        x = fid.variables['x'][:]
2658        y = fid.variables['y'][:]
2659
2660        #Check that first coordinate is correctly represented
2661        assert allclose(x[0], e-100000)
2662        assert allclose(y[0], n-200000)
2663
2664        fid.close()
2665
2666        #Cleanup
2667        os.remove(self.test_MOST_file + '.sww')
2668
2669
2670
2671    def test_sww_extent(self):
2672        """Not a test, rather a look at the sww format
2673        """
2674
2675        import time, os
2676        from Numeric import array, zeros, allclose, Float, concatenate
2677        from Scientific.IO.NetCDF import NetCDFFile
2678
2679        self.domain.filename = 'datatest' + str(id(self))
2680        self.domain.format = 'sww'
2681        self.domain.smooth = True
2682        self.domain.reduction = mean
2683        self.domain.set_datadir('.')
2684
2685
2686        sww = get_dataobject(self.domain)
2687        sww.store_connectivity()
2688        sww.store_timestep('stage')
2689        self.domain.time = 2.
2690
2691        #Modify stage at second timestep
2692        stage = self.domain.quantities['stage'].vertex_values
2693        self.domain.set_quantity('stage', stage/2)
2694
2695        sww.store_timestep('stage')
2696
2697        file_and_extension_name = self.domain.filename + ".sww"
2698        #print "file_and_extension_name",file_and_extension_name
2699        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2700               extent_sww(file_and_extension_name )
2701
2702        assert allclose(xmin, 0.0)
2703        assert allclose(xmax, 1.0)
2704        assert allclose(ymin, 0.0)
2705        assert allclose(ymax, 1.0)
2706        assert allclose(stagemin, -0.85)
2707        assert allclose(stagemax, 0.15)
2708
2709
2710        #Cleanup
2711        os.remove(sww.filename)
2712
2713
2714
2715    def test_sww2domain1(self):
2716        ################################################
2717        #Create a test domain, and evolve and save it.
2718        ################################################
2719        from mesh_factory import rectangular
2720        from Numeric import array
2721
2722        #Create basic mesh
2723
2724        yiel=0.01
2725        points, vertices, boundary = rectangular(10,10)
2726
2727        #Create shallow water domain
2728        domain = Domain(points, vertices, boundary)
2729        domain.geo_reference = Geo_reference(56,11,11)
2730        domain.smooth = False
2731        domain.visualise = False
2732        domain.store = True
2733        domain.filename = 'bedslope'
2734        domain.default_order=2
2735        #Bed-slope and friction
2736        domain.set_quantity('elevation', lambda x,y: -x/3)
2737        domain.set_quantity('friction', 0.1)
2738        # Boundary conditions
2739        from math import sin, pi
2740        Br = Reflective_boundary(domain)
2741        Bt = Transmissive_boundary(domain)
2742        Bd = Dirichlet_boundary([0.2,0.,0.])
2743        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2744
2745        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2746        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2747
2748        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2749        #Initial condition
2750        h = 0.05
2751        elevation = domain.quantities['elevation'].vertex_values
2752        domain.set_quantity('stage', elevation + h)
2753
2754        domain.check_integrity()
2755        #Evolution
2756        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2757            #domain.write_time()
2758            pass
2759
2760
2761        ##########################################
2762        #Import the example's file as a new domain
2763        ##########################################
2764        from data_manager import sww2domain
2765        from Numeric import allclose
2766        import os
2767
2768        filename = domain.datadir+os.sep+domain.filename+'.sww'
2769        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2770        #points, vertices, boundary = rectangular(15,15)
2771        #domain2.boundary = boundary
2772        ###################
2773        ##NOW TEST IT!!!
2774        ###################
2775
2776        #os.remove(domain.filename + '.sww')
2777        os.remove(filename)
2778
2779        bits = ['vertex_coordinates']
2780        for quantity in ['elevation']+domain.quantities_to_be_stored:
2781            bits.append('get_quantity("%s").get_integral()' %quantity)
2782            bits.append('get_quantity("%s").get_values()' %quantity)
2783
2784        for bit in bits:
2785            #print 'testing that domain.'+bit+' has been restored'
2786            #print bit
2787            #print 'done'
2788            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2789
2790        ######################################
2791        #Now evolve them both, just to be sure
2792        ######################################x
2793        visualise = False
2794        #visualise = True
2795        domain.visualise = visualise
2796        domain.time = 0.
2797        from time import sleep
2798
2799        final = .1
2800        domain.set_quantity('friction', 0.1)
2801        domain.store = False
2802        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2803
2804
2805        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2806            if visualise: sleep(1.)
2807            #domain.write_time()
2808            pass
2809
2810        final = final - (domain2.starttime-domain.starttime)
2811        #BUT since domain1 gets time hacked back to 0:
2812        final = final + (domain2.starttime-domain.starttime)
2813
2814        domain2.smooth = False
2815        domain2.visualise = visualise
2816        domain2.store = False
2817        domain2.default_order=2
2818        domain2.set_quantity('friction', 0.1)
2819        #Bed-slope and friction
2820        # Boundary conditions
2821        Bd2=Dirichlet_boundary([0.2,0.,0.])
2822        domain2.boundary = domain.boundary
2823        #print 'domain2.boundary'
2824        #print domain2.boundary
2825        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2826        #domain2.set_boundary({'exterior': Bd})
2827
2828        domain2.check_integrity()
2829
2830        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2831            if visualise: sleep(1.)
2832            #domain2.write_time()
2833            pass
2834
2835        ###################
2836        ##NOW TEST IT!!!
2837        ##################
2838
2839        bits = ['vertex_coordinates']
2840
2841        for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
2842            bits.append('get_quantity("%s").get_integral()' %quantity)
2843            bits.append('get_quantity("%s").get_values()' %quantity)
2844
2845        #print bits
2846        for bit in bits:
2847            #print bit
2848            #print eval('domain.'+bit)
2849            #print eval('domain2.'+bit)
2850           
2851            #print eval('domain.'+bit+'-domain2.'+bit)
2852            msg = 'Values in the two domains are different for ' + bit
2853            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),
2854                            rtol=1.e-5, atol=3.e-8), msg
2855
2856
2857    def test_sww2domain2(self):
2858        ##################################################################
2859        #Same as previous test, but this checks how NaNs are handled.
2860        ##################################################################
2861
2862
2863        from mesh_factory import rectangular
2864        from Numeric import array
2865
2866        #Create basic mesh
2867        points, vertices, boundary = rectangular(2,2)
2868
2869        #Create shallow water domain
2870        domain = Domain(points, vertices, boundary)
2871        domain.smooth = False
2872        domain.visualise = False
2873        domain.store = True
2874        domain.set_name('test_file')
2875        domain.set_datadir('.')
2876        domain.default_order=2
2877        domain.quantities_to_be_stored=['stage']
2878
2879        domain.set_quantity('elevation', lambda x,y: -x/3)
2880        domain.set_quantity('friction', 0.1)
2881
2882        from math import sin, pi
2883        Br = Reflective_boundary(domain)
2884        Bt = Transmissive_boundary(domain)
2885        Bd = Dirichlet_boundary([0.2,0.,0.])
2886        Bw = Time_boundary(domain=domain,
2887                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2888
2889        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2890
2891        h = 0.05
2892        elevation = domain.quantities['elevation'].vertex_values
2893        domain.set_quantity('stage', elevation + h)
2894
2895        domain.check_integrity()
2896
2897        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2898            pass
2899            #domain.write_time()
2900
2901
2902
2903        ##################################
2904        #Import the file as a new domain
2905        ##################################
2906        from data_manager import sww2domain
2907        from Numeric import allclose
2908        import os
2909
2910        filename = domain.datadir+os.sep+domain.filename+'.sww'
2911
2912        #Fail because NaNs are present
2913        try:
2914            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2915        except:
2916            #Now import it, filling NaNs to be 0
2917            filler = 0
2918            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2919
2920        #Clean up
2921        os.remove(filename)
2922
2923
2924        bits = [ 'geo_reference.get_xllcorner()',
2925                'geo_reference.get_yllcorner()',
2926                'vertex_coordinates']
2927
2928        for quantity in ['elevation']+domain.quantities_to_be_stored:
2929            bits.append('get_quantity("%s").get_integral()' %quantity)
2930            bits.append('get_quantity("%s").get_values()' %quantity)
2931
2932        for bit in bits:
2933        #    print 'testing that domain.'+bit+' has been restored'
2934            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2935
2936        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
2937        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
2938        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
2939        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
2940
2941
2942
2943    #def test_weed(self):
2944        from data_manager import weed
2945
2946        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2947        volumes1 = [[0,1,2],[3,4,5]]
2948        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2949        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2950
2951        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2952
2953        assert len(points2)==len(coordinates2)
2954        for i in range(len(coordinates2)):
2955            coordinate = tuple(coordinates2[i])
2956            assert points2.has_key(coordinate)
2957            points2[coordinate]=i
2958
2959        for triangle in volumes1:
2960            for coordinate in triangle:
2961                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2962                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2963
2964
2965    #FIXME This fails - smooth makes the comparism too hard for allclose
2966    def ztest_sww2domain3(self):
2967        ################################################
2968        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2969        ################################################
2970        from mesh_factory import rectangular
2971        from Numeric import array
2972        #Create basic mesh
2973
2974        yiel=0.01
2975        points, vertices, boundary = rectangular(10,10)
2976
2977        #Create shallow water domain
2978        domain = Domain(points, vertices, boundary)
2979        domain.geo_reference = Geo_reference(56,11,11)
2980        domain.smooth = True
2981        domain.visualise = False
2982        domain.store = True
2983        domain.filename = 'bedslope'
2984        domain.default_order=2
2985        #Bed-slope and friction
2986        domain.set_quantity('elevation', lambda x,y: -x/3)
2987        domain.set_quantity('friction', 0.1)
2988        # Boundary conditions
2989        from math import sin, pi
2990        Br = Reflective_boundary(domain)
2991        Bt = Transmissive_boundary(domain)
2992        Bd = Dirichlet_boundary([0.2,0.,0.])
2993        Bw = Time_boundary(domain=domain,
2994                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2995
2996        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2997
2998        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2999        #Initial condition
3000        h = 0.05
3001        elevation = domain.quantities['elevation'].vertex_values
3002        domain.set_quantity('stage', elevation + h)
3003
3004
3005        domain.check_integrity()
3006        #Evolution
3007        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3008        #    domain.write_time()
3009            pass
3010
3011
3012        ##########################################
3013        #Import the example's file as a new domain
3014        ##########################################
3015        from data_manager import sww2domain
3016        from Numeric import allclose
3017        import os
3018
3019        filename = domain.datadir+os.sep+domain.filename+'.sww'
3020        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
3021        #points, vertices, boundary = rectangular(15,15)
3022        #domain2.boundary = boundary
3023        ###################
3024        ##NOW TEST IT!!!
3025        ###################
3026
3027        os.remove(domain.filename + '.sww')
3028
3029        #FIXME smooth domain so that they can be compared
3030
3031
3032        bits = []#'vertex_coordinates']
3033        for quantity in ['elevation']+domain.quantities_to_be_stored:
3034            bits.append('quantities["%s"].get_integral()'%quantity)
3035
3036
3037        for bit in bits:
3038            #print 'testing that domain.'+bit+' has been restored'
3039            #print bit
3040            #print 'done'
3041            #print ('domain.'+bit), eval('domain.'+bit)
3042            #print ('domain2.'+bit), eval('domain2.'+bit)
3043            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
3044            pass
3045
3046        ######################################
3047        #Now evolve them both, just to be sure
3048        ######################################x
3049        visualise = False
3050        visualise = True
3051        domain.visualise = visualise
3052        domain.time = 0.
3053        from time import sleep
3054
3055        final = .5
3056        domain.set_quantity('friction', 0.1)
3057        domain.store = False
3058        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
3059
3060        for t in domain.evolve(yieldstep = yiel, finaltime = final):
3061            if visualise: sleep(.03)
3062            #domain.write_time()
3063            pass
3064
3065        domain2.smooth = True
3066        domain2.visualise = visualise
3067        domain2.store = False
3068        domain2.default_order=2
3069        domain2.set_quantity('friction', 0.1)
3070        #Bed-slope and friction
3071        # Boundary conditions
3072        Bd2=Dirichlet_boundary([0.2,0.,0.])
3073        Br2 = Reflective_boundary(domain2)
3074        domain2.boundary = domain.boundary
3075        #print 'domain2.boundary'
3076        #print domain2.boundary
3077        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
3078        #domain2.boundary = domain.boundary
3079        #domain2.set_boundary({'exterior': Bd})
3080
3081        domain2.check_integrity()
3082
3083        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3084            if visualise: sleep(.03)
3085            #domain2.write_time()
3086            pass
3087
3088        ###################
3089        ##NOW TEST IT!!!
3090        ##################
3091
3092        print '><><><><>>'
3093        bits = [ 'vertex_coordinates']
3094
3095        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
3096            #bits.append('quantities["%s"].get_integral()'%quantity)
3097            bits.append('get_quantity("%s").get_values()' %quantity)
3098
3099        for bit in bits:
3100            print bit
3101            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3102
3103
3104    def test_decimate_dem(self):
3105        """Test decimation of dem file
3106        """
3107
3108        import os
3109        from Numeric import ones, allclose, Float, arange
3110        from Scientific.IO.NetCDF import NetCDFFile
3111
3112        #Write test dem file
3113        root = 'decdemtest'
3114
3115        filename = root + '.dem'
3116        fid = NetCDFFile(filename, 'w')
3117
3118        fid.institution = 'Geoscience Australia'
3119        fid.description = 'NetCDF DEM format for compact and portable ' +\
3120                          'storage of spatial point data'
3121
3122        nrows = 15
3123        ncols = 18
3124
3125        fid.ncols = ncols
3126        fid.nrows = nrows
3127        fid.xllcorner = 2000.5
3128        fid.yllcorner = 3000.5
3129        fid.cellsize = 25
3130        fid.NODATA_value = -9999
3131
3132        fid.zone = 56
3133        fid.false_easting = 0.0
3134        fid.false_northing = 0.0
3135        fid.projection = 'UTM'
3136        fid.datum = 'WGS84'
3137        fid.units = 'METERS'
3138
3139        fid.createDimension('number_of_points', nrows*ncols)
3140
3141        fid.createVariable('elevation', Float, ('number_of_points',))
3142
3143        elevation = fid.variables['elevation']
3144
3145        elevation[:] = (arange(nrows*ncols))
3146
3147        fid.close()
3148
3149        #generate the elevation values expected in the decimated file
3150        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3151                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3152                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3153                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3154                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3155                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3156                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3157                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3158                         (144+145+146+162+163+164+180+181+182) / 9.0,
3159                         (148+149+150+166+167+168+184+185+186) / 9.0,
3160                         (152+153+154+170+171+172+188+189+190) / 9.0,
3161                         (156+157+158+174+175+176+192+193+194) / 9.0,
3162                         (216+217+218+234+235+236+252+253+254) / 9.0,
3163                         (220+221+222+238+239+240+256+257+258) / 9.0,
3164                         (224+225+226+242+243+244+260+261+262) / 9.0,
3165                         (228+229+230+246+247+248+264+265+266) / 9.0]
3166
3167        #generate a stencil for computing the decimated values
3168        stencil = ones((3,3), Float) / 9.0
3169
3170        decimate_dem(root, stencil=stencil, cellsize_new=100)
3171
3172        #Open decimated NetCDF file
3173        fid = NetCDFFile(root + '_100.dem', 'r')
3174
3175        # Get decimated elevation
3176        elevation = fid.variables['elevation']
3177
3178        #Check values
3179        assert allclose(elevation, ref_elevation)
3180
3181        #Cleanup
3182        fid.close()
3183
3184        os.remove(root + '.dem')
3185        os.remove(root + '_100.dem')
3186
3187    def test_decimate_dem_NODATA(self):
3188        """Test decimation of dem file that includes NODATA values
3189        """
3190
3191        import os
3192        from Numeric import ones, allclose, Float, arange, reshape
3193        from Scientific.IO.NetCDF import NetCDFFile
3194
3195        #Write test dem file
3196        root = 'decdemtest'
3197
3198        filename = root + '.dem'
3199        fid = NetCDFFile(filename, 'w')
3200
3201        fid.institution = 'Geoscience Australia'
3202        fid.description = 'NetCDF DEM format for compact and portable ' +\
3203                          'storage of spatial point data'
3204
3205        nrows = 15
3206        ncols = 18
3207        NODATA_value = -9999
3208
3209        fid.ncols = ncols
3210        fid.nrows = nrows
3211        fid.xllcorner = 2000.5
3212        fid.yllcorner = 3000.5
3213        fid.cellsize = 25
3214        fid.NODATA_value = NODATA_value
3215
3216        fid.zone = 56
3217        fid.false_easting = 0.0
3218        fid.false_northing = 0.0
3219        fid.projection = 'UTM'
3220        fid.datum = 'WGS84'
3221        fid.units = 'METERS'
3222
3223        fid.createDimension('number_of_points', nrows*ncols)
3224
3225        fid.createVariable('elevation', Float, ('number_of_points',))
3226
3227        elevation = fid.variables['elevation']
3228
3229        #generate initial elevation values
3230        elevation_tmp = (arange(nrows*ncols))
3231        #add some NODATA values
3232        elevation_tmp[0]   = NODATA_value
3233        elevation_tmp[95]  = NODATA_value
3234        elevation_tmp[188] = NODATA_value
3235        elevation_tmp[189] = NODATA_value
3236        elevation_tmp[190] = NODATA_value
3237        elevation_tmp[209] = NODATA_value
3238        elevation_tmp[252] = NODATA_value
3239
3240        elevation[:] = elevation_tmp
3241
3242        fid.close()
3243
3244        #generate the elevation values expected in the decimated file
3245        ref_elevation = [NODATA_value,
3246                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3247                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3248                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3249                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3250                         NODATA_value,
3251                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3252                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3253                         (144+145+146+162+163+164+180+181+182) / 9.0,
3254                         (148+149+150+166+167+168+184+185+186) / 9.0,
3255                         NODATA_value,
3256                         (156+157+158+174+175+176+192+193+194) / 9.0,
3257                         NODATA_value,
3258                         (220+221+222+238+239+240+256+257+258) / 9.0,
3259                         (224+225+226+242+243+244+260+261+262) / 9.0,
3260                         (228+229+230+246+247+248+264+265+266) / 9.0]
3261
3262        #generate a stencil for computing the decimated values
3263        stencil = ones((3,3), Float) / 9.0
3264
3265        decimate_dem(root, stencil=stencil, cellsize_new=100)
3266
3267        #Open decimated NetCDF file
3268        fid = NetCDFFile(root + '_100.dem', 'r')
3269
3270        # Get decimated elevation
3271        elevation = fid.variables['elevation']
3272
3273        #Check values
3274        assert allclose(elevation, ref_elevation)
3275
3276        #Cleanup
3277        fid.close()
3278
3279        os.remove(root + '.dem')
3280        os.remove(root + '_100.dem')
3281
3282    def xxxtestz_sww2ers_real(self):
3283        """Test that sww information can be converted correctly to asc/prj
3284        format readable by e.g. ArcView
3285        """
3286
3287        import time, os
3288        from Numeric import array, zeros, allclose, Float, concatenate
3289        from Scientific.IO.NetCDF import NetCDFFile
3290
3291        # the memory optimised least squares
3292        #  cellsize = 20,   # this one seems to hang
3293        #  cellsize = 200000, # Ran 1 test in 269.703s
3294                                #Ran 1 test in 267.344s
3295        #  cellsize = 20000,  # Ran 1 test in 460.922s
3296        #  cellsize = 2000   #Ran 1 test in 5340.250s
3297        #  cellsize = 200   #this one seems to hang, building matirx A
3298
3299        # not optimised
3300        # seems to hang
3301        #  cellsize = 2000   # Ran 1 test in 5334.563s
3302        #Export to ascii/prj files
3303        sww2dem('karratha_100m',
3304                quantity = 'depth',
3305                cellsize = 200000,
3306                verbose = True)
3307
3308    def test_read_asc(self):
3309        """Test conversion from dem in ascii format to native NetCDF xya format
3310        """
3311
3312        import time, os
3313        from Numeric import array, zeros, allclose, Float, concatenate
3314        from Scientific.IO.NetCDF import NetCDFFile
3315
3316        from data_manager import _read_asc
3317        #Write test asc file
3318        filename = tempfile.mktemp(".000")
3319        fid = open(filename, 'w')
3320        fid.write("""ncols         7
3321nrows         4
3322xllcorner     2000.5
3323yllcorner     3000.5
3324cellsize      25
3325NODATA_value  -9999
3326    97.921    99.285   125.588   180.830   258.645   342.872   415.836
3327   473.157   514.391   553.893   607.120   678.125   777.283   883.038
3328   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
3329   502.645   516.230   504.739   450.604   388.500   338.097   514.980
3330""")
3331        fid.close()
3332        bath_metadata, grid = _read_asc(filename, verbose=False)
3333        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
3334        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
3335        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
3336        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
3337        self.failUnless(grid[0][0]  == 97.921,  'Failed')
3338        self.failUnless(grid[3][6]  == 514.980,  'Failed')
3339
3340        os.remove(filename)
3341
3342    def test_asc_csiro2sww(self):
3343        import tempfile
3344
3345        bath_dir = tempfile.mkdtemp()
3346        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3347        #bath_dir = 'bath_data_manager_test'
3348        #print "os.getcwd( )",os.getcwd( )
3349        elevation_dir =  tempfile.mkdtemp()
3350        #elevation_dir = 'elev_expanded'
3351        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3352        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3353
3354        fid = open(bath_dir_filename, 'w')
3355        fid.write(""" ncols             3
3356 nrows             2
3357 xllcorner    148.00000
3358 yllcorner    -38.00000
3359 cellsize       0.25
3360 nodata_value   -9999.0
3361    9000.000    -1000.000    3000.0
3362   -1000.000    9000.000  -1000.000
3363""")
3364        fid.close()
3365
3366        fid = open(elevation_dir_filename1, 'w')
3367        fid.write(""" ncols             3
3368 nrows             2
3369 xllcorner    148.00000
3370 yllcorner    -38.00000
3371 cellsize       0.25
3372 nodata_value   -9999.0
3373    9000.000    0.000    3000.0
3374     0.000     9000.000     0.000
3375""")
3376        fid.close()
3377
3378        fid = open(elevation_dir_filename2, 'w')
3379        fid.write(""" ncols             3
3380 nrows             2
3381 xllcorner    148.00000
3382 yllcorner    -38.00000
3383 cellsize       0.25
3384 nodata_value   -9999.0
3385    9000.000    4000.000    4000.0
3386    4000.000    9000.000    4000.000
3387""")
3388        fid.close()
3389
3390        ucur_dir =  tempfile.mkdtemp()
3391        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3392        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3393
3394        fid = open(ucur_dir_filename1, 'w')
3395        fid.write(""" ncols             3
3396 nrows             2
3397 xllcorner    148.00000
3398 yllcorner    -38.00000
3399 cellsize       0.25
3400 nodata_value   -9999.0
3401    90.000    60.000    30.0
3402    10.000    10.000    10.000
3403""")
3404        fid.close()
3405        fid = open(ucur_dir_filename2, 'w')
3406        fid.write(""" ncols             3
3407 nrows             2
3408 xllcorner    148.00000
3409 yllcorner    -38.00000
3410 cellsize       0.25
3411 nodata_value   -9999.0
3412    90.000    60.000    30.0
3413    10.000    10.000    10.000
3414""")
3415        fid.close()
3416
3417        vcur_dir =  tempfile.mkdtemp()
3418        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3419        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3420
3421        fid = open(vcur_dir_filename1, 'w')
3422        fid.write(""" ncols             3
3423 nrows             2
3424 xllcorner    148.00000
3425 yllcorner    -38.00000
3426 cellsize       0.25
3427 nodata_value   -9999.0
3428    90.000    60.000    30.0
3429    10.000    10.000    10.000
3430""")
3431        fid.close()
3432        fid = open(vcur_dir_filename2, 'w')
3433        fid.write(""" ncols             3
3434 nrows             2
3435 xllcorner    148.00000
3436 yllcorner    -38.00000
3437 cellsize       0.25
3438 nodata_value   -9999.0
3439    90.000    60.000    30.0
3440    10.000    10.000    10.000
3441""")
3442        fid.close()
3443
3444        sww_file = 'a_test.sww'
3445        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
3446
3447        # check the sww file
3448
3449        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3450        x = fid.variables['x'][:]
3451        y = fid.variables['y'][:]
3452        z = fid.variables['z'][:]
3453        stage = fid.variables['stage'][:]
3454        xmomentum = fid.variables['xmomentum'][:]
3455        geo_ref = Geo_reference(NetCDFObject=fid)
3456        #print "geo_ref",geo_ref
3457        x_ref = geo_ref.get_xllcorner()
3458        y_ref = geo_ref.get_yllcorner()
3459        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3460        assert allclose(x_ref, 587798.418) # (-38, 148)
3461        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3462
3463        #Zone:   55
3464        #Easting:  588095.674  Northing: 5821451.722
3465        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3466        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3467
3468        #Zone:   55
3469        #Easting:  632145.632  Northing: 5820863.269
3470        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3471        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3472
3473        #Zone:   55
3474        #Easting:  609748.788  Northing: 5793447.860
3475        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3476        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3477
3478        assert allclose(z[0],9000.0 )
3479        assert allclose(stage[0][1],0.0 )
3480
3481        #(4000+1000)*60
3482        assert allclose(xmomentum[1][1],300000.0 )
3483
3484
3485        fid.close()
3486
3487        #tidy up
3488        os.remove(bath_dir_filename)
3489        os.rmdir(bath_dir)
3490
3491        os.remove(elevation_dir_filename1)
3492        os.remove(elevation_dir_filename2)
3493        os.rmdir(elevation_dir)
3494
3495        os.remove(ucur_dir_filename1)
3496        os.remove(ucur_dir_filename2)
3497        os.rmdir(ucur_dir)
3498
3499        os.remove(vcur_dir_filename1)
3500        os.remove(vcur_dir_filename2)
3501        os.rmdir(vcur_dir)
3502
3503
3504        # remove sww file
3505        os.remove(sww_file)
3506
3507    def test_asc_csiro2sww2(self):
3508        import tempfile
3509
3510        bath_dir = tempfile.mkdtemp()
3511        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3512        #bath_dir = 'bath_data_manager_test'
3513        #print "os.getcwd( )",os.getcwd( )
3514        elevation_dir =  tempfile.mkdtemp()
3515        #elevation_dir = 'elev_expanded'
3516        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3517        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3518
3519        fid = open(bath_dir_filename, 'w')
3520        fid.write(""" ncols             3
3521 nrows             2
3522 xllcorner    148.00000
3523 yllcorner    -38.00000
3524 cellsize       0.25
3525 nodata_value   -9999.0
3526    9000.000    -1000.000    3000.0
3527   -1000.000    9000.000  -1000.000
3528""")
3529        fid.close()
3530
3531        fid = open(elevation_dir_filename1, 'w')
3532        fid.write(""" ncols             3
3533 nrows             2
3534 xllcorner    148.00000
3535 yllcorner    -38.00000
3536 cellsize       0.25
3537 nodata_value   -9999.0
3538    9000.000    0.000    3000.0
3539     0.000     -9999.000     -9999.000
3540""")
3541        fid.close()
3542
3543        fid = open(elevation_dir_filename2, 'w')
3544        fid.write(""" ncols             3
3545 nrows             2
3546 xllcorner    148.00000
3547 yllcorner    -38.00000
3548 cellsize       0.25
3549 nodata_value   -9999.0
3550    9000.000    4000.000    4000.0
3551    4000.000    9000.000    4000.000
3552""")
3553        fid.close()
3554
3555        ucur_dir =  tempfile.mkdtemp()
3556        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3557        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3558
3559        fid = open(ucur_dir_filename1, 'w')
3560        fid.write(""" ncols             3
3561 nrows             2
3562 xllcorner    148.00000
3563 yllcorner    -38.00000
3564 cellsize       0.25
3565 nodata_value   -9999.0
3566    90.000    60.000    30.0
3567    10.000    10.000    10.000
3568""")
3569        fid.close()
3570        fid = open(ucur_dir_filename2, 'w')
3571        fid.write(""" ncols             3
3572 nrows             2
3573 xllcorner    148.00000
3574 yllcorner    -38.00000
3575 cellsize       0.25
3576 nodata_value   -9999.0
3577    90.000    60.000    30.0
3578    10.000    10.000    10.000
3579""")
3580        fid.close()
3581
3582        vcur_dir =  tempfile.mkdtemp()
3583        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3584        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3585
3586        fid = open(vcur_dir_filename1, 'w')
3587        fid.write(""" ncols             3
3588 nrows             2
3589 xllcorner    148.00000
3590 yllcorner    -38.00000
3591 cellsize       0.25
3592 nodata_value   -9999.0
3593    90.000    60.000    30.0
3594    10.000    10.000    10.000
3595""")
3596        fid.close()
3597        fid = open(vcur_dir_filename2, 'w')
3598        fid.write(""" ncols             3
3599 nrows             2
3600 xllcorner    148.00000
3601 yllcorner    -38.00000
3602 cellsize       0.25
3603 nodata_value   -9999.0
3604    90.000    60.000    30.0
3605    10.000    10.000    10.000
3606""")
3607        fid.close()
3608
3609        try:
3610            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3611                          vcur_dir, sww_file)
3612        except:
3613            #tidy up
3614            os.remove(bath_dir_filename)
3615            os.rmdir(bath_dir)
3616
3617            os.remove(elevation_dir_filename1)
3618            os.remove(elevation_dir_filename2)
3619            os.rmdir(elevation_dir)
3620
3621            os.remove(ucur_dir_filename1)
3622            os.remove(ucur_dir_filename2)
3623            os.rmdir(ucur_dir)
3624
3625            os.remove(vcur_dir_filename1)
3626            os.remove(vcur_dir_filename2)
3627            os.rmdir(vcur_dir)
3628        else:
3629            #tidy up
3630            os.remove(bath_dir_filename)
3631            os.rmdir(bath_dir)
3632
3633            os.remove(elevation_dir_filename1)
3634            os.remove(elevation_dir_filename2)
3635            os.rmdir(elevation_dir)
3636            raise 'Should raise exception'
3637
3638            os.remove(ucur_dir_filename1)
3639            os.remove(ucur_dir_filename2)
3640            os.rmdir(ucur_dir)
3641
3642            os.remove(vcur_dir_filename1)
3643            os.remove(vcur_dir_filename2)
3644            os.rmdir(vcur_dir)
3645
3646
3647
3648    def test_asc_csiro2sww3(self):
3649        import tempfile
3650
3651        bath_dir = tempfile.mkdtemp()
3652        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3653        #bath_dir = 'bath_data_manager_test'
3654        #print "os.getcwd( )",os.getcwd( )
3655        elevation_dir =  tempfile.mkdtemp()
3656        #elevation_dir = 'elev_expanded'
3657        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3658        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3659
3660        fid = open(bath_dir_filename, 'w')
3661        fid.write(""" ncols             3
3662 nrows             2
3663 xllcorner    148.00000
3664 yllcorner    -38.00000
3665 cellsize       0.25
3666 nodata_value   -9999.0
3667    9000.000    -1000.000    3000.0
3668   -1000.000    9000.000  -1000.000
3669""")
3670        fid.close()
3671
3672        fid = open(elevation_dir_filename1, 'w')
3673        fid.write(""" ncols             3
3674 nrows             2
3675 xllcorner    148.00000
3676 yllcorner    -38.00000
3677 cellsize       0.25
3678 nodata_value   -9999.0
3679    9000.000    0.000    3000.0
3680     0.000     -9999.000     -9999.000
3681""")
3682        fid.close()
3683
3684        fid = open(elevation_dir_filename2, 'w')
3685        fid.write(""" ncols             3
3686 nrows             2
3687 xllcorner    148.00000
3688 yllcorner    -38.00000
3689 cellsize       0.25
3690 nodata_value   -9999.0
3691    9000.000    4000.000    4000.0
3692    4000.000    9000.000    4000.000
3693""")
3694        fid.close()
3695
3696        ucur_dir =  tempfile.mkdtemp()
3697        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3698        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3699
3700        fid = open(ucur_dir_filename1, 'w')
3701        fid.write(""" ncols             3
3702 nrows             2
3703 xllcorner    148.00000
3704 yllcorner    -38.00000
3705 cellsize       0.25
3706 nodata_value   -9999.0
3707    90.000    60.000    30.0
3708    10.000    10.000    10.000
3709""")
3710        fid.close()
3711        fid = open(ucur_dir_filename2, 'w')
3712        fid.write(""" ncols             3
3713 nrows             2
3714 xllcorner    148.00000
3715 yllcorner    -38.00000
3716 cellsize       0.25
3717 nodata_value   -9999.0
3718    90.000    60.000    30.0
3719    10.000    10.000    10.000
3720""")
3721        fid.close()
3722
3723        vcur_dir =  tempfile.mkdtemp()
3724        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3725        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3726
3727        fid = open(vcur_dir_filename1, 'w')
3728        fid.write(""" ncols             3
3729 nrows             2
3730 xllcorner    148.00000
3731 yllcorner    -38.00000
3732 cellsize       0.25
3733 nodata_value   -9999.0
3734    90.000    60.000    30.0
3735    10.000    10.000    10.000
3736""")
3737        fid.close()
3738        fid = open(vcur_dir_filename2, 'w')
3739        fid.write(""" ncols             3
3740 nrows             2
3741 xllcorner    148.00000
3742 yllcorner    -38.00000
3743 cellsize       0.25
3744 nodata_value   -9999.0
3745    90.000    60.000    30.0
3746    10.000    10.000    10.000
3747""")
3748        fid.close()
3749
3750        sww_file = 'a_test.sww'
3751        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3752                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3753                      mean_stage = 100)
3754
3755        # check the sww file
3756
3757        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3758        x = fid.variables['x'][:]
3759        y = fid.variables['y'][:]
3760        z = fid.variables['z'][:]
3761        stage = fid.variables['stage'][:]
3762        xmomentum = fid.variables['xmomentum'][:]
3763        geo_ref = Geo_reference(NetCDFObject=fid)
3764        #print "geo_ref",geo_ref
3765        x_ref = geo_ref.get_xllcorner()
3766        y_ref = geo_ref.get_yllcorner()
3767        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3768        assert allclose(x_ref, 587798.418) # (-38, 148)
3769        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3770
3771        #Zone:   55
3772        #Easting:  588095.674  Northing: 5821451.722
3773        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3774        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3775
3776        #Zone:   55
3777        #Easting:  632145.632  Northing: 5820863.269
3778        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3779        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3780
3781        #Zone:   55
3782        #Easting:  609748.788  Northing: 5793447.860
3783        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3784        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3785
3786        assert allclose(z[0],9000.0 )
3787        assert allclose(stage[0][4],100.0 )
3788        assert allclose(stage[0][5],100.0 )
3789
3790        #(100.0 - 9000)*10
3791        assert allclose(xmomentum[0][4], -89000.0 )
3792
3793        #(100.0 - -1000.000)*10
3794        assert allclose(xmomentum[0][5], 11000.0 )
3795
3796        fid.close()
3797
3798        #tidy up
3799        os.remove(bath_dir_filename)
3800        os.rmdir(bath_dir)
3801
3802        os.remove(elevation_dir_filename1)
3803        os.remove(elevation_dir_filename2)
3804        os.rmdir(elevation_dir)
3805
3806        os.remove(ucur_dir_filename1)
3807        os.remove(ucur_dir_filename2)
3808        os.rmdir(ucur_dir)
3809
3810        os.remove(vcur_dir_filename1)
3811        os.remove(vcur_dir_filename2)
3812        os.rmdir(vcur_dir)
3813
3814        # remove sww file
3815        os.remove(sww_file)
3816
3817
3818    def test_asc_csiro2sww4(self):
3819        """
3820        Test specifying the extent
3821        """
3822
3823        import tempfile
3824
3825        bath_dir = tempfile.mkdtemp()
3826        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3827        #bath_dir = 'bath_data_manager_test'
3828        #print "os.getcwd( )",os.getcwd( )
3829        elevation_dir =  tempfile.mkdtemp()
3830        #elevation_dir = 'elev_expanded'
3831        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3832        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3833
3834        fid = open(bath_dir_filename, 'w')
3835        fid.write(""" ncols             4
3836 nrows             4
3837 xllcorner    148.00000
3838 yllcorner    -38.00000
3839 cellsize       0.25
3840 nodata_value   -9999.0
3841   -9000.000    -1000.000   -3000.0 -2000.000
3842   -1000.000    9000.000  -1000.000 -3000.000
3843   -4000.000    6000.000   2000.000 -5000.000
3844   -9000.000    -1000.000   -3000.0 -2000.000
3845""")
3846        fid.close()
3847
3848        fid = open(elevation_dir_filename1, 'w')
3849        fid.write(""" ncols             4
3850 nrows             4
3851 xllcorner    148.00000
3852 yllcorner    -38.00000
3853 cellsize       0.25
3854 nodata_value   -9999.0
3855   -900.000    -100.000   -300.0 -200.000
3856   -100.000    900.000  -100.000 -300.000
3857   -400.000    600.000   200.000 -500.000
3858   -900.000    -100.000   -300.0 -200.000
3859""")
3860        fid.close()
3861
3862        fid = open(elevation_dir_filename2, 'w')
3863        fid.write(""" ncols             4
3864 nrows             4
3865 xllcorner    148.00000
3866 yllcorner    -38.00000
3867 cellsize       0.25
3868 nodata_value   -9999.0
3869   -990.000    -110.000   -330.0 -220.000
3870   -110.000    990.000  -110.000 -330.000
3871   -440.000    660.000   220.000 -550.000
3872   -990.000    -110.000   -330.0 -220.000
3873""")
3874        fid.close()
3875
3876        ucur_dir =  tempfile.mkdtemp()
3877        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3878        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3879
3880        fid = open(ucur_dir_filename1, 'w')
3881        fid.write(""" ncols             4
3882 nrows             4
3883 xllcorner    148.00000
3884 yllcorner    -38.00000
3885 cellsize       0.25
3886 nodata_value   -9999.0
3887   -90.000    -10.000   -30.0 -20.000
3888   -10.000    90.000  -10.000 -30.000
3889   -40.000    60.000   20.000 -50.000
3890   -90.000    -10.000   -30.0 -20.000
3891""")
3892        fid.close()
3893        fid = open(ucur_dir_filename2, 'w')
3894        fid.write(""" ncols             4
3895 nrows             4
3896 xllcorner    148.00000
3897 yllcorner    -38.00000
3898 cellsize       0.25
3899 nodata_value   -9999.0
3900   -90.000    -10.000   -30.0 -20.000
3901   -10.000    99.000  -11.000 -30.000
3902   -40.000    66.000   22.000 -50.000
3903   -90.000    -10.000   -30.0 -20.000
3904""")
3905        fid.close()
3906
3907        vcur_dir =  tempfile.mkdtemp()
3908        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3909        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3910
3911        fid = open(vcur_dir_filename1, 'w')
3912        fid.write(""" ncols             4
3913 nrows             4
3914 xllcorner    148.00000
3915 yllcorner    -38.00000
3916 cellsize       0.25
3917 nodata_value   -9999.0
3918   -90.000    -10.000   -30.0 -20.000
3919   -10.000    80.000  -20.000 -30.000
3920   -40.000    50.000   10.000 -50.000
3921   -90.000    -10.000   -30.0 -20.000
3922""")
3923        fid.close()
3924        fid = open(vcur_dir_filename2, 'w')
3925        fid.write(""" ncols             4
3926 nrows             4
3927 xllcorner    148.00000
3928 yllcorner    -38.00000
3929 cellsize       0.25
3930 nodata_value   -9999.0
3931   -90.000    -10.000   -30.0 -20.000
3932   -10.000    88.000  -22.000 -30.000
3933   -40.000    55.000   11.000 -50.000
3934   -90.000    -10.000   -30.0 -20.000
3935""")
3936        fid.close()
3937
3938        sww_file = tempfile.mktemp(".sww")
3939        #sww_file = 'a_test.sww'
3940        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3941                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3942                      mean_stage = 100,
3943                       minlat = -37.6, maxlat = -37.6,
3944                  minlon = 148.3, maxlon = 148.3
3945                      #,verbose = True
3946                      )
3947
3948        # check the sww file
3949
3950        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3951        x = fid.variables['x'][:]
3952        y = fid.variables['y'][:]
3953        z = fid.variables['z'][:]
3954        stage = fid.variables['stage'][:]
3955        xmomentum = fid.variables['xmomentum'][:]
3956        ymomentum = fid.variables['ymomentum'][:]
3957        geo_ref = Geo_reference(NetCDFObject=fid)
3958        #print "geo_ref",geo_ref
3959        x_ref = geo_ref.get_xllcorner()
3960        y_ref = geo_ref.get_yllcorner()
3961        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3962
3963        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3964        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3965        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3966
3967        #Easting:  632145.632  Northing: 5820863.269
3968        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3969
3970        #print "x",x
3971        #print "y",y
3972        self.failUnless(len(x) == 4,'failed') # 2*2
3973        self.failUnless(len(x) == 4,'failed') # 2*2
3974
3975        #Zone:   55
3976        #Easting:  632145.632  Northing: 5820863.269
3977        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3978        # magic number - y is close enough for me.
3979        assert allclose(x[3], 632145.63 - x_ref)
3980        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3981
3982        assert allclose(z[0],9000.0 ) #z is elevation info
3983        #print "z",z
3984        # 2 time steps, 4 points
3985        self.failUnless(xmomentum.shape == (2,4), 'failed')
3986        self.failUnless(ymomentum.shape == (2,4), 'failed')
3987
3988        #(100.0 - -1000.000)*10
3989        #assert allclose(xmomentum[0][5], 11000.0 )
3990
3991        fid.close()
3992
3993        # is the sww file readable?
3994        #Lets see if we can convert it to a dem!
3995        #print "sww_file",sww_file
3996        #dem_file = tempfile.mktemp(".dem")
3997        domain = sww2domain(sww_file) ###, dem_file)
3998        domain.check_integrity()
3999
4000        #tidy up
4001        os.remove(bath_dir_filename)
4002        os.rmdir(bath_dir)
4003
4004        os.remove(elevation_dir_filename1)
4005        os.remove(elevation_dir_filename2)
4006        os.rmdir(elevation_dir)
4007
4008        os.remove(ucur_dir_filename1)
4009        os.remove(ucur_dir_filename2)
4010        os.rmdir(ucur_dir)
4011
4012        os.remove(vcur_dir_filename1)
4013        os.remove(vcur_dir_filename2)
4014        os.rmdir(vcur_dir)
4015
4016
4017
4018
4019        # remove sww file
4020        os.remove(sww_file)
4021
4022        # remove dem file
4023        #os.remove(dem_file)
4024
4025    def test_get_min_max_indexes(self):
4026        latitudes = [3,2,1,0]
4027        longitudes = [0,10,20,30]
4028
4029        # k - lat
4030        # l - lon
4031        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4032            latitudes,longitudes,
4033            -10,4,-10,31)
4034
4035        #print "kmin",kmin;print "kmax",kmax
4036        #print "lmin",lmin;print "lmax",lmax
4037        latitudes_new = latitudes[kmin:kmax]
4038        longitudes_news = longitudes[lmin:lmax]
4039        #print "latitudes_new", latitudes_new
4040        #print "longitudes_news",longitudes_news
4041        self.failUnless(latitudes == latitudes_new and \
4042                        longitudes == longitudes_news,
4043                         'failed')
4044
4045        ## 2nd test
4046        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4047            latitudes,longitudes,
4048            0.5,2.5,5,25)
4049        #print "kmin",kmin;print "kmax",kmax
4050        #print "lmin",lmin;print "lmax",lmax
4051        latitudes_new = latitudes[kmin:kmax]
4052        longitudes_news = longitudes[lmin:lmax]
4053        #print "latitudes_new", latitudes_new
4054        #print "longitudes_news",longitudes_news
4055
4056        self.failUnless(latitudes == latitudes_new and \
4057                        longitudes == longitudes_news,
4058                         'failed')
4059
4060        ## 3rd test
4061        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(\
4062            latitudes,
4063            longitudes,
4064            1.1,1.9,12,17)
4065        #print "kmin",kmin;print "kmax",kmax
4066        #print "lmin",lmin;print "lmax",lmax
4067        latitudes_new = latitudes[kmin:kmax]
4068        longitudes_news = longitudes[lmin:lmax]
4069        #print "latitudes_new", latitudes_new
4070        #print "longitudes_news",longitudes_news
4071
4072        self.failUnless(latitudes_new == [2, 1] and \
4073                        longitudes_news == [10, 20],
4074                         'failed')
4075
4076
4077        ## 4th test
4078        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4079            latitudes,longitudes,
4080                                                      -0.1,1.9,-2,17)
4081        #print "kmin",kmin;print "kmax",kmax
4082        #print "lmin",lmin;print "lmax",lmax
4083        latitudes_new = latitudes[kmin:kmax]
4084        longitudes_news = longitudes[lmin:lmax]
4085        #print "latitudes_new", latitudes_new
4086        #print "longitudes_news",longitudes_news
4087
4088        self.failUnless(latitudes_new == [2, 1, 0] and \
4089                        longitudes_news == [0, 10, 20],
4090                         'failed')
4091        ## 5th test
4092        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4093            latitudes,longitudes,
4094            0.1,1.9,2,17)
4095        #print "kmin",kmin;print "kmax",kmax
4096        #print "lmin",lmin;print "lmax",lmax
4097        latitudes_new = latitudes[kmin:kmax]
4098        longitudes_news = longitudes[lmin:lmax]
4099        #print "latitudes_new", latitudes_new
4100        #print "longitudes_news",longitudes_news
4101
4102        self.failUnless(latitudes_new == [2, 1, 0] and \
4103                        longitudes_news == [0, 10, 20],
4104                         'failed')
4105
4106        ## 6th test
4107
4108        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4109            latitudes,longitudes,
4110            1.5,4,18,32)
4111        #print "kmin",kmin;print "kmax",kmax
4112        #print "lmin",lmin;print "lmax",lmax
4113        latitudes_new = latitudes[kmin:kmax]
4114        longitudes_news = longitudes[lmin:lmax]
4115        #print "latitudes_new", latitudes_new
4116        #print "longitudes_news",longitudes_news
4117
4118        self.failUnless(latitudes_new == [3, 2, 1] and \
4119                        longitudes_news == [10, 20, 30],
4120                         'failed')
4121
4122
4123        ## 7th test
4124        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4125        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4126            latitudes,longitudes,
4127            1.5,1.5,15,15)
4128        #print "kmin",kmin;print "kmax",kmax
4129        #print "lmin",lmin;print "lmax",lmax
4130        latitudes_new = latitudes[kmin:kmax]
4131        longitudes_news = longitudes[lmin:lmax]
4132        m2d = m2d[kmin:kmax,lmin:lmax]
4133        #print "m2d", m2d
4134        #print "latitudes_new", latitudes_new
4135        #print "longitudes_news",longitudes_news
4136
4137        self.failUnless(latitudes_new == [2, 1] and \
4138                        longitudes_news == [10, 20],
4139                         'failed')
4140
4141        self.failUnless(m2d == [[5,6],[9,10]],
4142                         'failed')
4143
4144    def test_get_min_max_indexes2(self):
4145        latitudes = [-30,-35,-40,-45]
4146        longitudes = [148,149,150,151]
4147
4148        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4149
4150        # k - lat
4151        # l - lon
4152        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4153            latitudes,longitudes,
4154            -37,-27,147,149.5)
4155
4156        #print "kmin",kmin;print "kmax",kmax
4157        #print "lmin",lmin;print "lmax",lmax
4158        #print "m2d", m2d
4159        #print "latitudes", latitudes
4160        #print "longitudes",longitudes
4161        #print "latitudes[kmax]", latitudes[kmax]
4162        latitudes_new = latitudes[kmin:kmax]
4163        longitudes_new = longitudes[lmin:lmax]
4164        m2d = m2d[kmin:kmax,lmin:lmax]
4165        #print "m2d", m2d
4166        #print "latitudes_new", latitudes_new
4167        #print "longitudes_new",longitudes_new
4168
4169        self.failUnless(latitudes_new == [-30, -35, -40] and \
4170                        longitudes_new == [148, 149,150],
4171                         'failed')
4172        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
4173                         'failed')
4174
4175    def test_get_min_max_indexes3(self):
4176        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4177        longitudes = [148,149,150,151]
4178
4179        # k - lat
4180        # l - lon
4181        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4182            latitudes,longitudes,
4183            -43,-37,148.5,149.5)
4184
4185
4186        #print "kmin",kmin;print "kmax",kmax
4187        #print "lmin",lmin;print "lmax",lmax
4188        #print "latitudes", latitudes
4189        #print "longitudes",longitudes
4190        latitudes_new = latitudes[kmin:kmax]
4191        longitudes_news = longitudes[lmin:lmax]
4192        #print "latitudes_new", latitudes_new
4193        #print "longitudes_news",longitudes_news
4194
4195        self.failUnless(latitudes_new == [-35, -40, -45] and \
4196                        longitudes_news == [148, 149,150],
4197                         'failed')
4198
4199    def test_get_min_max_indexes4(self):
4200        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4201        longitudes = [148,149,150,151]
4202
4203        # k - lat
4204        # l - lon
4205        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4206            latitudes,longitudes)
4207
4208        #print "kmin",kmin;print "kmax",kmax
4209        #print "lmin",lmin;print "lmax",lmax
4210        #print "latitudes", latitudes
4211        #print "longitudes",longitudes
4212        latitudes_new = latitudes[kmin:kmax]
4213        longitudes_news = longitudes[lmin:lmax]
4214        #print "latitudes_new", latitudes_new
4215        #print "longitudes_news",longitudes_news
4216
4217        self.failUnless(latitudes_new == latitudes  and \
4218                        longitudes_news == longitudes,
4219                         'failed')
4220
4221    def test_tsh2sww(self):
4222        import os
4223        import tempfile
4224
4225        tsh_file = tempfile.mktemp(".tsh")
4226        file = open(tsh_file,"w")
4227        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
42280 0.0 0.0 0.0 0.0 0.01 \n \
42291 1.0 0.0 10.0 10.0 0.02  \n \
42302 0.0 1.0 0.0 10.0 0.03  \n \
42313 0.5 0.25 8.0 12.0 0.04  \n \
4232# Vert att title  \n \
4233elevation  \n \
4234stage  \n \
4235friction  \n \
42362 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
42370 0 3 2 -1  -1  1 dsg\n\
42381 0 1 3 -1  0 -1   ole nielsen\n\
42394 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
42400 1 0 2 \n\
42411 0 2 3 \n\
42422 2 3 \n\
42433 3 1 1 \n\
42443 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
42450 216.0 -86.0 \n \
42461 160.0 -167.0 \n \
42472 114.0 -91.0 \n \
42483 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
42490 0 1 0 \n \
42501 1 2 0 \n \
42512 2 0 0 \n \
42520 # <x> <y> ...Mesh Holes... \n \
42530 # <x> <y> <attribute>...Mesh Regions... \n \
42540 # <x> <y> <attribute>...Mesh Regions, area... \n\
4255#Geo reference \n \
425656 \n \
4257140 \n \
4258120 \n")
4259        file.close()
4260
4261        #sww_file = tempfile.mktemp(".sww")
4262        #print "sww_file",sww_file
4263        #print "sww_file",tsh_file
4264        tsh2sww(tsh_file)
4265
4266        os.remove(tsh_file)
4267        os.remove(tsh_file[:-4] + '.sww')
4268       
4269
4270
4271    def NOT_test_dem2pts(self):
4272        """Test conversion from dem in ascii format to native NetCDF xya format
4273        """
4274
4275        import time, os
4276        from Numeric import array, zeros, allclose, Float, concatenate
4277        from Scientific.IO.NetCDF import NetCDFFile
4278
4279        #Write test asc file
4280        root = 'demtest'
4281
4282        filename = root+'.asc'
4283        fid = open(filename, 'w')
4284        fid.write("""ncols         5
4285nrows         6
4286xllcorner     2000.5
4287yllcorner     3000.5
4288cellsize      25
4289NODATA_value  -9999
4290""")
4291        #Create linear function
4292
4293        ref_points = []
4294        ref_elevation = []
4295        for i in range(6):
4296            y = (6-i)*25.0
4297            for j in range(5):
4298                x = j*25.0
4299                z = x+2*y
4300
4301                ref_points.append( [x,y] )
4302                ref_elevation.append(z)
4303                fid.write('%f ' %z)
4304            fid.write('\n')
4305
4306        fid.close()
4307
4308        #Write prj file with metadata
4309        metafilename = root+'.prj'
4310        fid = open(metafilename, 'w')
4311
4312
4313        fid.write("""Projection UTM
4314Zone 56
4315Datum WGS84
4316Zunits NO
4317Units METERS
4318Spheroid WGS84
4319Xshift 0.0000000000
4320Yshift 10000000.0000000000
4321Parameters
4322""")
4323        fid.close()
4324
4325        #Convert to NetCDF pts
4326        convert_dem_from_ascii2netcdf(root)
4327        dem2pts(root)
4328
4329        #Check contents
4330        #Get NetCDF
4331        fid = NetCDFFile(root+'.pts', 'r')
4332
4333        # Get the variables
4334        #print fid.variables.keys()
4335        points = fid.variables['points']
4336        elevation = fid.variables['elevation']
4337
4338        #Check values
4339
4340        #print points[:]
4341        #print ref_points
4342        assert allclose(points, ref_points)
4343
4344        #print attributes[:]
4345        #print ref_elevation
4346        assert allclose(elevation, ref_elevation)
4347
4348        #Cleanup
4349        fid.close()
4350
4351
4352        os.remove(root + '.pts')
4353        os.remove(root + '.dem')
4354        os.remove(root + '.asc')
4355        os.remove(root + '.prj')
4356
4357
4358
4359    def NOT_test_dem2pts_bounding_box(self):
4360        """Test conversion from dem in ascii format to native NetCDF xya format
4361        """
4362
4363        import time, os
4364        from Numeric import array, zeros, allclose, Float, concatenate
4365        from Scientific.IO.NetCDF import NetCDFFile
4366
4367        #Write test asc file
4368        root = 'demtest'
4369
4370        filename = root+'.asc'
4371        fid = open(filename, 'w')
4372        fid.write("""ncols         5
4373nrows         6
4374xllcorner     2000.5
4375yllcorner     3000.5
4376cellsize      25
4377NODATA_value  -9999
4378""")
4379        #Create linear function
4380
4381        ref_points = []
4382        ref_elevation = []
4383        for i in range(6):
4384            y = (6-i)*25.0
4385            for j in range(5):
4386                x = j*25.0
4387                z = x+2*y
4388
4389                ref_points.append( [x,y] )
4390                ref_elevation.append(z)
4391                fid.write('%f ' %z)
4392            fid.write('\n')
4393
4394        fid.close()
4395
4396        #Write prj file with metadata
4397        metafilename = root+'.prj'
4398        fid = open(metafilename, 'w')
4399
4400
4401        fid.write("""Projection UTM
4402Zone 56
4403Datum WGS84
4404Zunits NO
4405Units METERS
4406Spheroid WGS84
4407Xshift 0.0000000000
4408Yshift 10000000.0000000000
4409Parameters
4410""")
4411        fid.close()
4412
4413        #Convert to NetCDF pts
4414        convert_dem_from_ascii2netcdf(root)
4415        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4416                northing_min=3035.0, northing_max=3125.5)
4417
4418        #Check contents
4419        #Get NetCDF
4420        fid = NetCDFFile(root+'.pts', 'r')
4421
4422        # Get the variables
4423        #print fid.variables.keys()
4424        points = fid.variables['points']
4425        elevation = fid.variables['elevation']
4426
4427        #Check values
4428        assert fid.xllcorner[0] == 2010.0
4429        assert fid.yllcorner[0] == 3035.0
4430
4431        #create new reference points
4432        ref_points = []
4433        ref_elevation = []
4434        for i in range(4):
4435            y = (4-i)*25.0 + 25.0
4436            y_new = y + 3000.5 - 3035.0
4437            for j in range(4):
4438                x = j*25.0 + 25.0
4439                x_new = x + 2000.5 - 2010.0
4440                z = x+2*y
4441
4442                ref_points.append( [x_new,y_new] )
4443                ref_elevation.append(z)
4444
4445        #print points[:]
4446        #print ref_points
4447        assert allclose(points, ref_points)
4448
4449        #print attributes[:]
4450        #print ref_elevation
4451        assert allclose(elevation, ref_elevation)
4452
4453        #Cleanup
4454        fid.close()
4455
4456
4457        os.remove(root + '.pts')
4458        os.remove(root + '.dem')
4459        os.remove(root + '.asc')
4460        os.remove(root + '.prj')
4461
4462
4463
4464    def NOT_test_dem2pts_remove_Nullvalues(self):
4465        """Test conversion from dem in ascii format to native NetCDF xya format
4466        """
4467
4468        import time, os
4469        from Numeric import array, zeros, allclose, Float, concatenate
4470        from Scientific.IO.NetCDF import NetCDFFile
4471
4472        #Write test asc file
4473        root = 'demtest'
4474
4475        filename = root+'.asc'
4476        fid = open(filename, 'w')
4477        fid.write("""ncols         5
4478nrows         6
4479xllcorner     2000.5
4480yllcorner     3000.5
4481cellsize      25
4482NODATA_value  -9999
4483""")
4484        #Create linear function
4485        # ref_ will write all the values
4486        # new_ref_ will write the values except for NODATA_values
4487        ref_points = []
4488        ref_elevation = []
4489        new_ref_pts = []
4490        new_ref_elev = []
4491        NODATA_value = -9999
4492        for i in range(6):
4493            y = (6-i)*25.0
4494            for j in range(5):
4495                x = j*25.0
4496                z = x+2*y
4497                if j == 4: z = NODATA_value # column
4498                if i == 2 and j == 2: z = NODATA_value # random
4499                if i == 5 and j == 1: z = NODATA_value
4500                if i == 1: z = NODATA_value # row
4501                if i == 3 and j == 1: z = NODATA_value # two pts/row
4502                if i == 3 and j == 3: z = NODATA_value
4503
4504
4505                if z <> NODATA_value:
4506                    new_ref_elev.append(z)
4507                    new_ref_pts.append( [x,y] )
4508
4509                ref_points.append( [x,y] )
4510                ref_elevation.append(z)
4511
4512                fid.write('%f ' %z)
4513            fid.write('\n')
4514
4515        fid.close()
4516
4517
4518        #Write prj file with metadata
4519        metafilename = root+'.prj'
4520        fid = open(metafilename, 'w')
4521
4522
4523        fid.write("""Projection UTM
4524Zone 56
4525Datum WGS84
4526Zunits NO
4527Units METERS
4528Spheroid WGS84
4529Xshift 0.0000000000
4530Yshift 10000000.0000000000
4531Parameters
4532""")
4533        fid.close()
4534
4535        #Convert to NetCDF pts
4536        convert_dem_from_ascii2netcdf(root)
4537        dem2pts(root)
4538
4539        #Check contents
4540        #Get NetCDF
4541        fid = NetCDFFile(root+'.pts', 'r')
4542
4543        # Get the variables
4544        #print fid.variables.keys()
4545        points = fid.variables['points']
4546        elevation = fid.variables['elevation']
4547
4548        #Check values
4549        #print 'points', points[:]
4550        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
4551        assert allclose(points, new_ref_pts), 'points do not align'
4552
4553        #print 'elevation', elevation[:]
4554        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
4555        assert allclose(elevation, new_ref_elev), 'elevations do not align'
4556
4557        #Cleanup
4558        fid.close()
4559
4560
4561        os.remove(root + '.pts')
4562        os.remove(root + '.dem')
4563        os.remove(root + '.asc')
4564        os.remove(root + '.prj')
4565
4566    def NOT_test_dem2pts_bounding_box_Nullvalues(self):
4567        """Test conversion from dem in ascii format to native NetCDF xya format
4568        """
4569
4570        import time, os
4571        from Numeric import array, zeros, allclose, Float, concatenate
4572        from Scientific.IO.NetCDF import NetCDFFile
4573
4574        #Write test asc file
4575        root = 'demtest'
4576
4577        filename = root+'.asc'
4578        fid = open(filename, 'w')
4579        fid.write("""ncols         5
4580nrows         6
4581xllcorner     2000.5
4582yllcorner     3000.5
4583cellsize      25
4584NODATA_value  -9999
4585""")
4586        #Create linear function
4587
4588        ref_points = []
4589        ref_elevation = []
4590        new_ref_pts1 = []
4591        new_ref_elev1 = []
4592        NODATA_value = -9999
4593        for i in range(6):
4594            y = (6-i)*25.0
4595            for j in range(5):
4596                x = j*25.0
4597                z = x+2*y
4598                if j == 4: z = NODATA_value # column
4599                if i == 2 and j == 2: z = NODATA_value # random
4600                if i == 5 and j == 1: z = NODATA_value
4601                if i == 1: z = NODATA_value # row
4602                if i == 3 and j == 1: z = NODATA_value # two pts/row
4603                if i == 3 and j == 3: z = NODATA_value
4604
4605                if z <> NODATA_value:
4606                    new_ref_elev1.append(z)
4607                    new_ref_pts1.append( [x,y] )
4608
4609                ref_points.append( [x,y] )
4610                ref_elevation.append(z)
4611                fid.write('%f ' %z)
4612            fid.write('\n')
4613
4614        fid.close()
4615
4616        #Write prj file with metadata
4617        metafilename = root+'.prj'
4618        fid = open(metafilename, 'w')
4619
4620
4621        fid.write("""Projection UTM
4622Zone 56
4623Datum WGS84
4624Zunits NO
4625Units METERS
4626Spheroid WGS84
4627Xshift 0.0000000000
4628Yshift 10000000.0000000000
4629Parameters
4630""")
4631        fid.close()
4632
4633        #Convert to NetCDF pts
4634        convert_dem_from_ascii2netcdf(root)
4635        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4636                northing_min=3035.0, northing_max=3125.5)
4637
4638        #Check contents
4639        #Get NetCDF
4640        fid = NetCDFFile(root+'.pts', 'r')
4641
4642        # Get the variables
4643        #print fid.variables.keys()
4644        points = fid.variables['points']
4645        elevation = fid.variables['elevation']
4646
4647        #Check values
4648        assert fid.xllcorner[0] == 2010.0
4649        assert fid.yllcorner[0] == 3035.0
4650
4651        #create new reference points
4652        ref_points = []
4653        ref_elevation = []
4654        new_ref_pts2 = []
4655        new_ref_elev2 = []
4656        for i in range(4):
4657            y = (4-i)*25.0 + 25.0
4658            y_new = y + 3000.5 - 3035.0
4659            for j in range(4):
4660                x = j*25.0 + 25.0
4661                x_new = x + 2000.5 - 2010.0
4662                z = x+2*y
4663
4664                if j == 3: z = NODATA_value # column
4665                if i == 1 and j == 1: z = NODATA_value # random
4666                if i == 4 and j == 0: z = NODATA_value
4667                if i == 0: z = NODATA_value # row
4668                if i == 2 and j == 0: z = NODATA_value # two pts/row
4669                if i == 2 and j == 2: z = NODATA_value
4670
4671                if z <> NODATA_value:
4672                    new_ref_elev2.append(z)
4673                    new_ref_pts2.append( [x_new,y_new] )
4674
4675
4676                ref_points.append( [x_new,y_new] )
4677                ref_elevation.append(z)
4678
4679        #print points[:]
4680        #print ref_points
4681        #assert allclose(points, ref_points)
4682
4683        #print attributes[:]
4684        #print ref_elevation
4685        #assert allclose(elevation, ref_elevation)
4686
4687
4688        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
4689        assert allclose(points, new_ref_pts2), 'points do not align'
4690
4691        #print 'elevation', elevation[:]
4692        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
4693        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
4694        #Cleanup
4695        fid.close()
4696
4697
4698        os.remove(root + '.pts')
4699        os.remove(root + '.dem')
4700        os.remove(root + '.asc')
4701        os.remove(root + '.prj')
4702
4703
4704########## testing nbed class ##################
4705    def test_exposure_csv_loading(self):
4706       
4707
4708        file_name = tempfile.mktemp(".xya")
4709        file = open(file_name,"w")
4710        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
4711115.0, -21.0, splat, 0.0\n\
4712114.0, -21.7, pow, 10.0\n\
4713114.5, -21.4, bang, 40.0\n")
4714        file.close()
4715        exposure = Exposure_csv(file_name)
4716        exposure.get_column("sound")
4717       
4718        self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
4719                        'FAILED!')
4720        self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
4721                        'FAILED!')
4722       
4723        os.remove(file_name)
4724       
4725    def test_exposure_csv_loading(self):
4726       
4727
4728        file_name = tempfile.mktemp(".xya")
4729        file = open(file_name,"w")
4730        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
4731115.0, -21.0, splat, 0.0\n\
4732114.0, -21.7, pow, 10.0\n\
4733114.5, -21.4, bang, 40.0\n")
4734        file.close()
4735        exposure = Exposure_csv(file_name)
4736        exposure.get_column("sound")
4737       
4738        self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
4739                        'FAILED!')
4740        self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
4741                        'FAILED!')
4742       
4743        os.remove(file_name)
4744
4745    def test_exposure_csv_cmp(self):
4746        file_name = tempfile.mktemp(".xya")
4747        file = open(file_name,"w")
4748        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
4749115.0, -21.0, splat, 0.0\n\
4750114.0, -21.7, pow, 10.0\n\
4751114.5, -21.4, bang, 40.0\n")
4752        file.close()
4753       
4754        e1 = Exposure_csv(file_name)
4755        e2 = Exposure_csv(file_name)
4756        os.remove(file_name)
4757
4758        self.failUnless(cmp(e1,e2)==0,
4759                        'FAILED!')
4760       
4761        self.failUnless(cmp(e1,"hey")==1,
4762                        'FAILED!')
4763       
4764        file_name = tempfile.mktemp(".xya")
4765        file = open(file_name,"w")
4766        # Note, this has less spaces in the title,
4767        # the instances will be the same.
4768        file.write("LATITUDE,LONGITUDE ,sound, speed \n\
4769115.0, -21.0, splat, 0.0\n\
4770114.0, -21.7, pow, 10.0\n\
4771114.5, -21.4, bang, 40.0\n")
4772        file.close()
4773        e3 = Exposure_csv(file_name)
4774        os.remove(file_name)
4775
4776        self.failUnless(cmp(e3,e2)==0,
4777                        'FAILED!')
4778       
4779        file_name = tempfile.mktemp(".xya")
4780        file = open(file_name,"w")
4781        # Note, 40 changed to 44 .
4782        file.write("LATITUDE,LONGITUDE ,sound, speed \n\
4783115.0, -21.0, splat, 0.0\n\
4784114.0, -21.7, pow, 10.0\n\
4785114.5, -21.4, bang, 44.0\n")
4786        file.close()
4787        e4 = Exposure_csv(file_name)
4788        os.remove(file_name)
4789        #print "e4",e4._attribute_dic
4790        #print "e2",e2._attribute_dic
4791        self.failUnless(cmp(e4,e2)<>0,
4792                        'FAILED!')
4793       
4794        file_name = tempfile.mktemp(".xya")
4795        file = open(file_name,"w")
4796        # Note, the first two columns are swapped.
4797        file.write("LONGITUDE,LATITUDE ,sound, speed \n\
4798 -21.0,115.0, splat, 0.0\n\
4799 -21.7,114.0, pow, 10.0\n\
4800 -21.4,114.5, bang, 40.0\n")
4801        file.close()
4802        e5 = Exposure_csv(file_name)
4803        os.remove(file_name)
4804
4805        self.failUnless(cmp(e3,e5)<>0,
4806                        'FAILED!')
4807       
4808    def test_exposure_csv_saving(self):
4809       
4810
4811        file_name = tempfile.mktemp(".xya")
4812        file = open(file_name,"w")
4813        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
4814115.0, -21.0, splat, 0.0\n\
4815114.0, -21.7, pow, 10.0\n\
4816114.5, -21.4, bang, 40.0\n")
4817        file.close()
4818        e1 = Exposure_csv(file_name)
4819       
4820        file_name2 = tempfile.mktemp(".xya")
4821        e1.save(file_name = file_name2)
4822        e2 = Exposure_csv(file_name2)
4823       
4824        self.failUnless(cmp(e1,e2)==0,
4825                        'FAILED!')
4826        os.remove(file_name)
4827        os.remove(file_name2)
4828
4829    def test_exposure_csv_get_location(self):
4830        file_name = tempfile.mktemp(".xya")
4831        file = open(file_name,"w")
4832        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
4833150.916666667, -34.5, splat, 0.0\n\
4834150.0, -34.0, pow, 10.0\n")
4835        file.close()
4836        e1 = Exposure_csv(file_name)
4837
4838        gsd = e1.get_location()
4839       
4840        points = gsd.get_data_points(absolute=True)
4841       
4842        assert allclose(points[0][0], 308728.009)
4843        assert allclose(points[0][1], 6180432.601)
4844        assert allclose(points[1][0],  222908.705)
4845        assert allclose(points[1][1], 6233785.284)
4846        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
4847                        'Bad zone error!')
4848
4849        os.remove(file_name)
4850       
4851    def test_exposure_csv_set_column_get_column(self):
4852        file_name = tempfile.mktemp(".xya")
4853        file = open(file_name,"w")
4854        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
4855150.916666667, -34.5, splat, 0.0\n\
4856150.0, -34.0, pow, 10.0\n")
4857        file.close()
4858        e1 = Exposure_csv(file_name)     
4859        os.remove(file_name)
4860
4861        new_title = "feast"
4862        new_values = ["chicken","soup"]
4863        e1.set_column(new_title, new_values)
4864        returned_values = e1.get_column(new_title)
4865        self.failUnless(returned_values == new_values,
4866                        ' Error!')
4867       
4868        file_name2 = tempfile.mktemp(".xya")
4869        e1.save(file_name = file_name2)
4870        e2 = Exposure_csv(file_name2)
4871        returned_values = e2.get_column(new_title)
4872        self.failUnless(returned_values == new_values,
4873                        ' Error!')       
4874        os.remove(file_name2)
4875
4876    def test_exposure_csv_set_column_get_column_error_checking(self):
4877        file_name = tempfile.mktemp(".xya")
4878        file = open(file_name,"w")
4879        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
4880150.916666667, -34.5, splat, 0.0\n\
4881150.0, -34.0, pow, 10.0\n")
4882        file.close()
4883        e1 = Exposure_csv(file_name)     
4884        os.remove(file_name)
4885
4886        new_title = "sound"
4887        new_values = [12.5,7.6]
4888        try:
4889            e1.set_column(new_title, new_values)
4890        except TitleValueError:
4891            pass
4892        else:
4893            self.failUnless(0 ==1,  'Error not thrown error!')
4894           
4895        e1.set_column(new_title, new_values, overwrite=True)
4896        returned_values = e1.get_column(new_title)
4897        self.failUnless(returned_values == new_values,
4898                        ' Error!')       
4899       
4900        new2_title = "short list"
4901        new2_values = [12.5]
4902        try:
4903            e1.set_column(new2_title, new2_values)
4904        except DataMissingValuesError:
4905            pass
4906        else:
4907            self.failUnless(0 ==1,  'Error not thrown error!')
4908           
4909        new2_title = "long list"
4910        new2_values = [12.5, 7,8]
4911        try:
4912            e1.set_column(new2_title, new2_values)
4913        except DataMissingValuesError:
4914            pass
4915        else:
4916            self.failUnless(0 ==1,  'Error not thrown error!')
4917        file_name2 = tempfile.mktemp(".xya")
4918        e1.save(file_name = file_name2)
4919        e2 = Exposure_csv(file_name2)
4920        returned_values = e2.get_column(new_title)
4921        for returned, new in map(None, returned_values, new_values):
4922            self.failUnless(returned == str(new), ' Error!')
4923        #self.failUnless(returned_values == new_values, ' Error!')       
4924        os.remove(file_name2)
4925       
4926        try:
4927            e1.get_column("toe jam")
4928        except TitleValueError:
4929            pass
4930        else:
4931            self.failUnless(0 ==1,  'Error not thrown error!')
4932           
4933    def test_exposure_csv_loading_x_y(self):
4934       
4935
4936        file_name = tempfile.mktemp(".xya")
4937        file = open(file_name,"w")
4938        file.write("x, y ,sound  , speed \n\
4939115.0, 7, splat, 0.0\n\
4940114.0, 8.0, pow, 10.0\n\
4941114.5, 9., bang, 40.0\n")
4942        file.close()
4943        e1 = Exposure_csv(file_name, is_x_y_locations=True)
4944        gsd = e1.get_location()
4945       
4946        points = gsd.get_data_points(absolute=True)
4947       
4948        assert allclose(points[0][0], 115)
4949        assert allclose(points[0][1], 7)
4950        assert allclose(points[1][0], 114)
4951        assert allclose(points[1][1], 8)
4952        assert allclose(points[2][0], 114.5)
4953        assert allclose(points[2][1], 9)
4954        self.failUnless(gsd.get_geo_reference().get_zone() == -1,
4955                        'Bad zone error!')
4956
4957        os.remove(file_name)
4958
4959           
4960    def test_exposure_csv_loading_x_y2(self):
4961       
4962        csv_file = tempfile.mktemp(".csv")
4963        fd = open(csv_file,'wb')
4964        writer = csv.writer(fd)
4965        writer.writerow(['x','y','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
4966        writer.writerow([5.5,0.5,'199770','130000','Metal','Timber',20])
4967        writer.writerow([4.5,1.0,'150000','76000','Metal','Double Brick',20])
4968        writer.writerow([4.5,1.5,'150000','76000','Metal','Brick Veneer',20])
4969        fd.close()
4970
4971        e1 = Exposure_csv(csv_file)
4972        gsd = e1.get_location()
4973       
4974        points = gsd.get_data_points(absolute=True)
4975        assert allclose(points[0][0], 5.5)
4976        assert allclose(points[0][1], 0.5)
4977        assert allclose(points[1][0], 4.5)
4978        assert allclose(points[1][1], 1.0)
4979        assert allclose(points[2][0], 4.5)
4980        assert allclose(points[2][1], 1.5)
4981        self.failUnless(gsd.get_geo_reference().get_zone() == -1,
4982                        'Bad zone error!')
4983
4984        os.remove(csv_file)
4985
4986    #### TESTS FOR URS 2 SWW  ###     
4987   
4988    def create_mux(self, points_num=None):
4989        # write all the mux stuff.
4990        time_step_count = 3
4991        time_step = 0.5
4992       
4993        longitudes = [150.66667, 150.83334, 151., 151.16667]
4994        latitudes = [-34.5, -34.33333, -34.16667, -34]
4995
4996        if points_num == None:
4997            points_num = len(longitudes) * len(latitudes)
4998
4999        lonlatdeps = []
5000        quantities = ['HA','UA','VA']
5001        mux_names = ['-z-mux','-e-mux','-n-mux']
5002        quantities_init = [[],[],[]]
5003        for i,lon in enumerate(longitudes):
5004            for j,lat in enumerate(latitudes):
5005                _ , e, n = redfearn(lat, lon)
5006                lonlatdeps.append([lon, lat, n])
5007                quantities_init[0].append(e) # HA
5008                quantities_init[1].append(n ) # UA
5009                quantities_init[2].append(e) # VA
5010        #print "lonlatdeps",lonlatdeps
5011        _,base_name = tempfile.mkstemp("")
5012        files = []       
5013        for i,q in enumerate(quantities): 
5014            quantities_init[i] = ensure_numeric(quantities_init[i])
5015            #print "HA_init", HA_init
5016            q_time = zeros((time_step_count, points_num), Float)
5017            for time in range(time_step_count):
5018                q_time[time,:] = quantities_init[i] #* time * 4
5019           
5020            #Write C files
5021            columns = 3 # long, lat , depth
5022            file = base_name + mux_names[i]
5023            f = open(file, 'wb')
5024            files.append(file)
5025            f.write(pack('i',points_num))
5026            f.write(pack('i',time_step_count))
5027            f.write(pack('f',time_step))
5028
5029            #write lat/long info
5030            for lonlatdep in lonlatdeps:
5031                for float in lonlatdep:
5032                    f.write(pack('f',float))
5033                   
5034            # Write quantity info
5035            for time in  range(time_step_count):
5036                for i in range(points_num):
5037                    f.write(pack('f',q_time[time,i]))
5038            f.close()
5039        return base_name, files
5040       
5041   
5042    def delete_mux(self, files):
5043        for file in files:
5044            os.remove(file)
5045           
5046    def test_urs2sww_test_fail(self):
5047        points_num = -100
5048        time_step_count = 45
5049        time_step = -7
5050        _,base_name = tempfile.mkstemp("")
5051        files = []
5052        quantities = ['HA','UA','VA']
5053        mux_names = ['-z-mux','-e-mux','-n-mux']
5054        for i,q in enumerate(quantities): 
5055            #Write C files
5056            columns = 3 # long, lat , depth
5057            file = base_name + mux_names[i]
5058            f = open(file, 'wb')
5059            files.append(file)
5060            f.write(pack('i',points_num))
5061            f.write(pack('i',time_step_count))
5062            f.write(pack('f',time_step))
5063
5064            f.close()   
5065        tide = 1
5066        try:
5067            urs2sww(base_name, remove_nc_files=True, mean_stage=tide)       
5068        except ANUGAError:
5069            pass
5070        else:
5071            self.delete_mux(files)
5072            msg = 'Should have raised exception'
5073            raise msg
5074        sww_file = base_name + '.sww'
5075        self.delete_mux(files)
5076       
5077    def test_urs2sww_test_fail2(self):
5078        base_name = 'Harry-high-pants'
5079        try:
5080            urs2sww(base_name)       
5081        except IOError:
5082            pass
5083        else:
5084            self.delete_mux(files)
5085            msg = 'Should have raised exception'
5086            raise msg
5087           
5088    def test_urs2sww(self):
5089        tide = 1
5090        base_name, files = self.create_mux()
5091        urs2sww(base_name, mean_stage=tide)
5092        sww_file = base_name + '.sww'
5093       
5094        #Let's interigate the sww file
5095        fid = NetCDFFile(sww_file)
5096
5097        x = fid.variables['x'][:]
5098        y = fid.variables['y'][:]
5099        geo_reference = Geo_reference(NetCDFObject=fid)
5100       
5101        #Check that first coordinate is correctly represented       
5102        #Work out the UTM coordinates for first point
5103        zone, e, n = redfearn(-34.5, 150.66667)       
5104       
5105        assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
5106
5107        #Check first value
5108        stage = fid.variables['stage'][:]
5109        xmomentum = fid.variables['xmomentum'][:]
5110        ymomentum = fid.variables['ymomentum'][:]
5111
5112        #print ymomentum
5113        assert allclose(stage[0,0], e +tide)  #Meters
5114
5115        #Check the momentums - ua
5116        #momentum = velocity*(stage-elevation)
5117        #momentum = velocity*(stage+elevation)
5118        # -(-elevation) since elevation is inverted in mux files
5119        # = n*(e+tide+n) based on how I'm writing these files
5120        answer = n*(e+tide+n)
5121        actual = xmomentum[0,0]
5122        #print "answer",answer
5123        #print "actual",actual
5124        assert allclose(answer, actual)  #Meters
5125       
5126        fid.close()
5127
5128
5129        self.delete_mux(files)
5130        os.remove(sww_file)
5131       
5132    def bad_test_assuming_theres_a_mux_file(self):
5133        # these mux files aren't in the repository, plus they are bad!
5134        base_name = 'o-z-mux' 
5135        base_name = 'o-e-mux'
5136        file_name = base_name + '.nc'
5137        lonlatdep_numeric, lon, lat = _binary_c2nc(base_name, file_name, 'HA')
5138       
5139        #os.remove(file_name)
5140       
5141    def test_lon_lat2grid(self):
5142        lonlatdep = [
5143            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
5144            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
5145            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
5146            [ 113.19999695  ,  -26.33329964 ,   0.        ] ]
5147           
5148        long, lat = lon_lat2grid(lonlatdep)
5149
5150        for i, result in enumerate(lat):
5151            assert lonlatdep [i][1] == result
5152        assert len(lat) == 2 
5153
5154        for i, result in enumerate(long):
5155            assert lonlatdep [i*2][0] == result
5156        assert len(long) == 2
5157
5158    def test_lon_lat2grid_bad(self):
5159        lonlatdep  = [
5160            [ -26.06669998,  113.06700134,    1.        ],
5161            [ -26.06669998 , 113.19999695 ,   2.        ],
5162            [ -26.06669998 , 113.33300018,    3.        ],
5163            [ -26.06669998 , 113.43299866   , 4.        ],
5164            [ -26.20000076 , 113.06700134,    5.        ],
5165            [ -26.20000076 , 113.19999695 ,   6.        ],
5166            [ -26.20000076 , 113.33300018  ,  7.        ],
5167            [ -26.20000076 , 113.43299866   , 8.        ],
5168            [ -26.33329964 , 113.06700134,    9.        ],
5169            [ -26.33329964 , 113.19999695 ,   10.        ],
5170            [ -26.33329964 , 113.33300018  ,  11.        ],
5171            [ -26.33329964 , 113.43299866 ,   12.        ],
5172            [ -26.43330002 , 113.06700134 ,   13        ],
5173            [ -26.43330002 , 113.19999695 ,   14.        ],
5174            [ -26.43330002 , 113.33300018,    15.        ],
5175            [ -26.43330002 , 113.43299866,    16.        ]]
5176        try:
5177            long, lat = lon_lat2grid(lonlatdep)
5178        except AssertionError:
5179            pass
5180        else:
5181            msg = 'Should have raised exception'
5182            raise msg
5183       
5184    def test_lon_lat2gridII(self):
5185        lonlatdep = [
5186            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
5187            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
5188            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
5189            [ 113.19999695  ,  -26.344329964 ,   0.        ] ]
5190        try:
5191            long, lat = lon_lat2grid(lonlatdep)
5192        except AssertionError:
5193            pass
5194        else:
5195            msg = 'Should have raised exception'
5196            raise msg
5197       
5198    def trial_loading(self):
5199        basename_in = 'karratha'
5200        basename_out = basename_in
5201        urs2sww(basename_in, basename_out)
5202    #### END TESTS FOR URS 2 SWW  ###
5203
5204       
5205#-------------------------------------------------------------
5206if __name__ == "__main__":
5207    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_test_fail')
5208    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww')
5209    suite = unittest.makeSuite(Test_Data_Manager,'test')
5210    runner = unittest.TextTestRunner()
5211    runner.run(suite)
5212
5213
Note: See TracBrowser for help on using the repository browser.