source: anuga_core/source/anuga/file_conversion/test_file_conversion.py @ 7743

Last change on this file since 7743 was 7743, checked in by hudson, 14 years ago

Removed more modules from data_handler: code to do with building destruction.
File size: 73.1 KB
Line 
1
2from anuga.shallow_water.shallow_water_domain import Domain
3from ferret2sww import ferret2sww
4from anuga.utilities.numerical_tools import ensure_numeric, mean
5from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
6     import Transmissive_boundary
7from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
8from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
9from anuga.shallow_water.sww_file import SWW_file
10from data_manager import extent_sww
11from anuga.config import netcdf_float, epsilon, g
12from Scientific.IO.NetCDF import NetCDFFile
13import sys
14import unittest
15import numpy as num
16import copy
17import os
18
19
20class Test_File_Conversion(unittest.TestCase):
21    """ A suite of tests to test file conversion functions.
22        These tests are quite coarse-grained: converting a file
23        and checking that its headers and some of its contents
24        are correct.
25    """
26    verbose = False
27
28    def set_verbose(self):
29        Test_File_Conversion.verbose = True
30       
31    def setUp(self):
32        import time
33       
34        self.verbose = Test_File_Conversion.verbose
35        # Create basic mesh
36        points, vertices, boundary = rectangular(2, 2)
37
38        # Create shallow water domain
39        domain = Domain(points, vertices, boundary)
40        domain.default_order = 2
41
42        # Set some field values
43        domain.set_quantity('elevation', lambda x,y: -x)
44        domain.set_quantity('friction', 0.03)
45
46
47        ######################
48        # Boundary conditions
49        B = Transmissive_boundary(domain)
50        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
51
52
53        ######################
54        #Initial condition - with jumps
55        bed = domain.quantities['elevation'].vertex_values
56        stage = num.zeros(bed.shape, num.float)
57
58        h = 0.3
59        for i in range(stage.shape[0]):
60            if i % 2 == 0:
61                stage[i,:] = bed[i,:] + h
62            else:
63                stage[i,:] = bed[i,:]
64
65        domain.set_quantity('stage', stage)
66
67
68        domain.distribute_to_vertices_and_edges()               
69        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
70
71
72        self.domain = domain
73
74        C = domain.get_vertex_coordinates()
75        self.X = C[:,0:6:2].copy()
76        self.Y = C[:,1:6:2].copy()
77
78        self.F = bed
79
80        #Write A testfile (not realistic. Values aren't realistic)
81        self.test_MOST_file = 'most_small'
82
83        longitudes = [150.66667, 150.83334, 151., 151.16667]
84        latitudes = [-34.5, -34.33333, -34.16667, -34]
85
86        long_name = 'LON'
87        lat_name = 'LAT'
88
89        nx = 4
90        ny = 4
91        six = 6
92
93
94        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
95            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
96
97            fid.createDimension(long_name,nx)
98            fid.createVariable(long_name,netcdf_float,(long_name,))
99            fid.variables[long_name].point_spacing='uneven'
100            fid.variables[long_name].units='degrees_east'
101            fid.variables[long_name].assignValue(longitudes)
102
103            fid.createDimension(lat_name,ny)
104            fid.createVariable(lat_name,netcdf_float,(lat_name,))
105            fid.variables[lat_name].point_spacing='uneven'
106            fid.variables[lat_name].units='degrees_north'
107            fid.variables[lat_name].assignValue(latitudes)
108
109            fid.createDimension('TIME',six)
110            fid.createVariable('TIME',netcdf_float,('TIME',))
111            fid.variables['TIME'].point_spacing='uneven'
112            fid.variables['TIME'].units='seconds'
113            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
114
115
116            name = ext[1:3].upper()
117            if name == 'E.': name = 'ELEVATION'
118            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
119            fid.variables[name].units='CENTIMETERS'
120            fid.variables[name].missing_value=-1.e+034
121
122            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
123                                              [-0.1214216, 0, 0, 0],
124                                              [0, 0, 0, 0],
125                                              [0, 0, 0, 0]],
126                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
127                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
128                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
129                                              [0, 2.291054e-005, 2.291054e-005, 0]],
130                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
131                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
132                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
133                                              [0, 0.0001374632, 0.0001374632, 0]],
134                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
135                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
136                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
137                                              [0, 0.0002520159, 0.0002520159, 0]],
138                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
139                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
140                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
141                                              [0, 0.0003665686, 0.0003665686, 0]],
142                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
143                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
144                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
145                                              [0, 0.0004811212, 0.0004811212, 0]]])
146
147
148            fid.close()
149
150
151
152
153    def tearDown(self):
154        import os
155        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
156            #print 'Trying to remove', self.test_MOST_file + ext
157            os.remove(self.test_MOST_file + ext)
158
159    def test_ferret2sww1(self):
160        """Test that georeferencing etc works when converting from
161        ferret format (lat/lon) to sww format (UTM)
162        """
163        from Scientific.IO.NetCDF import NetCDFFile
164        import os, sys
165
166        #The test file has
167        # LON = 150.66667, 150.83334, 151, 151.16667
168        # LAT = -34.5, -34.33333, -34.16667, -34 ;
169        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
170        #
171        # First value (index=0) in small_ha.nc is 0.3400644 cm,
172        # Fourth value (index==3) is -6.50198 cm
173
174
175
176        #Read
177        from anuga.coordinate_transforms.redfearn import redfearn
178        #fid = NetCDFFile(self.test_MOST_file)
179        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
180        first_value = fid.variables['HA'][:][0,0,0]
181        fourth_value = fid.variables['HA'][:][0,0,3]
182        fid.close()
183
184
185        #Call conversion (with zero origin)
186        #ferret2sww('small', verbose=False,
187        #           origin = (56, 0, 0))
188        ferret2sww(self.test_MOST_file, verbose=self.verbose,
189                   origin = (56, 0, 0))
190
191        #Work out the UTM coordinates for first point
192        zone, e, n = redfearn(-34.5, 150.66667)
193        #print zone, e, n
194
195        #Read output file 'small.sww'
196        #fid = NetCDFFile('small.sww')
197        fid = NetCDFFile(self.test_MOST_file + '.sww')
198
199        x = fid.variables['x'][:]
200        y = fid.variables['y'][:]
201
202        #Check that first coordinate is correctly represented
203        assert num.allclose(x[0], e)
204        assert num.allclose(y[0], n)
205
206        #Check first value
207        stage = fid.variables['stage'][:]
208        xmomentum = fid.variables['xmomentum'][:]
209        ymomentum = fid.variables['ymomentum'][:]
210
211        #print ymomentum
212
213        assert num.allclose(stage[0,0], first_value/100)  #Meters
214
215        #Check fourth value
216        assert num.allclose(stage[0,3], fourth_value/100)  #Meters
217
218        fid.close()
219
220        #Cleanup
221        import os
222        os.remove(self.test_MOST_file + '.sww')
223
224
225
226    def test_ferret2sww_zscale(self):
227        """Test that zscale workse
228        """
229        from Scientific.IO.NetCDF import NetCDFFile
230        import os, sys
231
232        #The test file has
233        # LON = 150.66667, 150.83334, 151, 151.16667
234        # LAT = -34.5, -34.33333, -34.16667, -34 ;
235        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
236        #
237        # First value (index=0) in small_ha.nc is 0.3400644 cm,
238        # Fourth value (index==3) is -6.50198 cm
239
240
241        #Read
242        from anuga.coordinate_transforms.redfearn import redfearn
243        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
244        first_value = fid.variables['HA'][:][0,0,0]
245        fourth_value = fid.variables['HA'][:][0,0,3]
246        fid.close()
247
248        #Call conversion (with no scaling)
249        ferret2sww(self.test_MOST_file, verbose=self.verbose,
250                   origin = (56, 0, 0))
251
252        #Work out the UTM coordinates for first point
253        fid = NetCDFFile(self.test_MOST_file + '.sww')
254
255        #Check values
256        stage_1 = fid.variables['stage'][:]
257        xmomentum_1 = fid.variables['xmomentum'][:]
258        ymomentum_1 = fid.variables['ymomentum'][:]
259
260        assert num.allclose(stage_1[0,0], first_value/100)  #Meters
261        assert num.allclose(stage_1[0,3], fourth_value/100)  #Meters
262
263        fid.close()
264
265        #Call conversion (with scaling)
266        ferret2sww(self.test_MOST_file,
267                   zscale = 5,
268                   verbose=self.verbose,
269                   origin = (56, 0, 0))
270
271        #Work out the UTM coordinates for first point
272        fid = NetCDFFile(self.test_MOST_file + '.sww')
273
274        #Check values
275        stage_5 = fid.variables['stage'][:]
276        xmomentum_5 = fid.variables['xmomentum'][:]
277        ymomentum_5 = fid.variables['ymomentum'][:]
278        elevation = fid.variables['elevation'][:]
279
280        assert num.allclose(stage_5[0,0], 5*first_value/100)  #Meters
281        assert num.allclose(stage_5[0,3], 5*fourth_value/100)  #Meters
282
283        assert num.allclose(5*stage_1, stage_5)
284
285        # Momentum will also be changed due to new depth
286
287        depth_1 = stage_1-elevation
288        depth_5 = stage_5-elevation
289
290
291        for i in range(stage_1.shape[0]):
292            for j in range(stage_1.shape[1]):           
293                if depth_1[i,j] > epsilon:
294
295                    scale = depth_5[i,j]/depth_1[i,j]
296                    ref_xmomentum = xmomentum_1[i,j] * scale
297                    ref_ymomentum = ymomentum_1[i,j] * scale
298                   
299                    #print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
300                   
301                    assert num.allclose(xmomentum_5[i,j], ref_xmomentum)
302                    assert num.allclose(ymomentum_5[i,j], ref_ymomentum)
303                   
304       
305
306        fid.close()
307
308
309        #Cleanup
310        import os
311        os.remove(self.test_MOST_file + '.sww')
312
313
314
315    def test_ferret2sww_2(self):
316        """Test that georeferencing etc works when converting from
317        ferret format (lat/lon) to sww format (UTM)
318        """
319        from Scientific.IO.NetCDF import NetCDFFile
320
321        #The test file has
322        # LON = 150.66667, 150.83334, 151, 151.16667
323        # LAT = -34.5, -34.33333, -34.16667, -34 ;
324        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
325        #
326        # First value (index=0) in small_ha.nc is 0.3400644 cm,
327        # Fourth value (index==3) is -6.50198 cm
328
329
330        from anuga.coordinate_transforms.redfearn import redfearn
331        #fid = NetCDFFile('small_ha.nc')
332        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
333
334        #Pick a coordinate and a value
335
336        time_index = 1
337        lat_index = 0
338        lon_index = 2
339
340        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
341        test_time = fid.variables['TIME'][:][time_index]
342        test_lat = fid.variables['LAT'][:][lat_index]
343        test_lon = fid.variables['LON'][:][lon_index]
344
345        linear_point_index = lat_index*4 + lon_index
346        fid.close()
347
348        #Call conversion (with zero origin)
349        ferret2sww(self.test_MOST_file, verbose=self.verbose,
350                   origin = (56, 0, 0))
351
352
353        #Work out the UTM coordinates for test point
354        zone, e, n = redfearn(test_lat, test_lon)
355
356        #Read output file 'small.sww'
357        fid = NetCDFFile(self.test_MOST_file + '.sww')
358
359        x = fid.variables['x'][:]
360        y = fid.variables['y'][:]
361
362        #Check that test coordinate is correctly represented
363        assert num.allclose(x[linear_point_index], e)
364        assert num.allclose(y[linear_point_index], n)
365
366        #Check test value
367        stage = fid.variables['stage'][:]
368
369        assert num.allclose(stage[time_index, linear_point_index], test_value/100)
370
371        fid.close()
372
373        #Cleanup
374        import os
375        os.remove(self.test_MOST_file + '.sww')
376
377
378    def test_ferret2sww_lat_long(self):
379        # Test that min lat long works
380
381        #The test file has
382        # LON = 150.66667, 150.83334, 151, 151.16667
383        # LAT = -34.5, -34.33333, -34.16667, -34 ;
384       
385        #Read
386        from anuga.coordinate_transforms.redfearn import redfearn
387        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
388        first_value = fid.variables['HA'][:][0,0,0]
389        fourth_value = fid.variables['HA'][:][0,0,3]
390        fid.close()
391
392
393        #Call conversion (with zero origin)
394        #ferret2sww('small', verbose=self.verbose,
395        #           origin = (56, 0, 0))
396        ferret2sww(self.test_MOST_file, verbose=self.verbose,
397                   origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
398
399        #Work out the UTM coordinates for first point
400        zone, e, n = redfearn(-34.5, 150.66667)
401        #print zone, e, n
402
403        #Read output file 'small.sww'
404        #fid = NetCDFFile('small.sww')
405        fid = NetCDFFile(self.test_MOST_file + '.sww')
406
407        x = fid.variables['x'][:]
408        y = fid.variables['y'][:]
409        #Check that first coordinate is correctly represented
410        assert 16 == len(x)
411
412        fid.close()
413
414        #Cleanup
415        import os
416        os.remove(self.test_MOST_file + '.sww')
417
418
419    def test_ferret2sww_lat_longII(self):
420        # Test that min lat long works
421
422        #The test file has
423        # LON = 150.66667, 150.83334, 151, 151.16667
424        # LAT = -34.5, -34.33333, -34.16667, -34 ;
425       
426        #Read
427        from anuga.coordinate_transforms.redfearn import redfearn
428        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
429        first_value = fid.variables['HA'][:][0,0,0]
430        fourth_value = fid.variables['HA'][:][0,0,3]
431        fid.close()
432
433
434        #Call conversion (with zero origin)
435        #ferret2sww('small', verbose=False,
436        #           origin = (56, 0, 0))
437        ferret2sww(self.test_MOST_file, verbose=False,
438                   origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
439
440        #Work out the UTM coordinates for first point
441        zone, e, n = redfearn(-34.5, 150.66667)
442        #print zone, e, n
443
444        #Read output file 'small.sww'
445        #fid = NetCDFFile('small.sww')
446        fid = NetCDFFile(self.test_MOST_file + '.sww')
447
448        x = fid.variables['x'][:]
449        y = fid.variables['y'][:]
450        #Check that first coordinate is correctly represented
451        assert 12 == len(x)
452
453        fid.close()
454
455        #Cleanup
456        import os
457        os.remove(self.test_MOST_file + '.sww')
458       
459    def test_ferret2sww3(self):
460        """Elevation included
461        """
462        from Scientific.IO.NetCDF import NetCDFFile
463
464        #The test file has
465        # LON = 150.66667, 150.83334, 151, 151.16667
466        # LAT = -34.5, -34.33333, -34.16667, -34 ;
467        # ELEVATION = [-1 -2 -3 -4
468        #             -5 -6 -7 -8
469        #              ...
470        #              ...      -16]
471        # where the top left corner is -1m,
472        # and the ll corner is -13.0m
473        #
474        # First value (index=0) in small_ha.nc is 0.3400644 cm,
475        # Fourth value (index==3) is -6.50198 cm
476
477        from anuga.coordinate_transforms.redfearn import redfearn
478        import os
479        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
480        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
481        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
482        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
483
484        h1_list = [150.66667,150.83334,151.]
485        h2_list = [-34.5,-34.33333]
486
487        long_name = 'LON'
488        lat_name = 'LAT'
489        time_name = 'TIME'
490
491        nx = 3
492        ny = 2
493
494        for fid in [fid1,fid2,fid3]:
495            fid.createDimension(long_name,nx)
496            fid.createVariable(long_name,netcdf_float,(long_name,))
497            fid.variables[long_name].point_spacing='uneven'
498            fid.variables[long_name].units='degrees_east'
499            fid.variables[long_name].assignValue(h1_list)
500
501            fid.createDimension(lat_name,ny)
502            fid.createVariable(lat_name,netcdf_float,(lat_name,))
503            fid.variables[lat_name].point_spacing='uneven'
504            fid.variables[lat_name].units='degrees_north'
505            fid.variables[lat_name].assignValue(h2_list)
506
507            fid.createDimension(time_name,2)
508            fid.createVariable(time_name,netcdf_float,(time_name,))
509            fid.variables[time_name].point_spacing='uneven'
510            fid.variables[time_name].units='seconds'
511            fid.variables[time_name].assignValue([0.,1.])
512            #if fid == fid3: break
513
514
515        for fid in [fid4]:
516            fid.createDimension(long_name,nx)
517            fid.createVariable(long_name,netcdf_float,(long_name,))
518            fid.variables[long_name].point_spacing='uneven'
519            fid.variables[long_name].units='degrees_east'
520            fid.variables[long_name].assignValue(h1_list)
521
522            fid.createDimension(lat_name,ny)
523            fid.createVariable(lat_name,netcdf_float,(lat_name,))
524            fid.variables[lat_name].point_spacing='uneven'
525            fid.variables[lat_name].units='degrees_north'
526            fid.variables[lat_name].assignValue(h2_list)
527
528        name = {}
529        name[fid1]='HA'
530        name[fid2]='UA'
531        name[fid3]='VA'
532        name[fid4]='ELEVATION'
533
534        units = {}
535        units[fid1]='cm'
536        units[fid2]='cm/s'
537        units[fid3]='cm/s'
538        units[fid4]='m'
539
540        values = {}
541        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
542        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
543        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
544        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
545
546        for fid in [fid1,fid2,fid3]:
547          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
548          fid.variables[name[fid]].point_spacing='uneven'
549          fid.variables[name[fid]].units=units[fid]
550          fid.variables[name[fid]].assignValue(values[fid])
551          fid.variables[name[fid]].missing_value = -99999999.
552          #if fid == fid3: break
553
554        for fid in [fid4]:
555            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
556            fid.variables[name[fid]].point_spacing='uneven'
557            fid.variables[name[fid]].units=units[fid]
558            fid.variables[name[fid]].assignValue(values[fid])
559            fid.variables[name[fid]].missing_value = -99999999.
560
561
562        fid1.sync(); fid1.close()
563        fid2.sync(); fid2.close()
564        fid3.sync(); fid3.close()
565        fid4.sync(); fid4.close()
566
567        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
568        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
569        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
570
571
572        first_amp = fid1.variables['HA'][:][0,0,0]
573        third_amp = fid1.variables['HA'][:][0,0,2]
574        first_elevation = fid2.variables['ELEVATION'][0,0]
575        third_elevation= fid2.variables['ELEVATION'][:][0,2]
576        first_speed = fid3.variables['VA'][0,0,0]
577        third_speed = fid3.variables['VA'][:][0,0,2]
578
579        fid1.close()
580        fid2.close()
581        fid3.close()
582
583        #Call conversion (with zero origin)
584        ferret2sww('test', verbose=self.verbose,
585                   origin = (56, 0, 0), inverted_bathymetry=False)
586
587        os.remove('test_va.nc')
588        os.remove('test_ua.nc')
589        os.remove('test_ha.nc')
590        os.remove('test_e.nc')
591
592        #Read output file 'test.sww'
593        fid = NetCDFFile('test.sww')
594
595
596        #Check first value
597        elevation = fid.variables['elevation'][:]
598        stage = fid.variables['stage'][:]
599        xmomentum = fid.variables['xmomentum'][:]
600        ymomentum = fid.variables['ymomentum'][:]
601
602        #print ymomentum
603        first_height = first_amp/100 - first_elevation
604        third_height = third_amp/100 - third_elevation
605        first_momentum=first_speed*first_height/100
606        third_momentum=third_speed*third_height/100
607
608        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
609        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
610
611        fid.close()
612
613        #Cleanup
614        os.remove('test.sww')
615
616
617
618    def test_ferret2sww4(self):
619        """Like previous but with augmented variable names as
620        in files produced by ferret as opposed to MOST
621        """
622        from Scientific.IO.NetCDF import NetCDFFile
623
624        #The test file has
625        # LON = 150.66667, 150.83334, 151, 151.16667
626        # LAT = -34.5, -34.33333, -34.16667, -34 ;
627        # ELEVATION = [-1 -2 -3 -4
628        #             -5 -6 -7 -8
629        #              ...
630        #              ...      -16]
631        # where the top left corner is -1m,
632        # and the ll corner is -13.0m
633        #
634        # First value (index=0) in small_ha.nc is 0.3400644 cm,
635        # Fourth value (index==3) is -6.50198 cm
636
637        from anuga.coordinate_transforms.redfearn import redfearn
638        import os
639        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
640        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
641        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
642        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
643
644        h1_list = [150.66667,150.83334,151.]
645        h2_list = [-34.5,-34.33333]
646
647#        long_name = 'LON961_1261'
648#        lat_name = 'LAT481_841'
649#        time_name = 'TIME1'
650
651        long_name = 'LON'
652        lat_name = 'LAT'
653        time_name = 'TIME'
654
655        nx = 3
656        ny = 2
657
658        for fid in [fid1,fid2,fid3]:
659            fid.createDimension(long_name,nx)
660            fid.createVariable(long_name,netcdf_float,(long_name,))
661            fid.variables[long_name].point_spacing='uneven'
662            fid.variables[long_name].units='degrees_east'
663            fid.variables[long_name].assignValue(h1_list)
664
665            fid.createDimension(lat_name,ny)
666            fid.createVariable(lat_name,netcdf_float,(lat_name,))
667            fid.variables[lat_name].point_spacing='uneven'
668            fid.variables[lat_name].units='degrees_north'
669            fid.variables[lat_name].assignValue(h2_list)
670
671            fid.createDimension(time_name,2)
672            fid.createVariable(time_name,netcdf_float,(time_name,))
673            fid.variables[time_name].point_spacing='uneven'
674            fid.variables[time_name].units='seconds'
675            fid.variables[time_name].assignValue([0.,1.])
676            #if fid == fid3: break
677
678
679        for fid in [fid4]:
680            fid.createDimension(long_name,nx)
681            fid.createVariable(long_name,netcdf_float,(long_name,))
682            fid.variables[long_name].point_spacing='uneven'
683            fid.variables[long_name].units='degrees_east'
684            fid.variables[long_name].assignValue(h1_list)
685
686            fid.createDimension(lat_name,ny)
687            fid.createVariable(lat_name,netcdf_float,(lat_name,))
688            fid.variables[lat_name].point_spacing='uneven'
689            fid.variables[lat_name].units='degrees_north'
690            fid.variables[lat_name].assignValue(h2_list)
691
692        name = {}
693        name[fid1]='HA'
694        name[fid2]='UA'
695        name[fid3]='VA'
696        name[fid4]='ELEVATION'
697
698        units = {}
699        units[fid1]='cm'
700        units[fid2]='cm/s'
701        units[fid3]='cm/s'
702        units[fid4]='m'
703
704        values = {}
705        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
706        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
707        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
708        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
709
710        for fid in [fid1,fid2,fid3]:
711          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
712          fid.variables[name[fid]].point_spacing='uneven'
713          fid.variables[name[fid]].units=units[fid]
714          fid.variables[name[fid]].assignValue(values[fid])
715          fid.variables[name[fid]].missing_value = -99999999.
716          #if fid == fid3: break
717
718        for fid in [fid4]:
719            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
720            fid.variables[name[fid]].point_spacing='uneven'
721            fid.variables[name[fid]].units=units[fid]
722            fid.variables[name[fid]].assignValue(values[fid])
723            fid.variables[name[fid]].missing_value = -99999999.
724
725
726        fid1.sync(); fid1.close()
727        fid2.sync(); fid2.close()
728        fid3.sync(); fid3.close()
729        fid4.sync(); fid4.close()
730
731        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
732        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
733        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
734
735
736        first_amp = fid1.variables['HA'][:][0,0,0]
737        third_amp = fid1.variables['HA'][:][0,0,2]
738        first_elevation = fid2.variables['ELEVATION'][0,0]
739        third_elevation= fid2.variables['ELEVATION'][:][0,2]
740        first_speed = fid3.variables['VA'][0,0,0]
741        third_speed = fid3.variables['VA'][:][0,0,2]
742
743        fid1.close()
744        fid2.close()
745        fid3.close()
746
747        #Call conversion (with zero origin)
748        ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
749                   , inverted_bathymetry=False)
750
751        os.remove('test_va.nc')
752        os.remove('test_ua.nc')
753        os.remove('test_ha.nc')
754        os.remove('test_e.nc')
755
756        #Read output file 'test.sww'
757        fid = NetCDFFile('test.sww')
758
759
760        #Check first value
761        elevation = fid.variables['elevation'][:]
762        stage = fid.variables['stage'][:]
763        xmomentum = fid.variables['xmomentum'][:]
764        ymomentum = fid.variables['ymomentum'][:]
765
766        #print ymomentum
767        first_height = first_amp/100 - first_elevation
768        third_height = third_amp/100 - third_elevation
769        first_momentum=first_speed*first_height/100
770        third_momentum=third_speed*third_height/100
771
772        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
773        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
774
775        fid.close()
776
777        #Cleanup
778        os.remove('test.sww')
779
780
781
782
783    def test_ferret2sww_nz_origin(self):
784        from Scientific.IO.NetCDF import NetCDFFile
785        from anuga.coordinate_transforms.redfearn import redfearn
786
787        #Call conversion (with nonzero origin)
788        ferret2sww(self.test_MOST_file, verbose=self.verbose,
789                   origin = (56, 100000, 200000))
790
791
792        #Work out the UTM coordinates for first point
793        zone, e, n = redfearn(-34.5, 150.66667)
794
795        #Read output file 'small.sww'
796        #fid = NetCDFFile('small.sww', netcdf_mode_r)
797        fid = NetCDFFile(self.test_MOST_file + '.sww')
798
799        x = fid.variables['x'][:]
800        y = fid.variables['y'][:]
801
802        #Check that first coordinate is correctly represented
803        assert num.allclose(x[0], e-100000)
804        assert num.allclose(y[0], n-200000)
805
806        fid.close()
807
808        #Cleanup
809        os.remove(self.test_MOST_file + '.sww')
810
811
812    def test_ferret2sww_lat_longII(self):
813        # Test that min lat long works
814
815        #The test file has
816        # LON = 150.66667, 150.83334, 151, 151.16667
817        # LAT = -34.5, -34.33333, -34.16667, -34 ;
818       
819        #Read
820        from anuga.coordinate_transforms.redfearn import redfearn
821        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
822        first_value = fid.variables['HA'][:][0,0,0]
823        fourth_value = fid.variables['HA'][:][0,0,3]
824        fid.close()
825
826
827        #Call conversion (with zero origin)
828        #ferret2sww('small', verbose=self.verbose,
829        #           origin = (56, 0, 0))
830        try:
831            ferret2sww(self.test_MOST_file, verbose=self.verbose,
832                   origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
833        except AssertionError:
834            pass
835        else:
836            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
837
838    def test_sww_extent(self):
839        """Not a test, rather a look at the sww format
840        """
841
842        import time, os
843        from Scientific.IO.NetCDF import NetCDFFile
844
845        self.domain.set_name('datatest' + str(id(self)))
846        self.domain.format = 'sww'
847        self.domain.smooth = True
848        self.domain.reduction = mean
849        self.domain.set_datadir('.')
850        #self.domain.tight_slope_limiters = 1       
851
852
853        sww = SWW_file(self.domain)
854        sww.store_connectivity()
855        sww.store_timestep()
856        self.domain.time = 2.
857
858        #Modify stage at second timestep
859        stage = self.domain.quantities['stage'].vertex_values
860        self.domain.set_quantity('stage', stage/2)
861
862        sww.store_timestep()
863
864        file_and_extension_name = self.domain.get_name() + ".sww"
865        #print "file_and_extension_name",file_and_extension_name
866        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
867               extent_sww(file_and_extension_name )
868
869        assert num.allclose(xmin, 0.0)
870        assert num.allclose(xmax, 1.0)
871        assert num.allclose(ymin, 0.0)
872        assert num.allclose(ymax, 1.0)
873
874        # FIXME (Ole): Revisit these numbers
875        #assert num.allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
876        #assert num.allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
877
878
879        #Cleanup
880        os.remove(sww.filename)
881
882
883
884#-------------------------------------------------------------
885
886if __name__ == "__main__":
887    suite = unittest.makeSuite(Test_File_Conversion, 'test_sww')
888   
889    # FIXME(Ole): When Ross has implemented logging, we can
890    # probably get rid of all this:
891    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
892        Test_File_Conversion.verbose=True
893        saveout = sys.stdout   
894        filename = ".temp_verbose"
895        fid = open(filename, 'w')
896        sys.stdout = fid
897    else:
898        pass
899    runner = unittest.TextTestRunner() #verbosity=2)
900    runner.run(suite)
901
902    # Cleaning up
903    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
904        sys.stdout = saveout
905        fid.close() 
906        os.remove(filename)
907
908
909    def test_asc_csiro2sww(self):
910        import tempfile
911
912        bath_dir = tempfile.mkdtemp()
913        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
914        #bath_dir = 'bath_data_manager_test'
915        #print "os.getcwd( )",os.getcwd( )
916        elevation_dir =  tempfile.mkdtemp()
917        #elevation_dir = 'elev_expanded'
918        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
919        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
920
921        fid = open(bath_dir_filename, 'w')
922        fid.write(""" ncols             3
923 nrows             2
924 xllcorner    148.00000
925 yllcorner    -38.00000
926 cellsize       0.25
927 nodata_value   -9999.0
928    9000.000    -1000.000    3000.0
929   -1000.000    9000.000  -1000.000
930""")
931        fid.close()
932
933        fid = open(elevation_dir_filename1, 'w')
934        fid.write(""" ncols             3
935 nrows             2
936 xllcorner    148.00000
937 yllcorner    -38.00000
938 cellsize       0.25
939 nodata_value   -9999.0
940    9000.000    0.000    3000.0
941     0.000     9000.000     0.000
942""")
943        fid.close()
944
945        fid = open(elevation_dir_filename2, 'w')
946        fid.write(""" ncols             3
947 nrows             2
948 xllcorner    148.00000
949 yllcorner    -38.00000
950 cellsize       0.25
951 nodata_value   -9999.0
952    9000.000    4000.000    4000.0
953    4000.000    9000.000    4000.000
954""")
955        fid.close()
956
957        ucur_dir =  tempfile.mkdtemp()
958        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
959        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
960
961        fid = open(ucur_dir_filename1, 'w')
962        fid.write(""" ncols             3
963 nrows             2
964 xllcorner    148.00000
965 yllcorner    -38.00000
966 cellsize       0.25
967 nodata_value   -9999.0
968    90.000    60.000    30.0
969    10.000    10.000    10.000
970""")
971        fid.close()
972        fid = open(ucur_dir_filename2, 'w')
973        fid.write(""" ncols             3
974 nrows             2
975 xllcorner    148.00000
976 yllcorner    -38.00000
977 cellsize       0.25
978 nodata_value   -9999.0
979    90.000    60.000    30.0
980    10.000    10.000    10.000
981""")
982        fid.close()
983
984        vcur_dir =  tempfile.mkdtemp()
985        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
986        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
987
988        fid = open(vcur_dir_filename1, 'w')
989        fid.write(""" ncols             3
990 nrows             2
991 xllcorner    148.00000
992 yllcorner    -38.00000
993 cellsize       0.25
994 nodata_value   -9999.0
995    90.000    60.000    30.0
996    10.000    10.000    10.000
997""")
998        fid.close()
999        fid = open(vcur_dir_filename2, 'w')
1000        fid.write(""" ncols             3
1001 nrows             2
1002 xllcorner    148.00000
1003 yllcorner    -38.00000
1004 cellsize       0.25
1005 nodata_value   -9999.0
1006    90.000    60.000    30.0
1007    10.000    10.000    10.000
1008""")
1009        fid.close()
1010
1011        sww_file = 'a_test.sww'
1012        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
1013                      verbose=self.verbose)
1014
1015        # check the sww file
1016
1017        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
1018        x = fid.variables['x'][:]
1019        y = fid.variables['y'][:]
1020        z = fid.variables['elevation'][:]
1021        stage = fid.variables['stage'][:]
1022        xmomentum = fid.variables['xmomentum'][:]
1023        geo_ref = Geo_reference(NetCDFObject=fid)
1024        #print "geo_ref",geo_ref
1025        x_ref = geo_ref.get_xllcorner()
1026        y_ref = geo_ref.get_yllcorner()
1027        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
1028        assert num.allclose(x_ref, 587798.418) # (-38, 148)
1029        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
1030
1031        #Zone:   55
1032        #Easting:  588095.674  Northing: 5821451.722
1033        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
1034        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
1035
1036        #Zone:   55
1037        #Easting:  632145.632  Northing: 5820863.269
1038        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
1039        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
1040
1041        #Zone:   55
1042        #Easting:  609748.788  Northing: 5793447.860
1043        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
1044        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
1045
1046        assert num.allclose(z[0],9000.0 )
1047        assert num.allclose(stage[0][1],0.0 )
1048
1049        #(4000+1000)*60
1050        assert num.allclose(xmomentum[1][1],300000.0 )
1051
1052
1053        fid.close()
1054
1055        #tidy up
1056        os.remove(bath_dir_filename)
1057        os.rmdir(bath_dir)
1058
1059        os.remove(elevation_dir_filename1)
1060        os.remove(elevation_dir_filename2)
1061        os.rmdir(elevation_dir)
1062
1063        os.remove(ucur_dir_filename1)
1064        os.remove(ucur_dir_filename2)
1065        os.rmdir(ucur_dir)
1066
1067        os.remove(vcur_dir_filename1)
1068        os.remove(vcur_dir_filename2)
1069        os.rmdir(vcur_dir)
1070
1071
1072        # remove sww file
1073        os.remove(sww_file)
1074
1075    def test_asc_csiro2sww2(self):
1076        import tempfile
1077
1078        bath_dir = tempfile.mkdtemp()
1079        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
1080        #bath_dir = 'bath_data_manager_test'
1081        #print "os.getcwd( )",os.getcwd( )
1082        elevation_dir =  tempfile.mkdtemp()
1083        #elevation_dir = 'elev_expanded'
1084        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
1085        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
1086
1087        fid = open(bath_dir_filename, 'w')
1088        fid.write(""" ncols             3
1089 nrows             2
1090 xllcorner    148.00000
1091 yllcorner    -38.00000
1092 cellsize       0.25
1093 nodata_value   -9999.0
1094    9000.000    -1000.000    3000.0
1095   -1000.000    9000.000  -1000.000
1096""")
1097        fid.close()
1098
1099        fid = open(elevation_dir_filename1, 'w')
1100        fid.write(""" ncols             3
1101 nrows             2
1102 xllcorner    148.00000
1103 yllcorner    -38.00000
1104 cellsize       0.25
1105 nodata_value   -9999.0
1106    9000.000    0.000    3000.0
1107     0.000     -9999.000     -9999.000
1108""")
1109        fid.close()
1110
1111        fid = open(elevation_dir_filename2, 'w')
1112        fid.write(""" ncols             3
1113 nrows             2
1114 xllcorner    148.00000
1115 yllcorner    -38.00000
1116 cellsize       0.25
1117 nodata_value   -9999.0
1118    9000.000    4000.000    4000.0
1119    4000.000    9000.000    4000.000
1120""")
1121        fid.close()
1122
1123        ucur_dir =  tempfile.mkdtemp()
1124        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
1125        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
1126
1127        fid = open(ucur_dir_filename1, 'w')
1128        fid.write(""" ncols             3
1129 nrows             2
1130 xllcorner    148.00000
1131 yllcorner    -38.00000
1132 cellsize       0.25
1133 nodata_value   -9999.0
1134    90.000    60.000    30.0
1135    10.000    10.000    10.000
1136""")
1137        fid.close()
1138        fid = open(ucur_dir_filename2, 'w')
1139        fid.write(""" ncols             3
1140 nrows             2
1141 xllcorner    148.00000
1142 yllcorner    -38.00000
1143 cellsize       0.25
1144 nodata_value   -9999.0
1145    90.000    60.000    30.0
1146    10.000    10.000    10.000
1147""")
1148        fid.close()
1149
1150        vcur_dir =  tempfile.mkdtemp()
1151        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
1152        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
1153
1154        fid = open(vcur_dir_filename1, 'w')
1155        fid.write(""" ncols             3
1156 nrows             2
1157 xllcorner    148.00000
1158 yllcorner    -38.00000
1159 cellsize       0.25
1160 nodata_value   -9999.0
1161    90.000    60.000    30.0
1162    10.000    10.000    10.000
1163""")
1164        fid.close()
1165        fid = open(vcur_dir_filename2, 'w')
1166        fid.write(""" ncols             3
1167 nrows             2
1168 xllcorner    148.00000
1169 yllcorner    -38.00000
1170 cellsize       0.25
1171 nodata_value   -9999.0
1172    90.000    60.000    30.0
1173    10.000    10.000    10.000
1174""")
1175        fid.close()
1176
1177        try:
1178            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
1179                          vcur_dir, sww_file,
1180                      verbose=self.verbose)
1181        except:
1182            #tidy up
1183            os.remove(bath_dir_filename)
1184            os.rmdir(bath_dir)
1185
1186            os.remove(elevation_dir_filename1)
1187            os.remove(elevation_dir_filename2)
1188            os.rmdir(elevation_dir)
1189
1190            os.remove(ucur_dir_filename1)
1191            os.remove(ucur_dir_filename2)
1192            os.rmdir(ucur_dir)
1193
1194            os.remove(vcur_dir_filename1)
1195            os.remove(vcur_dir_filename2)
1196            os.rmdir(vcur_dir)
1197        else:
1198            #tidy up
1199            os.remove(bath_dir_filename)
1200            os.rmdir(bath_dir)
1201
1202            os.remove(elevation_dir_filename1)
1203            os.remove(elevation_dir_filename2)
1204            os.rmdir(elevation_dir)
1205            raise 'Should raise exception'
1206
1207            os.remove(ucur_dir_filename1)
1208            os.remove(ucur_dir_filename2)
1209            os.rmdir(ucur_dir)
1210
1211            os.remove(vcur_dir_filename1)
1212            os.remove(vcur_dir_filename2)
1213            os.rmdir(vcur_dir)
1214
1215
1216
1217    def test_asc_csiro2sww3(self):
1218        import tempfile
1219
1220        bath_dir = tempfile.mkdtemp()
1221        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
1222        #bath_dir = 'bath_data_manager_test'
1223        #print "os.getcwd( )",os.getcwd( )
1224        elevation_dir =  tempfile.mkdtemp()
1225        #elevation_dir = 'elev_expanded'
1226        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
1227        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
1228
1229        fid = open(bath_dir_filename, 'w')
1230        fid.write(""" ncols             3
1231 nrows             2
1232 xllcorner    148.00000
1233 yllcorner    -38.00000
1234 cellsize       0.25
1235 nodata_value   -9999.0
1236    9000.000    -1000.000    3000.0
1237   -1000.000    9000.000  -1000.000
1238""")
1239        fid.close()
1240
1241        fid = open(elevation_dir_filename1, 'w')
1242        fid.write(""" ncols             3
1243 nrows             2
1244 xllcorner    148.00000
1245 yllcorner    -38.00000
1246 cellsize       0.25
1247 nodata_value   -9999.0
1248    9000.000    0.000    3000.0
1249     0.000     -9999.000     -9999.000
1250""")
1251        fid.close()
1252
1253        fid = open(elevation_dir_filename2, 'w')
1254        fid.write(""" ncols             3
1255 nrows             2
1256 xllcorner    148.00000
1257 yllcorner    -38.00000
1258 cellsize       0.25
1259 nodata_value   -9999.0
1260    9000.000    4000.000    4000.0
1261    4000.000    9000.000    4000.000
1262""")
1263        fid.close()
1264
1265        ucur_dir =  tempfile.mkdtemp()
1266        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
1267        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
1268
1269        fid = open(ucur_dir_filename1, 'w')
1270        fid.write(""" ncols             3
1271 nrows             2
1272 xllcorner    148.00000
1273 yllcorner    -38.00000
1274 cellsize       0.25
1275 nodata_value   -9999.0
1276    90.000    60.000    30.0
1277    10.000    10.000    10.000
1278""")
1279        fid.close()
1280        fid = open(ucur_dir_filename2, 'w')
1281        fid.write(""" ncols             3
1282 nrows             2
1283 xllcorner    148.00000
1284 yllcorner    -38.00000
1285 cellsize       0.25
1286 nodata_value   -9999.0
1287    90.000    60.000    30.0
1288    10.000    10.000    10.000
1289""")
1290        fid.close()
1291
1292        vcur_dir =  tempfile.mkdtemp()
1293        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
1294        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
1295
1296        fid = open(vcur_dir_filename1, 'w')
1297        fid.write(""" ncols             3
1298 nrows             2
1299 xllcorner    148.00000
1300 yllcorner    -38.00000
1301 cellsize       0.25
1302 nodata_value   -9999.0
1303    90.000    60.000    30.0
1304    10.000    10.000    10.000
1305""")
1306        fid.close()
1307        fid = open(vcur_dir_filename2, 'w')
1308        fid.write(""" ncols             3
1309 nrows             2
1310 xllcorner    148.00000
1311 yllcorner    -38.00000
1312 cellsize       0.25
1313 nodata_value   -9999.0
1314    90.000    60.000    30.0
1315    10.000    10.000    10.000
1316""")
1317        fid.close()
1318
1319        sww_file = 'a_test.sww'
1320        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
1321                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
1322                      mean_stage = 100,
1323                      verbose=self.verbose)
1324
1325        # check the sww file
1326
1327        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
1328        x = fid.variables['x'][:]
1329        y = fid.variables['y'][:]
1330        z = fid.variables['elevation'][:]
1331        stage = fid.variables['stage'][:]
1332        xmomentum = fid.variables['xmomentum'][:]
1333        geo_ref = Geo_reference(NetCDFObject=fid)
1334        #print "geo_ref",geo_ref
1335        x_ref = geo_ref.get_xllcorner()
1336        y_ref = geo_ref.get_yllcorner()
1337        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
1338        assert num.allclose(x_ref, 587798.418) # (-38, 148)
1339        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
1340
1341        #Zone:   55
1342        #Easting:  588095.674  Northing: 5821451.722
1343        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
1344        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
1345
1346        #Zone:   55
1347        #Easting:  632145.632  Northing: 5820863.269
1348        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
1349        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
1350
1351        #Zone:   55
1352        #Easting:  609748.788  Northing: 5793447.860
1353        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
1354        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
1355
1356        assert num.allclose(z[0],9000.0 )
1357        assert num.allclose(stage[0][4],100.0 )
1358        assert num.allclose(stage[0][5],100.0 )
1359
1360        #(100.0 - 9000)*10
1361        assert num.allclose(xmomentum[0][4], -89000.0 )
1362
1363        #(100.0 - -1000.000)*10
1364        assert num.allclose(xmomentum[0][5], 11000.0 )
1365
1366        fid.close()
1367
1368        #tidy up
1369        os.remove(bath_dir_filename)
1370        os.rmdir(bath_dir)
1371
1372        os.remove(elevation_dir_filename1)
1373        os.remove(elevation_dir_filename2)
1374        os.rmdir(elevation_dir)
1375
1376        os.remove(ucur_dir_filename1)
1377        os.remove(ucur_dir_filename2)
1378        os.rmdir(ucur_dir)
1379
1380        os.remove(vcur_dir_filename1)
1381        os.remove(vcur_dir_filename2)
1382        os.rmdir(vcur_dir)
1383
1384        # remove sww file
1385        os.remove(sww_file)
1386
1387
1388    def test_asc_csiro2sww4(self):
1389        """
1390        Test specifying the extent
1391        """
1392
1393        import tempfile
1394
1395        bath_dir = tempfile.mkdtemp()
1396        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
1397        #bath_dir = 'bath_data_manager_test'
1398        #print "os.getcwd( )",os.getcwd( )
1399        elevation_dir =  tempfile.mkdtemp()
1400        #elevation_dir = 'elev_expanded'
1401        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
1402        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
1403
1404        fid = open(bath_dir_filename, 'w')
1405        fid.write(""" ncols             4
1406 nrows             4
1407 xllcorner    148.00000
1408 yllcorner    -38.00000
1409 cellsize       0.25
1410 nodata_value   -9999.0
1411   -9000.000    -1000.000   -3000.0 -2000.000
1412   -1000.000    9000.000  -1000.000 -3000.000
1413   -4000.000    6000.000   2000.000 -5000.000
1414   -9000.000    -1000.000   -3000.0 -2000.000
1415""")
1416        fid.close()
1417
1418        fid = open(elevation_dir_filename1, 'w')
1419        fid.write(""" ncols             4
1420 nrows             4
1421 xllcorner    148.00000
1422 yllcorner    -38.00000
1423 cellsize       0.25
1424 nodata_value   -9999.0
1425   -900.000    -100.000   -300.0 -200.000
1426   -100.000    900.000  -100.000 -300.000
1427   -400.000    600.000   200.000 -500.000
1428   -900.000    -100.000   -300.0 -200.000
1429""")
1430        fid.close()
1431
1432        fid = open(elevation_dir_filename2, 'w')
1433        fid.write(""" ncols             4
1434 nrows             4
1435 xllcorner    148.00000
1436 yllcorner    -38.00000
1437 cellsize       0.25
1438 nodata_value   -9999.0
1439   -990.000    -110.000   -330.0 -220.000
1440   -110.000    990.000  -110.000 -330.000
1441   -440.000    660.000   220.000 -550.000
1442   -990.000    -110.000   -330.0 -220.000
1443""")
1444        fid.close()
1445
1446        ucur_dir =  tempfile.mkdtemp()
1447        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
1448        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
1449
1450        fid = open(ucur_dir_filename1, 'w')
1451        fid.write(""" ncols             4
1452 nrows             4
1453 xllcorner    148.00000
1454 yllcorner    -38.00000
1455 cellsize       0.25
1456 nodata_value   -9999.0
1457   -90.000    -10.000   -30.0 -20.000
1458   -10.000    90.000  -10.000 -30.000
1459   -40.000    60.000   20.000 -50.000
1460   -90.000    -10.000   -30.0 -20.000
1461""")
1462        fid.close()
1463        fid = open(ucur_dir_filename2, 'w')
1464        fid.write(""" ncols             4
1465 nrows             4
1466 xllcorner    148.00000
1467 yllcorner    -38.00000
1468 cellsize       0.25
1469 nodata_value   -9999.0
1470   -90.000    -10.000   -30.0 -20.000
1471   -10.000    99.000  -11.000 -30.000
1472   -40.000    66.000   22.000 -50.000
1473   -90.000    -10.000   -30.0 -20.000
1474""")
1475        fid.close()
1476
1477        vcur_dir =  tempfile.mkdtemp()
1478        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
1479        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
1480
1481        fid = open(vcur_dir_filename1, 'w')
1482        fid.write(""" ncols             4
1483 nrows             4
1484 xllcorner    148.00000
1485 yllcorner    -38.00000
1486 cellsize       0.25
1487 nodata_value   -9999.0
1488   -90.000    -10.000   -30.0 -20.000
1489   -10.000    80.000  -20.000 -30.000
1490   -40.000    50.000   10.000 -50.000
1491   -90.000    -10.000   -30.0 -20.000
1492""")
1493        fid.close()
1494        fid = open(vcur_dir_filename2, 'w')
1495        fid.write(""" ncols             4
1496 nrows             4
1497 xllcorner    148.00000
1498 yllcorner    -38.00000
1499 cellsize       0.25
1500 nodata_value   -9999.0
1501   -90.000    -10.000   -30.0 -20.000
1502   -10.000    88.000  -22.000 -30.000
1503   -40.000    55.000   11.000 -50.000
1504   -90.000    -10.000   -30.0 -20.000
1505""")
1506        fid.close()
1507
1508        sww_file = tempfile.mktemp(".sww")
1509        #sww_file = 'a_test.sww'
1510        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
1511                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
1512                      mean_stage = 100,
1513                       minlat = -37.6, maxlat = -37.6,
1514                  minlon = 148.3, maxlon = 148.3,
1515                      verbose=self.verbose
1516                      #,verbose = True
1517                      )
1518
1519        # check the sww file
1520
1521        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
1522        x = fid.variables['x'][:]
1523        y = fid.variables['y'][:]
1524        z = fid.variables['elevation'][:]
1525        stage = fid.variables['stage'][:]
1526        xmomentum = fid.variables['xmomentum'][:]
1527        ymomentum = fid.variables['ymomentum'][:]
1528        geo_ref = Geo_reference(NetCDFObject=fid)
1529        #print "geo_ref",geo_ref
1530        x_ref = geo_ref.get_xllcorner()
1531        y_ref = geo_ref.get_yllcorner()
1532        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
1533
1534        assert num.allclose(fid.starttime, 0.0) # (-37.45, 148.25)
1535        assert num.allclose(x_ref, 610120.388) # (-37.45, 148.25)
1536        assert num.allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
1537
1538        #Easting:  632145.632  Northing: 5820863.269
1539        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
1540
1541        #print "x",x
1542        #print "y",y
1543        self.failUnless(len(x) == 4,'failed') # 2*2
1544        self.failUnless(len(x) == 4,'failed') # 2*2
1545
1546        #Zone:   55
1547        #Easting:  632145.632  Northing: 5820863.269
1548        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
1549        # magic number - y is close enough for me.
1550        assert num.allclose(x[3], 632145.63 - x_ref)
1551        assert num.allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
1552
1553        assert num.allclose(z[0],9000.0 ) #z is elevation info
1554        #print "z",z
1555        # 2 time steps, 4 points
1556        self.failUnless(xmomentum.shape == (2,4), 'failed')
1557        self.failUnless(ymomentum.shape == (2,4), 'failed')
1558
1559        #(100.0 - -1000.000)*10
1560        #assert num.allclose(xmomentum[0][5], 11000.0 )
1561
1562        fid.close()
1563
1564        # is the sww file readable?
1565        #Lets see if we can convert it to a dem!
1566        # if you uncomment, remember to delete the file
1567        #print "sww_file",sww_file
1568        #dem_file = tempfile.mktemp(".dem")
1569        domain = sww2domain(sww_file) ###, dem_file)
1570        domain.check_integrity()
1571
1572        #tidy up
1573        os.remove(bath_dir_filename)
1574        os.rmdir(bath_dir)
1575
1576        os.remove(elevation_dir_filename1)
1577        os.remove(elevation_dir_filename2)
1578        os.rmdir(elevation_dir)
1579
1580        os.remove(ucur_dir_filename1)
1581        os.remove(ucur_dir_filename2)
1582        os.rmdir(ucur_dir)
1583
1584        os.remove(vcur_dir_filename1)
1585        os.remove(vcur_dir_filename2)
1586        os.rmdir(vcur_dir)
1587
1588        # remove sww file
1589        os.remove(sww_file)
1590
1591
1592
1593    def test_get_min_max_indexes(self):
1594        latitudes = [3,2,1,0]
1595        longitudes = [0,10,20,30]
1596
1597        # k - lat
1598        # l - lon
1599        kmin, kmax, lmin, lmax = get_min_max_indices(
1600            latitudes,longitudes,
1601            -10,4,-10,31)
1602
1603        #print "kmin",kmin;print "kmax",kmax
1604        #print "lmin",lmin;print "lmax",lmax
1605        latitudes_new = latitudes[kmin:kmax]
1606        longitudes_news = longitudes[lmin:lmax]
1607        #print "latitudes_new", latitudes_new
1608        #print "longitudes_news",longitudes_news
1609        self.failUnless(latitudes == latitudes_new and \
1610                        longitudes == longitudes_news,
1611                         'failed')
1612
1613        ## 2nd test
1614        kmin, kmax, lmin, lmax = get_min_max_indices(
1615            latitudes,longitudes,
1616            0.5,2.5,5,25)
1617        #print "kmin",kmin;print "kmax",kmax
1618        #print "lmin",lmin;print "lmax",lmax
1619        latitudes_new = latitudes[kmin:kmax]
1620        longitudes_news = longitudes[lmin:lmax]
1621        #print "latitudes_new", latitudes_new
1622        #print "longitudes_news",longitudes_news
1623
1624        self.failUnless(latitudes == latitudes_new and \
1625                        longitudes == longitudes_news,
1626                         'failed')
1627
1628        ## 3rd test
1629        kmin, kmax, lmin, lmax = get_min_max_indices(\
1630            latitudes,
1631            longitudes,
1632            1.1,1.9,12,17)
1633        #print "kmin",kmin;print "kmax",kmax
1634        #print "lmin",lmin;print "lmax",lmax
1635        latitudes_new = latitudes[kmin:kmax]
1636        longitudes_news = longitudes[lmin:lmax]
1637        #print "latitudes_new", latitudes_new
1638        #print "longitudes_news",longitudes_news
1639
1640        self.failUnless(latitudes_new == [2, 1] and \
1641                        longitudes_news == [10, 20],
1642                         'failed')
1643
1644
1645        ## 4th test
1646        kmin, kmax, lmin, lmax = get_min_max_indices(
1647            latitudes,longitudes,
1648                                                      -0.1,1.9,-2,17)
1649        #print "kmin",kmin;print "kmax",kmax
1650        #print "lmin",lmin;print "lmax",lmax
1651        latitudes_new = latitudes[kmin:kmax]
1652        longitudes_news = longitudes[lmin:lmax]
1653        #print "latitudes_new", latitudes_new
1654        #print "longitudes_news",longitudes_news
1655
1656        self.failUnless(latitudes_new == [2, 1, 0] and \
1657                        longitudes_news == [0, 10, 20],
1658                         'failed')
1659        ## 5th test
1660        kmin, kmax, lmin, lmax = get_min_max_indices(
1661            latitudes,longitudes,
1662            0.1,1.9,2,17)
1663        #print "kmin",kmin;print "kmax",kmax
1664        #print "lmin",lmin;print "lmax",lmax
1665        latitudes_new = latitudes[kmin:kmax]
1666        longitudes_news = longitudes[lmin:lmax]
1667        #print "latitudes_new", latitudes_new
1668        #print "longitudes_news",longitudes_news
1669
1670        self.failUnless(latitudes_new == [2, 1, 0] and \
1671                        longitudes_news == [0, 10, 20],
1672                         'failed')
1673
1674        ## 6th test
1675
1676        kmin, kmax, lmin, lmax = get_min_max_indices(
1677            latitudes,longitudes,
1678            1.5,4,18,32)
1679        #print "kmin",kmin;print "kmax",kmax
1680        #print "lmin",lmin;print "lmax",lmax
1681        latitudes_new = latitudes[kmin:kmax]
1682        longitudes_news = longitudes[lmin:lmax]
1683        #print "latitudes_new", latitudes_new
1684        #print "longitudes_news",longitudes_news
1685
1686        self.failUnless(latitudes_new == [3, 2, 1] and \
1687                        longitudes_news == [10, 20, 30],
1688                         'failed')
1689
1690
1691        ## 7th test
1692        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
1693        kmin, kmax, lmin, lmax = get_min_max_indices(
1694            latitudes,longitudes,
1695            1.5,1.5,15,15)
1696        #print "kmin",kmin;print "kmax",kmax
1697        #print "lmin",lmin;print "lmax",lmax
1698        latitudes_new = latitudes[kmin:kmax]
1699        longitudes_news = longitudes[lmin:lmax]
1700        m2d = m2d[kmin:kmax,lmin:lmax]
1701        #print "m2d", m2d
1702        #print "latitudes_new", latitudes_new
1703        #print "longitudes_news",longitudes_news
1704
1705        self.failUnless(num.alltrue(latitudes_new == [2, 1]) and 
1706                        num.alltrue(longitudes_news == [10, 20]),
1707                        'failed')
1708
1709        self.failUnless(num.alltrue(m2d == [[5,6],[9,10]]), 'failed')
1710
1711    def test_get_min_max_indexes_lat_ascending(self):
1712        latitudes = [0,1,2,3]
1713        longitudes = [0,10,20,30]
1714
1715        # k - lat
1716        # l - lon
1717        kmin, kmax, lmin, lmax = get_min_max_indices(
1718            latitudes,longitudes,
1719            -10,4,-10,31)
1720
1721        #print "kmin",kmin;print "kmax",kmax
1722        #print "lmin",lmin;print "lmax",lmax
1723        latitudes_new = latitudes[kmin:kmax]
1724        longitudes_news = longitudes[lmin:lmax]
1725        #print "latitudes_new", latitudes_new
1726        #print "longitudes_news",longitudes_news
1727        self.failUnless(latitudes == latitudes_new and \
1728                        longitudes == longitudes_news,
1729                         'failed')
1730
1731        ## 3rd test
1732        kmin, kmax, lmin, lmax = get_min_max_indices(\
1733            latitudes,
1734            longitudes,
1735            1.1,1.9,12,17)
1736        #print "kmin",kmin;print "kmax",kmax
1737        #print "lmin",lmin;print "lmax",lmax
1738        latitudes_new = latitudes[kmin:kmax]
1739        longitudes_news = longitudes[lmin:lmax]
1740        #print "latitudes_new", latitudes_new
1741        #print "longitudes_news",longitudes_news
1742
1743        self.failUnless(latitudes_new == [1, 2] and \
1744                        longitudes_news == [10, 20],
1745                         'failed')
1746
1747    def test_get_min_max_indexes2(self):
1748        latitudes = [-30,-35,-40,-45]
1749        longitudes = [148,149,150,151]
1750
1751        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
1752
1753        # k - lat
1754        # l - lon
1755        kmin, kmax, lmin, lmax = get_min_max_indices(
1756            latitudes,longitudes,
1757            -37,-27,147,149.5)
1758
1759        #print "kmin",kmin;print "kmax",kmax
1760        #print "lmin",lmin;print "lmax",lmax
1761        #print "m2d", m2d
1762        #print "latitudes", latitudes
1763        #print "longitudes",longitudes
1764        #print "latitudes[kmax]", latitudes[kmax]
1765        latitudes_new = latitudes[kmin:kmax]
1766        longitudes_new = longitudes[lmin:lmax]
1767        m2d = m2d[kmin:kmax,lmin:lmax]
1768        #print "m2d", m2d
1769        #print "latitudes_new", latitudes_new
1770        #print "longitudes_new",longitudes_new
1771
1772        self.failUnless(latitudes_new == [-30, -35, -40] and
1773                        longitudes_new == [148, 149,150],
1774                        'failed')
1775        self.failUnless(num.alltrue(m2d == [[0,1,2],[4,5,6],[8,9,10]]),
1776                        'failed')
1777
1778    def test_get_min_max_indexes3(self):
1779        latitudes = [-30,-35,-40,-45,-50,-55,-60]
1780        longitudes = [148,149,150,151]
1781
1782        # k - lat
1783        # l - lon
1784        kmin, kmax, lmin, lmax = get_min_max_indices(
1785            latitudes,longitudes,
1786            -43,-37,148.5,149.5)
1787
1788
1789        #print "kmin",kmin;print "kmax",kmax
1790        #print "lmin",lmin;print "lmax",lmax
1791        #print "latitudes", latitudes
1792        #print "longitudes",longitudes
1793        latitudes_new = latitudes[kmin:kmax]
1794        longitudes_news = longitudes[lmin:lmax]
1795        #print "latitudes_new", latitudes_new
1796        #print "longitudes_news",longitudes_news
1797
1798        self.failUnless(latitudes_new == [-35, -40, -45] and
1799                        longitudes_news == [148, 149,150],
1800                         'failed')
1801
1802    def test_get_min_max_indexes4(self):
1803        latitudes = [-30,-35,-40,-45,-50,-55,-60]
1804        longitudes = [148,149,150,151]
1805
1806        # k - lat
1807        # l - lon
1808        kmin, kmax, lmin, lmax = get_min_max_indices(
1809            latitudes,longitudes)
1810
1811        #print "kmin",kmin;print "kmax",kmax
1812        #print "lmin",lmin;print "lmax",lmax
1813        #print "latitudes", latitudes
1814        #print "longitudes",longitudes
1815        latitudes_new = latitudes[kmin:kmax]
1816        longitudes_news = longitudes[lmin:lmax]
1817        #print "latitudes_new", latitudes_new
1818        #print "longitudes_news",longitudes_news
1819
1820        self.failUnless(latitudes_new == latitudes  and \
1821                        longitudes_news == longitudes,
1822                         'failed')
1823
1824    def test_tsh2sww(self):
1825        import os
1826        import tempfile
1827
1828        tsh_file = tempfile.mktemp(".tsh")
1829        file = open(tsh_file,"w")
1830        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
18310 0.0 0.0 0.0 0.0 0.01 \n \
18321 1.0 0.0 10.0 10.0 0.02  \n \
18332 0.0 1.0 0.0 10.0 0.03  \n \
18343 0.5 0.25 8.0 12.0 0.04  \n \
1835# Vert att title  \n \
1836elevation  \n \
1837stage  \n \
1838friction  \n \
18392 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
18400 0 3 2 -1  -1  1 dsg\n\
18411 0 1 3 -1  0 -1   ole nielsen\n\
18424 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
18430 1 0 2 \n\
18441 0 2 3 \n\
18452 2 3 \n\
18463 3 1 1 \n\
18473 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
18480 216.0 -86.0 \n \
18491 160.0 -167.0 \n \
18502 114.0 -91.0 \n \
18513 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
18520 0 1 0 \n \
18531 1 2 0 \n \
18542 2 0 0 \n \
18550 # <x> <y> ...Mesh Holes... \n \
18560 # <x> <y> <attribute>...Mesh Regions... \n \
18570 # <x> <y> <attribute>...Mesh Regions, area... \n\
1858#Geo reference \n \
185956 \n \
1860140 \n \
1861120 \n")
1862        file.close()
1863
1864        #sww_file = tempfile.mktemp(".sww")
1865        #print "sww_file",sww_file
1866        #print "sww_file",tsh_file
1867        tsh2sww(tsh_file,
1868                verbose=self.verbose)
1869
1870        os.remove(tsh_file)
1871        os.remove(tsh_file[:-4] + '.sww')
1872       
1873         
1874    def test_urs2sww_test_fail(self):
1875        points_num = -100
1876        time_step_count = 45
1877        time_step = -7
1878        file_handle, base_name = tempfile.mkstemp("")       
1879        os.close(file_handle)
1880        os.remove(base_name)
1881       
1882        files = []
1883        quantities = ['HA','UA','VA']
1884       
1885        mux_names = [WAVEHEIGHT_MUX_LABEL,
1886                     EAST_VELOCITY_LABEL,
1887                     NORTH_VELOCITY_LABEL]
1888        for i,q in enumerate(quantities): 
1889            #Write C files
1890            columns = 3 # long, lat , depth
1891            file = base_name + mux_names[i]
1892            f = open(file, 'wb')
1893            files.append(file)
1894            f.write(pack('i',points_num))
1895            f.write(pack('i',time_step_count))
1896            f.write(pack('f',time_step))
1897
1898            f.close()   
1899        tide = 1
1900        try:
1901            urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
1902                      verbose=self.verbose)       
1903        except ANUGAError:
1904            pass
1905        else:
1906            self.delete_mux(files)
1907            msg = 'Should have raised exception'
1908            raise msg
1909        sww_file = base_name + '.sww'
1910        self.delete_mux(files)
1911       
1912    def test_urs2sww_test_fail2(self):
1913        base_name = 'Harry-high-pants'
1914        try:
1915            urs2sww(base_name)       
1916        except IOError:
1917            pass
1918        else:
1919            self.delete_mux(files)
1920            msg = 'Should have raised exception'
1921            raise msg
1922           
1923    def test_urs2sww(self):
1924        tide = 1
1925        base_name, files = self.create_mux()
1926        urs2sww(base_name
1927                #, origin=(0,0,0)
1928                , mean_stage=tide
1929                , remove_nc_files=True,
1930                      verbose=self.verbose
1931                )
1932        sww_file = base_name + '.sww'
1933       
1934        #Let's interigate the sww file
1935        # Note, the sww info is not gridded.  It is point data.
1936        fid = NetCDFFile(sww_file)
1937
1938        x = fid.variables['x'][:]
1939        y = fid.variables['y'][:]
1940        geo_reference = Geo_reference(NetCDFObject=fid)
1941
1942       
1943        #Check that first coordinate is correctly represented       
1944        #Work out the UTM coordinates for first point
1945        zone, e, n = redfearn(-34.5, 150.66667)
1946       
1947        assert num.allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
1948
1949        # Make x and y absolute
1950        points = geo_reference.get_absolute(map(None, x, y))
1951        points = ensure_numeric(points)
1952        x = points[:,0]
1953        y = points[:,1]
1954       
1955        #Check first value
1956        stage = fid.variables['stage'][:]
1957        xmomentum = fid.variables['xmomentum'][:]
1958        ymomentum = fid.variables['ymomentum'][:]
1959        elevation = fid.variables['elevation'][:]
1960        assert num.allclose(stage[0,0], e +tide)  #Meters
1961
1962        #Check the momentums - ua
1963        #momentum = velocity*(stage-elevation)
1964        # elevation = - depth
1965        #momentum = velocity_ua *(stage+depth)
1966        # = n*(e+tide+n) based on how I'm writing these files
1967        #
1968        answer_x = n*(e+tide+n)
1969        actual_x = xmomentum[0,0]
1970        #print "answer_x",answer_x
1971        #print "actual_x",actual_x
1972        assert num.allclose(answer_x, actual_x)  #Meters
1973
1974        #Check the momentums - va
1975        #momentum = velocity*(stage-elevation)
1976        # -(-elevation) since elevation is inverted in mux files
1977        #momentum = velocity_va *(stage+elevation)
1978        # = e*(e+tide+n) based on how I'm writing these files
1979        answer_y = e*(e+tide+n) * -1 # talking into account mux file format
1980        actual_y = ymomentum[0,0]
1981        #print "answer_y",answer_y
1982        #print "actual_y",actual_y
1983        assert num.allclose(answer_y, actual_y)  #Meters
1984       
1985        assert num.allclose(answer_x, actual_x)  #Meters
1986       
1987        # check the stage values, first time step.
1988        # These arrays are equal since the Easting values were used as
1989        # the stage
1990        assert num.allclose(stage[0], x +tide)  #Meters
1991
1992        # check the elevation values.
1993        # -ve since urs measures depth, sww meshers height,
1994        # these arrays are equal since the northing values were used as
1995        # the elevation
1996        assert num.allclose(-elevation, y)  #Meters
1997       
1998        fid.close()
1999        self.delete_mux(files)
2000        os.remove(sww_file)
2001       
2002    def test_urs2sww_momentum(self):
2003        tide = 1
2004        time_step_count = 3
2005        time_step = 2
2006        #lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
2007        # This is gridded
2008        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
2009        depth=20
2010        ha=2
2011        ua=5
2012        va=-10 #-ve added to take into account mux file format where south
2013               # is positive.
2014        base_name, files = self.write_mux(lat_long_points,
2015                                          time_step_count, time_step,
2016                                          depth=depth,
2017                                          ha=ha,
2018                                          ua=ua,
2019                                          va=va)
2020        # write_mux(self,lat_long_points, time_step_count, time_step,
2021        #          depth=None, ha=None, ua=None, va=None
2022        urs2sww(base_name
2023                #, origin=(0,0,0)
2024                , mean_stage=tide
2025                , remove_nc_files=True,
2026                      verbose=self.verbose
2027                )
2028        sww_file = base_name + '.sww'
2029       
2030        #Let's interigate the sww file
2031        # Note, the sww info is not gridded.  It is point data.
2032        fid = NetCDFFile(sww_file)
2033
2034        x = fid.variables['x'][:]
2035        y = fid.variables['y'][:]
2036        geo_reference = Geo_reference(NetCDFObject=fid)
2037       
2038        #Check first value
2039        stage = fid.variables['stage'][:]
2040        xmomentum = fid.variables['xmomentum'][:]
2041        ymomentum = fid.variables['ymomentum'][:]
2042        elevation = fid.variables['elevation'][:]
2043        #assert allclose(stage[0,0], e + tide)  #Meters
2044        #print "xmomentum", xmomentum
2045        #print "ymomentum", ymomentum
2046        #Check the momentums - ua
2047        #momentum = velocity*water height
2048        #water height = mux_depth + mux_height +tide
2049        #water height = mux_depth + mux_height +tide
2050        #momentum = velocity*(mux_depth + mux_height +tide)
2051        #
2052       
2053        answer = 115
2054        actual = xmomentum[0,0]
2055        assert num.allclose(answer, actual)  #Meters^2/ sec
2056        answer = 230
2057        actual = ymomentum[0,0]
2058        #print "answer",answer
2059        #print "actual",actual
2060        assert num.allclose(answer, actual)  #Meters^2/ sec
2061
2062        # check the stage values, first time step.
2063        # These arrays are equal since the Easting values were used as
2064        # the stage
2065
2066        #assert allclose(stage[0], x +tide)  #Meters
2067
2068        # check the elevation values.
2069        # -ve since urs measures depth, sww meshers height,
2070        # these arrays are equal since the northing values were used as
2071        # the elevation
2072        #assert allclose(-elevation, y)  #Meters
2073       
2074        fid.close()
2075        self.delete_mux(files)
2076        os.remove(sww_file)
2077       
2078 
2079    def test_urs2sww_origin(self):
2080        tide = 1
2081        base_name, files = self.create_mux()
2082        urs2sww(base_name
2083                , origin=(0,0,0)
2084                , mean_stage=tide
2085                , remove_nc_files=True,
2086                      verbose=self.verbose
2087                )
2088        sww_file = base_name + '.sww'
2089       
2090        #Let's interigate the sww file
2091        # Note, the sww info is not gridded.  It is point data.
2092        fid = NetCDFFile(sww_file)
2093
2094        #  x and y are absolute
2095        x = fid.variables['x'][:]
2096        y = fid.variables['y'][:]
2097        geo_reference = Geo_reference(NetCDFObject=fid)
2098
2099       
2100        time = fid.variables['time'][:]
2101        #print "time", time
2102        assert num.allclose([0.,0.5,1.], time)
2103        assert fid.starttime == 0.0
2104        #Check that first coordinate is correctly represented       
2105        #Work out the UTM coordinates for first point
2106        zone, e, n = redfearn(-34.5, 150.66667)       
2107       
2108        assert num.allclose([x[0],y[0]], [e,n])
2109
2110       
2111        #Check first value
2112        stage = fid.variables['stage'][:]
2113        xmomentum = fid.variables['xmomentum'][:]
2114        ymomentum = fid.variables['ymomentum'][:]
2115        elevation = fid.variables['elevation'][:]
2116        assert num.allclose(stage[0,0], e +tide)  #Meters
2117
2118        #Check the momentums - ua
2119        #momentum = velocity*(stage-elevation)
2120        #momentum = velocity*(stage+elevation)
2121        # -(-elevation) since elevation is inverted in mux files
2122        # = n*(e+tide+n) based on how I'm writing these files
2123        answer = n*(e+tide+n)
2124        actual = xmomentum[0,0]
2125        assert num.allclose(answer, actual)  #Meters
2126
2127        # check the stage values, first time step.
2128        # These arrays are equal since the Easting values were used as
2129        # the stage
2130        assert num.allclose(stage[0], x +tide)  #Meters
2131
2132        # check the elevation values.
2133        # -ve since urs measures depth, sww meshers height,
2134        # these arrays are equal since the northing values were used as
2135        # the elevation
2136        assert num.allclose(-elevation, y)  #Meters
2137       
2138        fid.close()
2139        self.delete_mux(files)
2140        os.remove(sww_file)
2141 
2142    def test_urs2sww_minmaxlatlong(self):
2143       
2144        #longitudes = [150.66667, 150.83334, 151., 151.16667]
2145        #latitudes = [-34.5, -34.33333, -34.16667, -34]
2146
2147        tide = 1
2148        base_name, files = self.create_mux()
2149        urs2sww(base_name,
2150                minlat=-34.5,
2151                maxlat=-34,
2152                minlon= 150.66667,
2153                maxlon= 151.16667,
2154                mean_stage=tide,
2155                remove_nc_files=True,
2156                      verbose=self.verbose
2157                )
2158        sww_file = base_name + '.sww'
2159       
2160        #Let's interigate the sww file
2161        # Note, the sww info is not gridded.  It is point data.
2162        fid = NetCDFFile(sww_file)
2163       
2164
2165        # Make x and y absolute
2166        x = fid.variables['x'][:]
2167        y = fid.variables['y'][:]
2168        geo_reference = Geo_reference(NetCDFObject=fid)
2169        points = geo_reference.get_absolute(map(None, x, y))
2170        points = ensure_numeric(points)
2171        x = points[:,0]
2172        y = points[:,1]
2173       
2174        #Check that first coordinate is correctly represented       
2175        #Work out the UTM coordinates for first point
2176        zone, e, n = redfearn(-34.5, 150.66667) 
2177        assert num.allclose([x[0],y[0]], [e,n])
2178
2179       
2180        #Check first value
2181        stage = fid.variables['stage'][:]
2182        xmomentum = fid.variables['xmomentum'][:]
2183        ymomentum = fid.variables['ymomentum'][:]
2184        elevation = fid.variables['elevation'][:]
2185        assert num.allclose(stage[0,0], e +tide)  #Meters
2186
2187        #Check the momentums - ua
2188        #momentum = velocity*(stage-elevation)
2189        #momentum = velocity*(stage+elevation)
2190        # -(-elevation) since elevation is inverted in mux files
2191        # = n*(e+tide+n) based on how I'm writing these files
2192        answer = n*(e+tide+n)
2193        actual = xmomentum[0,0]
2194        assert num.allclose(answer, actual)  #Meters
2195
2196        # check the stage values, first time step.
2197        # These arrays are equal since the Easting values were used as
2198        # the stage
2199        assert num.allclose(stage[0], x +tide)  #Meters
2200
2201        # check the elevation values.
2202        # -ve since urs measures depth, sww meshers height,
2203        # these arrays are equal since the northing values were used as
2204        # the elevation
2205        assert num.allclose(-elevation, y)  #Meters
2206       
2207        fid.close()
2208        self.delete_mux(files)
2209        os.remove(sww_file)
2210       
2211    def test_urs2sww_minmaxmintmaxt(self):
2212       
2213        #longitudes = [150.66667, 150.83334, 151., 151.16667]
2214        #latitudes = [-34.5, -34.33333, -34.16667, -34]
2215
2216        tide = 1
2217        base_name, files = self.create_mux()
2218       
2219        urs2sww(base_name,
2220                mint=0.25,
2221                maxt=0.75,
2222                mean_stage=tide,
2223                remove_nc_files=True,
2224                verbose=self.verbose)
2225        sww_file = base_name + '.sww'
2226       
2227        #Let's interigate the sww file
2228        # Note, the sww info is not gridded.  It is point data.
2229        fid = NetCDFFile(sww_file)
2230
2231       
2232        time = fid.variables['time'][:]
2233        assert num.allclose(time, [0.0]) # the time is relative
2234        assert fid.starttime == 0.5
2235       
2236        fid.close()
2237        self.delete_mux(files)
2238        #print "sww_file", sww_file
2239        os.remove(sww_file)
2240 
2241
2242
2243#-------------------------------------------------------------
2244
2245if __name__ == "__main__":
2246    suite = unittest.makeSuite(Test_File_Conversion,'test')
2247    runner = unittest.TextTestRunner()
2248    runner.run(suite)
2249
Note: See TracBrowser for help on using the repository browser.