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

Last change on this file since 4701 was 4701, checked in by ole, 16 years ago

Work towards ticket:192

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