source: anuga_core/source/anuga/pyvolution/test_data_manager.py @ 3529

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

minimum_allowed_depth is now an attribute of domain

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