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

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

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