source: trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py @ 7776

Last change on this file since 7776 was 7776, checked in by hudson, 15 years ago

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