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

Last change on this file since 7780 was 7780, checked in by hudson, 13 years ago

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