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

Last change on this file since 5253 was 5253, checked in by duncan, 16 years ago

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