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

Last change on this file since 3975 was 3975, checked in by duncan, 17 years ago

extend tests of urs2sww.

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