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

Last change on this file since 4805 was 4805, checked in by ole, 17 years ago

Work towards making tight_slope_limiters the default. Next step is to run a large simulation to make sure timestepping isn't being adversely affected.
File size: 235.1 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
234        self.domain.tight_slope_limiters = 0 # Backwards compatibility
235       
236        sww = get_dataobject(self.domain)       
237
238        for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
239            pass
240           
241        # Get NetCDF
242        fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
243
244        # Get the variables
245        range = fid.variables['stage_range'][:]
246        assert allclose(range,[-0.93519, 0.15]) or\
247               allclose(range,[-0.9352743, 0.15]) # Old slope limiters
248       
249        range = fid.variables['xmomentum_range'][:]
250        assert allclose(range,[0,0.4695096]) or\
251               allclose(range,[0,0.47790655]) # Old slope limiters
252       
253        range = fid.variables['ymomentum_range'][:]
254        #assert allclose(range,[0,0.02174380])
255        assert allclose(range,[0,0.02174439]) or\
256               allclose(range,[0,0.02283983]) # Old slope limiters
257       
258        fid.close()
259        os.remove(sww.filename)
260
261    def test_sww_extrema(self):
262        """Test that extrema of quantities can be retrieved at every vertex
263        Extrema are updated at every *internal* timestep
264        """
265
266        domain = self.domain
267       
268        domain.set_name('extrema_test' + str(id(self)))
269        domain.format = 'sww'
270        domain.smooth = True
271
272        assert domain.quantities_to_be_monitored is None
273        assert domain.monitor_polygon is None
274        assert domain.monitor_time_interval is None       
275       
276        domain.set_quantities_to_be_monitored(['xmomentum',
277                                               'ymomentum',
278                                               'stage-elevation'])
279
280        assert domain.monitor_polygon is None
281        assert domain.monitor_time_interval is None
282
283
284        domain.set_quantities_to_be_monitored(['xmomentum',
285                                               'ymomentum',
286                                               'stage-elevation'],
287                                              polygon=domain.get_boundary_polygon(),
288                                              time_interval=[0,1])
289       
290       
291        assert len(domain.quantities_to_be_monitored) == 3
292        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
293        assert domain.quantities_to_be_monitored.has_key('xmomentum')               
294        assert domain.quantities_to_be_monitored.has_key('ymomentum')       
295
296       
297        #domain.protect_against_isolated_degenerate_timesteps = True
298        #domain.tight_slope_limiters = 1
299        domain.tight_slope_limiters = 0 # Backwards compatibility
300       
301        sww = get_dataobject(domain)
302
303        for t in domain.evolve(yieldstep = 1, finaltime = 1):
304            pass
305            #print domain.timestepping_statistics()
306            domain.quantity_statistics(precision = '%.8f') # Silent
307
308           
309        # Get NetCDF
310        fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
311
312        # Get the variables
313        extrema = fid.variables['stage-elevation.extrema'][:]
314        assert allclose(extrema, [0.00, 0.30])
315
316        loc = fid.variables['stage-elevation.min_location'][:]
317        assert allclose(loc, [0.16666667, 0.33333333])
318
319        loc = fid.variables['stage-elevation.max_location'][:]       
320        assert allclose(loc, [0.8333333, 0.16666667])       
321
322        time = fid.variables['stage-elevation.max_time'][:]
323        assert allclose(time, 0.0)               
324
325        extrema = fid.variables['xmomentum.extrema'][:]
326        assert allclose(extrema,[-0.06062178, 0.47886313])
327       
328        extrema = fid.variables['ymomentum.extrema'][:]
329        assert allclose(extrema,[0.00, 0.06241221])       
330
331        time_interval = fid.variables['extrema.time_interval'][:]
332        assert allclose(time_interval, [0,1])
333       
334        polygon = fid.variables['extrema.polygon'][:]       
335        assert allclose(polygon, domain.get_boundary_polygon())
336       
337        fid.close()
338        #print "sww.filename", sww.filename
339        os.remove(sww.filename)
340
341       
342       
343    def test_sww_constant_smooth(self):
344        """Test that constant sww information can be written correctly
345        (non smooth)
346        """
347        self.domain.set_name('datatest' + str(id(self)))
348        self.domain.format = 'sww'
349        self.domain.smooth = True
350
351        sww = get_dataobject(self.domain)
352        sww.store_connectivity()
353
354        #Check contents
355        #Get NetCDF
356        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
357
358        # Get the variables
359        X = fid.variables['x'][:]
360        Y = fid.variables['y'][:]
361        Z = fid.variables['elevation'][:]
362        V = fid.variables['volumes']
363
364        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
365        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
366        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
367        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
368        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
369
370        assert Z[4] == -0.5
371
372        assert V[2,0] == 4
373        assert V[2,1] == 5
374        assert V[2,2] == 1
375        assert V[4,0] == 6
376        assert V[4,1] == 7
377        assert V[4,2] == 3
378
379        fid.close()
380        os.remove(sww.filename)
381       
382
383    def test_sww_variable(self):
384        """Test that sww information can be written correctly
385        """
386        self.domain.set_name('datatest' + str(id(self)))
387        self.domain.format = 'sww'
388        self.domain.smooth = True
389        self.domain.reduction = mean
390
391        sww = get_dataobject(self.domain)
392        sww.store_connectivity()
393        sww.store_timestep('stage')
394
395        #Check contents
396        #Get NetCDF
397        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
398
399
400        # Get the variables
401        time = fid.variables['time']
402        stage = fid.variables['stage']
403
404        Q = self.domain.quantities['stage']
405        Q0 = Q.vertex_values[:,0]
406        Q1 = Q.vertex_values[:,1]
407        Q2 = Q.vertex_values[:,2]
408
409        A = stage[0,:]
410        #print A[0], (Q2[0,0] + Q1[1,0])/2
411        assert allclose(A[0], (Q2[0] + Q1[1])/2)
412        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
413        assert allclose(A[2], Q0[3])
414        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
415
416        #Center point
417        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
418                                 Q0[5] + Q2[6] + Q1[7])/6)
419       
420        fid.close()
421        os.remove(sww.filename)
422
423
424    def test_sww_variable2(self):
425        """Test that sww information can be written correctly
426        multiple timesteps. Use average as reduction operator
427        """
428
429        import time, os
430        from Numeric import array, zeros, allclose, Float, concatenate
431        from Scientific.IO.NetCDF import NetCDFFile
432
433        self.domain.set_name('datatest' + str(id(self)))
434        self.domain.format = 'sww'
435        self.domain.smooth = True
436
437        self.domain.reduction = mean
438
439        sww = get_dataobject(self.domain)
440        sww.store_connectivity()
441        sww.store_timestep('stage')
442        #self.domain.tight_slope_limiters = 1
443        self.domain.evolve_to_end(finaltime = 0.01)
444        sww.store_timestep('stage')
445
446
447        #Check contents
448        #Get NetCDF
449        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
450
451        # Get the variables
452        x = fid.variables['x']
453        y = fid.variables['y']
454        z = fid.variables['elevation']
455        time = fid.variables['time']
456        stage = fid.variables['stage']
457
458        #Check values
459        Q = self.domain.quantities['stage']
460        Q0 = Q.vertex_values[:,0]
461        Q1 = Q.vertex_values[:,1]
462        Q2 = Q.vertex_values[:,2]
463
464        A = stage[1,:]
465        assert allclose(A[0], (Q2[0] + Q1[1])/2)
466        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
467        assert allclose(A[2], Q0[3])
468        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
469
470        #Center point
471        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
472                                 Q0[5] + Q2[6] + Q1[7])/6)
473
474
475        fid.close()
476
477        #Cleanup
478        os.remove(sww.filename)
479
480    def no_test_sww_variable3(self):
481        """Test that sww information can be written correctly
482        multiple timesteps using a different reduction operator (min)
483        """
484
485        # Different reduction in sww files has been made obsolete.
486       
487        import time, os
488        from Numeric import array, zeros, allclose, Float, concatenate
489        from Scientific.IO.NetCDF import NetCDFFile
490
491        self.domain.set_name('datatest' + str(id(self)))
492        self.domain.format = 'sww'
493        self.domain.smooth = True
494        self.domain.reduction = min
495
496        sww = get_dataobject(self.domain)
497        sww.store_connectivity()
498        sww.store_timestep('stage')
499        #self.domain.tight_slope_limiters = 1
500        self.domain.evolve_to_end(finaltime = 0.01)
501        sww.store_timestep('stage')
502
503
504        #Check contents
505        #Get NetCDF
506        fid = NetCDFFile(sww.filename, 'r')
507
508        # Get the variables
509        x = fid.variables['x']
510        y = fid.variables['y']
511        z = fid.variables['elevation']
512        time = fid.variables['time']
513        stage = fid.variables['stage']
514
515        #Check values
516        Q = self.domain.quantities['stage']
517        Q0 = Q.vertex_values[:,0]
518        Q1 = Q.vertex_values[:,1]
519        Q2 = Q.vertex_values[:,2]
520
521        A = stage[1,:]
522        assert allclose(A[0], min(Q2[0], Q1[1]))
523        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
524        assert allclose(A[2], Q0[3])
525        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
526
527        #Center point
528        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
529                                  Q0[5], Q2[6], Q1[7]))
530
531
532        fid.close()
533
534        #Cleanup
535        os.remove(sww.filename)
536
537
538    def test_sync(self):
539        """test_sync - Test info stored at each timestep is as expected (incl initial condition)
540        """
541
542        import time, os, config
543        from Numeric import array, zeros, allclose, Float, concatenate
544        from Scientific.IO.NetCDF import NetCDFFile
545
546        self.domain.set_name('synctest')
547        self.domain.format = 'sww'
548        self.domain.smooth = False
549        self.domain.store = True
550        self.domain.beta_h = 0
551
552        # In this case tight_slope_limiters as default
553        # in conjunction with protection
554        # against isolated degenerate timesteps works.
555        #self.domain.tight_slope_limiters = 1
556        #self.domain.protect_against_isolated_degenerate_timesteps = True
557
558        #print 'tight_sl', self.domain.tight_slope_limiters
559       
560
561        #Evolution
562        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
563           
564            #########self.domain.write_time(track_speeds=True)
565            stage = self.domain.quantities['stage'].vertex_values
566
567            #Get NetCDF
568            fid = NetCDFFile(self.domain.writer.filename, 'r')
569            stage_file = fid.variables['stage']
570           
571            if t == 0.0:
572                assert allclose(stage, self.initial_stage)
573                assert allclose(stage_file[:], stage.flat)
574            else:
575                assert not allclose(stage, self.initial_stage)
576                assert not allclose(stage_file[:], stage.flat)
577
578            fid.close()
579
580        os.remove(self.domain.writer.filename)
581
582
583    def test_sww_minimum_storable_height(self):
584        """Test that sww information can be written correctly
585        multiple timesteps using a different reduction operator (min)
586        """
587
588        import time, os
589        from Numeric import array, zeros, allclose, Float, concatenate
590        from Scientific.IO.NetCDF import NetCDFFile
591
592        self.domain.set_name('datatest' + str(id(self)))
593        self.domain.format = 'sww'
594        self.domain.smooth = True
595        self.domain.reduction = min
596        self.domain.minimum_storable_height = 100
597
598        sww = get_dataobject(self.domain)
599        sww.store_connectivity()
600        sww.store_timestep('stage')
601
602        #self.domain.tight_slope_limiters = 1
603        self.domain.evolve_to_end(finaltime = 0.01)
604        sww.store_timestep('stage')
605
606
607        #Check contents
608        #Get NetCDF
609        fid = NetCDFFile(sww.filename, 'r')
610
611
612        # Get the variables
613        x = fid.variables['x']
614        y = fid.variables['y']
615        z = fid.variables['elevation']
616        time = fid.variables['time']
617        stage = fid.variables['stage']
618
619        #Check values
620        Q = self.domain.quantities['stage']
621        Q0 = Q.vertex_values[:,0]
622        Q1 = Q.vertex_values[:,1]
623        Q2 = Q.vertex_values[:,2]
624
625        A = stage[1,:]
626        assert allclose(stage[1,:], z[:])
627        fid.close()
628
629        #Cleanup
630        os.remove(sww.filename)
631
632
633    def Not_a_test_sww_DSG(self):
634        """Not a test, rather a look at the sww format
635        """
636
637        import time, os
638        from Numeric import array, zeros, allclose, Float, concatenate
639        from Scientific.IO.NetCDF import NetCDFFile
640
641        self.domain.set_name('datatest' + str(id(self)))
642        self.domain.format = 'sww'
643        self.domain.smooth = True
644        self.domain.reduction = mean
645
646        sww = get_dataobject(self.domain)
647        sww.store_connectivity()
648        sww.store_timestep('stage')
649
650        #Check contents
651        #Get NetCDF
652        fid = NetCDFFile(sww.filename, 'r')
653
654        # Get the variables
655        x = fid.variables['x']
656        y = fid.variables['y']
657        z = fid.variables['elevation']
658
659        volumes = fid.variables['volumes']
660        time = fid.variables['time']
661
662        # 2D
663        stage = fid.variables['stage']
664
665        X = x[:]
666        Y = y[:]
667        Z = z[:]
668        V = volumes[:]
669        T = time[:]
670        S = stage[:,:]
671
672#         print "****************************"
673#         print "X ",X
674#         print "****************************"
675#         print "Y ",Y
676#         print "****************************"
677#         print "Z ",Z
678#         print "****************************"
679#         print "V ",V
680#         print "****************************"
681#         print "Time ",T
682#         print "****************************"
683#         print "Stage ",S
684#         print "****************************"
685
686
687        fid.close()
688
689        #Cleanup
690        os.remove(sww.filename)
691
692
693
694    def test_dem2pts_bounding_box_v2(self):
695        """Test conversion from dem in ascii format to native NetCDF format
696        """
697
698        import time, os
699        from Numeric import array, zeros, allclose, Float, concatenate, ones
700        from Scientific.IO.NetCDF import NetCDFFile
701
702        #Write test asc file
703        root = 'demtest'
704
705        filename = root+'.asc'
706        fid = open(filename, 'w')
707        fid.write("""ncols         10
708nrows         10
709xllcorner     2000
710yllcorner     3000
711cellsize      1
712NODATA_value  -9999
713""")
714        #Create linear function
715        ref_points = []
716        ref_elevation = []
717        x0 = 2000
718        y = 3010
719        yvec = range(10)
720        xvec = range(10)
721        z = -1
722        for i in range(10):
723            y = y - 1
724            for j in range(10):
725                x = x0 + xvec[j]
726                z += 1
727                ref_points.append ([x,y])
728                ref_elevation.append(z)
729                fid.write('%f ' %z)
730            fid.write('\n')
731
732        fid.close()
733
734        #print 'sending pts', ref_points
735        #print 'sending elev', ref_elevation
736
737        #Write prj file with metadata
738        metafilename = root+'.prj'
739        fid = open(metafilename, 'w')
740
741
742        fid.write("""Projection UTM
743Zone 56
744Datum WGS84
745Zunits NO
746Units METERS
747Spheroid WGS84
748Xshift 0.0000000000
749Yshift 10000000.0000000000
750Parameters
751""")
752        fid.close()
753
754        #Convert to NetCDF pts
755        convert_dem_from_ascii2netcdf(root)
756        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
757                northing_min=3003.0, northing_max=3006.0,
758                verbose=self.verbose)
759
760        #Check contents
761        #Get NetCDF
762        fid = NetCDFFile(root+'.pts', 'r')
763
764        # Get the variables
765        #print fid.variables.keys()
766        points = fid.variables['points']
767        elevation = fid.variables['elevation']
768
769        #Check values
770        assert fid.xllcorner[0] == 2002.0
771        assert fid.yllcorner[0] == 3003.0
772
773        #create new reference points
774        newz = []
775        newz[0:5] = ref_elevation[32:38]
776        newz[6:11] = ref_elevation[42:48]
777        newz[12:17] = ref_elevation[52:58]
778        newz[18:23] = ref_elevation[62:68]
779        ref_elevation = []
780        ref_elevation = newz
781        ref_points = []
782        x0 = 2002
783        y = 3007
784        yvec = range(4)
785        xvec = range(6)
786        for i in range(4):
787            y = y - 1
788            ynew = y - 3003.0
789            for j in range(6):
790                x = x0 + xvec[j]
791                xnew = x - 2002.0
792                ref_points.append ([xnew,ynew]) #Relative point values
793
794        assert allclose(points, ref_points)
795
796        assert allclose(elevation, ref_elevation)
797
798        #Cleanup
799        fid.close()
800
801
802        os.remove(root + '.pts')
803        os.remove(root + '.dem')
804        os.remove(root + '.asc')
805        os.remove(root + '.prj')
806
807
808    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
809        """Test conversion from dem in ascii format to native NetCDF format
810        """
811
812        import time, os
813        from Numeric import array, zeros, allclose, Float, concatenate, ones
814        from Scientific.IO.NetCDF import NetCDFFile
815
816        #Write test asc file
817        root = 'demtest'
818
819        filename = root+'.asc'
820        fid = open(filename, 'w')
821        fid.write("""ncols         10
822nrows         10
823xllcorner     2000
824yllcorner     3000
825cellsize      1
826NODATA_value  -9999
827""")
828        #Create linear function
829        ref_points = []
830        ref_elevation = []
831        x0 = 2000
832        y = 3010
833        yvec = range(10)
834        xvec = range(10)
835        #z = range(100)
836        z = zeros(100)
837        NODATA_value = -9999
838        count = -1
839        for i in range(10):
840            y = y - 1
841            for j in range(10):
842                x = x0 + xvec[j]
843                ref_points.append ([x,y])
844                count += 1
845                z[count] = (4*i - 3*j)%13
846                if j == 4: z[count] = NODATA_value #column inside clipping region
847                if j == 8: z[count] = NODATA_value #column outside clipping region
848                if i == 9: z[count] = NODATA_value #row outside clipping region
849                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
850                ref_elevation.append( z[count] )
851                fid.write('%f ' %z[count])
852            fid.write('\n')
853
854        fid.close()
855
856        #print 'sending elev', ref_elevation
857
858        #Write prj file with metadata
859        metafilename = root+'.prj'
860        fid = open(metafilename, 'w')
861
862
863        fid.write("""Projection UTM
864Zone 56
865Datum WGS84
866Zunits NO
867Units METERS
868Spheroid WGS84
869Xshift 0.0000000000
870Yshift 10000000.0000000000
871Parameters
872""")
873        fid.close()
874
875        #Convert to NetCDF pts
876        convert_dem_from_ascii2netcdf(root)
877        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
878                northing_min=3003.0, northing_max=3006.0,
879                verbose=self.verbose)
880
881        #Check contents
882        #Get NetCDF
883        fid = NetCDFFile(root+'.pts', 'r')
884
885        # Get the variables
886        #print fid.variables.keys()
887        points = fid.variables['points']
888        elevation = fid.variables['elevation']
889
890        #Check values
891        assert fid.xllcorner[0] == 2002.0
892        assert fid.yllcorner[0] == 3003.0
893
894        #create new reference points
895        newz = zeros(19)
896        newz[0:2] = ref_elevation[32:34]
897        newz[2:5] = ref_elevation[35:38]
898        newz[5:7] = ref_elevation[42:44]
899        newz[7] = ref_elevation[45]
900        newz[8] = ref_elevation[47]
901        newz[9:11] = ref_elevation[52:54]
902        newz[11:14] = ref_elevation[55:58]
903        newz[14:16] = ref_elevation[62:64]
904        newz[16:19] = ref_elevation[65:68]
905
906
907        ref_elevation = newz
908        ref_points = []
909        new_ref_points = []
910        x0 = 2002
911        y = 3007
912        yvec = range(4)
913        xvec = range(6)
914        for i in range(4):
915            y = y - 1
916            ynew = y - 3003.0
917            for j in range(6):
918                x = x0 + xvec[j]
919                xnew = x - 2002.0
920                if j <> 2 and (i<>1 or j<>4):
921                    ref_points.append([x,y])
922                    new_ref_points.append ([xnew,ynew])
923
924
925        assert allclose(points, new_ref_points)
926        assert allclose(elevation, ref_elevation)
927
928        #Cleanup
929        fid.close()
930
931
932        os.remove(root + '.pts')
933        os.remove(root + '.dem')
934        os.remove(root + '.asc')
935        os.remove(root + '.prj')
936
937
938    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
939        """Test conversion from dem in ascii format to native NetCDF format
940        Check missing values on clipping boundary
941        """
942
943        import time, os
944        from Numeric import array, zeros, allclose, Float, concatenate, ones
945        from Scientific.IO.NetCDF import NetCDFFile
946
947        #Write test asc file
948        root = 'demtest'
949
950        filename = root+'.asc'
951        fid = open(filename, 'w')
952        fid.write("""ncols         10
953nrows         10
954xllcorner     2000
955yllcorner     3000
956cellsize      1
957NODATA_value  -9999
958""")
959        #Create linear function
960        ref_points = []
961        ref_elevation = []
962        x0 = 2000
963        y = 3010
964        yvec = range(10)
965        xvec = range(10)
966        #z = range(100)
967        z = zeros(100)
968        NODATA_value = -9999
969        count = -1
970        for i in range(10):
971            y = y - 1
972            for j in range(10):
973                x = x0 + xvec[j]
974                ref_points.append ([x,y])
975                count += 1
976                z[count] = (4*i - 3*j)%13
977                if j == 4: z[count] = NODATA_value #column inside clipping region
978                if j == 8: z[count] = NODATA_value #column outside clipping region
979                if i == 6: z[count] = NODATA_value #row on clipping boundary
980                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
981                ref_elevation.append( z[count] )
982                fid.write('%f ' %z[count])
983            fid.write('\n')
984
985        fid.close()
986
987        #print 'sending elev', ref_elevation
988
989        #Write prj file with metadata
990        metafilename = root+'.prj'
991        fid = open(metafilename, 'w')
992
993
994        fid.write("""Projection UTM
995Zone 56
996Datum WGS84
997Zunits NO
998Units METERS
999Spheroid WGS84
1000Xshift 0.0000000000
1001Yshift 10000000.0000000000
1002Parameters
1003""")
1004        fid.close()
1005
1006        #Convert to NetCDF pts
1007        convert_dem_from_ascii2netcdf(root)
1008        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
1009                northing_min=3003.0, northing_max=3006.0,
1010                verbose=self.verbose)
1011
1012        #Check contents
1013        #Get NetCDF
1014        fid = NetCDFFile(root+'.pts', 'r')
1015
1016        # Get the variables
1017        #print fid.variables.keys()
1018        points = fid.variables['points']
1019        elevation = fid.variables['elevation']
1020
1021        #Check values
1022        assert fid.xllcorner[0] == 2002.0
1023        assert fid.yllcorner[0] == 3003.0
1024
1025        #create new reference points
1026        newz = zeros(14)
1027        newz[0:2] = ref_elevation[32:34]
1028        newz[2:5] = ref_elevation[35:38]
1029        newz[5:7] = ref_elevation[42:44]
1030        newz[7] = ref_elevation[45]
1031        newz[8] = ref_elevation[47]
1032        newz[9:11] = ref_elevation[52:54]
1033        newz[11:14] = ref_elevation[55:58]
1034
1035
1036
1037        ref_elevation = newz
1038        ref_points = []
1039        new_ref_points = []
1040        x0 = 2002
1041        y = 3007
1042        yvec = range(4)
1043        xvec = range(6)
1044        for i in range(4):
1045            y = y - 1
1046            ynew = y - 3003.0
1047            for j in range(6):
1048                x = x0 + xvec[j]
1049                xnew = x - 2002.0
1050                if j <> 2 and (i<>1 or j<>4) and i<>3:
1051                    ref_points.append([x,y])
1052                    new_ref_points.append ([xnew,ynew])
1053
1054
1055        #print points[:],points[:].shape
1056        #print new_ref_points, len(new_ref_points)
1057
1058        assert allclose(elevation, ref_elevation)
1059        assert allclose(points, new_ref_points)
1060
1061
1062        #Cleanup
1063        fid.close()
1064
1065
1066        os.remove(root + '.pts')
1067        os.remove(root + '.dem')
1068        os.remove(root + '.asc')
1069        os.remove(root + '.prj')
1070
1071
1072    def test_hecras_cross_sections2pts(self):
1073        """Test conversion from HECRAS cross sections in ascii format
1074        to native NetCDF pts format
1075        """
1076
1077        import time, os
1078        from Numeric import array, zeros, allclose, Float, concatenate
1079        from Scientific.IO.NetCDF import NetCDFFile
1080
1081        #Write test asc file
1082        root = 'hecrastest'
1083
1084        filename = root+'.sdf'
1085        fid = open(filename, 'w')
1086        fid.write("""
1087# RAS export file created on Mon 15Aug2005 11:42
1088# by HEC-RAS Version 3.1.1
1089
1090BEGIN HEADER:
1091  UNITS: METRIC
1092  DTM TYPE: TIN
1093  DTM: v:\1\cit\perth_topo\river_tin
1094  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
1095  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
1096  MAP PROJECTION: UTM
1097  PROJECTION ZONE: 50
1098  DATUM: AGD66
1099  VERTICAL DATUM:
1100  NUMBER OF REACHES:  19
1101  NUMBER OF CROSS-SECTIONS:  2
1102END HEADER:
1103
1104
1105BEGIN CROSS-SECTIONS:
1106
1107  CROSS-SECTION:
1108    STREAM ID:Southern-Wungong
1109    REACH ID:Southern-Wungong
1110    STATION:21410
1111    CUT LINE:
1112      407546.08 , 6437277.542
1113      407329.32 , 6437489.482
1114      407283.11 , 6437541.232
1115    SURFACE LINE:
1116     407546.08,   6437277.54,   52.14
1117     407538.88,   6437284.58,   51.07
1118     407531.68,   6437291.62,   50.56
1119     407524.48,   6437298.66,   49.58
1120     407517.28,   6437305.70,   49.09
1121     407510.08,   6437312.74,   48.76
1122  END:
1123
1124  CROSS-SECTION:
1125    STREAM ID:Swan River
1126    REACH ID:Swan Mouth
1127    STATION:840.*
1128    CUT LINE:
1129      381178.0855 , 6452559.0685
1130      380485.4755 , 6453169.272
1131    SURFACE LINE:
1132     381178.09,   6452559.07,   4.17
1133     381169.49,   6452566.64,   4.26
1134     381157.78,   6452576.96,   4.34
1135     381155.97,   6452578.56,   4.35
1136     381143.72,   6452589.35,   4.43
1137     381136.69,   6452595.54,   4.58
1138     381114.74,   6452614.88,   4.41
1139     381075.53,   6452649.43,   4.17
1140     381071.47,   6452653.00,   3.99
1141     381063.46,   6452660.06,   3.67
1142     381054.41,   6452668.03,   3.67
1143  END:
1144END CROSS-SECTIONS:
1145""")
1146
1147        fid.close()
1148
1149
1150        #Convert to NetCDF pts
1151        hecras_cross_sections2pts(root)
1152
1153        #Check contents
1154        #Get NetCDF
1155        fid = NetCDFFile(root+'.pts', 'r')
1156
1157        # Get the variables
1158        #print fid.variables.keys()
1159        points = fid.variables['points']
1160        elevation = fid.variables['elevation']
1161
1162        #Check values
1163        ref_points = [[407546.08, 6437277.54],
1164                      [407538.88, 6437284.58],
1165                      [407531.68, 6437291.62],
1166                      [407524.48, 6437298.66],
1167                      [407517.28, 6437305.70],
1168                      [407510.08, 6437312.74]]
1169
1170        ref_points += [[381178.09, 6452559.07],
1171                       [381169.49, 6452566.64],
1172                       [381157.78, 6452576.96],
1173                       [381155.97, 6452578.56],
1174                       [381143.72, 6452589.35],
1175                       [381136.69, 6452595.54],
1176                       [381114.74, 6452614.88],
1177                       [381075.53, 6452649.43],
1178                       [381071.47, 6452653.00],
1179                       [381063.46, 6452660.06],
1180                       [381054.41, 6452668.03]]
1181
1182
1183        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
1184        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
1185
1186        #print points[:]
1187        #print ref_points
1188        assert allclose(points, ref_points)
1189
1190        #print attributes[:]
1191        #print ref_elevation
1192        assert allclose(elevation, ref_elevation)
1193
1194        #Cleanup
1195        fid.close()
1196
1197
1198        os.remove(root + '.sdf')
1199        os.remove(root + '.pts')
1200
1201
1202    def test_sww2dem_asc_elevation_depth(self):
1203        """
1204        test_sww2dem_asc_elevation_depth(self):
1205        Test that sww information can be converted correctly to asc/prj
1206        format readable by e.g. ArcView
1207        """
1208
1209        import time, os
1210        from Numeric import array, zeros, allclose, Float, concatenate
1211        from Scientific.IO.NetCDF import NetCDFFile
1212
1213        #Setup
1214        self.domain.set_name('datatest')
1215
1216        prjfile = self.domain.get_name() + '_elevation.prj'
1217        ascfile = self.domain.get_name() + '_elevation.asc'
1218        swwfile = self.domain.get_name() + '.sww'
1219
1220        self.domain.set_datadir('.')
1221        self.domain.format = 'sww'
1222        self.domain.smooth = True
1223        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1224        self.domain.set_quantity('stage', 1.0)
1225
1226        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1227
1228        sww = get_dataobject(self.domain)
1229        sww.store_connectivity()
1230        sww.store_timestep('stage')
1231
1232        #self.domain.tight_slope_limiters = 1
1233
1234        self.domain.evolve_to_end(finaltime = 0.01)
1235        sww.store_timestep('stage')
1236
1237        cellsize = 0.25
1238        #Check contents
1239        #Get NetCDF
1240
1241        fid = NetCDFFile(sww.filename, 'r')
1242
1243        # Get the variables
1244        x = fid.variables['x'][:]
1245        y = fid.variables['y'][:]
1246        z = fid.variables['elevation'][:]
1247        time = fid.variables['time'][:]
1248        stage = fid.variables['stage'][:]
1249
1250        fid.close()
1251
1252        #Export to ascii/prj files
1253        sww2dem(self.domain.get_name(),
1254                quantity = 'elevation',
1255                cellsize = cellsize,
1256                verbose = self.verbose,
1257                format = 'asc')
1258
1259        #Check prj (meta data)
1260        prjid = open(prjfile)
1261        lines = prjid.readlines()
1262        prjid.close()
1263
1264        L = lines[0].strip().split()
1265        assert L[0].strip().lower() == 'projection'
1266        assert L[1].strip().lower() == 'utm'
1267
1268        L = lines[1].strip().split()
1269        assert L[0].strip().lower() == 'zone'
1270        assert L[1].strip().lower() == '56'
1271
1272        L = lines[2].strip().split()
1273        assert L[0].strip().lower() == 'datum'
1274        assert L[1].strip().lower() == 'wgs84'
1275
1276        L = lines[3].strip().split()
1277        assert L[0].strip().lower() == 'zunits'
1278        assert L[1].strip().lower() == 'no'
1279
1280        L = lines[4].strip().split()
1281        assert L[0].strip().lower() == 'units'
1282        assert L[1].strip().lower() == 'meters'
1283
1284        L = lines[5].strip().split()
1285        assert L[0].strip().lower() == 'spheroid'
1286        assert L[1].strip().lower() == 'wgs84'
1287
1288        L = lines[6].strip().split()
1289        assert L[0].strip().lower() == 'xshift'
1290        assert L[1].strip().lower() == '500000'
1291
1292        L = lines[7].strip().split()
1293        assert L[0].strip().lower() == 'yshift'
1294        assert L[1].strip().lower() == '10000000'
1295
1296        L = lines[8].strip().split()
1297        assert L[0].strip().lower() == 'parameters'
1298
1299
1300        #Check asc file
1301        ascid = open(ascfile)
1302        lines = ascid.readlines()
1303        ascid.close()
1304
1305        L = lines[0].strip().split()
1306        assert L[0].strip().lower() == 'ncols'
1307        assert L[1].strip().lower() == '5'
1308
1309        L = lines[1].strip().split()
1310        assert L[0].strip().lower() == 'nrows'
1311        assert L[1].strip().lower() == '5'
1312
1313        L = lines[2].strip().split()
1314        assert L[0].strip().lower() == 'xllcorner'
1315        assert allclose(float(L[1].strip().lower()), 308500)
1316
1317        L = lines[3].strip().split()
1318        assert L[0].strip().lower() == 'yllcorner'
1319        assert allclose(float(L[1].strip().lower()), 6189000)
1320
1321        L = lines[4].strip().split()
1322        assert L[0].strip().lower() == 'cellsize'
1323        assert allclose(float(L[1].strip().lower()), cellsize)
1324
1325        L = lines[5].strip().split()
1326        assert L[0].strip() == 'NODATA_value'
1327        assert L[1].strip().lower() == '-9999'
1328
1329        #Check grid values
1330        for j in range(5):
1331            L = lines[6+j].strip().split()
1332            y = (4-j) * cellsize
1333            for i in range(5):
1334                assert allclose(float(L[i]), -i*cellsize - y)
1335               
1336        #Cleanup
1337        os.remove(prjfile)
1338        os.remove(ascfile)
1339
1340        #Export to ascii/prj files
1341        sww2dem(self.domain.get_name(),
1342                quantity = 'depth',
1343                cellsize = cellsize,
1344                verbose = self.verbose,
1345                format = 'asc')
1346       
1347        #Check asc file
1348        ascfile = self.domain.get_name() + '_depth.asc'
1349        prjfile = self.domain.get_name() + '_depth.prj'
1350        ascid = open(ascfile)
1351        lines = ascid.readlines()
1352        ascid.close()
1353
1354        L = lines[0].strip().split()
1355        assert L[0].strip().lower() == 'ncols'
1356        assert L[1].strip().lower() == '5'
1357
1358        L = lines[1].strip().split()
1359        assert L[0].strip().lower() == 'nrows'
1360        assert L[1].strip().lower() == '5'
1361
1362        L = lines[2].strip().split()
1363        assert L[0].strip().lower() == 'xllcorner'
1364        assert allclose(float(L[1].strip().lower()), 308500)
1365
1366        L = lines[3].strip().split()
1367        assert L[0].strip().lower() == 'yllcorner'
1368        assert allclose(float(L[1].strip().lower()), 6189000)
1369
1370        L = lines[4].strip().split()
1371        assert L[0].strip().lower() == 'cellsize'
1372        assert allclose(float(L[1].strip().lower()), cellsize)
1373
1374        L = lines[5].strip().split()
1375        assert L[0].strip() == 'NODATA_value'
1376        assert L[1].strip().lower() == '-9999'
1377
1378        #Check grid values
1379        for j in range(5):
1380            L = lines[6+j].strip().split()
1381            y = (4-j) * cellsize
1382            for i in range(5):
1383                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1384
1385
1386        #Cleanup
1387        os.remove(prjfile)
1388        os.remove(ascfile)
1389        os.remove(swwfile)
1390
1391
1392    def test_export_grid(self):
1393        """
1394        test_export_grid(self):
1395        Test that sww information can be converted correctly to asc/prj
1396        format readable by e.g. ArcView
1397        """
1398
1399        import time, os
1400        from Numeric import array, zeros, allclose, Float, concatenate
1401        from Scientific.IO.NetCDF import NetCDFFile
1402
1403        try:
1404            os.remove('teg*.sww')
1405        except:
1406            pass
1407
1408        #Setup
1409        self.domain.set_name('teg')
1410
1411        prjfile = self.domain.get_name() + '_elevation.prj'
1412        ascfile = self.domain.get_name() + '_elevation.asc'
1413        swwfile = self.domain.get_name() + '.sww'
1414
1415        self.domain.set_datadir('.')
1416        self.domain.format = 'sww'
1417        self.domain.smooth = True
1418        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1419        self.domain.set_quantity('stage', 1.0)
1420
1421        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1422
1423        sww = get_dataobject(self.domain)
1424        sww.store_connectivity()
1425        sww.store_timestep('stage')
1426        self.domain.evolve_to_end(finaltime = 0.01)
1427        sww.store_timestep('stage')
1428
1429        cellsize = 0.25
1430        #Check contents
1431        #Get NetCDF
1432
1433        fid = NetCDFFile(sww.filename, 'r')
1434
1435        # Get the variables
1436        x = fid.variables['x'][:]
1437        y = fid.variables['y'][:]
1438        z = fid.variables['elevation'][:]
1439        time = fid.variables['time'][:]
1440        stage = fid.variables['stage'][:]
1441
1442        fid.close()
1443
1444        #Export to ascii/prj files
1445        export_grid(self.domain.get_name(),
1446                quantities = 'elevation',
1447                cellsize = cellsize,
1448                verbose = self.verbose,
1449                format = 'asc')
1450
1451        #Check asc file
1452        ascid = open(ascfile)
1453        lines = ascid.readlines()
1454        ascid.close()
1455
1456        L = lines[2].strip().split()
1457        assert L[0].strip().lower() == 'xllcorner'
1458        assert allclose(float(L[1].strip().lower()), 308500)
1459
1460        L = lines[3].strip().split()
1461        assert L[0].strip().lower() == 'yllcorner'
1462        assert allclose(float(L[1].strip().lower()), 6189000)
1463
1464        #Check grid values
1465        for j in range(5):
1466            L = lines[6+j].strip().split()
1467            y = (4-j) * cellsize
1468            for i in range(5):
1469                assert allclose(float(L[i]), -i*cellsize - y)
1470               
1471        #Cleanup
1472        os.remove(prjfile)
1473        os.remove(ascfile)
1474        os.remove(swwfile)
1475
1476    def test_export_gridII(self):
1477        """
1478        test_export_gridII(self):
1479        Test that sww information can be converted correctly to asc/prj
1480        format readable by e.g. ArcView
1481        """
1482
1483        import time, os
1484        from Numeric import array, zeros, allclose, Float, concatenate
1485        from Scientific.IO.NetCDF import NetCDFFile
1486
1487        try:
1488            os.remove('teg*.sww')
1489        except:
1490            pass
1491
1492        #Setup
1493        self.domain.set_name('tegII')
1494
1495        swwfile = self.domain.get_name() + '.sww'
1496
1497        self.domain.set_datadir('.')
1498        self.domain.format = 'sww'
1499        self.domain.smooth = True
1500        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1501        self.domain.set_quantity('stage', 1.0)
1502
1503        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1504
1505        sww = get_dataobject(self.domain)
1506        sww.store_connectivity()
1507        sww.store_timestep('stage')
1508        self.domain.evolve_to_end(finaltime = 0.01)
1509        sww.store_timestep('stage')
1510
1511        cellsize = 0.25
1512        #Check contents
1513        #Get NetCDF
1514
1515        fid = NetCDFFile(sww.filename, 'r')
1516
1517        # Get the variables
1518        x = fid.variables['x'][:]
1519        y = fid.variables['y'][:]
1520        z = fid.variables['elevation'][:]
1521        time = fid.variables['time'][:]
1522        stage = fid.variables['stage'][:]
1523
1524        fid.close()
1525
1526        #Export to ascii/prj files
1527        if True:
1528            export_grid(self.domain.get_name(),
1529                        quantities = ['elevation', 'depth'],
1530                        cellsize = cellsize,
1531                        verbose = self.verbose,
1532                        format = 'asc')
1533
1534        else:
1535            export_grid(self.domain.get_name(),
1536                quantities = ['depth'],
1537                cellsize = cellsize,
1538                verbose = self.verbose,
1539                format = 'asc')
1540
1541
1542            export_grid(self.domain.get_name(),
1543                quantities = ['elevation'],
1544                cellsize = cellsize,
1545                verbose = self.verbose,
1546                format = 'asc')
1547
1548        prjfile = self.domain.get_name() + '_elevation.prj'
1549        ascfile = self.domain.get_name() + '_elevation.asc'
1550       
1551        #Check asc file
1552        ascid = open(ascfile)
1553        lines = ascid.readlines()
1554        ascid.close()
1555
1556        L = lines[2].strip().split()
1557        assert L[0].strip().lower() == 'xllcorner'
1558        assert allclose(float(L[1].strip().lower()), 308500)
1559
1560        L = lines[3].strip().split()
1561        assert L[0].strip().lower() == 'yllcorner'
1562        assert allclose(float(L[1].strip().lower()), 6189000)
1563
1564        #print "ascfile", ascfile
1565        #Check grid values
1566        for j in range(5):
1567            L = lines[6+j].strip().split()
1568            y = (4-j) * cellsize
1569            for i in range(5):
1570                #print " -i*cellsize - y",  -i*cellsize - y
1571                #print "float(L[i])", float(L[i])
1572                assert allclose(float(L[i]), -i*cellsize - y)
1573               
1574        #Cleanup
1575        os.remove(prjfile)
1576        os.remove(ascfile)
1577       
1578        #Check asc file
1579        ascfile = self.domain.get_name() + '_depth.asc'
1580        prjfile = self.domain.get_name() + '_depth.prj'
1581        ascid = open(ascfile)
1582        lines = ascid.readlines()
1583        ascid.close()
1584
1585        L = lines[2].strip().split()
1586        assert L[0].strip().lower() == 'xllcorner'
1587        assert allclose(float(L[1].strip().lower()), 308500)
1588
1589        L = lines[3].strip().split()
1590        assert L[0].strip().lower() == 'yllcorner'
1591        assert allclose(float(L[1].strip().lower()), 6189000)
1592
1593        #Check grid values
1594        for j in range(5):
1595            L = lines[6+j].strip().split()
1596            y = (4-j) * cellsize
1597            for i in range(5):
1598                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1599
1600        #Cleanup
1601        os.remove(prjfile)
1602        os.remove(ascfile)
1603        os.remove(swwfile)
1604
1605
1606    def test_export_gridIII(self):
1607        """
1608        test_export_gridIII
1609        Test that sww information can be converted correctly to asc/prj
1610        format readable by e.g. ArcView
1611        """
1612
1613        import time, os
1614        from Numeric import array, zeros, allclose, Float, concatenate
1615        from Scientific.IO.NetCDF import NetCDFFile
1616
1617        try:
1618            os.remove('teg*.sww')
1619        except:
1620            pass
1621
1622        #Setup
1623       
1624        self.domain.set_name('tegIII')
1625
1626        swwfile = self.domain.get_name() + '.sww'
1627
1628        self.domain.set_datadir('.')
1629        self.domain.format = 'sww'
1630        self.domain.smooth = True
1631        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1632        self.domain.set_quantity('stage', 1.0)
1633
1634        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1635       
1636        sww = get_dataobject(self.domain)
1637        sww.store_connectivity()
1638        sww.store_timestep('stage')
1639        self.domain.evolve_to_end(finaltime = 0.01)
1640        sww.store_timestep('stage')
1641
1642        cellsize = 0.25
1643        #Check contents
1644        #Get NetCDF
1645
1646        fid = NetCDFFile(sww.filename, 'r')
1647
1648        # Get the variables
1649        x = fid.variables['x'][:]
1650        y = fid.variables['y'][:]
1651        z = fid.variables['elevation'][:]
1652        time = fid.variables['time'][:]
1653        stage = fid.variables['stage'][:]
1654
1655        fid.close()
1656
1657        #Export to ascii/prj files
1658        extra_name_out = 'yeah'
1659        if True:
1660            export_grid(self.domain.get_name(),
1661                        quantities = ['elevation', 'depth'],
1662                        extra_name_out = extra_name_out,
1663                        cellsize = cellsize,
1664                        verbose = self.verbose,
1665                        format = 'asc')
1666
1667        else:
1668            export_grid(self.domain.get_name(),
1669                quantities = ['depth'],
1670                cellsize = cellsize,
1671                verbose = self.verbose,
1672                format = 'asc')
1673
1674
1675            export_grid(self.domain.get_name(),
1676                quantities = ['elevation'],
1677                cellsize = cellsize,
1678                verbose = self.verbose,
1679                format = 'asc')
1680
1681        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
1682        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
1683       
1684        #Check asc file
1685        ascid = open(ascfile)
1686        lines = ascid.readlines()
1687        ascid.close()
1688
1689        L = lines[2].strip().split()
1690        assert L[0].strip().lower() == 'xllcorner'
1691        assert allclose(float(L[1].strip().lower()), 308500)
1692
1693        L = lines[3].strip().split()
1694        assert L[0].strip().lower() == 'yllcorner'
1695        assert allclose(float(L[1].strip().lower()), 6189000)
1696
1697        #print "ascfile", ascfile
1698        #Check grid values
1699        for j in range(5):
1700            L = lines[6+j].strip().split()
1701            y = (4-j) * cellsize
1702            for i in range(5):
1703                #print " -i*cellsize - y",  -i*cellsize - y
1704                #print "float(L[i])", float(L[i])
1705                assert allclose(float(L[i]), -i*cellsize - y)
1706               
1707        #Cleanup
1708        os.remove(prjfile)
1709        os.remove(ascfile)
1710       
1711        #Check asc file
1712        ascfile = self.domain.get_name() + '_depth_yeah.asc'
1713        prjfile = self.domain.get_name() + '_depth_yeah.prj'
1714        ascid = open(ascfile)
1715        lines = ascid.readlines()
1716        ascid.close()
1717
1718        L = lines[2].strip().split()
1719        assert L[0].strip().lower() == 'xllcorner'
1720        assert allclose(float(L[1].strip().lower()), 308500)
1721
1722        L = lines[3].strip().split()
1723        assert L[0].strip().lower() == 'yllcorner'
1724        assert allclose(float(L[1].strip().lower()), 6189000)
1725
1726        #Check grid values
1727        for j in range(5):
1728            L = lines[6+j].strip().split()
1729            y = (4-j) * cellsize
1730            for i in range(5):
1731                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1732
1733        #Cleanup
1734        os.remove(prjfile)
1735        os.remove(ascfile)
1736        os.remove(swwfile)
1737
1738    def test_export_grid_bad(self):
1739        """Test that sww information can be converted correctly to asc/prj
1740        format readable by e.g. ArcView
1741        """
1742
1743        try:
1744            export_grid('a_small_round-egg',
1745                        quantities = ['elevation', 'depth'],
1746                        cellsize = 99,
1747                        verbose = self.verbose,
1748                        format = 'asc')
1749        except IOError:
1750            pass
1751        else:
1752            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
1753
1754    def test_export_grid_parallel(self):
1755        """Test that sww information can be converted correctly to asc/prj
1756        format readable by e.g. ArcView
1757        """
1758
1759        import time, os
1760        from Numeric import array, zeros, allclose, Float, concatenate
1761        from Scientific.IO.NetCDF import NetCDFFile
1762
1763        base_name = 'tegp'
1764        #Setup
1765        self.domain.set_name(base_name+'_P0_8')
1766        swwfile = self.domain.get_name() + '.sww'
1767
1768        self.domain.set_datadir('.')
1769        self.domain.format = 'sww'
1770        self.domain.smooth = True
1771        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1772        self.domain.set_quantity('stage', 1.0)
1773
1774        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1775
1776        sww = get_dataobject(self.domain)
1777        sww.store_connectivity()
1778        sww.store_timestep('stage')
1779        self.domain.evolve_to_end(finaltime = 0.0001)
1780        #Setup
1781        self.domain.set_name(base_name+'_P1_8')
1782        swwfile2 = self.domain.get_name() + '.sww'
1783        sww = get_dataobject(self.domain)
1784        sww.store_connectivity()
1785        sww.store_timestep('stage')
1786        self.domain.evolve_to_end(finaltime = 0.0002)
1787        sww.store_timestep('stage')
1788
1789        cellsize = 0.25
1790        #Check contents
1791        #Get NetCDF
1792
1793        fid = NetCDFFile(sww.filename, 'r')
1794
1795        # Get the variables
1796        x = fid.variables['x'][:]
1797        y = fid.variables['y'][:]
1798        z = fid.variables['elevation'][:]
1799        time = fid.variables['time'][:]
1800        stage = fid.variables['stage'][:]
1801
1802        fid.close()
1803
1804        #Export to ascii/prj files
1805        extra_name_out = 'yeah'
1806        export_grid(base_name,
1807                    quantities = ['elevation', 'depth'],
1808                    extra_name_out = extra_name_out,
1809                    cellsize = cellsize,
1810                    verbose = self.verbose,
1811                    format = 'asc')
1812
1813        prjfile = base_name + '_P0_8_elevation_yeah.prj'
1814        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
1815        #Check asc file
1816        ascid = open(ascfile)
1817        lines = ascid.readlines()
1818        ascid.close()
1819        #Check grid values
1820        for j in range(5):
1821            L = lines[6+j].strip().split()
1822            y = (4-j) * cellsize
1823            for i in range(5):
1824                #print " -i*cellsize - y",  -i*cellsize - y
1825                #print "float(L[i])", float(L[i])
1826                assert allclose(float(L[i]), -i*cellsize - y)               
1827        #Cleanup
1828        os.remove(prjfile)
1829        os.remove(ascfile)
1830
1831        prjfile = base_name + '_P1_8_elevation_yeah.prj'
1832        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
1833        #Check asc file
1834        ascid = open(ascfile)
1835        lines = ascid.readlines()
1836        ascid.close()
1837        #Check grid values
1838        for j in range(5):
1839            L = lines[6+j].strip().split()
1840            y = (4-j) * cellsize
1841            for i in range(5):
1842                #print " -i*cellsize - y",  -i*cellsize - y
1843                #print "float(L[i])", float(L[i])
1844                assert allclose(float(L[i]), -i*cellsize - y)               
1845        #Cleanup
1846        os.remove(prjfile)
1847        os.remove(ascfile)
1848        os.remove(swwfile)
1849
1850        #Check asc file
1851        ascfile = base_name + '_P0_8_depth_yeah.asc'
1852        prjfile = base_name + '_P0_8_depth_yeah.prj'
1853        ascid = open(ascfile)
1854        lines = ascid.readlines()
1855        ascid.close()
1856        #Check grid values
1857        for j in range(5):
1858            L = lines[6+j].strip().split()
1859            y = (4-j) * cellsize
1860            for i in range(5):
1861                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1862        #Cleanup
1863        os.remove(prjfile)
1864        os.remove(ascfile)
1865
1866        #Check asc file
1867        ascfile = base_name + '_P1_8_depth_yeah.asc'
1868        prjfile = base_name + '_P1_8_depth_yeah.prj'
1869        ascid = open(ascfile)
1870        lines = ascid.readlines()
1871        ascid.close()
1872        #Check grid values
1873        for j in range(5):
1874            L = lines[6+j].strip().split()
1875            y = (4-j) * cellsize
1876            for i in range(5):
1877                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1878        #Cleanup
1879        os.remove(prjfile)
1880        os.remove(ascfile)
1881        os.remove(swwfile2)
1882
1883    def test_sww2dem_larger(self):
1884        """Test that sww information can be converted correctly to asc/prj
1885        format readable by e.g. ArcView. Here:
1886
1887        ncols         11
1888        nrows         11
1889        xllcorner     308500
1890        yllcorner     6189000
1891        cellsize      10.000000
1892        NODATA_value  -9999
1893        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1894         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1895         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1896         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
1897         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
1898         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
1899         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
1900         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
1901         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
1902         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
1903           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
1904
1905        """
1906
1907        import time, os
1908        from Numeric import array, zeros, allclose, Float, concatenate
1909        from Scientific.IO.NetCDF import NetCDFFile
1910
1911        #Setup
1912
1913        from mesh_factory import rectangular
1914
1915        #Create basic mesh (100m x 100m)
1916        points, vertices, boundary = rectangular(2, 2, 100, 100)
1917
1918        #Create shallow water domain
1919        domain = Domain(points, vertices, boundary)
1920        domain.default_order = 2
1921
1922        domain.set_name('datatest')
1923
1924        prjfile = domain.get_name() + '_elevation.prj'
1925        ascfile = domain.get_name() + '_elevation.asc'
1926        swwfile = domain.get_name() + '.sww'
1927
1928        domain.set_datadir('.')
1929        domain.format = 'sww'
1930        domain.smooth = True
1931        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1932
1933        #
1934        domain.set_quantity('elevation', lambda x,y: -x-y)
1935        domain.set_quantity('stage', 0)
1936
1937        B = Transmissive_boundary(domain)
1938        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1939
1940
1941        #
1942        sww = get_dataobject(domain)
1943        sww.store_connectivity()
1944        sww.store_timestep('stage')
1945       
1946        domain.tight_slope_limiters = 1
1947        domain.evolve_to_end(finaltime = 0.01)
1948        sww.store_timestep('stage')
1949
1950        cellsize = 10  #10m grid
1951
1952
1953        #Check contents
1954        #Get NetCDF
1955
1956        fid = NetCDFFile(sww.filename, 'r')
1957
1958        # Get the variables
1959        x = fid.variables['x'][:]
1960        y = fid.variables['y'][:]
1961        z = fid.variables['elevation'][:]
1962        time = fid.variables['time'][:]
1963        stage = fid.variables['stage'][:]
1964
1965
1966        #Export to ascii/prj files
1967        sww2dem(domain.get_name(),
1968                quantity = 'elevation',
1969                cellsize = cellsize,
1970                verbose = self.verbose,
1971                format = 'asc')
1972
1973
1974        #Check prj (meta data)
1975        prjid = open(prjfile)
1976        lines = prjid.readlines()
1977        prjid.close()
1978
1979        L = lines[0].strip().split()
1980        assert L[0].strip().lower() == 'projection'
1981        assert L[1].strip().lower() == 'utm'
1982
1983        L = lines[1].strip().split()
1984        assert L[0].strip().lower() == 'zone'
1985        assert L[1].strip().lower() == '56'
1986
1987        L = lines[2].strip().split()
1988        assert L[0].strip().lower() == 'datum'
1989        assert L[1].strip().lower() == 'wgs84'
1990
1991        L = lines[3].strip().split()
1992        assert L[0].strip().lower() == 'zunits'
1993        assert L[1].strip().lower() == 'no'
1994
1995        L = lines[4].strip().split()
1996        assert L[0].strip().lower() == 'units'
1997        assert L[1].strip().lower() == 'meters'
1998
1999        L = lines[5].strip().split()
2000        assert L[0].strip().lower() == 'spheroid'
2001        assert L[1].strip().lower() == 'wgs84'
2002
2003        L = lines[6].strip().split()
2004        assert L[0].strip().lower() == 'xshift'
2005        assert L[1].strip().lower() == '500000'
2006
2007        L = lines[7].strip().split()
2008        assert L[0].strip().lower() == 'yshift'
2009        assert L[1].strip().lower() == '10000000'
2010
2011        L = lines[8].strip().split()
2012        assert L[0].strip().lower() == 'parameters'
2013
2014
2015        #Check asc file
2016        ascid = open(ascfile)
2017        lines = ascid.readlines()
2018        ascid.close()
2019
2020        L = lines[0].strip().split()
2021        assert L[0].strip().lower() == 'ncols'
2022        assert L[1].strip().lower() == '11'
2023
2024        L = lines[1].strip().split()
2025        assert L[0].strip().lower() == 'nrows'
2026        assert L[1].strip().lower() == '11'
2027
2028        L = lines[2].strip().split()
2029        assert L[0].strip().lower() == 'xllcorner'
2030        assert allclose(float(L[1].strip().lower()), 308500)
2031
2032        L = lines[3].strip().split()
2033        assert L[0].strip().lower() == 'yllcorner'
2034        assert allclose(float(L[1].strip().lower()), 6189000)
2035
2036        L = lines[4].strip().split()
2037        assert L[0].strip().lower() == 'cellsize'
2038        assert allclose(float(L[1].strip().lower()), cellsize)
2039
2040        L = lines[5].strip().split()
2041        assert L[0].strip() == 'NODATA_value'
2042        assert L[1].strip().lower() == '-9999'
2043
2044        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
2045        for i, line in enumerate(lines[6:]):
2046            for j, value in enumerate( line.split() ):
2047                #assert allclose(float(value), -(10-i+j)*cellsize)
2048                assert float(value) == -(10-i+j)*cellsize
2049
2050
2051        fid.close()
2052
2053        #Cleanup
2054        os.remove(prjfile)
2055        os.remove(ascfile)
2056        os.remove(swwfile)
2057
2058
2059
2060
2061    def test_sww2dem_boundingbox(self):
2062        """Test that sww information can be converted correctly to asc/prj
2063        format readable by e.g. ArcView.
2064        This will test that mesh can be restricted by bounding box
2065
2066        Original extent is 100m x 100m:
2067
2068        Eastings:   308500 -  308600
2069        Northings: 6189000 - 6189100
2070
2071        Bounding box changes this to the 50m x 50m square defined by
2072
2073        Eastings:   308530 -  308570
2074        Northings: 6189050 - 6189100
2075
2076        The cropped values should be
2077
2078         -130 -140 -150 -160 -170
2079         -120 -130 -140 -150 -160
2080         -110 -120 -130 -140 -150
2081         -100 -110 -120 -130 -140
2082          -90 -100 -110 -120 -130
2083          -80  -90 -100 -110 -120
2084
2085        and the new lower reference point should be
2086        Eastings:   308530
2087        Northings: 6189050
2088
2089        Original dataset is the same as in test_sww2dem_larger()
2090
2091        """
2092
2093        import time, os
2094        from Numeric import array, zeros, allclose, Float, concatenate
2095        from Scientific.IO.NetCDF import NetCDFFile
2096
2097        #Setup
2098
2099        from mesh_factory import rectangular
2100
2101        #Create basic mesh (100m x 100m)
2102        points, vertices, boundary = rectangular(2, 2, 100, 100)
2103
2104        #Create shallow water domain
2105        domain = Domain(points, vertices, boundary)
2106        domain.default_order = 2
2107
2108        domain.set_name('datatest')
2109
2110        prjfile = domain.get_name() + '_elevation.prj'
2111        ascfile = domain.get_name() + '_elevation.asc'
2112        swwfile = domain.get_name() + '.sww'
2113
2114        domain.set_datadir('.')
2115        domain.format = 'sww'
2116        domain.smooth = True
2117        domain.geo_reference = Geo_reference(56, 308500, 6189000)
2118
2119        #
2120        domain.set_quantity('elevation', lambda x,y: -x-y)
2121        domain.set_quantity('stage', 0)
2122
2123        B = Transmissive_boundary(domain)
2124        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
2125
2126
2127        #
2128        sww = get_dataobject(domain)
2129        sww.store_connectivity()
2130        sww.store_timestep('stage')
2131
2132        #domain.tight_slope_limiters = 1
2133        domain.evolve_to_end(finaltime = 0.01)
2134        sww.store_timestep('stage')
2135
2136        cellsize = 10  #10m grid
2137
2138
2139        #Check contents
2140        #Get NetCDF
2141
2142        fid = NetCDFFile(sww.filename, 'r')
2143
2144        # Get the variables
2145        x = fid.variables['x'][:]
2146        y = fid.variables['y'][:]
2147        z = fid.variables['elevation'][:]
2148        time = fid.variables['time'][:]
2149        stage = fid.variables['stage'][:]
2150
2151
2152        #Export to ascii/prj files
2153        sww2dem(domain.get_name(),
2154                quantity = 'elevation',
2155                cellsize = cellsize,
2156                easting_min = 308530,
2157                easting_max = 308570,
2158                northing_min = 6189050,
2159                northing_max = 6189100,
2160                verbose = self.verbose,
2161                format = 'asc')
2162
2163        fid.close()
2164
2165
2166        #Check prj (meta data)
2167        prjid = open(prjfile)
2168        lines = prjid.readlines()
2169        prjid.close()
2170
2171        L = lines[0].strip().split()
2172        assert L[0].strip().lower() == 'projection'
2173        assert L[1].strip().lower() == 'utm'
2174
2175        L = lines[1].strip().split()
2176        assert L[0].strip().lower() == 'zone'
2177        assert L[1].strip().lower() == '56'
2178
2179        L = lines[2].strip().split()
2180        assert L[0].strip().lower() == 'datum'
2181        assert L[1].strip().lower() == 'wgs84'
2182
2183        L = lines[3].strip().split()
2184        assert L[0].strip().lower() == 'zunits'
2185        assert L[1].strip().lower() == 'no'
2186
2187        L = lines[4].strip().split()
2188        assert L[0].strip().lower() == 'units'
2189        assert L[1].strip().lower() == 'meters'
2190
2191        L = lines[5].strip().split()
2192        assert L[0].strip().lower() == 'spheroid'
2193        assert L[1].strip().lower() == 'wgs84'
2194
2195        L = lines[6].strip().split()
2196        assert L[0].strip().lower() == 'xshift'
2197        assert L[1].strip().lower() == '500000'
2198
2199        L = lines[7].strip().split()
2200        assert L[0].strip().lower() == 'yshift'
2201        assert L[1].strip().lower() == '10000000'
2202
2203        L = lines[8].strip().split()
2204        assert L[0].strip().lower() == 'parameters'
2205
2206
2207        #Check asc file
2208        ascid = open(ascfile)
2209        lines = ascid.readlines()
2210        ascid.close()
2211
2212        L = lines[0].strip().split()
2213        assert L[0].strip().lower() == 'ncols'
2214        assert L[1].strip().lower() == '5'
2215
2216        L = lines[1].strip().split()
2217        assert L[0].strip().lower() == 'nrows'
2218        assert L[1].strip().lower() == '6'
2219
2220        L = lines[2].strip().split()
2221        assert L[0].strip().lower() == 'xllcorner'
2222        assert allclose(float(L[1].strip().lower()), 308530)
2223
2224        L = lines[3].strip().split()
2225        assert L[0].strip().lower() == 'yllcorner'
2226        assert allclose(float(L[1].strip().lower()), 6189050)
2227
2228        L = lines[4].strip().split()
2229        assert L[0].strip().lower() == 'cellsize'
2230        assert allclose(float(L[1].strip().lower()), cellsize)
2231
2232        L = lines[5].strip().split()
2233        assert L[0].strip() == 'NODATA_value'
2234        assert L[1].strip().lower() == '-9999'
2235
2236        #Check grid values
2237        for i, line in enumerate(lines[6:]):
2238            for j, value in enumerate( line.split() ):
2239                #assert float(value) == -(10-i+j)*cellsize
2240                assert float(value) == -(10-i+j+3)*cellsize
2241
2242
2243
2244        #Cleanup
2245        os.remove(prjfile)
2246        os.remove(ascfile)
2247        os.remove(swwfile)
2248
2249
2250
2251    def test_sww2dem_asc_stage_reduction(self):
2252        """Test that sww information can be converted correctly to asc/prj
2253        format readable by e.g. ArcView
2254
2255        This tests the reduction of quantity stage using min
2256        """
2257
2258        import time, os
2259        from Numeric import array, zeros, allclose, Float, concatenate
2260        from Scientific.IO.NetCDF import NetCDFFile
2261
2262        #Setup
2263        self.domain.set_name('datatest')
2264
2265        prjfile = self.domain.get_name() + '_stage.prj'
2266        ascfile = self.domain.get_name() + '_stage.asc'
2267        swwfile = self.domain.get_name() + '.sww'
2268
2269        self.domain.set_datadir('.')
2270        self.domain.format = 'sww'
2271        self.domain.smooth = True
2272        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2273
2274        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2275
2276
2277        sww = get_dataobject(self.domain)
2278        sww.store_connectivity()
2279        sww.store_timestep('stage')
2280
2281        #self.domain.tight_slope_limiters = 1
2282        self.domain.evolve_to_end(finaltime = 0.01)
2283        sww.store_timestep('stage')
2284
2285        cellsize = 0.25
2286        #Check contents
2287        #Get NetCDF
2288
2289        fid = NetCDFFile(sww.filename, 'r')
2290
2291        # Get the variables
2292        x = fid.variables['x'][:]
2293        y = fid.variables['y'][:]
2294        z = fid.variables['elevation'][:]
2295        time = fid.variables['time'][:]
2296        stage = fid.variables['stage'][:]
2297
2298
2299        #Export to ascii/prj files
2300        sww2dem(self.domain.get_name(),
2301                quantity = 'stage',
2302                cellsize = cellsize,
2303                reduction = min,
2304                format = 'asc',
2305                verbose=self.verbose)
2306
2307
2308        #Check asc file
2309        ascid = open(ascfile)
2310        lines = ascid.readlines()
2311        ascid.close()
2312
2313        L = lines[0].strip().split()
2314        assert L[0].strip().lower() == 'ncols'
2315        assert L[1].strip().lower() == '5'
2316
2317        L = lines[1].strip().split()
2318        assert L[0].strip().lower() == 'nrows'
2319        assert L[1].strip().lower() == '5'
2320
2321        L = lines[2].strip().split()
2322        assert L[0].strip().lower() == 'xllcorner'
2323        assert allclose(float(L[1].strip().lower()), 308500)
2324
2325        L = lines[3].strip().split()
2326        assert L[0].strip().lower() == 'yllcorner'
2327        assert allclose(float(L[1].strip().lower()), 6189000)
2328
2329        L = lines[4].strip().split()
2330        assert L[0].strip().lower() == 'cellsize'
2331        assert allclose(float(L[1].strip().lower()), cellsize)
2332
2333        L = lines[5].strip().split()
2334        assert L[0].strip() == 'NODATA_value'
2335        assert L[1].strip().lower() == '-9999'
2336
2337
2338        #Check grid values (where applicable)
2339        for j in range(5):
2340            if j%2 == 0:
2341                L = lines[6+j].strip().split()
2342                jj = 4-j
2343                for i in range(5):
2344                    if i%2 == 0:
2345                        index = jj/2 + i/2*3
2346                        val0 = stage[0,index]
2347                        val1 = stage[1,index]
2348
2349                        #print i, j, index, ':', L[i], val0, val1
2350                        assert allclose(float(L[i]), min(val0, val1))
2351
2352
2353        fid.close()
2354
2355        #Cleanup
2356        os.remove(prjfile)
2357        os.remove(ascfile)
2358        os.remove(swwfile)
2359
2360
2361
2362    def test_sww2dem_asc_derived_quantity(self):
2363        """Test that sww information can be converted correctly to asc/prj
2364        format readable by e.g. ArcView
2365
2366        This tests the use of derived quantities
2367        """
2368
2369        import time, os
2370        from Numeric import array, zeros, allclose, Float, concatenate
2371        from Scientific.IO.NetCDF import NetCDFFile
2372
2373        #Setup
2374        self.domain.set_name('datatest')
2375
2376        prjfile = self.domain.get_name() + '_depth.prj'
2377        ascfile = self.domain.get_name() + '_depth.asc'
2378        swwfile = self.domain.get_name() + '.sww'
2379
2380        self.domain.set_datadir('.')
2381        self.domain.format = 'sww'
2382        self.domain.smooth = True
2383        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2384        self.domain.set_quantity('stage', 0.0)
2385
2386        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2387
2388
2389        sww = get_dataobject(self.domain)
2390        sww.store_connectivity()
2391        sww.store_timestep('stage')
2392
2393        #self.domain.tight_slope_limiters = 1
2394        self.domain.evolve_to_end(finaltime = 0.01)
2395        sww.store_timestep('stage')
2396
2397        cellsize = 0.25
2398        #Check contents
2399        #Get NetCDF
2400
2401        fid = NetCDFFile(sww.filename, 'r')
2402
2403        # Get the variables
2404        x = fid.variables['x'][:]
2405        y = fid.variables['y'][:]
2406        z = fid.variables['elevation'][:]
2407        time = fid.variables['time'][:]
2408        stage = fid.variables['stage'][:]
2409
2410
2411        #Export to ascii/prj files
2412        sww2dem(self.domain.get_name(),
2413                basename_out = 'datatest_depth',
2414                quantity = 'stage - elevation',
2415                cellsize = cellsize,
2416                reduction = min,
2417                format = 'asc',
2418                verbose = self.verbose)
2419
2420
2421        #Check asc file
2422        ascid = open(ascfile)
2423        lines = ascid.readlines()
2424        ascid.close()
2425
2426        L = lines[0].strip().split()
2427        assert L[0].strip().lower() == 'ncols'
2428        assert L[1].strip().lower() == '5'
2429
2430        L = lines[1].strip().split()
2431        assert L[0].strip().lower() == 'nrows'
2432        assert L[1].strip().lower() == '5'
2433
2434        L = lines[2].strip().split()
2435        assert L[0].strip().lower() == 'xllcorner'
2436        assert allclose(float(L[1].strip().lower()), 308500)
2437
2438        L = lines[3].strip().split()
2439        assert L[0].strip().lower() == 'yllcorner'
2440        assert allclose(float(L[1].strip().lower()), 6189000)
2441
2442        L = lines[4].strip().split()
2443        assert L[0].strip().lower() == 'cellsize'
2444        assert allclose(float(L[1].strip().lower()), cellsize)
2445
2446        L = lines[5].strip().split()
2447        assert L[0].strip() == 'NODATA_value'
2448        assert L[1].strip().lower() == '-9999'
2449
2450
2451        #Check grid values (where applicable)
2452        for j in range(5):
2453            if j%2 == 0:
2454                L = lines[6+j].strip().split()
2455                jj = 4-j
2456                for i in range(5):
2457                    if i%2 == 0:
2458                        index = jj/2 + i/2*3
2459                        val0 = stage[0,index] - z[index]
2460                        val1 = stage[1,index] - z[index]
2461
2462                        #print i, j, index, ':', L[i], val0, val1
2463                        assert allclose(float(L[i]), min(val0, val1))
2464
2465
2466        fid.close()
2467
2468        #Cleanup
2469        os.remove(prjfile)
2470        os.remove(ascfile)
2471        os.remove(swwfile)
2472
2473
2474
2475
2476
2477    def test_sww2dem_asc_missing_points(self):
2478        """Test that sww information can be converted correctly to asc/prj
2479        format readable by e.g. ArcView
2480
2481        This test includes the writing of missing values
2482        """
2483
2484        import time, os
2485        from Numeric import array, zeros, allclose, Float, concatenate
2486        from Scientific.IO.NetCDF import NetCDFFile
2487
2488        #Setup mesh not coinciding with rectangle.
2489        #This will cause missing values to occur in gridded data
2490
2491
2492        points = [                        [1.0, 1.0],
2493                              [0.5, 0.5], [1.0, 0.5],
2494                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
2495
2496        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
2497
2498        #Create shallow water domain
2499        domain = Domain(points, vertices)
2500        domain.default_order=2
2501
2502
2503        #Set some field values
2504        domain.set_quantity('elevation', lambda x,y: -x-y)
2505        domain.set_quantity('friction', 0.03)
2506
2507
2508        ######################
2509        # Boundary conditions
2510        B = Transmissive_boundary(domain)
2511        domain.set_boundary( {'exterior': B} )
2512
2513
2514        ######################
2515        #Initial condition - with jumps
2516
2517        bed = domain.quantities['elevation'].vertex_values
2518        stage = zeros(bed.shape, Float)
2519
2520        h = 0.3
2521        for i in range(stage.shape[0]):
2522            if i % 2 == 0:
2523                stage[i,:] = bed[i,:] + h
2524            else:
2525                stage[i,:] = bed[i,:]
2526
2527        domain.set_quantity('stage', stage)
2528        domain.distribute_to_vertices_and_edges()
2529
2530        domain.set_name('datatest')
2531
2532        prjfile = domain.get_name() + '_elevation.prj'
2533        ascfile = domain.get_name() + '_elevation.asc'
2534        swwfile = domain.get_name() + '.sww'
2535
2536        domain.set_datadir('.')
2537        domain.format = 'sww'
2538        domain.smooth = True
2539
2540        domain.geo_reference = Geo_reference(56,308500,6189000)
2541
2542        sww = get_dataobject(domain)
2543        sww.store_connectivity()
2544        sww.store_timestep('stage')
2545
2546        cellsize = 0.25
2547        #Check contents
2548        #Get NetCDF
2549
2550        fid = NetCDFFile(swwfile, 'r')
2551
2552        # Get the variables
2553        x = fid.variables['x'][:]
2554        y = fid.variables['y'][:]
2555        z = fid.variables['elevation'][:]
2556        time = fid.variables['time'][:]
2557
2558        try:
2559            geo_reference = Geo_reference(NetCDFObject=fid)
2560        except AttributeError, e:
2561            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
2562
2563        #Export to ascii/prj files
2564        sww2dem(domain.get_name(),
2565                quantity = 'elevation',
2566                cellsize = cellsize,
2567                verbose = self.verbose,
2568                format = 'asc')
2569
2570
2571        #Check asc file
2572        ascid = open(ascfile)
2573        lines = ascid.readlines()
2574        ascid.close()
2575
2576        L = lines[0].strip().split()
2577        assert L[0].strip().lower() == 'ncols'
2578        assert L[1].strip().lower() == '5'
2579
2580        L = lines[1].strip().split()
2581        assert L[0].strip().lower() == 'nrows'
2582        assert L[1].strip().lower() == '5'
2583
2584        L = lines[2].strip().split()
2585        assert L[0].strip().lower() == 'xllcorner'
2586        assert allclose(float(L[1].strip().lower()), 308500)
2587
2588        L = lines[3].strip().split()
2589        assert L[0].strip().lower() == 'yllcorner'
2590        assert allclose(float(L[1].strip().lower()), 6189000)
2591
2592        L = lines[4].strip().split()
2593        assert L[0].strip().lower() == 'cellsize'
2594        assert allclose(float(L[1].strip().lower()), cellsize)
2595
2596        L = lines[5].strip().split()
2597        assert L[0].strip() == 'NODATA_value'
2598        assert L[1].strip().lower() == '-9999'
2599
2600        #Check grid values
2601        for j in range(5):
2602            L = lines[6+j].strip().split()
2603            assert len(L) == 5
2604            y = (4-j) * cellsize
2605
2606            for i in range(5):
2607                #print i
2608                if i+j >= 4:
2609                    assert allclose(float(L[i]), -i*cellsize - y)
2610                else:
2611                    #Missing values
2612                    assert allclose(float(L[i]), -9999)
2613
2614
2615
2616        fid.close()
2617
2618        #Cleanup
2619        os.remove(prjfile)
2620        os.remove(ascfile)
2621        os.remove(swwfile)
2622
2623    def test_sww2ers_simple(self):
2624        """Test that sww information can be converted correctly to asc/prj
2625        format readable by e.g. ArcView
2626        """
2627
2628        import time, os
2629        from Numeric import array, zeros, allclose, Float, concatenate
2630        from Scientific.IO.NetCDF import NetCDFFile
2631
2632
2633        NODATA_value = 1758323
2634
2635        #Setup
2636        self.domain.set_name('datatest')
2637
2638        headerfile = self.domain.get_name() + '.ers'
2639        swwfile = self.domain.get_name() + '.sww'
2640
2641        self.domain.set_datadir('.')
2642        self.domain.format = 'sww'
2643        self.domain.smooth = True
2644        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2645
2646        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2647
2648        sww = get_dataobject(self.domain)
2649        sww.store_connectivity()
2650        sww.store_timestep('stage')
2651
2652        #self.domain.tight_slope_limiters = 1
2653        self.domain.evolve_to_end(finaltime = 0.01)
2654        sww.store_timestep('stage')
2655
2656        cellsize = 0.25
2657        #Check contents
2658        #Get NetCDF
2659
2660        fid = NetCDFFile(sww.filename, 'r')
2661
2662        # Get the variables
2663        x = fid.variables['x'][:]
2664        y = fid.variables['y'][:]
2665        z = fid.variables['elevation'][:]
2666        time = fid.variables['time'][:]
2667        stage = fid.variables['stage'][:]
2668
2669
2670        #Export to ers files
2671        sww2dem(self.domain.get_name(),
2672                quantity = 'elevation',
2673                cellsize = cellsize,
2674                NODATA_value = NODATA_value,
2675                verbose = self.verbose,
2676                format = 'ers')
2677
2678        #Check header data
2679        from ermapper_grids import read_ermapper_header, read_ermapper_data
2680
2681        header = read_ermapper_header(self.domain.get_name() + '_elevation.ers')
2682        #print header
2683        assert header['projection'].lower() == '"utm-56"'
2684        assert header['datum'].lower() == '"wgs84"'
2685        assert header['units'].lower() == '"meters"'
2686        assert header['value'].lower() == '"elevation"'
2687        assert header['xdimension'] == '0.25'
2688        assert header['ydimension'] == '0.25'
2689        assert float(header['eastings']) == 308500.0   #xllcorner
2690        assert float(header['northings']) == 6189000.0 #yllcorner
2691        assert int(header['nroflines']) == 5
2692        assert int(header['nrofcellsperline']) == 5
2693        assert int(header['nullcellvalue']) == NODATA_value
2694        #FIXME - there is more in the header
2695
2696
2697        #Check grid data
2698        grid = read_ermapper_data(self.domain.get_name() + '_elevation')
2699
2700        #FIXME (Ole): Why is this the desired reference grid for -x-y?
2701        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
2702                    -1,    -1.25, -1.5,  -1.75, -2.0,
2703                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
2704                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
2705                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
2706
2707
2708        #print grid
2709        assert allclose(grid, ref_grid)
2710
2711        fid.close()
2712
2713        #Cleanup
2714        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
2715        #Done (Ole) - it was because sww2ers didn't close it's sww file
2716        os.remove(sww.filename)
2717        os.remove(self.domain.get_name() + '_elevation')
2718        os.remove(self.domain.get_name() + '_elevation.ers')
2719
2720
2721
2722    def test_sww2pts_centroids(self):
2723        """Test that sww information can be converted correctly to pts data at specified coordinates
2724        - in this case, the centroids.
2725        """
2726
2727        import time, os
2728        from Numeric import array, zeros, allclose, Float, concatenate, NewAxis
2729        from Scientific.IO.NetCDF import NetCDFFile
2730        from anuga.geospatial_data.geospatial_data import Geospatial_data
2731
2732        # Used for points that lie outside mesh
2733        NODATA_value = 1758323
2734
2735        # Setup
2736        self.domain.set_name('datatest')
2737
2738        ptsfile = self.domain.get_name() + '_elevation.pts'
2739        swwfile = self.domain.get_name() + '.sww'
2740
2741        self.domain.set_datadir('.')
2742        self.domain.format = 'sww'
2743        self.smooth = True #self.set_store_vertices_uniquely(False)
2744        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2745
2746        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2747
2748        sww = get_dataobject(self.domain)
2749        sww.store_connectivity()
2750        sww.store_timestep('stage')
2751
2752        #self.domain.tight_slope_limiters = 1
2753        self.domain.evolve_to_end(finaltime = 0.01)
2754        sww.store_timestep('stage')
2755
2756        # Check contents in NetCDF
2757        fid = NetCDFFile(sww.filename, 'r')
2758
2759        # Get the variables
2760        x = fid.variables['x'][:]
2761        y = fid.variables['y'][:]
2762        elevation = fid.variables['elevation'][:]
2763        time = fid.variables['time'][:]
2764        stage = fid.variables['stage'][:]
2765
2766        volumes = fid.variables['volumes'][:]
2767
2768
2769        # Invoke interpolation for vertex points       
2770        points = concatenate( (x[:,NewAxis],y[:,NewAxis]), axis=1 )
2771        sww2pts(self.domain.get_name(),
2772                quantity = 'elevation',
2773                data_points = points,
2774                NODATA_value = NODATA_value,
2775                verbose = self.verbose)
2776        ref_point_values = elevation
2777        point_values = Geospatial_data(ptsfile).get_attributes()
2778        #print 'P', point_values
2779        #print 'Ref', ref_point_values       
2780        assert allclose(point_values, ref_point_values)       
2781
2782
2783
2784        # Invoke interpolation for centroids
2785        points = self.domain.get_centroid_coordinates()
2786        #print points
2787        sww2pts(self.domain.get_name(),
2788                quantity = 'elevation',
2789                data_points = points,
2790                NODATA_value = NODATA_value,
2791                verbose = self.verbose)
2792        ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids
2793
2794       
2795        point_values = Geospatial_data(ptsfile).get_attributes()
2796        #print 'P', point_values
2797        #print 'Ref', ref_point_values       
2798        assert allclose(point_values, ref_point_values)       
2799
2800
2801
2802        fid.close()
2803
2804        #Cleanup
2805        os.remove(sww.filename)
2806        os.remove(ptsfile)
2807
2808
2809
2810
2811    def test_ferret2sww1(self):
2812        """Test that georeferencing etc works when converting from
2813        ferret format (lat/lon) to sww format (UTM)
2814        """
2815        from Scientific.IO.NetCDF import NetCDFFile
2816        import os, sys
2817
2818        #The test file has
2819        # LON = 150.66667, 150.83334, 151, 151.16667
2820        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2821        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2822        #
2823        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2824        # Fourth value (index==3) is -6.50198 cm
2825
2826
2827
2828        #Read
2829        from anuga.coordinate_transforms.redfearn import redfearn
2830        #fid = NetCDFFile(self.test_MOST_file)
2831        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2832        first_value = fid.variables['HA'][:][0,0,0]
2833        fourth_value = fid.variables['HA'][:][0,0,3]
2834        fid.close()
2835
2836
2837        #Call conversion (with zero origin)
2838        #ferret2sww('small', verbose=False,
2839        #           origin = (56, 0, 0))
2840        ferret2sww(self.test_MOST_file, verbose=self.verbose,
2841                   origin = (56, 0, 0))
2842
2843        #Work out the UTM coordinates for first point
2844        zone, e, n = redfearn(-34.5, 150.66667)
2845        #print zone, e, n
2846
2847        #Read output file 'small.sww'
2848        #fid = NetCDFFile('small.sww')
2849        fid = NetCDFFile(self.test_MOST_file + '.sww')
2850
2851        x = fid.variables['x'][:]
2852        y = fid.variables['y'][:]
2853
2854        #Check that first coordinate is correctly represented
2855        assert allclose(x[0], e)
2856        assert allclose(y[0], n)
2857
2858        #Check first value
2859        stage = fid.variables['stage'][:]
2860        xmomentum = fid.variables['xmomentum'][:]
2861        ymomentum = fid.variables['ymomentum'][:]
2862
2863        #print ymomentum
2864
2865        assert allclose(stage[0,0], first_value/100)  #Meters
2866
2867        #Check fourth value
2868        assert allclose(stage[0,3], fourth_value/100)  #Meters
2869
2870        fid.close()
2871
2872        #Cleanup
2873        import os
2874        os.remove(self.test_MOST_file + '.sww')
2875
2876
2877
2878    def test_ferret2sww_zscale(self):
2879        """Test that zscale workse
2880        """
2881        from Scientific.IO.NetCDF import NetCDFFile
2882        import os, sys
2883
2884        #The test file has
2885        # LON = 150.66667, 150.83334, 151, 151.16667
2886        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2887        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2888        #
2889        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2890        # Fourth value (index==3) is -6.50198 cm
2891
2892
2893        #Read
2894        from anuga.coordinate_transforms.redfearn import redfearn
2895        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2896        first_value = fid.variables['HA'][:][0,0,0]
2897        fourth_value = fid.variables['HA'][:][0,0,3]
2898        fid.close()
2899
2900        #Call conversion (with no scaling)
2901        ferret2sww(self.test_MOST_file, verbose=self.verbose,
2902                   origin = (56, 0, 0))
2903
2904        #Work out the UTM coordinates for first point
2905        fid = NetCDFFile(self.test_MOST_file + '.sww')
2906
2907        #Check values
2908        stage_1 = fid.variables['stage'][:]
2909        xmomentum_1 = fid.variables['xmomentum'][:]
2910        ymomentum_1 = fid.variables['ymomentum'][:]
2911
2912        assert allclose(stage_1[0,0], first_value/100)  #Meters
2913        assert allclose(stage_1[0,3], fourth_value/100)  #Meters
2914
2915        fid.close()
2916
2917        #Call conversion (with scaling)
2918        ferret2sww(self.test_MOST_file,
2919                   zscale = 5,
2920                   verbose=self.verbose,
2921                   origin = (56, 0, 0))
2922
2923        #Work out the UTM coordinates for first point
2924        fid = NetCDFFile(self.test_MOST_file + '.sww')
2925
2926        #Check values
2927        stage_5 = fid.variables['stage'][:]
2928        xmomentum_5 = fid.variables['xmomentum'][:]
2929        ymomentum_5 = fid.variables['ymomentum'][:]
2930        elevation = fid.variables['elevation'][:]
2931
2932        assert allclose(stage_5[0,0], 5*first_value/100)  #Meters
2933        assert allclose(stage_5[0,3], 5*fourth_value/100)  #Meters
2934
2935        assert allclose(5*stage_1, stage_5)
2936
2937        # Momentum will also be changed due to new depth
2938
2939        depth_1 = stage_1-elevation
2940        depth_5 = stage_5-elevation
2941
2942
2943        for i in range(stage_1.shape[0]):
2944            for j in range(stage_1.shape[1]):           
2945                if depth_1[i,j] > epsilon:
2946
2947                    scale = depth_5[i,j]/depth_1[i,j]
2948                    ref_xmomentum = xmomentum_1[i,j] * scale
2949                    ref_ymomentum = ymomentum_1[i,j] * scale
2950                   
2951                    #print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
2952                   
2953                    assert allclose(xmomentum_5[i,j], ref_xmomentum)
2954                    assert allclose(ymomentum_5[i,j], ref_ymomentum)
2955                   
2956       
2957
2958        fid.close()
2959
2960
2961        #Cleanup
2962        import os
2963        os.remove(self.test_MOST_file + '.sww')
2964
2965
2966
2967    def test_ferret2sww_2(self):
2968        """Test that georeferencing etc works when converting from
2969        ferret format (lat/lon) to sww format (UTM)
2970        """
2971        from Scientific.IO.NetCDF import NetCDFFile
2972
2973        #The test file has
2974        # LON = 150.66667, 150.83334, 151, 151.16667
2975        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2976        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2977        #
2978        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2979        # Fourth value (index==3) is -6.50198 cm
2980
2981
2982        from anuga.coordinate_transforms.redfearn import redfearn
2983
2984        #fid = NetCDFFile('small_ha.nc')
2985        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2986
2987        #Pick a coordinate and a value
2988
2989        time_index = 1
2990        lat_index = 0
2991        lon_index = 2
2992
2993        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
2994        test_time = fid.variables['TIME'][:][time_index]
2995        test_lat = fid.variables['LAT'][:][lat_index]
2996        test_lon = fid.variables['LON'][:][lon_index]
2997
2998        linear_point_index = lat_index*4 + lon_index
2999        fid.close()
3000
3001        #Call conversion (with zero origin)
3002        ferret2sww(self.test_MOST_file, verbose=self.verbose,
3003                   origin = (56, 0, 0))
3004
3005
3006        #Work out the UTM coordinates for test point
3007        zone, e, n = redfearn(test_lat, test_lon)
3008
3009        #Read output file 'small.sww'
3010        fid = NetCDFFile(self.test_MOST_file + '.sww')
3011
3012        x = fid.variables['x'][:]
3013        y = fid.variables['y'][:]
3014
3015        #Check that test coordinate is correctly represented
3016        assert allclose(x[linear_point_index], e)
3017        assert allclose(y[linear_point_index], n)
3018
3019        #Check test value
3020        stage = fid.variables['stage'][:]
3021
3022        assert allclose(stage[time_index, linear_point_index], test_value/100)
3023
3024        fid.close()
3025
3026        #Cleanup
3027        import os
3028        os.remove(self.test_MOST_file + '.sww')
3029
3030
3031    def test_ferret2sww_lat_long(self):
3032        # Test that min lat long works
3033
3034        #The test file has
3035        # LON = 150.66667, 150.83334, 151, 151.16667
3036        # LAT = -34.5, -34.33333, -34.16667, -34 ;
3037       
3038        #Read
3039        from anuga.coordinate_transforms.redfearn import redfearn
3040        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3041        first_value = fid.variables['HA'][:][0,0,0]
3042        fourth_value = fid.variables['HA'][:][0,0,3]
3043        fid.close()
3044
3045
3046        #Call conversion (with zero origin)
3047        #ferret2sww('small', verbose=self.verbose,
3048        #           origin = (56, 0, 0))
3049        ferret2sww(self.test_MOST_file, verbose=self.verbose,
3050                   origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
3051
3052        #Work out the UTM coordinates for first point
3053        zone, e, n = redfearn(-34.5, 150.66667)
3054        #print zone, e, n
3055
3056        #Read output file 'small.sww'
3057        #fid = NetCDFFile('small.sww')
3058        fid = NetCDFFile(self.test_MOST_file + '.sww')
3059
3060        x = fid.variables['x'][:]
3061        y = fid.variables['y'][:]
3062        #Check that first coordinate is correctly represented
3063        assert 16 == len(x)
3064
3065        fid.close()
3066
3067        #Cleanup
3068        import os
3069        os.remove(self.test_MOST_file + '.sww')
3070
3071
3072    def test_ferret2sww_lat_longII(self):
3073        # Test that min lat long works
3074
3075        #The test file has
3076        # LON = 150.66667, 150.83334, 151, 151.16667
3077        # LAT = -34.5, -34.33333, -34.16667, -34 ;
3078       
3079        #Read
3080        from anuga.coordinate_transforms.redfearn import redfearn
3081        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3082        first_value = fid.variables['HA'][:][0,0,0]
3083        fourth_value = fid.variables['HA'][:][0,0,3]
3084        fid.close()
3085
3086
3087        #Call conversion (with zero origin)
3088        #ferret2sww('small', verbose=False,
3089        #           origin = (56, 0, 0))
3090        ferret2sww(self.test_MOST_file, verbose=False,
3091                   origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
3092
3093        #Work out the UTM coordinates for first point
3094        zone, e, n = redfearn(-34.5, 150.66667)
3095        #print zone, e, n
3096
3097        #Read output file 'small.sww'
3098        #fid = NetCDFFile('small.sww')
3099        fid = NetCDFFile(self.test_MOST_file + '.sww')
3100
3101        x = fid.variables['x'][:]
3102        y = fid.variables['y'][:]
3103        #Check that first coordinate is correctly represented
3104        assert 12 == len(x)
3105
3106        fid.close()
3107
3108        #Cleanup
3109        import os
3110        os.remove(self.test_MOST_file + '.sww')
3111       
3112    def test_ferret2sww3(self):
3113        """Elevation included
3114        """
3115        from Scientific.IO.NetCDF import NetCDFFile
3116
3117        #The test file has
3118        # LON = 150.66667, 150.83334, 151, 151.16667
3119        # LAT = -34.5, -34.33333, -34.16667, -34 ;
3120        # ELEVATION = [-1 -2 -3 -4
3121        #             -5 -6 -7 -8
3122        #              ...
3123        #              ...      -16]
3124        # where the top left corner is -1m,
3125        # and the ll corner is -13.0m
3126        #
3127        # First value (index=0) in small_ha.nc is 0.3400644 cm,
3128        # Fourth value (index==3) is -6.50198 cm
3129
3130        from anuga.coordinate_transforms.redfearn import redfearn
3131        import os
3132        fid1 = NetCDFFile('test_ha.nc','w')
3133        fid2 = NetCDFFile('test_ua.nc','w')
3134        fid3 = NetCDFFile('test_va.nc','w')
3135        fid4 = NetCDFFile('test_e.nc','w')
3136
3137        h1_list = [150.66667,150.83334,151.]
3138        h2_list = [-34.5,-34.33333]
3139
3140        long_name = 'LON'
3141        lat_name = 'LAT'
3142        time_name = 'TIME'
3143
3144        nx = 3
3145        ny = 2
3146
3147        for fid in [fid1,fid2,fid3]:
3148            fid.createDimension(long_name,nx)
3149            fid.createVariable(long_name,'d',(long_name,))
3150            fid.variables[long_name].point_spacing='uneven'
3151            fid.variables[long_name].units='degrees_east'
3152            fid.variables[long_name].assignValue(h1_list)
3153
3154            fid.createDimension(lat_name,ny)
3155            fid.createVariable(lat_name,'d',(lat_name,))
3156            fid.variables[lat_name].point_spacing='uneven'
3157            fid.variables[lat_name].units='degrees_north'
3158            fid.variables[lat_name].assignValue(h2_list)
3159
3160            fid.createDimension(time_name,2)
3161            fid.createVariable(time_name,'d',(time_name,))
3162            fid.variables[time_name].point_spacing='uneven'
3163            fid.variables[time_name].units='seconds'
3164            fid.variables[time_name].assignValue([0.,1.])
3165            if fid == fid3: break
3166
3167
3168        for fid in [fid4]:
3169            fid.createDimension(long_name,nx)
3170            fid.createVariable(long_name,'d',(long_name,))
3171            fid.variables[long_name].point_spacing='uneven'
3172            fid.variables[long_name].units='degrees_east'
3173            fid.variables[long_name].assignValue(h1_list)
3174
3175            fid.createDimension(lat_name,ny)
3176            fid.createVariable(lat_name,'d',(lat_name,))
3177            fid.variables[lat_name].point_spacing='uneven'
3178            fid.variables[lat_name].units='degrees_north'
3179            fid.variables[lat_name].assignValue(h2_list)
3180
3181        name = {}
3182        name[fid1]='HA'
3183        name[fid2]='UA'
3184        name[fid3]='VA'
3185        name[fid4]='ELEVATION'
3186
3187        units = {}
3188        units[fid1]='cm'
3189        units[fid2]='cm/s'
3190        units[fid3]='cm/s'
3191        units[fid4]='m'
3192
3193        values = {}
3194        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
3195        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
3196        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
3197        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
3198
3199        for fid in [fid1,fid2,fid3]:
3200          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
3201          fid.variables[name[fid]].point_spacing='uneven'
3202          fid.variables[name[fid]].units=units[fid]
3203          fid.variables[name[fid]].assignValue(values[fid])
3204          fid.variables[name[fid]].missing_value = -99999999.
3205          if fid == fid3: break
3206
3207        for fid in [fid4]:
3208            fid.createVariable(name[fid],'d',(lat_name,long_name))
3209            fid.variables[name[fid]].point_spacing='uneven'
3210            fid.variables[name[fid]].units=units[fid]
3211            fid.variables[name[fid]].assignValue(values[fid])
3212            fid.variables[name[fid]].missing_value = -99999999.
3213
3214
3215        fid1.sync(); fid1.close()
3216        fid2.sync(); fid2.close()
3217        fid3.sync(); fid3.close()
3218        fid4.sync(); fid4.close()
3219
3220        fid1 = NetCDFFile('test_ha.nc','r')
3221        fid2 = NetCDFFile('test_e.nc','r')
3222        fid3 = NetCDFFile('test_va.nc','r')
3223
3224
3225        first_amp = fid1.variables['HA'][:][0,0,0]
3226        third_amp = fid1.variables['HA'][:][0,0,2]
3227        first_elevation = fid2.variables['ELEVATION'][0,0]
3228        third_elevation= fid2.variables['ELEVATION'][:][0,2]
3229        first_speed = fid3.variables['VA'][0,0,0]
3230        third_speed = fid3.variables['VA'][:][0,0,2]
3231
3232        fid1.close()
3233        fid2.close()
3234        fid3.close()
3235
3236        #Call conversion (with zero origin)
3237        ferret2sww('test', verbose=self.verbose,
3238                   origin = (56, 0, 0), inverted_bathymetry=False)
3239
3240        os.remove('test_va.nc')
3241        os.remove('test_ua.nc')
3242        os.remove('test_ha.nc')
3243        os.remove('test_e.nc')
3244
3245        #Read output file 'test.sww'
3246        fid = NetCDFFile('test.sww')
3247
3248
3249        #Check first value
3250        elevation = fid.variables['elevation'][:]
3251        stage = fid.variables['stage'][:]
3252        xmomentum = fid.variables['xmomentum'][:]
3253        ymomentum = fid.variables['ymomentum'][:]
3254
3255        #print ymomentum
3256        first_height = first_amp/100 - first_elevation
3257        third_height = third_amp/100 - third_elevation
3258        first_momentum=first_speed*first_height/100
3259        third_momentum=third_speed*third_height/100
3260
3261        assert allclose(ymomentum[0][0],first_momentum)  #Meters
3262        assert allclose(ymomentum[0][2],third_momentum)  #Meters
3263
3264        fid.close()
3265
3266        #Cleanup
3267        os.remove('test.sww')
3268
3269
3270
3271    def test_ferret2sww4(self):
3272        """Like previous but with augmented variable names as
3273        in files produced by ferret as opposed to MOST
3274        """
3275        from Scientific.IO.NetCDF import NetCDFFile
3276
3277        #The test file has
3278        # LON = 150.66667, 150.83334, 151, 151.16667
3279        # LAT = -34.5, -34.33333, -34.16667, -34 ;
3280        # ELEVATION = [-1 -2 -3 -4
3281        #             -5 -6 -7 -8
3282        #              ...
3283        #              ...      -16]
3284        # where the top left corner is -1m,
3285        # and the ll corner is -13.0m
3286        #
3287        # First value (index=0) in small_ha.nc is 0.3400644 cm,
3288        # Fourth value (index==3) is -6.50198 cm
3289
3290        from anuga.coordinate_transforms.redfearn import redfearn
3291        import os
3292        fid1 = NetCDFFile('test_ha.nc','w')
3293        fid2 = NetCDFFile('test_ua.nc','w')
3294        fid3 = NetCDFFile('test_va.nc','w')
3295        fid4 = NetCDFFile('test_e.nc','w')
3296
3297        h1_list = [150.66667,150.83334,151.]
3298        h2_list = [-34.5,-34.33333]
3299
3300#        long_name = 'LON961_1261'
3301#        lat_name = 'LAT481_841'
3302#        time_name = 'TIME1'
3303
3304        long_name = 'LON'
3305        lat_name = 'LAT'
3306        time_name = 'TIME'
3307
3308        nx = 3
3309        ny = 2
3310
3311        for fid in [fid1,fid2,fid3]:
3312            fid.createDimension(long_name,nx)
3313            fid.createVariable(long_name,'d',(long_name,))
3314            fid.variables[long_name].point_spacing='uneven'
3315            fid.variables[long_name].units='degrees_east'
3316            fid.variables[long_name].assignValue(h1_list)
3317
3318            fid.createDimension(lat_name,ny)
3319            fid.createVariable(lat_name,'d',(lat_name,))
3320            fid.variables[lat_name].point_spacing='uneven'
3321            fid.variables[lat_name].units='degrees_north'
3322            fid.variables[lat_name].assignValue(h2_list)
3323
3324            fid.createDimension(time_name,2)
3325            fid.createVariable(time_name,'d',(time_name,))
3326            fid.variables[time_name].point_spacing='uneven'
3327            fid.variables[time_name].units='seconds'
3328            fid.variables[time_name].assignValue([0.,1.])
3329            if fid == fid3: break
3330
3331
3332        for fid in [fid4]:
3333            fid.createDimension(long_name,nx)
3334            fid.createVariable(long_name,'d',(long_name,))
3335            fid.variables[long_name].point_spacing='uneven'
3336            fid.variables[long_name].units='degrees_east'
3337            fid.variables[long_name].assignValue(h1_list)
3338
3339            fid.createDimension(lat_name,ny)
3340            fid.createVariable(lat_name,'d',(lat_name,))
3341            fid.variables[lat_name].point_spacing='uneven'
3342            fid.variables[lat_name].units='degrees_north'
3343            fid.variables[lat_name].assignValue(h2_list)
3344
3345        name = {}
3346        name[fid1]='HA'
3347        name[fid2]='UA'
3348        name[fid3]='VA'
3349        name[fid4]='ELEVATION'
3350
3351        units = {}
3352        units[fid1]='cm'
3353        units[fid2]='cm/s'
3354        units[fid3]='cm/s'
3355        units[fid4]='m'
3356
3357        values = {}
3358        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
3359        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
3360        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
3361        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
3362
3363        for fid in [fid1,fid2,fid3]:
3364          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
3365          fid.variables[name[fid]].point_spacing='uneven'
3366          fid.variables[name[fid]].units=units[fid]
3367          fid.variables[name[fid]].assignValue(values[fid])
3368          fid.variables[name[fid]].missing_value = -99999999.
3369          if fid == fid3: break
3370
3371        for fid in [fid4]:
3372            fid.createVariable(name[fid],'d',(lat_name,long_name))
3373            fid.variables[name[fid]].point_spacing='uneven'
3374            fid.variables[name[fid]].units=units[fid]
3375            fid.variables[name[fid]].assignValue(values[fid])
3376            fid.variables[name[fid]].missing_value = -99999999.
3377
3378
3379        fid1.sync(); fid1.close()
3380        fid2.sync(); fid2.close()
3381        fid3.sync(); fid3.close()
3382        fid4.sync(); fid4.close()
3383
3384        fid1 = NetCDFFile('test_ha.nc','r')
3385        fid2 = NetCDFFile('test_e.nc','r')
3386        fid3 = NetCDFFile('test_va.nc','r')
3387
3388
3389        first_amp = fid1.variables['HA'][:][0,0,0]
3390        third_amp = fid1.variables['HA'][:][0,0,2]
3391        first_elevation = fid2.variables['ELEVATION'][0,0]
3392        third_elevation= fid2.variables['ELEVATION'][:][0,2]
3393        first_speed = fid3.variables['VA'][0,0,0]
3394        third_speed = fid3.variables['VA'][:][0,0,2]
3395
3396        fid1.close()
3397        fid2.close()
3398        fid3.close()
3399
3400        #Call conversion (with zero origin)
3401        ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
3402                   , inverted_bathymetry=False)
3403
3404        os.remove('test_va.nc')
3405        os.remove('test_ua.nc')
3406        os.remove('test_ha.nc')
3407        os.remove('test_e.nc')
3408
3409        #Read output file 'test.sww'
3410        fid = NetCDFFile('test.sww')
3411
3412
3413        #Check first value
3414        elevation = fid.variables['elevation'][:]
3415        stage = fid.variables['stage'][:]
3416        xmomentum = fid.variables['xmomentum'][:]
3417        ymomentum = fid.variables['ymomentum'][:]
3418
3419        #print ymomentum
3420        first_height = first_amp/100 - first_elevation
3421        third_height = third_amp/100 - third_elevation
3422        first_momentum=first_speed*first_height/100
3423        third_momentum=third_speed*third_height/100
3424
3425        assert allclose(ymomentum[0][0],first_momentum)  #Meters
3426        assert allclose(ymomentum[0][2],third_momentum)  #Meters
3427
3428        fid.close()
3429
3430        #Cleanup
3431        os.remove('test.sww')
3432
3433
3434
3435
3436    def test_ferret2sww_nz_origin(self):
3437        from Scientific.IO.NetCDF import NetCDFFile
3438        from anuga.coordinate_transforms.redfearn import redfearn
3439
3440        #Call conversion (with nonzero origin)
3441        ferret2sww(self.test_MOST_file, verbose=self.verbose,
3442                   origin = (56, 100000, 200000))
3443
3444
3445        #Work out the UTM coordinates for first point
3446        zone, e, n = redfearn(-34.5, 150.66667)
3447
3448        #Read output file 'small.sww'
3449        #fid = NetCDFFile('small.sww', 'r')
3450        fid = NetCDFFile(self.test_MOST_file + '.sww')
3451
3452        x = fid.variables['x'][:]
3453        y = fid.variables['y'][:]
3454
3455        #Check that first coordinate is correctly represented
3456        assert allclose(x[0], e-100000)
3457        assert allclose(y[0], n-200000)
3458
3459        fid.close()
3460
3461        #Cleanup
3462        os.remove(self.test_MOST_file + '.sww')
3463
3464
3465    def test_ferret2sww_lat_longII(self):
3466        # Test that min lat long works
3467
3468        #The test file has
3469        # LON = 150.66667, 150.83334, 151, 151.16667
3470        # LAT = -34.5, -34.33333, -34.16667, -34 ;
3471       
3472        #Read
3473        from anuga.coordinate_transforms.redfearn import redfearn
3474        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3475        first_value = fid.variables['HA'][:][0,0,0]
3476        fourth_value = fid.variables['HA'][:][0,0,3]
3477        fid.close()
3478
3479
3480        #Call conversion (with zero origin)
3481        #ferret2sww('small', verbose=self.verbose,
3482        #           origin = (56, 0, 0))
3483        try:
3484            ferret2sww(self.test_MOST_file, verbose=self.verbose,
3485                   origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
3486        except AssertionError:
3487            pass
3488        else:
3489            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
3490
3491    def test_sww_extent(self):
3492        """Not a test, rather a look at the sww format
3493        """
3494
3495        import time, os
3496        from Numeric import array, zeros, allclose, Float, concatenate
3497        from Scientific.IO.NetCDF import NetCDFFile
3498
3499        self.domain.set_name('datatest' + str(id(self)))
3500        self.domain.format = 'sww'
3501        self.domain.smooth = True
3502        self.domain.reduction = mean
3503        self.domain.set_datadir('.')
3504        #self.domain.tight_slope_limiters = 1       
3505
3506
3507        sww = get_dataobject(self.domain)
3508        sww.store_connectivity()
3509        sww.store_timestep('stage')
3510        self.domain.time = 2.
3511
3512        #Modify stage at second timestep
3513        stage = self.domain.quantities['stage'].vertex_values
3514        self.domain.set_quantity('stage', stage/2)
3515
3516        sww.store_timestep('stage')
3517
3518        file_and_extension_name = self.domain.get_name() + ".sww"
3519        #print "file_and_extension_name",file_and_extension_name
3520        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
3521               extent_sww(file_and_extension_name )
3522
3523        assert allclose(xmin, 0.0)
3524        assert allclose(xmax, 1.0)
3525        assert allclose(ymin, 0.0)
3526        assert allclose(ymax, 1.0)
3527
3528        # FIXME (Ole): Revisit these numbers
3529        #assert allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
3530        #assert allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
3531
3532
3533        #Cleanup
3534        os.remove(sww.filename)
3535
3536
3537
3538    def test_sww2domain1(self):
3539        ################################################
3540        #Create a test domain, and evolve and save it.
3541        ################################################
3542        from mesh_factory import rectangular
3543        from Numeric import array
3544
3545        #Create basic mesh
3546
3547        yiel=0.01
3548        points, vertices, boundary = rectangular(10,10)
3549
3550        #Create shallow water domain
3551        domain = Domain(points, vertices, boundary)
3552        domain.geo_reference = Geo_reference(56,11,11)
3553        domain.smooth = False
3554        domain.store = True
3555        domain.set_name('bedslope')
3556        domain.default_order=2
3557        #Bed-slope and friction
3558        domain.set_quantity('elevation', lambda x,y: -x/3)
3559        domain.set_quantity('friction', 0.1)
3560        # Boundary conditions
3561        from math import sin, pi
3562        Br = Reflective_boundary(domain)
3563        Bt = Transmissive_boundary(domain)
3564        Bd = Dirichlet_boundary([0.2,0.,0.])
3565        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
3566
3567        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3568        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3569
3570        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
3571        #Initial condition
3572        h = 0.05
3573        elevation = domain.quantities['elevation'].vertex_values
3574        domain.set_quantity('stage', elevation + h)
3575
3576        domain.check_integrity()
3577        #Evolution
3578        #domain.tight_slope_limiters = 1
3579        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3580            #domain.write_time()
3581            pass
3582
3583
3584        ##########################################
3585        #Import the example's file as a new domain
3586        ##########################################
3587        from data_manager import sww2domain
3588        from Numeric import allclose
3589        import os
3590
3591        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3592        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
3593        #points, vertices, boundary = rectangular(15,15)
3594        #domain2.boundary = boundary
3595        ###################
3596        ##NOW TEST IT!!!
3597        ###################
3598
3599        os.remove(filename)
3600
3601        bits = ['vertex_coordinates']
3602        for quantity in ['elevation']+domain.quantities_to_be_stored:
3603            bits.append('get_quantity("%s").get_integral()' %quantity)
3604            bits.append('get_quantity("%s").get_values()' %quantity)
3605
3606        for bit in bits:
3607            #print 'testing that domain.'+bit+' has been restored'
3608            #print bit
3609            #print 'done'
3610            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3611
3612        ######################################
3613        #Now evolve them both, just to be sure
3614        ######################################x
3615        domain.time = 0.
3616        from time import sleep
3617
3618        final = .1
3619        domain.set_quantity('friction', 0.1)
3620        domain.store = False
3621        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3622
3623
3624        for t in domain.evolve(yieldstep = yiel, finaltime = final):
3625            #domain.write_time()
3626            pass
3627
3628        final = final - (domain2.starttime-domain.starttime)
3629        #BUT since domain1 gets time hacked back to 0:
3630        final = final + (domain2.starttime-domain.starttime)
3631
3632        domain2.smooth = False
3633        domain2.store = False
3634        domain2.default_order=2
3635        domain2.set_quantity('friction', 0.1)
3636        #Bed-slope and friction
3637        # Boundary conditions
3638        Bd2=Dirichlet_boundary([0.2,0.,0.])
3639        domain2.boundary = domain.boundary
3640        #print 'domain2.boundary'
3641        #print domain2.boundary
3642        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3643        #domain2.set_boundary({'exterior': Bd})
3644
3645        domain2.check_integrity()
3646
3647        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3648            #domain2.write_time()
3649            pass
3650
3651        ###################
3652        ##NOW TEST IT!!!
3653        ##################
3654
3655        bits = ['vertex_coordinates']
3656
3657        for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
3658            bits.append('get_quantity("%s").get_integral()' %quantity)
3659            bits.append('get_quantity("%s").get_values()' %quantity)
3660
3661        #print bits
3662        for bit in bits:
3663            #print bit
3664            #print eval('domain.'+bit)
3665            #print eval('domain2.'+bit)
3666           
3667            #print eval('domain.'+bit+'-domain2.'+bit)
3668            msg = 'Values in the two domains are different for ' + bit
3669            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),
3670                            rtol=1.e-5, atol=3.e-8), msg
3671
3672
3673    def test_sww2domain2(self):
3674        ##################################################################
3675        #Same as previous test, but this checks how NaNs are handled.
3676        ##################################################################
3677
3678
3679        from mesh_factory import rectangular
3680        from Numeric import array
3681
3682        #Create basic mesh
3683        points, vertices, boundary = rectangular(2,2)
3684
3685        #Create shallow water domain
3686        domain = Domain(points, vertices, boundary)
3687        domain.smooth = False
3688        domain.store = True
3689        domain.set_name('test_file')
3690        domain.set_datadir('.')
3691        domain.default_order=2
3692        domain.quantities_to_be_stored=['stage']
3693
3694        domain.set_quantity('elevation', lambda x,y: -x/3)
3695        domain.set_quantity('friction', 0.1)
3696
3697        from math import sin, pi
3698        Br = Reflective_boundary(domain)
3699        Bt = Transmissive_boundary(domain)
3700        Bd = Dirichlet_boundary([0.2,0.,0.])
3701        Bw = Time_boundary(domain=domain,
3702                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
3703
3704        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3705
3706        h = 0.05
3707        elevation = domain.quantities['elevation'].vertex_values
3708        domain.set_quantity('stage', elevation + h)
3709
3710        domain.check_integrity()
3711
3712        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
3713            pass
3714            #domain.write_time()
3715
3716
3717
3718        ##################################
3719        #Import the file as a new domain
3720        ##################################
3721        from data_manager import sww2domain
3722        from Numeric import allclose
3723        import os
3724
3725        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3726
3727        #Fail because NaNs are present
3728        try:
3729            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=self.verbose)
3730        except:
3731            #Now import it, filling NaNs to be 0
3732            filler = 0
3733            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=self.verbose)
3734
3735        #Clean up
3736        os.remove(filename)
3737
3738
3739        bits = [ 'geo_reference.get_xllcorner()',
3740                'geo_reference.get_yllcorner()',
3741                'vertex_coordinates']
3742
3743        for quantity in ['elevation']+domain.quantities_to_be_stored:
3744            bits.append('get_quantity("%s").get_integral()' %quantity)
3745            bits.append('get_quantity("%s").get_values()' %quantity)
3746
3747        for bit in bits:
3748        #    print 'testing that domain.'+bit+' has been restored'
3749            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3750
3751        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
3752        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
3753        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
3754        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
3755
3756
3757
3758    #def test_weed(self):
3759        from data_manager import weed
3760
3761        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
3762        volumes1 = [[0,1,2],[3,4,5]]
3763        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
3764        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
3765
3766        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
3767
3768        assert len(points2)==len(coordinates2)
3769        for i in range(len(coordinates2)):
3770            coordinate = tuple(coordinates2[i])
3771            assert points2.has_key(coordinate)
3772            points2[coordinate]=i
3773
3774        for triangle in volumes1:
3775            for coordinate in triangle:
3776                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
3777                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
3778
3779
3780    #FIXME This fails - smooth makes the comparism too hard for allclose
3781    def ztest_sww2domain3(self):
3782        ################################################
3783        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
3784        ################################################
3785        from mesh_factory import rectangular
3786        from Numeric import array
3787        #Create basic mesh
3788
3789        yiel=0.01
3790        points, vertices, boundary = rectangular(10,10)
3791
3792        #Create shallow water domain
3793        domain = Domain(points, vertices, boundary)
3794        domain.geo_reference = Geo_reference(56,11,11)
3795        domain.smooth = True
3796        domain.store = True
3797        domain.set_name('bedslope')
3798        domain.default_order=2
3799        #Bed-slope and friction
3800        domain.set_quantity('elevation', lambda x,y: -x/3)
3801        domain.set_quantity('friction', 0.1)
3802        # Boundary conditions
3803        from math import sin, pi
3804        Br = Reflective_boundary(domain)
3805        Bt = Transmissive_boundary(domain)
3806        Bd = Dirichlet_boundary([0.2,0.,0.])
3807        Bw = Time_boundary(domain=domain,
3808                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
3809
3810        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3811
3812        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
3813        #Initial condition
3814        h = 0.05
3815        elevation = domain.quantities['elevation'].vertex_values
3816        domain.set_quantity('stage', elevation + h)
3817
3818
3819        domain.check_integrity()
3820        #Evolution
3821        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3822        #    domain.write_time()
3823            pass
3824
3825
3826        ##########################################
3827        #Import the example's file as a new domain
3828        ##########################################
3829        from data_manager import sww2domain
3830        from Numeric import allclose
3831        import os
3832
3833        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3834        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
3835        #points, vertices, boundary = rectangular(15,15)
3836        #domain2.boundary = boundary
3837        ###################
3838        ##NOW TEST IT!!!
3839        ###################
3840
3841        os.remove(domain.get_name() + '.sww')
3842
3843        #FIXME smooth domain so that they can be compared
3844
3845
3846        bits = []#'vertex_coordinates']
3847        for quantity in ['elevation']+domain.quantities_to_be_stored:
3848            bits.append('quantities["%s"].get_integral()'%quantity)
3849
3850
3851        for bit in bits:
3852            #print 'testing that domain.'+bit+' has been restored'
3853            #print bit
3854            #print 'done'
3855            #print ('domain.'+bit), eval('domain.'+bit)
3856            #print ('domain2.'+bit), eval('domain2.'+bit)
3857            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
3858            pass
3859
3860        ######################################
3861        #Now evolve them both, just to be sure
3862        ######################################x
3863        domain.time = 0.
3864        from time import sleep
3865
3866        final = .5
3867        domain.set_quantity('friction', 0.1)
3868        domain.store = False
3869        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
3870
3871        for t in domain.evolve(yieldstep = yiel, finaltime = final):
3872            #domain.write_time()
3873            pass
3874
3875        domain2.smooth = True
3876        domain2.store = False
3877        domain2.default_order=2
3878        domain2.set_quantity('friction', 0.1)
3879        #Bed-slope and friction
3880        # Boundary conditions
3881        Bd2=Dirichlet_boundary([0.2,0.,0.])
3882        Br2 = Reflective_boundary(domain2)
3883        domain2.boundary = domain.boundary
3884        #print 'domain2.boundary'
3885        #print domain2.boundary
3886        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
3887        #domain2.boundary = domain.boundary
3888        #domain2.set_boundary({'exterior': Bd})
3889
3890        domain2.check_integrity()
3891
3892        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3893            #domain2.write_time()
3894            pass
3895
3896        ###################
3897        ##NOW TEST IT!!!
3898        ##################
3899
3900        print '><><><><>>'
3901        bits = [ 'vertex_coordinates']
3902
3903        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
3904            #bits.append('quantities["%s"].get_integral()'%quantity)
3905            bits.append('get_quantity("%s").get_values()' %quantity)
3906
3907        for bit in bits:
3908            print bit
3909            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3910
3911
3912    def test_decimate_dem(self):
3913        """Test decimation of dem file
3914        """
3915
3916        import os
3917        from Numeric import ones, allclose, Float, arange
3918        from Scientific.IO.NetCDF import NetCDFFile
3919
3920        #Write test dem file
3921        root = 'decdemtest'
3922
3923        filename = root + '.dem'
3924        fid = NetCDFFile(filename, 'w')
3925
3926        fid.institution = 'Geoscience Australia'
3927        fid.description = 'NetCDF DEM format for compact and portable ' +\
3928                          'storage of spatial point data'
3929
3930        nrows = 15
3931        ncols = 18
3932
3933        fid.ncols = ncols
3934        fid.nrows = nrows
3935        fid.xllcorner = 2000.5
3936        fid.yllcorner = 3000.5
3937        fid.cellsize = 25
3938        fid.NODATA_value = -9999
3939
3940        fid.zone = 56
3941        fid.false_easting = 0.0
3942        fid.false_northing = 0.0
3943        fid.projection = 'UTM'
3944        fid.datum = 'WGS84'
3945        fid.units = 'METERS'
3946
3947        fid.createDimension('number_of_points', nrows*ncols)
3948
3949        fid.createVariable('elevation', Float, ('number_of_points',))
3950
3951        elevation = fid.variables['elevation']
3952
3953        elevation[:] = (arange(nrows*ncols))
3954
3955        fid.close()
3956
3957        #generate the elevation values expected in the decimated file
3958        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3959                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3960                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3961                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3962                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3963                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3964                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3965                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3966                         (144+145+146+162+163+164+180+181+182) / 9.0,
3967                         (148+149+150+166+167+168+184+185+186) / 9.0,
3968                         (152+153+154+170+171+172+188+189+190) / 9.0,
3969                         (156+157+158+174+175+176+192+193+194) / 9.0,
3970                         (216+217+218+234+235+236+252+253+254) / 9.0,
3971                         (220+221+222+238+239+240+256+257+258) / 9.0,
3972                         (224+225+226+242+243+244+260+261+262) / 9.0,
3973                         (228+229+230+246+247+248+264+265+266) / 9.0]
3974
3975        #generate a stencil for computing the decimated values
3976        stencil = ones((3,3), Float) / 9.0
3977
3978        decimate_dem(root, stencil=stencil, cellsize_new=100)
3979
3980        #Open decimated NetCDF file
3981        fid = NetCDFFile(root + '_100.dem', 'r')
3982
3983        # Get decimated elevation
3984        elevation = fid.variables['elevation']
3985
3986        #Check values
3987        assert allclose(elevation, ref_elevation)
3988
3989        #Cleanup
3990        fid.close()
3991
3992        os.remove(root + '.dem')
3993        os.remove(root + '_100.dem')
3994
3995    def test_decimate_dem_NODATA(self):
3996        """Test decimation of dem file that includes NODATA values
3997        """
3998
3999        import os
4000        from Numeric import ones, allclose, Float, arange, reshape
4001        from Scientific.IO.NetCDF import NetCDFFile
4002
4003        #Write test dem file
4004        root = 'decdemtest'
4005
4006        filename = root + '.dem'
4007        fid = NetCDFFile(filename, 'w')
4008
4009        fid.institution = 'Geoscience Australia'
4010        fid.description = 'NetCDF DEM format for compact and portable ' +\
4011                          'storage of spatial point data'
4012
4013        nrows = 15
4014        ncols = 18
4015        NODATA_value = -9999
4016
4017        fid.ncols = ncols
4018        fid.nrows = nrows
4019        fid.xllcorner = 2000.5
4020        fid.yllcorner = 3000.5
4021        fid.cellsize = 25
4022        fid.NODATA_value = NODATA_value
4023
4024        fid.zone = 56
4025        fid.false_easting = 0.0
4026        fid.false_northing = 0.0
4027        fid.projection = 'UTM'
4028        fid.datum = 'WGS84'
4029        fid.units = 'METERS'
4030
4031        fid.createDimension('number_of_points', nrows*ncols)
4032
4033        fid.createVariable('elevation', Float, ('number_of_points',))
4034
4035        elevation = fid.variables['elevation']
4036
4037        #generate initial elevation values
4038        elevation_tmp = (arange(nrows*ncols))
4039        #add some NODATA values
4040        elevation_tmp[0]   = NODATA_value
4041        elevation_tmp[95]  = NODATA_value
4042        elevation_tmp[188] = NODATA_value
4043        elevation_tmp[189] = NODATA_value
4044        elevation_tmp[190] = NODATA_value
4045        elevation_tmp[209] = NODATA_value
4046        elevation_tmp[252] = NODATA_value
4047
4048        elevation[:] = elevation_tmp
4049
4050        fid.close()
4051
4052        #generate the elevation values expected in the decimated file
4053        ref_elevation = [NODATA_value,
4054                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
4055                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
4056                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
4057                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
4058                         NODATA_value,
4059                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
4060                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
4061                         (144+145+146+162+163+164+180+181+182) / 9.0,
4062                         (148+149+150+166+167+168+184+185+186) / 9.0,
4063                         NODATA_value,
4064                         (156+157+158+174+175+176+192+193+194) / 9.0,
4065                         NODATA_value,
4066                         (220+221+222+238+239+240+256+257+258) / 9.0,
4067                         (224+225+226+242+243+244+260+261+262) / 9.0,
4068                         (228+229+230+246+247+248+264+265+266) / 9.0]
4069
4070        #generate a stencil for computing the decimated values
4071        stencil = ones((3,3), Float) / 9.0
4072
4073        decimate_dem(root, stencil=stencil, cellsize_new=100)
4074
4075        #Open decimated NetCDF file
4076        fid = NetCDFFile(root + '_100.dem', 'r')
4077
4078        # Get decimated elevation
4079        elevation = fid.variables['elevation']
4080
4081        #Check values
4082        assert allclose(elevation, ref_elevation)
4083
4084        #Cleanup
4085        fid.close()
4086
4087        os.remove(root + '.dem')
4088        os.remove(root + '_100.dem')
4089
4090    def xxxtestz_sww2ers_real(self):
4091        """Test that sww information can be converted correctly to asc/prj
4092        format readable by e.g. ArcView
4093        """
4094
4095        import time, os
4096        from Numeric import array, zeros, allclose, Float, concatenate
4097        from Scientific.IO.NetCDF import NetCDFFile
4098
4099        # the memory optimised least squares
4100        #  cellsize = 20,   # this one seems to hang
4101        #  cellsize = 200000, # Ran 1 test in 269.703s
4102                                #Ran 1 test in 267.344s
4103        #  cellsize = 20000,  # Ran 1 test in 460.922s
4104        #  cellsize = 2000   #Ran 1 test in 5340.250s
4105        #  cellsize = 200   #this one seems to hang, building matirx A
4106
4107        # not optimised
4108        # seems to hang
4109        #  cellsize = 2000   # Ran 1 test in 5334.563s
4110        #Export to ascii/prj files
4111        sww2dem('karratha_100m',
4112                quantity = 'depth',
4113                cellsize = 200000,
4114                verbose = True)
4115
4116    def test_read_asc(self):
4117        """Test conversion from dem in ascii format to native NetCDF format
4118        """
4119
4120        import time, os
4121        from Numeric import array, zeros, allclose, Float, concatenate
4122        from Scientific.IO.NetCDF import NetCDFFile
4123
4124        from data_manager import _read_asc
4125        #Write test asc file
4126        filename = tempfile.mktemp(".000")
4127        fid = open(filename, 'w')
4128        fid.write("""ncols         7
4129nrows         4
4130xllcorner     2000.5
4131yllcorner     3000.5
4132cellsize      25
4133NODATA_value  -9999
4134    97.921    99.285   125.588   180.830   258.645   342.872   415.836
4135   473.157   514.391   553.893   607.120   678.125   777.283   883.038
4136   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
4137   502.645   516.230   504.739   450.604   388.500   338.097   514.980
4138""")
4139        fid.close()
4140        bath_metadata, grid = _read_asc(filename, verbose=self.verbose)
4141        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
4142        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
4143        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
4144        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
4145        self.failUnless(grid[0][0]  == 97.921,  'Failed')
4146        self.failUnless(grid[3][6]  == 514.980,  'Failed')
4147
4148        os.remove(filename)
4149
4150    def test_asc_csiro2sww(self):
4151        import tempfile
4152
4153        bath_dir = tempfile.mkdtemp()
4154        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4155        #bath_dir = 'bath_data_manager_test'
4156        #print "os.getcwd( )",os.getcwd( )
4157        elevation_dir =  tempfile.mkdtemp()
4158        #elevation_dir = 'elev_expanded'
4159        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4160        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4161
4162        fid = open(bath_dir_filename, 'w')
4163        fid.write(""" ncols             3
4164 nrows             2
4165 xllcorner    148.00000
4166 yllcorner    -38.00000
4167 cellsize       0.25
4168 nodata_value   -9999.0
4169    9000.000    -1000.000    3000.0
4170   -1000.000    9000.000  -1000.000
4171""")
4172        fid.close()
4173
4174        fid = open(elevation_dir_filename1, 'w')
4175        fid.write(""" ncols             3
4176 nrows             2
4177 xllcorner    148.00000
4178 yllcorner    -38.00000
4179 cellsize       0.25
4180 nodata_value   -9999.0
4181    9000.000    0.000    3000.0
4182     0.000     9000.000     0.000
4183""")
4184        fid.close()
4185
4186        fid = open(elevation_dir_filename2, 'w')
4187        fid.write(""" ncols             3
4188 nrows             2
4189 xllcorner    148.00000
4190 yllcorner    -38.00000
4191 cellsize       0.25
4192 nodata_value   -9999.0
4193    9000.000    4000.000    4000.0
4194    4000.000    9000.000    4000.000
4195""")
4196        fid.close()
4197
4198        ucur_dir =  tempfile.mkdtemp()
4199        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4200        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4201
4202        fid = open(ucur_dir_filename1, 'w')
4203        fid.write(""" ncols             3
4204 nrows             2
4205 xllcorner    148.00000
4206 yllcorner    -38.00000
4207 cellsize       0.25
4208 nodata_value   -9999.0
4209    90.000    60.000    30.0
4210    10.000    10.000    10.000
4211""")
4212        fid.close()
4213        fid = open(ucur_dir_filename2, 'w')
4214        fid.write(""" ncols             3
4215 nrows             2
4216 xllcorner    148.00000
4217 yllcorner    -38.00000
4218 cellsize       0.25
4219 nodata_value   -9999.0
4220    90.000    60.000    30.0
4221    10.000    10.000    10.000
4222""")
4223        fid.close()
4224
4225        vcur_dir =  tempfile.mkdtemp()
4226        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4227        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4228
4229        fid = open(vcur_dir_filename1, 'w')
4230        fid.write(""" ncols             3
4231 nrows             2
4232 xllcorner    148.00000
4233 yllcorner    -38.00000
4234 cellsize       0.25
4235 nodata_value   -9999.0
4236    90.000    60.000    30.0
4237    10.000    10.000    10.000
4238""")
4239        fid.close()
4240        fid = open(vcur_dir_filename2, 'w')
4241        fid.write(""" ncols             3
4242 nrows             2
4243 xllcorner    148.00000
4244 yllcorner    -38.00000
4245 cellsize       0.25
4246 nodata_value   -9999.0
4247    90.000    60.000    30.0
4248    10.000    10.000    10.000
4249""")
4250        fid.close()
4251
4252        sww_file = 'a_test.sww'
4253        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
4254                      verbose=self.verbose)
4255
4256        # check the sww file
4257
4258        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
4259        x = fid.variables['x'][:]
4260        y = fid.variables['y'][:]
4261        z = fid.variables['z'][:]
4262        stage = fid.variables['stage'][:]
4263        xmomentum = fid.variables['xmomentum'][:]
4264        geo_ref = Geo_reference(NetCDFObject=fid)
4265        #print "geo_ref",geo_ref
4266        x_ref = geo_ref.get_xllcorner()
4267        y_ref = geo_ref.get_yllcorner()
4268        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
4269        assert allclose(x_ref, 587798.418) # (-38, 148)
4270        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
4271
4272        #Zone:   55
4273        #Easting:  588095.674  Northing: 5821451.722
4274        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
4275        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
4276
4277        #Zone:   55
4278        #Easting:  632145.632  Northing: 5820863.269
4279        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
4280        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
4281
4282        #Zone:   55
4283        #Easting:  609748.788  Northing: 5793447.860
4284        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
4285        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
4286
4287        assert allclose(z[0],9000.0 )
4288        assert allclose(stage[0][1],0.0 )
4289
4290        #(4000+1000)*60
4291        assert allclose(xmomentum[1][1],300000.0 )
4292
4293
4294        fid.close()
4295
4296        #tidy up
4297        os.remove(bath_dir_filename)
4298        os.rmdir(bath_dir)
4299
4300        os.remove(elevation_dir_filename1)
4301        os.remove(elevation_dir_filename2)
4302        os.rmdir(elevation_dir)
4303
4304        os.remove(ucur_dir_filename1)
4305        os.remove(ucur_dir_filename2)
4306        os.rmdir(ucur_dir)
4307
4308        os.remove(vcur_dir_filename1)
4309        os.remove(vcur_dir_filename2)
4310        os.rmdir(vcur_dir)
4311
4312
4313        # remove sww file
4314        os.remove(sww_file)
4315
4316    def test_asc_csiro2sww2(self):
4317        import tempfile
4318
4319        bath_dir = tempfile.mkdtemp()
4320        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4321        #bath_dir = 'bath_data_manager_test'
4322        #print "os.getcwd( )",os.getcwd( )
4323        elevation_dir =  tempfile.mkdtemp()
4324        #elevation_dir = 'elev_expanded'
4325        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4326        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4327
4328        fid = open(bath_dir_filename, 'w')
4329        fid.write(""" ncols             3
4330 nrows             2
4331 xllcorner    148.00000
4332 yllcorner    -38.00000
4333 cellsize       0.25
4334 nodata_value   -9999.0
4335    9000.000    -1000.000    3000.0
4336   -1000.000    9000.000  -1000.000
4337""")
4338        fid.close()
4339
4340        fid = open(elevation_dir_filename1, 'w')
4341        fid.write(""" ncols             3
4342 nrows             2
4343 xllcorner    148.00000
4344 yllcorner    -38.00000
4345 cellsize       0.25
4346 nodata_value   -9999.0
4347    9000.000    0.000    3000.0
4348     0.000     -9999.000     -9999.000
4349""")
4350        fid.close()
4351
4352        fid = open(elevation_dir_filename2, 'w')
4353        fid.write(""" ncols             3
4354 nrows             2
4355 xllcorner    148.00000
4356 yllcorner    -38.00000
4357 cellsize       0.25
4358 nodata_value   -9999.0
4359    9000.000    4000.000    4000.0
4360    4000.000    9000.000    4000.000
4361""")
4362        fid.close()
4363
4364        ucur_dir =  tempfile.mkdtemp()
4365        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4366        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4367
4368        fid = open(ucur_dir_filename1, 'w')
4369        fid.write(""" ncols             3
4370 nrows             2
4371 xllcorner    148.00000
4372 yllcorner    -38.00000
4373 cellsize       0.25
4374 nodata_value   -9999.0
4375    90.000    60.000    30.0
4376    10.000    10.000    10.000
4377""")
4378        fid.close()
4379        fid = open(ucur_dir_filename2, 'w')
4380        fid.write(""" ncols             3
4381 nrows             2
4382 xllcorner    148.00000
4383 yllcorner    -38.00000
4384 cellsize       0.25
4385 nodata_value   -9999.0
4386    90.000    60.000    30.0
4387    10.000    10.000    10.000
4388""")
4389        fid.close()
4390
4391        vcur_dir =  tempfile.mkdtemp()
4392        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4393        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4394
4395        fid = open(vcur_dir_filename1, 'w')
4396        fid.write(""" ncols             3
4397 nrows             2
4398 xllcorner    148.00000
4399 yllcorner    -38.00000
4400 cellsize       0.25
4401 nodata_value   -9999.0
4402    90.000    60.000    30.0
4403    10.000    10.000    10.000
4404""")
4405        fid.close()
4406        fid = open(vcur_dir_filename2, 'w')
4407        fid.write(""" ncols             3
4408 nrows             2
4409 xllcorner    148.00000
4410 yllcorner    -38.00000
4411 cellsize       0.25
4412 nodata_value   -9999.0
4413    90.000    60.000    30.0
4414    10.000    10.000    10.000
4415""")
4416        fid.close()
4417
4418        try:
4419            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
4420                          vcur_dir, sww_file,
4421                      verbose=self.verbose)
4422        except:
4423            #tidy up
4424            os.remove(bath_dir_filename)
4425            os.rmdir(bath_dir)
4426
4427            os.remove(elevation_dir_filename1)
4428            os.remove(elevation_dir_filename2)
4429            os.rmdir(elevation_dir)
4430
4431            os.remove(ucur_dir_filename1)
4432            os.remove(ucur_dir_filename2)
4433            os.rmdir(ucur_dir)
4434
4435            os.remove(vcur_dir_filename1)
4436            os.remove(vcur_dir_filename2)
4437            os.rmdir(vcur_dir)
4438        else:
4439            #tidy up
4440            os.remove(bath_dir_filename)
4441            os.rmdir(bath_dir)
4442
4443            os.remove(elevation_dir_filename1)
4444            os.remove(elevation_dir_filename2)
4445            os.rmdir(elevation_dir)
4446            raise 'Should raise exception'
4447
4448            os.remove(ucur_dir_filename1)
4449            os.remove(ucur_dir_filename2)
4450            os.rmdir(ucur_dir)
4451
4452            os.remove(vcur_dir_filename1)
4453            os.remove(vcur_dir_filename2)
4454            os.rmdir(vcur_dir)
4455
4456
4457
4458    def test_asc_csiro2sww3(self):
4459        import tempfile
4460
4461        bath_dir = tempfile.mkdtemp()
4462        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4463        #bath_dir = 'bath_data_manager_test'
4464        #print "os.getcwd( )",os.getcwd( )
4465        elevation_dir =  tempfile.mkdtemp()
4466        #elevation_dir = 'elev_expanded'
4467        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4468        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4469
4470        fid = open(bath_dir_filename, 'w')
4471        fid.write(""" ncols             3
4472 nrows             2
4473 xllcorner    148.00000
4474 yllcorner    -38.00000
4475 cellsize       0.25
4476 nodata_value   -9999.0
4477    9000.000    -1000.000    3000.0
4478   -1000.000    9000.000  -1000.000
4479""")
4480        fid.close()
4481
4482        fid = open(elevation_dir_filename1, 'w')
4483        fid.write(""" ncols             3
4484 nrows             2
4485 xllcorner    148.00000
4486 yllcorner    -38.00000
4487 cellsize       0.25
4488 nodata_value   -9999.0
4489    9000.000    0.000    3000.0
4490     0.000     -9999.000     -9999.000
4491""")
4492        fid.close()
4493
4494        fid = open(elevation_dir_filename2, 'w')
4495        fid.write(""" ncols             3
4496 nrows             2
4497 xllcorner    148.00000
4498 yllcorner    -38.00000
4499 cellsize       0.25
4500 nodata_value   -9999.0
4501    9000.000    4000.000    4000.0
4502    4000.000    9000.000    4000.000
4503""")
4504        fid.close()
4505
4506        ucur_dir =  tempfile.mkdtemp()
4507        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4508        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4509
4510        fid = open(ucur_dir_filename1, 'w')
4511        fid.write(""" ncols             3
4512 nrows             2
4513 xllcorner    148.00000
4514 yllcorner    -38.00000
4515 cellsize       0.25
4516 nodata_value   -9999.0
4517    90.000    60.000    30.0
4518    10.000    10.000    10.000
4519""")
4520        fid.close()
4521        fid = open(ucur_dir_filename2, 'w')
4522        fid.write(""" ncols             3
4523 nrows             2
4524 xllcorner    148.00000
4525 yllcorner    -38.00000
4526 cellsize       0.25
4527 nodata_value   -9999.0
4528    90.000    60.000    30.0
4529    10.000    10.000    10.000
4530""")
4531        fid.close()
4532
4533        vcur_dir =  tempfile.mkdtemp()
4534        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4535        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4536
4537        fid = open(vcur_dir_filename1, 'w')
4538        fid.write(""" ncols             3
4539 nrows             2
4540 xllcorner    148.00000
4541 yllcorner    -38.00000
4542 cellsize       0.25
4543 nodata_value   -9999.0
4544    90.000    60.000    30.0
4545    10.000    10.000    10.000
4546""")
4547        fid.close()
4548        fid = open(vcur_dir_filename2, 'w')
4549        fid.write(""" ncols             3
4550 nrows             2
4551 xllcorner    148.00000
4552 yllcorner    -38.00000
4553 cellsize       0.25
4554 nodata_value   -9999.0
4555    90.000    60.000    30.0
4556    10.000    10.000    10.000
4557""")
4558        fid.close()
4559
4560        sww_file = 'a_test.sww'
4561        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
4562                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
4563                      mean_stage = 100,
4564                      verbose=self.verbose)
4565
4566        # check the sww file
4567
4568        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
4569        x = fid.variables['x'][:]
4570        y = fid.variables['y'][:]
4571        z = fid.variables['z'][:]
4572        stage = fid.variables['stage'][:]
4573        xmomentum = fid.variables['xmomentum'][:]
4574        geo_ref = Geo_reference(NetCDFObject=fid)
4575        #print "geo_ref",geo_ref
4576        x_ref = geo_ref.get_xllcorner()
4577        y_ref = geo_ref.get_yllcorner()
4578        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
4579        assert allclose(x_ref, 587798.418) # (-38, 148)
4580        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
4581
4582        #Zone:   55
4583        #Easting:  588095.674  Northing: 5821451.722
4584        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
4585        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
4586
4587        #Zone:   55
4588        #Easting:  632145.632  Northing: 5820863.269
4589        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
4590        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
4591
4592        #Zone:   55
4593        #Easting:  609748.788  Northing: 5793447.860
4594        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
4595        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
4596
4597        assert allclose(z[0],9000.0 )
4598        assert allclose(stage[0][4],100.0 )
4599        assert allclose(stage[0][5],100.0 )
4600
4601        #(100.0 - 9000)*10
4602        assert allclose(xmomentum[0][4], -89000.0 )
4603
4604        #(100.0 - -1000.000)*10
4605        assert allclose(xmomentum[0][5], 11000.0 )
4606
4607        fid.close()
4608
4609        #tidy up
4610        os.remove(bath_dir_filename)
4611        os.rmdir(bath_dir)
4612
4613        os.remove(elevation_dir_filename1)
4614        os.remove(elevation_dir_filename2)
4615        os.rmdir(elevation_dir)
4616
4617        os.remove(ucur_dir_filename1)
4618        os.remove(ucur_dir_filename2)
4619        os.rmdir(ucur_dir)
4620
4621        os.remove(vcur_dir_filename1)
4622        os.remove(vcur_dir_filename2)
4623        os.rmdir(vcur_dir)
4624
4625        # remove sww file
4626        os.remove(sww_file)
4627
4628
4629    def test_asc_csiro2sww4(self):
4630        """
4631        Test specifying the extent
4632        """
4633
4634        import tempfile
4635
4636        bath_dir = tempfile.mkdtemp()
4637        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4638        #bath_dir = 'bath_data_manager_test'
4639        #print "os.getcwd( )",os.getcwd( )
4640        elevation_dir =  tempfile.mkdtemp()
4641        #elevation_dir = 'elev_expanded'
4642        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4643        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4644
4645        fid = open(bath_dir_filename, 'w')
4646        fid.write(""" ncols             4
4647 nrows             4
4648 xllcorner    148.00000
4649 yllcorner    -38.00000
4650 cellsize       0.25
4651 nodata_value   -9999.0
4652   -9000.000    -1000.000   -3000.0 -2000.000
4653   -1000.000    9000.000  -1000.000 -3000.000
4654   -4000.000    6000.000   2000.000 -5000.000
4655   -9000.000    -1000.000   -3000.0 -2000.000
4656""")
4657        fid.close()
4658
4659        fid = open(elevation_dir_filename1, '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   -900.000    -100.000   -300.0 -200.000
4667   -100.000    900.000  -100.000 -300.000
4668   -400.000    600.000   200.000 -500.000
4669   -900.000    -100.000   -300.0 -200.000
4670""")
4671        fid.close()
4672
4673        fid = open(elevation_dir_filename2, 'w')
4674        fid.write(""" ncols             4
4675 nrows             4
4676 xllcorner    148.00000
4677 yllcorner    -38.00000
4678 cellsize       0.25
4679 nodata_value   -9999.0
4680   -990.000    -110.000   -330.0 -220.000
4681   -110.000    990.000  -110.000 -330.000
4682   -440.000    660.000   220.000 -550.000
4683   -990.000    -110.000   -330.0 -220.000
4684""")
4685        fid.close()
4686
4687        ucur_dir =  tempfile.mkdtemp()
4688        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4689        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4690
4691        fid = open(ucur_dir_filename1, 'w')
4692        fid.write(""" ncols             4
4693 nrows             4
4694 xllcorner    148.00000
4695 yllcorner    -38.00000
4696 cellsize       0.25
4697 nodata_value   -9999.0
4698   -90.000    -10.000   -30.0 -20.000
4699   -10.000    90.000  -10.000 -30.000
4700   -40.000    60.000   20.000 -50.000
4701   -90.000    -10.000   -30.0 -20.000
4702""")
4703        fid.close()
4704        fid = open(ucur_dir_filename2, 'w')
4705        fid.write(""" ncols             4
4706 nrows             4
4707 xllcorner    148.00000
4708 yllcorner    -38.00000
4709 cellsize       0.25
4710 nodata_value   -9999.0
4711   -90.000    -10.000   -30.0 -20.000
4712   -10.000    99.000  -11.000 -30.000
4713   -40.000    66.000   22.000 -50.000
4714   -90.000    -10.000   -30.0 -20.000
4715""")
4716        fid.close()
4717
4718        vcur_dir =  tempfile.mkdtemp()
4719        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4720        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4721
4722        fid = open(vcur_dir_filename1, 'w')
4723        fid.write(""" ncols             4
4724 nrows             4
4725 xllcorner    148.00000
4726 yllcorner    -38.00000
4727 cellsize       0.25
4728 nodata_value   -9999.0
4729   -90.000    -10.000   -30.0 -20.000
4730   -10.000    80.000  -20.000 -30.000
4731   -40.000    50.000   10.000 -50.000
4732   -90.000    -10.000   -30.0 -20.000
4733""")
4734        fid.close()
4735        fid = open(vcur_dir_filename2, 'w')
4736        fid.write(""" ncols             4
4737 nrows             4
4738 xllcorner    148.00000
4739 yllcorner    -38.00000
4740 cellsize       0.25
4741 nodata_value   -9999.0
4742   -90.000    -10.000   -30.0 -20.000
4743   -10.000    88.000  -22.000 -30.000
4744   -40.000    55.000   11.000 -50.000
4745   -90.000    -10.000   -30.0 -20.000
4746""")
4747        fid.close()
4748
4749        sww_file = tempfile.mktemp(".sww")
4750        #sww_file = 'a_test.sww'
4751        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
4752                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
4753                      mean_stage = 100,
4754                       minlat = -37.6, maxlat = -37.6,
4755                  minlon = 148.3, maxlon = 148.3,
4756                      verbose=self.verbose
4757                      #,verbose = True
4758                      )
4759
4760        # check the sww file
4761
4762        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
4763        x = fid.variables['x'][:]
4764        y = fid.variables['y'][:]
4765        z = fid.variables['z'][:]
4766        stage = fid.variables['stage'][:]
4767        xmomentum = fid.variables['xmomentum'][:]
4768        ymomentum = fid.variables['ymomentum'][:]
4769        geo_ref = Geo_reference(NetCDFObject=fid)
4770        #print "geo_ref",geo_ref
4771        x_ref = geo_ref.get_xllcorner()
4772        y_ref = geo_ref.get_yllcorner()
4773        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
4774
4775        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
4776        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
4777        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
4778
4779        #Easting:  632145.632  Northing: 5820863.269
4780        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
4781
4782        #print "x",x
4783        #print "y",y
4784        self.failUnless(len(x) == 4,'failed') # 2*2
4785        self.failUnless(len(x) == 4,'failed') # 2*2
4786
4787        #Zone:   55
4788        #Easting:  632145.632  Northing: 5820863.269
4789        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
4790        # magic number - y is close enough for me.
4791        assert allclose(x[3], 632145.63 - x_ref)
4792        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
4793
4794        assert allclose(z[0],9000.0 ) #z is elevation info
4795        #print "z",z
4796        # 2 time steps, 4 points
4797        self.failUnless(xmomentum.shape == (2,4), 'failed')
4798        self.failUnless(ymomentum.shape == (2,4), 'failed')
4799
4800        #(100.0 - -1000.000)*10
4801        #assert allclose(xmomentum[0][5], 11000.0 )
4802
4803        fid.close()
4804
4805        # is the sww file readable?
4806        #Lets see if we can convert it to a dem!
4807        # if you uncomment, remember to delete the file
4808        #print "sww_file",sww_file
4809        #dem_file = tempfile.mktemp(".dem")
4810        domain = sww2domain(sww_file) ###, dem_file)
4811        domain.check_integrity()
4812
4813        #tidy up
4814        os.remove(bath_dir_filename)
4815        os.rmdir(bath_dir)
4816
4817        os.remove(elevation_dir_filename1)
4818        os.remove(elevation_dir_filename2)
4819        os.rmdir(elevation_dir)
4820
4821        os.remove(ucur_dir_filename1)
4822        os.remove(ucur_dir_filename2)
4823        os.rmdir(ucur_dir)
4824
4825        os.remove(vcur_dir_filename1)
4826        os.remove(vcur_dir_filename2)
4827        os.rmdir(vcur_dir)
4828
4829
4830
4831
4832        # remove sww file
4833        os.remove(sww_file)
4834
4835
4836    def test_get_min_max_indexes(self):
4837        latitudes = [3,2,1,0]
4838        longitudes = [0,10,20,30]
4839
4840        # k - lat
4841        # l - lon
4842        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4843            latitudes,longitudes,
4844            -10,4,-10,31)
4845
4846        #print "kmin",kmin;print "kmax",kmax
4847        #print "lmin",lmin;print "lmax",lmax
4848        latitudes_new = latitudes[kmin:kmax]
4849        longitudes_news = longitudes[lmin:lmax]
4850        #print "latitudes_new", latitudes_new
4851        #print "longitudes_news",longitudes_news
4852        self.failUnless(latitudes == latitudes_new and \
4853                        longitudes == longitudes_news,
4854                         'failed')
4855
4856        ## 2nd test
4857        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4858            latitudes,longitudes,
4859            0.5,2.5,5,25)
4860        #print "kmin",kmin;print "kmax",kmax
4861        #print "lmin",lmin;print "lmax",lmax
4862        latitudes_new = latitudes[kmin:kmax]
4863        longitudes_news = longitudes[lmin:lmax]
4864        #print "latitudes_new", latitudes_new
4865        #print "longitudes_news",longitudes_news
4866
4867        self.failUnless(latitudes == latitudes_new and \
4868                        longitudes == longitudes_news,
4869                         'failed')
4870
4871        ## 3rd test
4872        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(\
4873            latitudes,
4874            longitudes,
4875            1.1,1.9,12,17)
4876        #print "kmin",kmin;print "kmax",kmax
4877        #print "lmin",lmin;print "lmax",lmax
4878        latitudes_new = latitudes[kmin:kmax]
4879        longitudes_news = longitudes[lmin:lmax]
4880        #print "latitudes_new", latitudes_new
4881        #print "longitudes_news",longitudes_news
4882
4883        self.failUnless(latitudes_new == [2, 1] and \
4884                        longitudes_news == [10, 20],
4885                         'failed')
4886
4887
4888        ## 4th test
4889        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4890            latitudes,longitudes,
4891                                                      -0.1,1.9,-2,17)
4892        #print "kmin",kmin;print "kmax",kmax
4893        #print "lmin",lmin;print "lmax",lmax
4894        latitudes_new = latitudes[kmin:kmax]
4895        longitudes_news = longitudes[lmin:lmax]
4896        #print "latitudes_new", latitudes_new
4897        #print "longitudes_news",longitudes_news
4898
4899        self.failUnless(latitudes_new == [2, 1, 0] and \
4900                        longitudes_news == [0, 10, 20],
4901                         'failed')
4902        ## 5th test
4903        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4904            latitudes,longitudes,
4905            0.1,1.9,2,17)
4906        #print "kmin",kmin;print "kmax",kmax
4907        #print "lmin",lmin;print "lmax",lmax
4908        latitudes_new = latitudes[kmin:kmax]
4909        longitudes_news = longitudes[lmin:lmax]
4910        #print "latitudes_new", latitudes_new
4911        #print "longitudes_news",longitudes_news
4912
4913        self.failUnless(latitudes_new == [2, 1, 0] and \
4914                        longitudes_news == [0, 10, 20],
4915                         'failed')
4916
4917        ## 6th test
4918
4919        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4920            latitudes,longitudes,
4921            1.5,4,18,32)
4922        #print "kmin",kmin;print "kmax",kmax
4923        #print "lmin",lmin;print "lmax",lmax
4924        latitudes_new = latitudes[kmin:kmax]
4925        longitudes_news = longitudes[lmin:lmax]
4926        #print "latitudes_new", latitudes_new
4927        #print "longitudes_news",longitudes_news
4928
4929        self.failUnless(latitudes_new == [3, 2, 1] and \
4930                        longitudes_news == [10, 20, 30],
4931                         'failed')
4932
4933
4934        ## 7th test
4935        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4936        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4937            latitudes,longitudes,
4938            1.5,1.5,15,15)
4939        #print "kmin",kmin;print "kmax",kmax
4940        #print "lmin",lmin;print "lmax",lmax
4941        latitudes_new = latitudes[kmin:kmax]
4942        longitudes_news = longitudes[lmin:lmax]
4943        m2d = m2d[kmin:kmax,lmin:lmax]
4944        #print "m2d", m2d
4945        #print "latitudes_new", latitudes_new
4946        #print "longitudes_news",longitudes_news
4947
4948        self.failUnless(latitudes_new == [2, 1] and \
4949                        longitudes_news == [10, 20],
4950                         'failed')
4951
4952        self.failUnless(m2d == [[5,6],[9,10]],
4953                         'failed')
4954
4955    def test_get_min_max_indexes_lat_ascending(self):
4956        latitudes = [0,1,2,3]
4957        longitudes = [0,10,20,30]
4958
4959        # k - lat
4960        # l - lon
4961        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4962            latitudes,longitudes,
4963            -10,4,-10,31)
4964
4965        #print "kmin",kmin;print "kmax",kmax
4966        #print "lmin",lmin;print "lmax",lmax
4967        latitudes_new = latitudes[kmin:kmax]
4968        longitudes_news = longitudes[lmin:lmax]
4969        #print "latitudes_new", latitudes_new
4970        #print "longitudes_news",longitudes_news
4971        self.failUnless(latitudes == latitudes_new and \
4972                        longitudes == longitudes_news,
4973                         'failed')
4974
4975        ## 3rd test
4976        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(\
4977            latitudes,
4978            longitudes,
4979            1.1,1.9,12,17)
4980        #print "kmin",kmin;print "kmax",kmax
4981        #print "lmin",lmin;print "lmax",lmax
4982        latitudes_new = latitudes[kmin:kmax]
4983        longitudes_news = longitudes[lmin:lmax]
4984        #print "latitudes_new", latitudes_new
4985        #print "longitudes_news",longitudes_news
4986
4987        self.failUnless(latitudes_new == [1, 2] and \
4988                        longitudes_news == [10, 20],
4989                         'failed')
4990
4991    def test_get_min_max_indexes2(self):
4992        latitudes = [-30,-35,-40,-45]
4993        longitudes = [148,149,150,151]
4994
4995        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4996
4997        # k - lat
4998        # l - lon
4999        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5000            latitudes,longitudes,
5001            -37,-27,147,149.5)
5002
5003        #print "kmin",kmin;print "kmax",kmax
5004        #print "lmin",lmin;print "lmax",lmax
5005        #print "m2d", m2d
5006        #print "latitudes", latitudes
5007        #print "longitudes",longitudes
5008        #print "latitudes[kmax]", latitudes[kmax]
5009        latitudes_new = latitudes[kmin:kmax]
5010        longitudes_new = longitudes[lmin:lmax]
5011        m2d = m2d[kmin:kmax,lmin:lmax]
5012        #print "m2d", m2d
5013        #print "latitudes_new", latitudes_new
5014        #print "longitudes_new",longitudes_new
5015
5016        self.failUnless(latitudes_new == [-30, -35, -40] and \
5017                        longitudes_new == [148, 149,150],
5018                         'failed')
5019        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
5020                         'failed')
5021
5022    def test_get_min_max_indexes3(self):
5023        latitudes = [-30,-35,-40,-45,-50,-55,-60]
5024        longitudes = [148,149,150,151]
5025
5026        # k - lat
5027        # l - lon
5028        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5029            latitudes,longitudes,
5030            -43,-37,148.5,149.5)
5031
5032
5033        #print "kmin",kmin;print "kmax",kmax
5034        #print "lmin",lmin;print "lmax",lmax
5035        #print "latitudes", latitudes
5036        #print "longitudes",longitudes
5037        latitudes_new = latitudes[kmin:kmax]
5038        longitudes_news = longitudes[lmin:lmax]
5039        #print "latitudes_new", latitudes_new
5040        #print "longitudes_news",longitudes_news
5041
5042        self.failUnless(latitudes_new == [-35, -40, -45] and \
5043                        longitudes_news == [148, 149,150],
5044                         'failed')
5045
5046    def test_get_min_max_indexes4(self):
5047        latitudes = [-30,-35,-40,-45,-50,-55,-60]
5048        longitudes = [148,149,150,151]
5049
5050        # k - lat
5051        # l - lon
5052        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5053            latitudes,longitudes)
5054
5055        #print "kmin",kmin;print "kmax",kmax
5056        #print "lmin",lmin;print "lmax",lmax
5057        #print "latitudes", latitudes
5058        #print "longitudes",longitudes
5059        latitudes_new = latitudes[kmin:kmax]
5060        longitudes_news = longitudes[lmin:lmax]
5061        #print "latitudes_new", latitudes_new
5062        #print "longitudes_news",longitudes_news
5063
5064        self.failUnless(latitudes_new == latitudes  and \
5065                        longitudes_news == longitudes,
5066                         'failed')
5067
5068    def test_tsh2sww(self):
5069        import os
5070        import tempfile
5071
5072        tsh_file = tempfile.mktemp(".tsh")
5073        file = open(tsh_file,"w")
5074        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
50750 0.0 0.0 0.0 0.0 0.01 \n \
50761 1.0 0.0 10.0 10.0 0.02  \n \
50772 0.0 1.0 0.0 10.0 0.03  \n \
50783 0.5 0.25 8.0 12.0 0.04  \n \
5079# Vert att title  \n \
5080elevation  \n \
5081stage  \n \
5082friction  \n \
50832 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
50840 0 3 2 -1  -1  1 dsg\n\
50851 0 1 3 -1  0 -1   ole nielsen\n\
50864 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
50870 1 0 2 \n\
50881 0 2 3 \n\
50892 2 3 \n\
50903 3 1 1 \n\
50913 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
50920 216.0 -86.0 \n \
50931 160.0 -167.0 \n \
50942 114.0 -91.0 \n \
50953 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
50960 0 1 0 \n \
50971 1 2 0 \n \
50982 2 0 0 \n \
50990 # <x> <y> ...Mesh Holes... \n \
51000 # <x> <y> <attribute>...Mesh Regions... \n \
51010 # <x> <y> <attribute>...Mesh Regions, area... \n\
5102#Geo reference \n \
510356 \n \
5104140 \n \
5105120 \n")
5106        file.close()
5107
5108        #sww_file = tempfile.mktemp(".sww")
5109        #print "sww_file",sww_file
5110        #print "sww_file",tsh_file
5111        tsh2sww(tsh_file,
5112                      verbose=self.verbose)
5113
5114        os.remove(tsh_file)
5115        os.remove(tsh_file[:-4] + '.sww')
5116       
5117
5118
5119
5120########## testing nbed class ##################
5121    def test_exposure_csv_loading(self):
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        exposure = Exposure_csv(file_name, title_check_list = ['speed','sound'])
5130        exposure.get_column("sound")
5131       
5132        self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
5133                        'FAILED!')
5134        self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
5135                        'FAILED!')
5136       
5137        os.remove(file_name)
5138       
5139    def test_exposure_csv_loadingII(self):
5140       
5141
5142        file_name = tempfile.mktemp(".txt")
5143        file = open(file_name,"w")
5144        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
5145115.0, -21.0, splat, 0.0\n\
5146114.0, -21.7, pow, 10.0\n\
5147114.5, -21.4, bang, 40.0\n")
5148        file.close()
5149        exposure = Exposure_csv(file_name)
5150        exposure.get_column("sound")
5151       
5152        self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
5153                        'FAILED!')
5154        self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
5155                        'FAILED!')
5156       
5157        os.remove(file_name)
5158       
5159    def test_exposure_csv_loading_title_check_list(self):
5160
5161        # I can't get cvs.reader to close the exposure file
5162        # The hacks below are to get around this.       
5163        if sys.platform == 'win32':
5164            file_name = tempfile.gettempdir() + \
5165                    "test_exposure_csv_loading_title_check_list.csv"
5166        else:
5167            file_name = tempfile.mktemp(".csv")
5168        file = open(file_name,"w")
5169        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
5170115.0, -21.0, splat, 0.0\n\
5171114.0, -21.7, pow, 10.0\n\
5172114.5, -21.4, bang, 40.0\n")
5173        file.close()
5174        try:
5175            exposure = Exposure_csv(file_name, title_check_list = ['SOUND'])
5176        except IOError:
5177            pass
5178        else:
5179            self.failUnless(0 ==1,  'Assertion not thrown error!')
5180           
5181        if not sys.platform == 'win32':
5182            os.remove(file_name)
5183       
5184    def test_exposure_csv_cmp(self):
5185        file_name = tempfile.mktemp(".csv")
5186        file = open(file_name,"w")
5187        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
5188115.0, -21.0, splat, 0.0\n\
5189114.0, -21.7, pow, 10.0\n\
5190114.5, -21.4, bang, 40.0\n")
5191        file.close()
5192       
5193        e1 = Exposure_csv(file_name)
5194        e2 = Exposure_csv(file_name)
5195        os.remove(file_name)
5196
5197        self.failUnless(cmp(e1,e2)==0,
5198                        'FAILED!')
5199       
5200        self.failUnless(cmp(e1,"hey")==1,
5201                        'FAILED!')
5202       
5203        file_name = tempfile.mktemp(".csv")
5204        file = open(file_name,"w")
5205        # Note, this has less spaces in the title,
5206        # the instances will be the same.
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        e3 = Exposure_csv(file_name)
5213        os.remove(file_name)
5214
5215        self.failUnless(cmp(e3,e2)==0,
5216                        'FAILED!')
5217       
5218        file_name = tempfile.mktemp(".csv")
5219        file = open(file_name,"w")
5220        # Note, 40 changed to 44 .
5221        file.write("LATITUDE,LONGITUDE ,sound, speed \n\
5222115.0, -21.0, splat, 0.0\n\
5223114.0, -21.7, pow, 10.0\n\
5224114.5, -21.4, bang, 44.0\n")
5225        file.close()
5226        e4 = Exposure_csv(file_name)
5227        os.remove(file_name)
5228        #print "e4",e4._attribute_dic
5229        #print "e2",e2._attribute_dic
5230        self.failUnless(cmp(e4,e2)<>0,
5231                        'FAILED!')
5232       
5233        file_name = tempfile.mktemp(".csv")
5234        file = open(file_name,"w")
5235        # Note, the first two columns are swapped.
5236        file.write("LONGITUDE,LATITUDE ,sound, speed \n\
5237 -21.0,115.0, splat, 0.0\n\
5238 -21.7,114.0, pow, 10.0\n\
5239 -21.4,114.5, bang, 40.0\n")
5240        file.close()
5241        e5 = Exposure_csv(file_name)
5242        os.remove(file_name)
5243
5244        self.failUnless(cmp(e3,e5)<>0,
5245                        'FAILED!')
5246       
5247    def test_exposure_csv_saving(self):
5248       
5249
5250        file_name = tempfile.mktemp(".csv")
5251        file = open(file_name,"w")
5252        file.write("LATITUDE, LONGITUDE ,sound  , speed \n\
5253115.0, -21.0, splat, 0.0\n\
5254114.0, -21.7, pow, 10.0\n\
5255114.5, -21.4, bang, 40.0\n")
5256        file.close()
5257        e1 = Exposure_csv(file_name)
5258       
5259        file_name2 = tempfile.mktemp(".csv")
5260        e1.save(file_name = file_name2)
5261        e2 = Exposure_csv(file_name2)
5262       
5263        self.failUnless(cmp(e1,e2)==0,
5264                        'FAILED!')
5265        os.remove(file_name)
5266        os.remove(file_name2)
5267
5268    def test_exposure_csv_get_location(self):
5269        file_name = tempfile.mktemp(".csv")
5270        file = open(file_name,"w")
5271        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
5272150.916666667, -34.5, splat, 0.0\n\
5273150.0, -34.0, pow, 10.0\n")
5274        file.close()
5275        e1 = Exposure_csv(file_name)
5276
5277        gsd = e1.get_location()
5278       
5279        points = gsd.get_data_points(absolute=True)
5280       
5281        assert allclose(points[0][0], 308728.009)
5282        assert allclose(points[0][1], 6180432.601)
5283        assert allclose(points[1][0],  222908.705)
5284        assert allclose(points[1][1], 6233785.284)
5285        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
5286                        'Bad zone error!')
5287
5288        os.remove(file_name)
5289       
5290    def test_exposure_csv_set_column_get_column(self):
5291        file_name = tempfile.mktemp(".csv")
5292        file = open(file_name,"w")
5293        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
5294150.916666667, -34.5, splat, 0.0\n\
5295150.0, -34.0, pow, 10.0\n")
5296        file.close()
5297        e1 = Exposure_csv(file_name)     
5298        os.remove(file_name)
5299
5300        new_title = "feast"
5301        new_values = ["chicken","soup"]
5302        e1.set_column(new_title, new_values)
5303        returned_values = e1.get_column(new_title)
5304        self.failUnless(returned_values == new_values,
5305                        ' Error!')
5306       
5307        file_name2 = tempfile.mktemp(".csv")
5308        e1.save(file_name = file_name2)
5309        e2 = Exposure_csv(file_name2)
5310        returned_values = e2.get_column(new_title)
5311        self.failUnless(returned_values == new_values,
5312                        ' Error!')       
5313        os.remove(file_name2)
5314
5315    def test_exposure_csv_set_column_get_column_error_checking(self):
5316        file_name = tempfile.mktemp(".csv")
5317        file = open(file_name,"w")
5318        file.write("LONGITUDE , LATITUDE, sound  , speed \n\
5319150.916666667, -34.5, splat, 0.0\n\
5320150.0, -34.0, pow, 10.0\n")
5321        file.close()
5322        e1 = Exposure_csv(file_name)     
5323        os.remove(file_name)
5324
5325        new_title = "sound"
5326        new_values = [12.5,7.6]
5327        try:
5328            e1.set_column(new_title, new_values)
5329        except TitleValueError:
5330            pass
5331        else:
5332            self.failUnless(0 ==1,  'Error not thrown error!')
5333           
5334        e1.set_column(new_title, new_values, overwrite=True)
5335        returned_values = e1.get_column(new_title)
5336        self.failUnless(returned_values == new_values,
5337                        ' Error!')       
5338       
5339        new2_title = "short list"
5340        new2_values = [12.5]
5341        try:
5342            e1.set_column(new2_title, new2_values)
5343        except DataMissingValuesError:
5344            pass
5345        else:
5346            self.failUnless(0 ==1,  'Error not thrown error!')
5347           
5348        new2_title = "long list"
5349        new2_values = [12.5, 7,8]
5350        try:
5351            e1.set_column(new2_title, new2_values)
5352        except DataMissingValuesError:
5353            pass
5354        else:
5355            self.failUnless(0 ==1,  'Error not thrown error!')
5356        file_name2 = tempfile.mktemp(".csv")
5357        e1.save(file_name = file_name2)
5358        e2 = Exposure_csv(file_name2)
5359        returned_values = e2.get_column(new_title)
5360        for returned, new in map(None, returned_values, new_values):
5361            self.failUnless(returned == str(new), ' Error!')
5362        #self.failUnless(returned_values == new_values, ' Error!')       
5363        os.remove(file_name2)
5364       
5365        try:
5366            e1.get_column("toe jam")
5367        except TitleValueError:
5368            pass
5369        else:
5370            self.failUnless(0 ==1,  'Error not thrown error!')
5371           
5372    def test_exposure_csv_loading_x_y(self):
5373       
5374
5375        file_name = tempfile.mktemp(".csv")
5376        file = open(file_name,"w")
5377        file.write("x, y ,sound  , speed \n\
5378115.0, 7, splat, 0.0\n\
5379114.0, 8.0, pow, 10.0\n\
5380114.5, 9., bang, 40.0\n")
5381        file.close()
5382        e1 = Exposure_csv(file_name, is_x_y_locations=True)
5383        gsd = e1.get_location()
5384       
5385        points = gsd.get_data_points(absolute=True)
5386       
5387        assert allclose(points[0][0], 115)
5388        assert allclose(points[0][1], 7)
5389        assert allclose(points[1][0], 114)
5390        assert allclose(points[1][1], 8)
5391        assert allclose(points[2][0], 114.5)
5392        assert allclose(points[2][1], 9)
5393        self.failUnless(gsd.get_geo_reference().get_zone() == -1,
5394                        'Bad zone error!')
5395
5396        os.remove(file_name)
5397
5398           
5399    def test_exposure_csv_loading_x_y2(self):
5400       
5401        csv_file = tempfile.mktemp(".csv")
5402        fd = open(csv_file,'wb')
5403        writer = csv.writer(fd)
5404        writer.writerow(['x','y','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
5405        writer.writerow([5.5,0.5,'199770','130000','Metal','Timber',20])
5406        writer.writerow([4.5,1.0,'150000','76000','Metal','Double Brick',20])
5407        writer.writerow([4.5,1.5,'150000','76000','Metal','Brick Veneer',20])
5408        fd.close()
5409
5410        e1 = Exposure_csv(csv_file)
5411        gsd = e1.get_location()
5412       
5413        points = gsd.get_data_points(absolute=True)
5414        assert allclose(points[0][0], 5.5)
5415        assert allclose(points[0][1], 0.5)
5416        assert allclose(points[1][0], 4.5)
5417        assert allclose(points[1][1], 1.0)
5418        assert allclose(points[2][0], 4.5)
5419        assert allclose(points[2][1], 1.5)
5420        self.failUnless(gsd.get_geo_reference().get_zone() == -1,
5421                        'Bad zone error!')
5422
5423        os.remove(csv_file)
5424
5425    #### TESTS FOR URS 2 SWW  ###     
5426   
5427    def create_mux(self, points_num=None):
5428        # write all the mux stuff.
5429        time_step_count = 3
5430        time_step = 0.5
5431       
5432        longitudes = [150.66667, 150.83334, 151., 151.16667]
5433        latitudes = [-34.5, -34.33333, -34.16667, -34]
5434
5435        if points_num == None:
5436            points_num = len(longitudes) * len(latitudes)
5437
5438        lonlatdeps = []
5439        quantities = ['HA','UA','VA']
5440        mux_names = [WAVEHEIGHT_MUX_LABEL,
5441                     EAST_VELOCITY_LABEL,
5442                     NORTH_VELOCITY_LABEL]
5443        quantities_init = [[],[],[]]
5444        # urs binary is latitude fastest
5445        for i,lon in enumerate(longitudes):
5446            for j,lat in enumerate(latitudes):
5447                _ , e, n = redfearn(lat, lon)
5448                lonlatdeps.append([lon, lat, n])
5449                quantities_init[0].append(e) # HA
5450                quantities_init[1].append(n ) # UA
5451                quantities_init[2].append(e) # VA
5452        #print "lonlatdeps",lonlatdeps
5453
5454        file_handle, base_name = tempfile.mkstemp("")
5455        os.close(file_handle)
5456        os.remove(base_name)
5457       
5458        files = []       
5459        for i,q in enumerate(quantities): 
5460            quantities_init[i] = ensure_numeric(quantities_init[i])
5461            #print "HA_init", HA_init
5462            q_time = zeros((time_step_count, points_num), Float)
5463            for time in range(time_step_count):
5464                q_time[time,:] = quantities_init[i] #* time * 4
5465           
5466            #Write C files
5467            columns = 3 # long, lat , depth
5468            file = base_name + mux_names[i]
5469            #print "base_name file",file
5470            f = open(file, 'wb')
5471            files.append(file)
5472            f.write(pack('i',points_num))
5473            f.write(pack('i',time_step_count))
5474            f.write(pack('f',time_step))
5475
5476            #write lat/long info
5477            for lonlatdep in lonlatdeps:
5478                for float in lonlatdep:
5479                    f.write(pack('f',float))
5480                   
5481            # Write quantity info
5482            for time in  range(time_step_count):
5483                for point_i in range(points_num):
5484                    f.write(pack('f',q_time[time,point_i]))
5485                    #print " mux_names[i]", mux_names[i]
5486                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
5487            f.close()
5488        return base_name, files
5489   
5490    def write_mux(self,lat_long_points, time_step_count, time_step,
5491                  depth=None, ha=None, ua=None, va=None
5492                  ):
5493        """
5494        This will write 3 non-gridded mux files, for testing.
5495        If no quantities are passed in,
5496        na and va quantities will be the Easting values.
5497        Depth and ua will be the Northing value.
5498        """
5499        #print "lat_long_points", lat_long_points
5500        #print "time_step_count",time_step_count
5501        #print "time_step",
5502       
5503        points_num = len(lat_long_points)
5504        lonlatdeps = []
5505        quantities = ['HA','UA','VA']
5506       
5507        mux_names = [WAVEHEIGHT_MUX_LABEL,
5508                     EAST_VELOCITY_LABEL,
5509                     NORTH_VELOCITY_LABEL]
5510        quantities_init = [[],[],[]]
5511        # urs binary is latitude fastest
5512        for point in lat_long_points:
5513            lat = point[0]
5514            lon = point[1]
5515            _ , e, n = redfearn(lat, lon)
5516            if depth is None:
5517                this_depth = n
5518            else:
5519                this_depth = depth
5520            if ha is None:
5521                this_ha = e
5522            else:
5523                this_ha = ha
5524            if ua is None:
5525                this_ua = n
5526            else:
5527                this_ua = ua
5528            if va is None:
5529                this_va = e   
5530            else:
5531                this_va = va         
5532            lonlatdeps.append([lon, lat, this_depth])
5533            quantities_init[0].append(this_ha) # HA
5534            quantities_init[1].append(this_ua) # UA
5535            quantities_init[2].append(this_va) # VA
5536               
5537        file_handle, base_name = tempfile.mkstemp("")
5538        os.close(file_handle)
5539        os.remove(base_name)
5540       
5541        files = []       
5542        for i,q in enumerate(quantities): 
5543            quantities_init[i] = ensure_numeric(quantities_init[i])
5544            #print "HA_init", HA_init
5545            q_time = zeros((time_step_count, points_num), Float)
5546            for time in range(time_step_count):
5547                q_time[time,:] = quantities_init[i] #* time * 4
5548           
5549            #Write C files
5550            columns = 3 # long, lat , depth
5551            file = base_name + mux_names[i]
5552            #print "base_name file",file
5553            f = open(file, 'wb')
5554            files.append(file)
5555            f.write(pack('i',points_num))
5556            f.write(pack('i',time_step_count))
5557            f.write(pack('f',time_step))
5558
5559            #write lat/long info
5560            for lonlatdep in lonlatdeps:
5561                for float in lonlatdep:
5562                    f.write(pack('f',float))
5563                   
5564            # Write quantity info
5565            for time in  range(time_step_count):
5566                for point_i in range(points_num):
5567                    f.write(pack('f',q_time[time,point_i]))
5568                    #print " mux_names[i]", mux_names[i]
5569                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
5570            f.close()
5571        return base_name, files
5572       
5573   
5574    def delete_mux(self, files):
5575        for file in files:
5576            os.remove(file)
5577           
5578    def test_urs2sww_test_fail(self):
5579        points_num = -100
5580        time_step_count = 45
5581        time_step = -7
5582        file_handle, base_name = tempfile.mkstemp("")       
5583        os.close(file_handle)
5584        os.remove(base_name)
5585       
5586        files = []
5587        quantities = ['HA','UA','VA']
5588       
5589        mux_names = [WAVEHEIGHT_MUX_LABEL,
5590                     EAST_VELOCITY_LABEL,
5591                     NORTH_VELOCITY_LABEL]
5592        for i,q in enumerate(quantities): 
5593            #Write C files
5594            columns = 3 # long, lat , depth
5595            file = base_name + mux_names[i]
5596            f = open(file, 'wb')
5597            files.append(file)
5598            f.write(pack('i',points_num))
5599            f.write(pack('i',time_step_count))
5600            f.write(pack('f',time_step))
5601
5602            f.close()   
5603        tide = 1
5604        try:
5605            urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
5606                      verbose=self.verbose)       
5607        except ANUGAError:
5608            pass
5609        else:
5610            self.delete_mux(files)
5611            msg = 'Should have raised exception'
5612            raise msg
5613        sww_file = base_name + '.sww'
5614        self.delete_mux(files)
5615       
5616    def test_urs2sww_test_fail2(self):
5617        base_name = 'Harry-high-pants'
5618        try:
5619            urs2sww(base_name)       
5620        except IOError:
5621            pass
5622        else:
5623            self.delete_mux(files)
5624            msg = 'Should have raised exception'
5625            raise msg
5626           
5627    def test_urs2sww(self):
5628        tide = 1
5629        base_name, files = self.create_mux()
5630        urs2sww(base_name
5631                #, origin=(0,0,0)
5632                , mean_stage=tide
5633                , remove_nc_files=True,
5634                      verbose=self.verbose
5635                )
5636        sww_file = base_name + '.sww'
5637       
5638        #Let's interigate the sww file
5639        # Note, the sww info is not gridded.  It is point data.
5640        fid = NetCDFFile(sww_file)
5641
5642        x = fid.variables['x'][:]
5643        y = fid.variables['y'][:]
5644        geo_reference = Geo_reference(NetCDFObject=fid)
5645
5646       
5647        #Check that first coordinate is correctly represented       
5648        #Work out the UTM coordinates for first point
5649        zone, e, n = redfearn(-34.5, 150.66667)
5650       
5651        assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
5652
5653        # Make x and y absolute
5654        points = geo_reference.get_absolute(map(None, x, y))
5655        points = ensure_numeric(points)
5656        x = points[:,0]
5657        y = points[:,1]
5658       
5659        #Check first value
5660        stage = fid.variables['stage'][:]
5661        xmomentum = fid.variables['xmomentum'][:]
5662        ymomentum = fid.variables['ymomentum'][:]
5663        elevation = fid.variables['elevation'][:]
5664        assert allclose(stage[0,0], e +tide)  #Meters
5665
5666        #Check the momentums - ua
5667        #momentum = velocity*(stage-elevation)
5668        # elevation = - depth
5669        #momentum = velocity_ua *(stage+depth)
5670        # = n*(e+tide+n) based on how I'm writing these files
5671        #
5672        answer_x = n*(e+tide+n)
5673        actual_x = xmomentum[0,0]
5674        #print "answer_x",answer_x
5675        #print "actual_x",actual_x
5676        assert allclose(answer_x, actual_x)  #Meters
5677
5678        #Check the momentums - va
5679        #momentum = velocity*(stage-elevation)
5680        # -(-elevation) since elevation is inverted in mux files
5681        #momentum = velocity_va *(stage+elevation)
5682        # = e*(e+tide+n) based on how I'm writing these files
5683        answer_y = e*(e+tide+n) * -1 # talking into account mux file format
5684        actual_y = ymomentum[0,0]
5685        #print "answer_y",answer_y
5686        #print "actual_y",actual_y
5687        assert allclose(answer_y, actual_y)  #Meters
5688       
5689        assert allclose(answer_x, actual_x)  #Meters
5690       
5691        # check the stage values, first time step.
5692        # These arrays are equal since the Easting values were used as
5693        # the stage
5694        assert allclose(stage[0], x +tide)  #Meters
5695
5696        # check the elevation values.
5697        # -ve since urs measures depth, sww meshers height,
5698        # these arrays are equal since the northing values were used as
5699        # the elevation
5700        assert allclose(-elevation, y)  #Meters
5701       
5702        fid.close()
5703        self.delete_mux(files)
5704        os.remove(sww_file)
5705       
5706    def test_urs2sww_momentum(self):
5707        tide = 1
5708        time_step_count = 3
5709        time_step = 2
5710        #lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
5711        # This is gridded
5712        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
5713        depth=20
5714        ha=2
5715        ua=5
5716        va=-10 #-ve added to take into account mux file format where south
5717               # is positive.
5718        base_name, files = self.write_mux(lat_long_points,
5719                                          time_step_count, time_step,
5720                                          depth=depth,
5721                                          ha=ha,
5722                                          ua=ua,
5723                                          va=va)
5724        # write_mux(self,lat_long_points, time_step_count, time_step,
5725        #          depth=None, ha=None, ua=None, va=None
5726        urs2sww(base_name
5727                #, origin=(0,0,0)
5728                , mean_stage=tide
5729                , remove_nc_files=True,
5730                      verbose=self.verbose
5731                )
5732        sww_file = base_name + '.sww'
5733       
5734        #Let's interigate the sww file
5735        # Note, the sww info is not gridded.  It is point data.
5736        fid = NetCDFFile(sww_file)
5737
5738        x = fid.variables['x'][:]
5739        y = fid.variables['y'][:]
5740        geo_reference = Geo_reference(NetCDFObject=fid)
5741       
5742        #Check first value
5743        stage = fid.variables['stage'][:]
5744        xmomentum = fid.variables['xmomentum'][:]
5745        ymomentum = fid.variables['ymomentum'][:]
5746        elevation = fid.variables['elevation'][:]
5747        #assert allclose(stage[0,0], e + tide)  #Meters
5748        #print "xmomentum", xmomentum
5749        #print "ymomentum", ymomentum
5750        #Check the momentums - ua
5751        #momentum = velocity*water height
5752        #water height = mux_depth + mux_height +tide
5753        #water height = mux_depth + mux_height +tide
5754        #momentum = velocity*(mux_depth + mux_height +tide)
5755        #
5756       
5757        answer = 115
5758        actual = xmomentum[0,0]
5759        assert allclose(answer, actual)  #Meters^2/ sec
5760        answer = 230
5761        actual = ymomentum[0,0]
5762        #print "answer",answer
5763        #print "actual",actual
5764        assert allclose(answer, actual)  #Meters^2/ sec
5765
5766        # check the stage values, first time step.
5767        # These arrays are equal since the Easting values were used as
5768        # the stage
5769
5770        #assert allclose(stage[0], x +tide)  #Meters
5771
5772        # check the elevation values.
5773        # -ve since urs measures depth, sww meshers height,
5774        # these arrays are equal since the northing values were used as
5775        # the elevation
5776        #assert allclose(-elevation, y)  #Meters
5777       
5778        fid.close()
5779        self.delete_mux(files)
5780        os.remove(sww_file)
5781       
5782 
5783    def test_urs2sww_origin(self):
5784        tide = 1
5785        base_name, files = self.create_mux()
5786        urs2sww(base_name
5787                , origin=(0,0,0)
5788                , mean_stage=tide
5789                , remove_nc_files=True,
5790                      verbose=self.verbose
5791                )
5792        sww_file = base_name + '.sww'
5793       
5794        #Let's interigate the sww file
5795        # Note, the sww info is not gridded.  It is point data.
5796        fid = NetCDFFile(sww_file)
5797
5798        #  x and y are absolute
5799        x = fid.variables['x'][:]
5800        y = fid.variables['y'][:]
5801        geo_reference = Geo_reference(NetCDFObject=fid)
5802
5803       
5804        time = fid.variables['time'][:]
5805        #print "time", time
5806        assert allclose([0.,0.5,1.], time)
5807        assert fid.starttime == 0.0
5808        #Check that first coordinate is correctly represented       
5809        #Work out the UTM coordinates for first point
5810        zone, e, n = redfearn(-34.5, 150.66667)       
5811       
5812        assert allclose([x[0],y[0]], [e,n])
5813
5814       
5815        #Check first value
5816        stage = fid.variables['stage'][:]
5817        xmomentum = fid.variables['xmomentum'][:]
5818        ymomentum = fid.variables['ymomentum'][:]
5819        elevation = fid.variables['elevation'][:]
5820        assert allclose(stage[0,0], e +tide)  #Meters
5821
5822        #Check the momentums - ua
5823        #momentum = velocity*(stage-elevation)
5824        #momentum = velocity*(stage+elevation)
5825        # -(-elevation) since elevation is inverted in mux files
5826        # = n*(e+tide+n) based on how I'm writing these files
5827        answer = n*(e+tide+n)
5828        actual = xmomentum[0,0]
5829        assert allclose(answer, actual)  #Meters
5830
5831        # check the stage values, first time step.
5832        # These arrays are equal since the Easting values were used as
5833        # the stage
5834        assert allclose(stage[0], x +tide)  #Meters
5835
5836        # check the elevation values.
5837        # -ve since urs measures depth, sww meshers height,
5838        # these arrays are equal since the northing values were used as
5839        # the elevation
5840        assert allclose(-elevation, y)  #Meters
5841       
5842        fid.close()
5843        self.delete_mux(files)
5844        os.remove(sww_file)
5845 
5846    def test_urs2sww_minmaxlatlong(self):
5847       
5848        #longitudes = [150.66667, 150.83334, 151., 151.16667]
5849        #latitudes = [-34.5, -34.33333, -34.16667, -34]
5850
5851        tide = 1
5852        base_name, files = self.create_mux()
5853        urs2sww(base_name,
5854                minlat=-34.5,
5855                maxlat=-34,
5856                minlon= 150.66667,
5857                maxlon= 151.16667,
5858                mean_stage=tide,
5859                remove_nc_files=True,
5860                      verbose=self.verbose
5861                )
5862        sww_file = base_name + '.sww'
5863       
5864        #Let's interigate the sww file
5865        # Note, the sww info is not gridded.  It is point data.
5866        fid = NetCDFFile(sww_file)
5867       
5868
5869        # Make x and y absolute
5870        x = fid.variables['x'][:]
5871        y = fid.variables['y'][:]
5872        geo_reference = Geo_reference(NetCDFObject=fid)
5873        points = geo_reference.get_absolute(map(None, x, y))
5874        points = ensure_numeric(points)
5875        x = points[:,0]
5876        y = points[:,1]
5877       
5878        #Check that first coordinate is correctly represented       
5879        #Work out the UTM coordinates for first point
5880        zone, e, n = redfearn(-34.5, 150.66667) 
5881        assert allclose([x[0],y[0]], [e,n])
5882
5883       
5884        #Check first value
5885        stage = fid.variables['stage'][:]
5886        xmomentum = fid.variables['xmomentum'][:]
5887        ymomentum = fid.variables['ymomentum'][:]
5888        elevation = fid.variables['elevation'][:]
5889        assert allclose(stage[0,0], e +tide)  #Meters
5890
5891        #Check the momentums - ua
5892        #momentum = velocity*(stage-elevation)
5893        #momentum = velocity*(stage+elevation)
5894        # -(-elevation) since elevation is inverted in mux files
5895        # = n*(e+tide+n) based on how I'm writing these files
5896        answer = n*(e+tide+n)
5897        actual = xmomentum[0,0]
5898        assert allclose(answer, actual)  #Meters
5899
5900        # check the stage values, first time step.
5901        # These arrays are equal since the Easting values were used as
5902        # the stage
5903        assert allclose(stage[0], x +tide)  #Meters
5904
5905        # check the elevation values.
5906        # -ve since urs measures depth, sww meshers height,
5907        # these arrays are equal since the northing values were used as
5908        # the elevation
5909        assert allclose(-elevation, y)  #Meters
5910       
5911        fid.close()
5912        self.delete_mux(files)
5913        os.remove(sww_file)
5914       
5915    def test_urs2sww_minmaxmintmaxt(self):
5916       
5917        #longitudes = [150.66667, 150.83334, 151., 151.16667]
5918        #latitudes = [-34.5, -34.33333, -34.16667, -34]
5919
5920        tide = 1
5921        base_name, files = self.create_mux()
5922        urs2sww(base_name,
5923                mint=0.25,
5924                maxt=0.75,
5925                mean_stage=tide,
5926                remove_nc_files=True,
5927                      verbose=self.verbose
5928                )
5929        sww_file = base_name + '.sww'
5930       
5931        #Let's interigate the sww file
5932        # Note, the sww info is not gridded.  It is point data.
5933        fid = NetCDFFile(sww_file)
5934
5935       
5936        time = fid.variables['time'][:]
5937        assert allclose(time, [0.0]) # the time is relative
5938        assert fid.starttime == 0.5
5939       
5940        fid.close()
5941        self.delete_mux(files)
5942        #print "sww_file", sww_file
5943        os.remove(sww_file)
5944       
5945    def test_lon_lat2grid(self):
5946        lonlatdep = [
5947            [ 113.06700134  ,  -26.06669998 ,   1.        ] ,
5948            [ 113.06700134  ,  -26.33329964 ,   3.        ] ,
5949            [ 113.19999695  ,  -26.06669998 ,   2.        ] ,
5950            [ 113.19999695  ,  -26.33329964 ,   4.        ] ]
5951           
5952        long, lat, quantity = lon_lat2grid(lonlatdep)
5953
5954        for i, result in enumerate(lat):
5955            assert lonlatdep [i][1] == result
5956        assert len(lat) == 2 
5957
5958        for i, result in enumerate(long):
5959            assert lonlatdep [i*2][0] == result
5960        assert len(long) == 2
5961
5962        for i,q in enumerate(quantity):
5963            assert q == i+1
5964           
5965    def test_lon_lat2grid_bad(self):
5966        lonlatdep  = [
5967            [ -26.06669998,  113.06700134,    1.        ],
5968            [ -26.06669998 , 113.19999695 ,   2.        ],
5969            [ -26.06669998 , 113.33300018,    3.        ],
5970            [ -26.06669998 , 113.43299866   , 4.        ],
5971            [ -26.20000076 , 113.06700134,    5.        ],
5972            [ -26.20000076 , 113.19999695 ,   6.        ],
5973            [ -26.20000076 , 113.33300018  ,  7.        ],
5974            [ -26.20000076 , 113.43299866   , 8.        ],
5975            [ -26.33329964 , 113.06700134,    9.        ],
5976            [ -26.33329964 , 113.19999695 ,   10.        ],
5977            [ -26.33329964 , 113.33300018  ,  11.        ],
5978            [ -26.33329964 , 113.43299866 ,   12.        ],
5979            [ -26.43330002 , 113.06700134 ,   13        ],
5980            [ -26.43330002 , 113.19999695 ,   14.        ],
5981            [ -26.43330002 , 113.33300018,    15.        ],
5982            [ -26.43330002 , 113.43299866,    16.        ]]
5983        try:
5984            long, lat, quantity = lon_lat2grid(lonlatdep)
5985        except AssertionError:
5986            pass
5987        else:
5988            msg = 'Should have raised exception'
5989            raise msg
5990       
5991    def test_lon_lat2gridII(self):
5992        lonlatdep = [
5993            [ 113.06700134  ,  -26.06669998 ,   1.        ] ,
5994            [ 113.06700134  ,  -26.33329964 ,   2.        ] ,
5995            [ 113.19999695  ,  -26.06669998 ,   3.        ] ,
5996            [ 113.19999695  ,  -26.344329964 ,   4.        ] ]
5997        try:
5998            long, lat, quantity = lon_lat2grid(lonlatdep)
5999        except AssertionError:
6000            pass
6001        else:
6002            msg = 'Should have raised exception'
6003            raise msg
6004       
6005    #### END TESTS FOR URS 2 SWW  ###
6006
6007    #### TESTS URS UNGRIDDED 2 SWW ###
6008    def test_URS_points_needed(self):
6009       
6010        ll_lat = -21.5
6011        ll_long = 114.5
6012        grid_spacing = 1./60.
6013        lat_amount = 30
6014        long_amount = 30
6015        zone = 50
6016
6017        boundary_polygon = [[250000,7660000],[280000,7660000],
6018                             [280000,7630000],[250000,7630000]]
6019        geo=URS_points_needed(boundary_polygon, zone, 
6020                              ll_lat, ll_long, grid_spacing, 
6021                              lat_amount, long_amount,
6022                              verbose=self.verbose)
6023        # to test this geo, can info from it be transfered onto the boundary
6024        # poly?
6025        #Maybe see if I can fit the data to the polygon - have to represent
6026        # the poly as points though.
6027        #geo.export_points_file("results.txt", as_lat_long=True)
6028        results = ImmutableSet(geo.get_data_points(as_lat_long=True))
6029        #print 'results',results
6030
6031        # These are a set of points that have to be in results
6032        points = []
6033        for i in range(18):
6034            lat = -21.0 - 8./60 - grid_spacing * i
6035            points.append((lat,degminsec2decimal_degrees(114,35,0))) 
6036            points.append((lat,degminsec2decimal_degrees(114,36,0))) 
6037            points.append((lat,degminsec2decimal_degrees(114,52,0))) 
6038            points.append((lat,degminsec2decimal_degrees(114,53,0)))
6039        geo_answer = Geospatial_data(data_points=points,
6040                                     points_are_lats_longs=True) 
6041        #geo_answer.export_points_file("answer.txt", as_lat_long=True) 
6042        answer = ImmutableSet(points)
6043       
6044        outs = answer.difference(results)
6045        #print "outs", outs
6046        # This doesn't work.  Though visualising the results shows that
6047        # it is correct.
6048        #assert answer.issubset(results)
6049        # this is why;
6050        #point (-21.133333333333333, 114.58333333333333)
6051        #result (-21.133333332232368, 114.58333333300342)
6052       
6053        for point in points:
6054            found = False
6055            for result in results:
6056                if allclose(point, result):
6057                    found = True
6058                    break
6059            if not found:
6060                assert False
6061       
6062   
6063    def dave_test_URS_points_needed(self):
6064        ll_lat = -21.51667
6065        ll_long = 114.51667
6066        grid_spacing = 2./60.
6067        lat_amount = 15
6068        long_amount = 15
6069
6070       
6071        boundary_polygon = [[250000,7660000],[280000,7660000],
6072                             [280000,7630000],[250000,7630000]]
6073        URS_points_needed_to_file('a_test_example',boundary_polygon,
6074                                  ll_lat, ll_long, grid_spacing, 
6075                                  lat_amount, long_amount,
6076                                  verbose=self.verbose)
6077       
6078    def X_test_URS_points_neededII(self):
6079        ll_lat = -21.5
6080        ll_long = 114.5
6081        grid_spacing = 1./60.
6082        lat_amount = 30
6083        long_amount = 30
6084
6085        # change this so lats and longs are inputed, then converted
6086       
6087        #boundary_polygon = [[7660000,250000],[7660000,280000],
6088        #                     [7630000,280000],[7630000,250000]]
6089        URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, 
6090                          lat_amount, long_amount,
6091                          verbose=self.verbose)
6092       
6093    #### END TESTS URS UNGRIDDED 2 SWW ###
6094    def test_Urs_points(self):
6095        time_step_count = 3
6096        time_step = 2
6097        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
6098        base_name, files = self.write_mux(lat_long_points,
6099                                          time_step_count, time_step)
6100        for file in files:
6101            urs = Urs_points(file)
6102            assert time_step_count == urs.time_step_count
6103            assert time_step == urs.time_step
6104
6105            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
6106                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
6107                    assert allclose(n, dep[2])
6108                       
6109            count = 0
6110            for slice in urs:
6111                count += 1
6112                #print slice
6113                for lat_lon, quantity in map(None, lat_long_points, slice):
6114                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
6115                    #print "quantity", quantity
6116                    #print "e", e
6117                    #print "n", n
6118                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
6119                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
6120                        assert allclose(e, quantity)
6121                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
6122                        assert allclose(n, quantity)
6123            assert count == time_step_count
6124                     
6125        self.delete_mux(files)
6126
6127    def test_urs_ungridded2sww (self):
6128       
6129        #Zone:   50   
6130        #Easting:  240992.578  Northing: 7620442.472
6131        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6132        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6133        time_step_count = 2
6134        time_step = 400
6135        tide = 9000000
6136        base_name, files = self.write_mux(lat_long,
6137                                          time_step_count, time_step)
6138        urs_ungridded2sww(base_name, mean_stage=tide,
6139                          verbose=self.verbose)
6140       
6141        # now I want to check the sww file ...
6142        sww_file = base_name + '.sww'
6143       
6144        #Let's interigate the sww file
6145        # Note, the sww info is not gridded.  It is point data.
6146        fid = NetCDFFile(sww_file)
6147       
6148        # Make x and y absolute
6149        x = fid.variables['x'][:]
6150        y = fid.variables['y'][:]
6151        geo_reference = Geo_reference(NetCDFObject=fid)
6152        points = geo_reference.get_absolute(map(None, x, y))
6153        points = ensure_numeric(points)
6154        x = points[:,0]
6155        y = points[:,1]
6156       
6157        #Check that first coordinate is correctly represented       
6158        #Work out the UTM coordinates for first point
6159        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6160        assert allclose([x[0],y[0]], [e,n])
6161
6162        #Check the time vector
6163        times = fid.variables['time'][:]
6164       
6165        times_actual = []
6166        for i in range(time_step_count):
6167            times_actual.append(time_step * i)
6168       
6169        assert allclose(ensure_numeric(times),
6170                        ensure_numeric(times_actual))
6171       
6172        #Check first value
6173        stage = fid.variables['stage'][:]
6174        xmomentum = fid.variables['xmomentum'][:]
6175        ymomentum = fid.variables['ymomentum'][:]
6176        elevation = fid.variables['elevation'][:]
6177        assert allclose(stage[0,0], e +tide)  #Meters
6178
6179
6180        #Check the momentums - ua
6181        #momentum = velocity*(stage-elevation)
6182        # elevation = - depth
6183        #momentum = velocity_ua *(stage+depth)
6184        # = n*(e+tide+n) based on how I'm writing these files
6185        #
6186        answer_x = n*(e+tide+n)
6187        actual_x = xmomentum[0,0]
6188        #print "answer_x",answer_x
6189        #print "actual_x",actual_x
6190        assert allclose(answer_x, actual_x)  #Meters
6191       
6192        #Check the momentums - va
6193        #momentum = velocity*(stage-elevation)
6194        # elevation = - depth
6195        #momentum = velocity_va *(stage+depth)
6196        # = e*(e+tide+n) based on how I'm writing these files
6197        #
6198        answer_y = -1*e*(e+tide+n)
6199        actual_y = ymomentum[0,0]
6200        #print "answer_y",answer_y
6201        #print "actual_y",actual_y
6202        assert allclose(answer_y, actual_y)  #Meters
6203
6204        # check the stage values, first time step.
6205        # These arrays are equal since the Easting values were used as
6206        # the stage
6207        assert allclose(stage[0], x +tide)  #Meters
6208        # check the elevation values.
6209        # -ve since urs measures depth, sww meshers height,
6210        # these arrays are equal since the northing values were used as
6211        # the elevation
6212        assert allclose(-elevation, y)  #Meters
6213       
6214        fid.close()
6215        self.delete_mux(files)
6216        os.remove(sww_file)
6217 
6218    def test_urs_ungridded2swwII (self):
6219       
6220        #Zone:   50   
6221        #Easting:  240992.578  Northing: 7620442.472
6222        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6223        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6224        time_step_count = 2
6225        time_step = 400
6226        tide = 9000000
6227        geo_reference = Geo_reference(50, 3434543,34534543)
6228        base_name, files = self.write_mux(lat_long,
6229                                          time_step_count, time_step)
6230        urs_ungridded2sww(base_name, mean_stage=tide, origin = geo_reference,
6231                          verbose=self.verbose)
6232       
6233        # now I want to check the sww file ...
6234        sww_file = base_name + '.sww'
6235       
6236        #Let's interigate the sww file
6237        # Note, the sww info is not gridded.  It is point data.
6238        fid = NetCDFFile(sww_file)
6239       
6240        # Make x and y absolute
6241        x = fid.variables['x'][:]
6242        y = fid.variables['y'][:]
6243        geo_reference = Geo_reference(NetCDFObject=fid)
6244        points = geo_reference.get_absolute(map(None, x, y))
6245        points = ensure_numeric(points)
6246        x = points[:,0]
6247        y = points[:,1]
6248       
6249        #Check that first coordinate is correctly represented       
6250        #Work out the UTM coordinates for first point
6251        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6252        assert allclose([x[0],y[0]], [e,n])
6253
6254        #Check the time vector
6255        times = fid.variables['time'][:]
6256       
6257        times_actual = []
6258        for i in range(time_step_count):
6259            times_actual.append(time_step * i)
6260       
6261        assert allclose(ensure_numeric(times),
6262                        ensure_numeric(times_actual))
6263       
6264        #Check first value
6265        stage = fid.variables['stage'][:]
6266        xmomentum = fid.variables['xmomentum'][:]
6267        ymomentum = fid.variables['ymomentum'][:]
6268        elevation = fid.variables['elevation'][:]
6269        assert allclose(stage[0,0], e +tide)  #Meters
6270
6271        #Check the momentums - ua
6272        #momentum = velocity*(stage-elevation)
6273        # elevation = - depth
6274        #momentum = velocity_ua *(stage+depth)
6275        # = n*(e+tide+n) based on how I'm writing these files
6276        #
6277        answer_x = n*(e+tide+n)
6278        actual_x = xmomentum[0,0]
6279        #print "answer_x",answer_x
6280        #print "actual_x",actual_x
6281        assert allclose(answer_x, actual_x)  #Meters
6282       
6283        #Check the momentums - va
6284        #momentum = velocity*(stage-elevation)
6285        # elevation = - depth
6286        #momentum = velocity_va *(stage+depth)
6287        # = e*(e+tide+n) based on how I'm writing these files
6288        #
6289        answer_y = -1*e*(e+tide+n)
6290        actual_y = ymomentum[0,0]
6291        #print "answer_y",answer_y
6292        #print "actual_y",actual_y
6293        assert allclose(answer_y, actual_y)  #Meters
6294
6295        # check the stage values, first time step.
6296        # These arrays are equal since the Easting values were used as
6297        # the stage
6298        assert allclose(stage[0], x +tide)  #Meters
6299        # check the elevation values.
6300        # -ve since urs measures depth, sww meshers height,
6301        # these arrays are equal since the northing values were used as
6302        # the elevation
6303        assert allclose(-elevation, y)  #Meters
6304       
6305        fid.close()
6306        self.delete_mux(files)
6307        os.remove(sww_file)
6308 
6309    def test_urs_ungridded2swwIII (self):
6310       
6311        #Zone:   50   
6312        #Easting:  240992.578  Northing: 7620442.472
6313        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6314        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6315        time_step_count = 2
6316        time_step = 400
6317        tide = 9000000
6318        base_name, files = self.write_mux(lat_long,
6319                                          time_step_count, time_step)
6320        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6321                          verbose=self.verbose)
6322       
6323        # now I want to check the sww file ...
6324        sww_file = base_name + '.sww'
6325       
6326        #Let's interigate the sww file
6327        # Note, the sww info is not gridded.  It is point data.
6328        fid = NetCDFFile(sww_file)
6329       
6330        # Make x and y absolute
6331        x = fid.variables['x'][:]
6332        y = fid.variables['y'][:]
6333        geo_reference = Geo_reference(NetCDFObject=fid)
6334        points = geo_reference.get_absolute(map(None, x, y))
6335        points = ensure_numeric(points)
6336        x = points[:,0]
6337        y = points[:,1]
6338       
6339        #Check that first coordinate is correctly represented       
6340        #Work out the UTM coordinates for first point
6341        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6342        assert allclose([x[0],y[0]], [e,n])
6343
6344        #Check the time vector
6345        times = fid.variables['time'][:]
6346       
6347        times_actual = []
6348        for i in range(time_step_count):
6349            times_actual.append(time_step * i)
6350       
6351        assert allclose(ensure_numeric(times),
6352                        ensure_numeric(times_actual))
6353       
6354        #Check first value
6355        stage = fid.variables['stage'][:]
6356        xmomentum = fid.variables['xmomentum'][:]
6357        ymomentum = fid.variables['ymomentum'][:]
6358        elevation = fid.variables['elevation'][:]
6359        assert allclose(stage[0,0], e +tide)  #Meters
6360
6361        #Check the momentums - ua
6362        #momentum = velocity*(stage-elevation)
6363        # elevation = - depth
6364        #momentum = velocity_ua *(stage+depth)
6365        # = n*(e+tide+n) based on how I'm writing these files
6366        #
6367        answer_x = n*(e+tide+n)
6368        actual_x = xmomentum[0,0]
6369        #print "answer_x",answer_x
6370        #print "actual_x",actual_x
6371        assert allclose(answer_x, actual_x)  #Meters
6372       
6373        #Check the momentums - va
6374        #momentum = velocity*(stage-elevation)
6375        # elevation = - depth
6376        #momentum = velocity_va *(stage+depth)
6377        # = e*(e+tide+n) based on how I'm writing these files
6378        #
6379        answer_y = -1*e*(e+tide+n)
6380        actual_y = ymomentum[0,0]
6381        #print "answer_y",answer_y
6382        #print "actual_y",actual_y
6383        assert allclose(answer_y, actual_y)  #Meters
6384
6385        # check the stage values, first time step.
6386        # These arrays are equal since the Easting values were used as
6387        # the stage
6388        assert allclose(stage[0], x +tide)  #Meters
6389        # check the elevation values.
6390        # -ve since urs measures depth, sww meshers height,
6391        # these arrays are equal since the northing values were used as
6392        # the elevation
6393        assert allclose(-elevation, y)  #Meters
6394       
6395        fid.close()
6396        self.delete_mux(files)
6397        os.remove(sww_file)
6398
6399       
6400    def test_urs_ungridded_hole (self):
6401       
6402        #Zone:   50   
6403        #Easting:  240992.578  Northing: 7620442.472
6404        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6405        lat_long = [[-20.5, 114.5],
6406                    [-20.6, 114.6],
6407                    [-20.5, 115.],
6408                    [-20.6, 115.],
6409                    [-20.5, 115.5],
6410                    [-20.6, 115.4],
6411                   
6412                    [-21., 114.5],
6413                    [-21., 114.6],
6414                    [-21., 115.5],
6415                    [-21., 115.4],
6416                   
6417                    [-21.5, 114.5],
6418                    [-21.4, 114.6],
6419                    [-21.5, 115.],
6420                    [-21.4, 115.],
6421                    [-21.5, 115.5],
6422                    [-21.4, 115.4]
6423                    ]
6424        time_step_count = 6
6425        time_step = 100
6426        tide = 9000000
6427        base_name, files = self.write_mux(lat_long,
6428                                          time_step_count, time_step)
6429        #Easting:  292110.784  Northing: 7676551.710
6430        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
6431
6432        urs_ungridded2sww(base_name, mean_stage=-240992.0,
6433                          hole_points_UTM=[ 292110.784, 7676551.710 ],
6434                          verbose=self.verbose)
6435       
6436        # now I want to check the sww file ...
6437        sww_file = base_name + '.sww'
6438       
6439        #Let's interigate the sww file
6440        # Note, the sww info is not gridded.  It is point data.
6441        fid = NetCDFFile(sww_file)
6442       
6443        number_of_volumes = fid.variables['volumes']
6444        #print "number_of_volumes",len(number_of_volumes)
6445        assert allclose(16, len(number_of_volumes))
6446       
6447        fid.close()
6448        self.delete_mux(files)
6449        #print "sww_file", sww_file
6450        os.remove(sww_file)
6451       
6452    def test_urs_ungridded_holeII(self):
6453
6454        # Check that if using a hole that returns no triangles,
6455        # urs_ungridded2sww removes the hole label.
6456       
6457        lat_long = [[-20.5, 114.5],
6458                    [-20.6, 114.6],
6459                    [-20.5, 115.5],
6460                    [-20.6, 115.4],
6461                   
6462                   
6463                    [-21.5, 114.5],
6464                    [-21.4, 114.6],
6465                    [-21.5, 115.5],
6466                    [-21.4, 115.4]
6467                    ]
6468        time_step_count = 6
6469        time_step = 100
6470        tide = 9000000
6471        base_name, files = self.write_mux(lat_long,
6472                                          time_step_count, time_step)
6473        #Easting:  292110.784  Northing: 7676551.710
6474        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
6475
6476        urs_ungridded2sww(base_name, mean_stage=-240992.0,
6477                          hole_points_UTM=[ 292110.784, 7676551.710 ],
6478                          verbose=self.verbose)
6479       
6480        # now I want to check the sww file ...
6481        sww_file = base_name + '.sww'
6482        fid = NetCDFFile(sww_file)
6483       
6484        volumes = fid.variables['volumes']
6485        #print "number_of_volumes",len(volumes)
6486
6487        fid.close()
6488        os.remove(sww_file)
6489       
6490        urs_ungridded2sww(base_name, mean_stage=-240992.0)
6491       
6492        # now I want to check the sww file ...
6493        sww_file = base_name + '.sww'
6494        fid = NetCDFFile(sww_file)
6495       
6496        volumes_again = fid.variables['volumes']
6497        #print "number_of_volumes",len(volumes_again)
6498        assert allclose(len(volumes_again),
6499                        len(volumes))
6500        fid.close()
6501        os.remove(sww_file)
6502        self.delete_mux(files) 
6503       
6504    def test_urs_ungridded2sww_mint_maxt (self):
6505       
6506        #Zone:   50   
6507        #Easting:  240992.578  Northing: 7620442.472
6508        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6509        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6510        time_step_count = 6
6511        time_step = 100
6512        tide = 9000000
6513        base_name, files = self.write_mux(lat_long,
6514                                          time_step_count, time_step)
6515        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6516                          mint=101, maxt=500,
6517                          verbose=self.verbose)
6518       
6519        # now I want to check the sww file ...
6520        sww_file = base_name + '.sww'
6521       
6522        #Let's interigate the sww file
6523        # Note, the sww info is not gridded.  It is point data.
6524        fid = NetCDFFile(sww_file)
6525       
6526        # Make x and y absolute
6527        x = fid.variables['x'][:]
6528        y = fid.variables['y'][:]
6529        geo_reference = Geo_reference(NetCDFObject=fid)
6530        points = geo_reference.get_absolute(map(None, x, y))
6531        points = ensure_numeric(points)
6532        x = points[:,0]
6533        y = points[:,1]
6534       
6535        #Check that first coordinate is correctly represented       
6536        #Work out the UTM coordinates for first point
6537        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6538        assert allclose([x[0],y[0]], [e,n])
6539
6540        #Check the time vector
6541        times = fid.variables['time'][:]
6542       
6543        times_actual = [0,100,200,300]
6544       
6545        assert allclose(ensure_numeric(times),
6546                        ensure_numeric(times_actual))
6547       
6548               #Check first value
6549        stage = fid.variables['stage'][:]
6550        xmomentum = fid.variables['xmomentum'][:]
6551        ymomentum = fid.variables['ymomentum'][:]
6552        elevation = fid.variables['elevation'][:]
6553        assert allclose(stage[0,0], e +tide)  #Meters
6554
6555        #Check the momentums - ua
6556        #momentum = velocity*(stage-elevation)
6557        # elevation = - depth
6558        #momentum = velocity_ua *(stage+depth)
6559        # = n*(e+tide+n) based on how I'm writing these files
6560        #
6561        answer_x = n*(e+tide+n)
6562        actual_x = xmomentum[0,0]
6563        #print "answer_x",answer_x
6564        #print "actual_x",actual_x
6565        assert allclose(answer_x, actual_x)  #Meters
6566       
6567        #Check the momentums - va
6568        #momentum = velocity*(stage-elevation)
6569        # elevation = - depth
6570        #momentum = velocity_va *(stage+depth)
6571        # = e*(e+tide+n) based on how I'm writing these files
6572        #
6573        answer_y = -1*e*(e+tide+n)
6574        actual_y = ymomentum[0,0]
6575        #print "answer_y",answer_y
6576        #print "actual_y",actual_y
6577        assert allclose(answer_y, actual_y)  #Meters
6578
6579        # check the stage values, first time step.
6580        # These arrays are equal since the Easting values were used as
6581        # the stage
6582        assert allclose(stage[0], x +tide)  #Meters
6583        # check the elevation values.
6584        # -ve since urs measures depth, sww meshers height,
6585        # these arrays are equal since the northing values were used as
6586        # the elevation
6587        assert allclose(-elevation, y)  #Meters
6588       
6589        fid.close()
6590        self.delete_mux(files)
6591        os.remove(sww_file)
6592       
6593    def test_urs_ungridded2sww_mint_maxtII (self):
6594       
6595        #Zone:   50   
6596        #Easting:  240992.578  Northing: 7620442.472
6597        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6598        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6599        time_step_count = 6
6600        time_step = 100
6601        tide = 9000000
6602        base_name, files = self.write_mux(lat_long,
6603                                          time_step_count, time_step)
6604        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6605                          mint=0, maxt=100000,
6606                          verbose=self.verbose)
6607       
6608        # now I want to check the sww file ...
6609        sww_file = base_name + '.sww'
6610       
6611        #Let's interigate the sww file
6612        # Note, the sww info is not gridded.  It is point data.
6613        fid = NetCDFFile(sww_file)
6614       
6615        # Make x and y absolute
6616        geo_reference = Geo_reference(NetCDFObject=fid)
6617        points = geo_reference.get_absolute(map(None, fid.variables['x'][:],
6618                                                fid.variables['y'][:]))
6619        points = ensure_numeric(points)
6620        x = points[:,0]
6621       
6622        #Check the time vector
6623        times = fid.variables['time'][:]
6624       
6625        times_actual = [0,100,200,300,400,500]
6626        assert allclose(ensure_numeric(times),
6627                        ensure_numeric(times_actual))
6628       
6629        #Check first value
6630        stage = fid.variables['stage'][:]
6631        assert allclose(stage[0], x +tide)
6632       
6633        fid.close()
6634        self.delete_mux(files)
6635        os.remove(sww_file)
6636       
6637    def test_urs_ungridded2sww_mint_maxtIII (self):
6638       
6639        #Zone:   50   
6640        #Easting:  240992.578  Northing: 7620442.472
6641        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6642        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6643        time_step_count = 6
6644        time_step = 100
6645        tide = 9000000
6646        base_name, files = self.write_mux(lat_long,
6647                                          time_step_count, time_step)
6648        try:
6649            urs_ungridded2sww(base_name, mean_stage=tide,
6650                          origin =(50,23432,4343),
6651                          mint=301, maxt=399,
6652                              verbose=self.verbose)
6653        except: 
6654            pass
6655        else:
6656            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
6657
6658        self.delete_mux(files)
6659       
6660    def test_urs_ungridded2sww_mint_maxt_bad (self):       
6661        #Zone:   50   
6662        #Easting:  240992.578  Northing: 7620442.472
6663        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6664        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6665        time_step_count = 6
6666        time_step = 100
6667        tide = 9000000
6668        base_name, files = self.write_mux(lat_long,
6669                                          time_step_count, time_step)
6670        try:
6671            urs_ungridded2sww(base_name, mean_stage=tide,
6672                          origin =(50,23432,4343),
6673                          mint=301, maxt=301,
6674                              verbose=self.verbose)
6675        except: 
6676            pass
6677        else:
6678            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
6679
6680        self.delete_mux(files)
6681
6682       
6683    def test_URS_points_needed_and_urs_ungridded2sww(self):
6684        # This doesn't actually check anything
6685        # 
6686        ll_lat = -21.5
6687        ll_long = 114.5
6688        grid_spacing = 1./60.
6689        lat_amount = 30
6690        long_amount = 30
6691        time_step_count = 2
6692        time_step = 400
6693        tide = -200000
6694        zone = 50
6695
6696        boundary_polygon = [[250000,7660000],[280000,7660000],
6697                             [280000,7630000],[250000,7630000]]
6698        geo=URS_points_needed(boundary_polygon, zone,
6699                              ll_lat, ll_long, grid_spacing, 
6700                              lat_amount, long_amount,
6701                              verbose=self.verbose)
6702        lat_long = geo.get_data_points(as_lat_long=True)
6703        base_name, files = self.write_mux(lat_long,
6704                                          time_step_count, time_step)
6705        urs_ungridded2sww(base_name, mean_stage=tide,
6706                          verbose=self.verbose)
6707        self.delete_mux(files)
6708        os.remove( base_name + '.sww')
6709   
6710    def cache_test_URS_points_needed_and_urs_ungridded2sww(self):
6711       
6712        ll_lat = -21.5
6713        ll_long = 114.5
6714        grid_spacing = 1./60.
6715        lat_amount = 30
6716        long_amount = 30
6717        time_step_count = 2
6718        time_step = 400
6719        tide = -200000
6720        zone = 50
6721
6722        boundary_polygon = [[250000,7660000],[270000,7650000],
6723                             [280000,7630000],[250000,7630000]]
6724        geo=URS_points_needed(boundary_polygon, zone,
6725                              ll_lat, ll_long, grid_spacing, 
6726                              lat_amount, long_amount, use_cache=True,
6727                              verbose=True)
6728       
6729    def visual_test_URS_points_needed_and_urs_ungridded2sww(self):
6730       
6731        ll_lat = -21.5
6732        ll_long = 114.5
6733        grid_spacing = 1./60.
6734        lat_amount = 30
6735        long_amount = 30
6736        time_step_count = 2
6737        time_step = 400
6738        tide = -200000
6739        zone = 50
6740
6741        boundary_polygon = [[250000,7660000],[270000,7650000],
6742                             [280000,7630000],[250000,7630000]]
6743        geo=URS_points_needed(boundary_polygon, zone,
6744                              ll_lat, ll_long, grid_spacing, 
6745                              lat_amount, long_amount)
6746        lat_long = geo.get_data_points(as_lat_long=True)
6747        base_name, files = self.write_mux(lat_long,
6748                                          time_step_count, time_step)
6749        urs_ungridded2sww(base_name, mean_stage=tide)
6750        self.delete_mux(files)
6751        os.remove( base_name + '.sww')
6752        # extend this so it interpolates onto the boundary.
6753        # have it fail if there is NaN
6754
6755    def test_triangulation(self):
6756        #
6757        # 
6758       
6759        filename = tempfile.mktemp("_data_manager.sww")
6760        outfile = NetCDFFile(filename, "w")
6761        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6762        volumes = (0,1,2)
6763        elevation = [0,1,2]
6764        new_origin = None
6765        new_origin = Geo_reference(56, 0, 0)
6766        times = [0, 10]
6767        number_of_volumes = len(volumes)
6768        number_of_points = len(points_utm)
6769        sww = Write_sww()
6770        sww.store_header(outfile, times, number_of_volumes,
6771                         number_of_points, description='fully sick testing',
6772                         verbose=self.verbose)
6773        sww.store_triangulation(outfile, points_utm, volumes,
6774                                elevation,  new_origin=new_origin,
6775                                verbose=self.verbose)       
6776        outfile.close()
6777        fid = NetCDFFile(filename)
6778
6779        x = fid.variables['x'][:]
6780        y = fid.variables['y'][:]
6781        fid.close()
6782
6783        assert allclose(array(map(None, x,y)), points_utm)
6784        os.remove(filename)
6785
6786       
6787    def test_triangulationII(self):
6788        #
6789        # 
6790       
6791        filename = tempfile.mktemp("_data_manager.sww")
6792        outfile = NetCDFFile(filename, "w")
6793        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6794        volumes = (0,1,2)
6795        elevation = [0,1,2]
6796        new_origin = None
6797        #new_origin = Geo_reference(56, 0, 0)
6798        times = [0, 10]
6799        number_of_volumes = len(volumes)
6800        number_of_points = len(points_utm)
6801        sww = Write_sww()
6802        sww.store_header(outfile, times, number_of_volumes,
6803                         number_of_points, description='fully sick testing',
6804                         verbose=self.verbose)
6805        sww.store_triangulation(outfile, points_utm, volumes,
6806                                elevation,  new_origin=new_origin,
6807                                verbose=self.verbose)       
6808        outfile.close()
6809        fid = NetCDFFile(filename)
6810
6811        x = fid.variables['x'][:]
6812        y = fid.variables['y'][:]
6813        results_georef = Geo_reference()
6814        results_georef.read_NetCDF(fid)
6815        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
6816        fid.close()
6817
6818        assert allclose(array(map(None, x,y)), points_utm)
6819        os.remove(filename)
6820
6821       
6822    def test_triangulation_new_origin(self):
6823        #
6824        # 
6825       
6826        filename = tempfile.mktemp("_data_manager.sww")
6827        outfile = NetCDFFile(filename, "w")
6828        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6829        volumes = (0,1,2)
6830        elevation = [0,1,2]
6831        new_origin = None
6832        new_origin = Geo_reference(56, 1, 554354)
6833        points_utm = new_origin.change_points_geo_ref(points_utm)
6834        times = [0, 10]
6835        number_of_volumes = len(volumes)
6836        number_of_points = len(points_utm)
6837        sww = Write_sww()
6838        sww.store_header(outfile, times, number_of_volumes,
6839                         number_of_points, description='fully sick testing',
6840                         verbose=self.verbose)
6841        sww.store_triangulation(outfile, points_utm, volumes,
6842                                elevation,  new_origin=new_origin,
6843                                verbose=self.verbose)       
6844        outfile.close()
6845        fid = NetCDFFile(filename)
6846
6847        x = fid.variables['x'][:]
6848        y = fid.variables['y'][:]
6849        results_georef = Geo_reference()
6850        results_georef.read_NetCDF(fid)
6851        assert results_georef == new_origin
6852        fid.close()
6853
6854        absolute = Geo_reference(56, 0,0)
6855        assert allclose(array( \
6856            absolute.change_points_geo_ref(map(None, x,y),
6857                                           new_origin)),points_utm)
6858       
6859        os.remove(filename)
6860       
6861    def test_triangulation_points_georeference(self):
6862        #
6863        # 
6864       
6865        filename = tempfile.mktemp("_data_manager.sww")
6866        outfile = NetCDFFile(filename, "w")
6867        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6868        volumes = (0,1,2)
6869        elevation = [0,1,2]
6870        new_origin = None
6871        points_georeference = Geo_reference(56, 1, 554354)
6872        points_utm = points_georeference.change_points_geo_ref(points_utm)
6873        times = [0, 10]
6874        number_of_volumes = len(volumes)
6875        number_of_points = len(points_utm)
6876        sww = Write_sww()
6877        sww.store_header(outfile, times, number_of_volumes,
6878                         number_of_points, description='fully sick testing',
6879                         verbose=self.verbose)
6880        sww.store_triangulation(outfile, points_utm, volumes,
6881                                elevation,  new_origin=new_origin,
6882                                points_georeference=points_georeference,
6883                                verbose=self.verbose)       
6884        outfile.close()
6885        fid = NetCDFFile(filename)
6886
6887        x = fid.variables['x'][:]
6888        y = fid.variables['y'][:]
6889        results_georef = Geo_reference()
6890        results_georef.read_NetCDF(fid)
6891        assert results_georef == points_georeference
6892        fid.close()
6893
6894        assert allclose(array(map(None, x,y)), points_utm)
6895        os.remove(filename)
6896       
6897    def test_triangulation_2_geo_refs(self):
6898        #
6899        # 
6900       
6901        filename = tempfile.mktemp("_data_manager.sww")
6902        outfile = NetCDFFile(filename, "w")
6903        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6904        volumes = (0,1,2)
6905        elevation = [0,1,2]
6906        new_origin = Geo_reference(56, 1, 1)
6907        points_georeference = Geo_reference(56, 0, 0)
6908        points_utm = points_georeference.change_points_geo_ref(points_utm)
6909        times = [0, 10]
6910        number_of_volumes = len(volumes)
6911        number_of_points = len(points_utm)
6912        sww = Write_sww()
6913        sww.store_header(outfile, times, number_of_volumes,
6914                         number_of_points, description='fully sick testing',
6915                         verbose=self.verbose)
6916        sww.store_triangulation(outfile, points_utm, volumes,
6917                                elevation,  new_origin=new_origin,
6918                                points_georeference=points_georeference,
6919                                verbose=self.verbose)       
6920        outfile.close()
6921        fid = NetCDFFile(filename)
6922
6923        x = fid.variables['x'][:]
6924        y = fid.variables['y'][:]
6925        results_georef = Geo_reference()
6926        results_georef.read_NetCDF(fid)
6927        assert results_georef == new_origin
6928        fid.close()
6929
6930
6931        absolute = Geo_reference(56, 0,0)
6932        assert allclose(array( \
6933            absolute.change_points_geo_ref(map(None, x,y),
6934                                           new_origin)),points_utm)
6935        os.remove(filename)
6936       
6937    def test_get_data_from_file(self):
6938#    from anuga.abstract_2d_finite_volumes.util import get_data_from_file
6939       
6940        import os
6941       
6942        fileName = tempfile.mktemp(".txt")
6943#        print"filename",fileName
6944        file = open(fileName,"w")
6945        file.write("elevation, stage\n\
69461.0, 3  \n\
69470.0, 4 \n\
69484.0, 3 \n\
69491.0, 6 \n")
6950        file.close()
6951       
6952        header,x = get_data_from_file(fileName)
6953#        print 'x',x
6954        os.remove(fileName)
6955       
6956        assert allclose(x[:,0], [1.0, 0.0,4.0, 1.0])
6957       
6958    def test_get_data_from_file1(self):
6959#    from anuga.abstract_2d_finite_volumes.util import get_data_from_file
6960       
6961        import os
6962       
6963        fileName = tempfile.mktemp(".txt")
6964#        print"filename",fileName
6965        file = open(fileName,"w")
6966        file.write("elevation stage\n\
69671.3 3  \n\
69680.0 4 \n\
69694.5 3.5 \n\
69701.0 6 \n")
6971        file.close()
6972       
6973        header, x = get_data_from_file(fileName,separator_value=' ')
6974        os.remove(fileName)
6975#        x = get_data_from_file(fileName)
6976#        print '1x',x[:,0]
6977       
6978        assert allclose(x[:,0], [1.3, 0.0,4.5, 1.0])
6979       
6980    def test_store_parameters(self):
6981        """tests store temporary file
6982        """
6983       
6984        from os import sep, getenv
6985       
6986        output_dir=''
6987        file_name='details.csv'
6988       
6989        kwargs = {'file_name':'new2.txt',
6990                  'output_dir':output_dir,
6991                  'file_name':file_name,
6992                  'who':'me',
6993                  'what':'detail',
6994                  'how':2,
6995                  'why':241,
6996#                  'completed':345
6997                  }
6998        store_parameters(verbose=False,**kwargs)
6999
7000        temp='detail_temp.csv'
7001        fid = open(temp)
7002        file_header = fid.readline()
7003        file_line = fid.readline()
7004        fid.close()
7005       
7006       
7007        keys = kwargs.keys()
7008        keys.sort()
7009        line=''
7010        header=''
7011        count=0
7012        #used the sorted keys to create the header and line data
7013        for k in keys:
7014#            print "%s = %s" %(k, kwargs[k])
7015            header = header+str(k)
7016            line = line+str(kwargs[k])
7017            count+=1
7018            if count <len(kwargs):
7019                header = header+','
7020                line = line+','
7021        header+='\n'
7022        line+='\n'
7023       
7024       
7025        #file exists
7026        assert access(temp,F_OK)
7027        assert header == file_header
7028        assert line == file_line
7029       
7030        os.remove(temp)
7031       
7032    def test_store_parameters1(self):
7033        """tests store in temporary file and other file
7034        """
7035       
7036        from os import sep, getenv
7037       
7038        output_dir=''
7039        file_name='details.csv'
7040       
7041        kwargs = {'file_name':'new2.txt',
7042                  'output_dir':output_dir,
7043                  'file_name':file_name,
7044                  'who':'me',
7045                  'what':'detail',
7046                  'how':2,
7047                  'why':241,
7048#                  'completed':345
7049                  }
7050        store_parameters(verbose=False,**kwargs)
7051       
7052        kwargs['how']=55
7053        kwargs['completed']=345
7054
7055        keys = kwargs.keys()
7056        keys.sort()
7057        line=''
7058        header=''
7059        count=0
7060        #used the sorted keys to create the header and line data
7061        for k in keys:
7062#            print "%s = %s" %(k, kwargs[k])
7063            header = header+str(k)
7064            line = line+str(kwargs[k])
7065            count+=1
7066            if count <len(kwargs):
7067                header = header+','
7068                line = line+','
7069        header+='\n'
7070        line+='\n'
7071       
7072        kwargs['how']=55
7073        kwargs['completed']=345
7074       
7075        store_parameters(verbose=False,**kwargs)
7076       
7077#        temp='detail_temp.csv'
7078        fid = open(file_name)
7079        file_header = fid.readline()
7080        file_line1 = fid.readline()
7081        file_line2 = fid.readline()
7082        fid.close()
7083       
7084       
7085        #file exists
7086#        print 'header',header,'line',line
7087#        print 'file_header',file_header,'file_line1',file_line1,'file_line2',file_line2
7088        assert access(file_name,F_OK)
7089        assert header == file_header
7090        assert line == file_line1
7091       
7092        temp='detail_temp.csv'
7093        os.remove(temp)
7094        os.remove(file_name)       
7095       
7096    def test_store_parameters2(self):
7097        """tests appending the data to the end of an existing file
7098        """
7099       
7100        from os import sep, getenv
7101       
7102        output_dir=''
7103        file_name='details.csv'
7104       
7105        kwargs = {'file_name':'new2.txt',
7106                  'output_dir':output_dir,
7107                  'file_name':file_name,
7108                  'who':'me',
7109                  'what':'detail',
7110                  'how':2,
7111                  'why':241,
7112                  'completed':345
7113                  }
7114        store_parameters(verbose=False,**kwargs)
7115       
7116        kwargs['how']=55
7117        kwargs['completed']=23.54532
7118       
7119        store_parameters(verbose=False,**kwargs)
7120       
7121        keys = kwargs.keys()
7122        keys.sort()
7123        line=''
7124        header=''
7125        count=0
7126        #used the sorted keys to create the header and line data
7127        for k in keys:
7128#            print "%s = %s" %(k, kwargs[k])
7129            header = header+str(k)
7130            line = line+str(kwargs[k])
7131            count+=1
7132            if count <len(kwargs):
7133                header = header+','
7134                line = line+','
7135        header+='\n'
7136        line+='\n'
7137       
7138        fid = open(file_name)
7139        file_header = fid.readline()
7140        file_line1 = fid.readline()
7141        file_line2 = fid.readline()
7142        fid.close()
7143       
7144        assert access(file_name,F_OK)
7145        assert header == file_header
7146        assert line == file_line2
7147       
7148        os.remove(file_name)       
7149       
7150
7151    def test_get_maximum_inundation(self):
7152        """Test that sww information can be converted correctly to maximum
7153        runup elevation and location (without and with georeferencing)
7154
7155        This test creates a slope and a runup which is maximal (~11m) at around 10s
7156        and levels out to the boundary condition (1m) at about 30s.
7157        """
7158
7159        import time, os
7160        from Numeric import array, zeros, allclose, Float, concatenate
7161        from Scientific.IO.NetCDF import NetCDFFile
7162
7163        #Setup
7164
7165        from mesh_factory import rectangular
7166
7167        # Create basic mesh (100m x 100m)
7168        points, vertices, boundary = rectangular(20, 5, 100, 50)
7169
7170        # Create shallow water domain
7171        domain = Domain(points, vertices, boundary)
7172        domain.default_order = 2
7173        domain.set_minimum_storable_height(0.01)
7174
7175        domain.set_name('runuptest')
7176        swwfile = domain.get_name() + '.sww'
7177
7178        domain.set_datadir('.')
7179        domain.format = 'sww'
7180        domain.smooth = True
7181
7182        # FIXME (Ole): Backwards compatibility
7183        # Look at sww file and see what happens when
7184        # domain.tight_slope_limiters = 1
7185        domain.tight_slope_limiters = 0 
7186       
7187        Br = Reflective_boundary(domain)
7188        Bd = Dirichlet_boundary([1.0,0,0])
7189
7190
7191        #---------- First run without geo referencing
7192       
7193        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
7194        domain.set_quantity('stage', -6)
7195        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
7196
7197        for t in domain.evolve(yieldstep=1, finaltime = 50):
7198            pass
7199
7200
7201        # Check maximal runup
7202        runup = get_maximum_inundation_elevation(swwfile)
7203        location = get_maximum_inundation_location(swwfile)
7204        #print 'Runup, location', runup, location
7205        assert allclose(runup, 11) or allclose(runup, 12) # old limiters
7206        assert allclose(location[0], 15) or allclose(location[0], 10)
7207
7208        # Check final runup
7209        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
7210        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
7211        #print 'Runup, location:',runup, location       
7212        assert allclose(runup, 1)
7213        assert allclose(location[0], 65)
7214
7215        # Check runup restricted to a polygon
7216        p = [[50,1], [99,1], [99,49], [50,49]]
7217        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
7218        location = get_maximum_inundation_location(swwfile, polygon=p)
7219        #print runup, location       
7220        assert allclose(runup, 4)
7221        assert allclose(location[0], 50)               
7222
7223        #Cleanup
7224        os.remove(swwfile)
7225       
7226
7227
7228        #------------- Now the same with georeferencing
7229
7230        domain.time=0.0
7231        E = 308500
7232        N = 6189000
7233        #E = N = 0
7234        domain.geo_reference = Geo_reference(56, E, N)
7235
7236        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
7237        domain.set_quantity('stage', -6)
7238        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
7239
7240        for t in domain.evolve(yieldstep=1, finaltime = 50):
7241            pass
7242
7243        # Check maximal runup
7244        runup = get_maximum_inundation_elevation(swwfile)
7245        location = get_maximum_inundation_location(swwfile)
7246        assert allclose(runup, 11) or allclose(runup, 12) # old limiters
7247        assert allclose(location[0], 15+E) or allclose(location[0], 10+E)
7248
7249        # Check final runup
7250        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
7251        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
7252        assert allclose(runup, 1)
7253        assert allclose(location[0], 65+E)
7254
7255        # Check runup restricted to a polygon
7256        p = array([[50,1], [99,1], [99,49], [50,49]]) + array([E, N])
7257
7258        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
7259        location = get_maximum_inundation_location(swwfile, polygon=p)
7260        assert allclose(runup, 4)
7261        assert allclose(location[0], 50+E)               
7262
7263
7264        #Cleanup
7265        os.remove(swwfile)
7266       
7267    def test_get_all_swwfiles(self):
7268        try:
7269            swwfiles = get_all_swwfiles('','test.txt')  #Invalid
7270        except IOError:
7271            pass
7272        else:
7273            raise 'Should have raised exception' 
7274       
7275    def test_get_all_swwfiles1(self):
7276       
7277        temp_dir = tempfile.mkdtemp('','sww_test')
7278        filename0 = tempfile.mktemp('.sww','test',temp_dir)
7279        filename1 = tempfile.mktemp('.sww','test',temp_dir)
7280        filename2 = tempfile.mktemp('.sww','test',temp_dir)
7281        filename3 = tempfile.mktemp('.sww','test',temp_dir)
7282       
7283        #print'filename', filename0,filename1,filename2,filename3
7284       
7285        fid0 = open(filename0, 'w')
7286        fid1 = open(filename1, 'w')
7287        fid2 = open(filename2, 'w')
7288        fid3 = open(filename3, 'w')
7289
7290        fid0.write('hello')
7291        fid1.write('hello')
7292        fid2.write('hello')
7293        fid3.write('hello')
7294       
7295        fid0.close()
7296        fid1.close()
7297        fid2.close()
7298        fid3.close()
7299       
7300       
7301        dir, name0 = os.path.split(filename0)
7302        #print 'dir',dir,name0
7303       
7304        iterate=get_all_swwfiles(dir,'test')
7305       
7306        del_dir(temp_dir)
7307#        removeall(temp_dir)
7308
7309        _, name0 = os.path.split(filename0) 
7310        #print'name0',name0[:-4],iterate[0]   
7311        _, name1 = os.path.split(filename1)       
7312        _, name2 = os.path.split(filename2)       
7313        _, name3 = os.path.split(filename3)       
7314
7315        assert name0[:-4] in iterate
7316        assert name1[:-4] in iterate
7317        assert name2[:-4] in iterate
7318        assert name3[:-4] in iterate
7319       
7320        assert len(iterate)==4
7321
7322       
7323       
7324
7325
7326
7327#-------------------------------------------------------------
7328if __name__ == "__main__":
7329    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_maximum_inundation')
7330    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_header')
7331    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid_parallel')
7332    suite = unittest.makeSuite(Test_Data_Manager,'test')
7333    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_extrema')
7334
7335   
7336    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
7337        Test_Data_Manager.verbose=True
7338        saveout = sys.stdout   
7339        filename = ".temp_verbose"
7340        fid = open(filename, 'w')
7341        sys.stdout = fid
7342    else:
7343        pass
7344    runner = unittest.TextTestRunner() #verbosity=2)
7345    runner.run(suite)
7346
7347    # Cleaning up
7348    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
7349        sys.stdout = saveout
7350        fid.close() 
7351        os.remove(filename)
7352
7353
Note: See TracBrowser for help on using the repository browser.