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

Last change on this file since 5245 was 5245, checked in by sexton, 16 years ago

(1) update of production processes document (2) test for generating points for URS output (3) update of graduate proposal

File size: 248.9 KB
Line 
1#!/usr/bin/env python
2#
3
4
5import unittest
6import copy
7from Numeric import zeros, array, allclose, Float
8from anuga.utilities.numerical_tools import mean
9import tempfile
10import os
11from Scientific.IO.NetCDF import NetCDFFile
12from struct import pack
13from sets import ImmutableSet
14
15from anuga.shallow_water import *
16from anuga.shallow_water.data_manager import *
17from anuga.config import epsilon
18from anuga.utilities.anuga_exceptions import ANUGAError
19from anuga.utilities.numerical_tools import ensure_numeric
20from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
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):
6121        # domain in southern hemisphere zone 51
6122        grid_spacing = 2./60.
6123        poly1 = [[296361.89, 8091928.62],[429495.07,8028278.82], [447230.56,8000674.05],
6124                     [429661.2,7982177.6],[236945.9,7897453.16],[183493.44,7942782.27],
6125                     [226583.04,8008058.96]]
6126        URS_points_needed_to_file('test_example_poly1',poly1,51,
6127                                  grid_spacing,verbose=self.verbose) 
6128
6129    def test_URS_points_needed_poly2(self):
6130        # domain in northern hemisphere zone 47
6131        grid_spacing = 2./60.
6132        poly2 = [[419336.424,810100.845],[342405.0281,711455.8026],[274649.9152,723352.9603],
6133                  [272089.092,972393.0131],[347633.3754,968551.7784],
6134                        [427979.2022,885965.2313],[427659.0993,875721.9386],
6135                        [429259.6138,861317.3083],[436301.8775,840830.723]]
6136        URS_points_needed_to_file('test_example_poly2',poly2,47,
6137                                  grid_spacing,verbose=self.verbose) 
6138       
6139    #### END TESTS URS UNGRIDDED 2 SWW ###
6140    def test_Urs_points(self):
6141        time_step_count = 3
6142        time_step = 2
6143        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
6144        base_name, files = self.write_mux(lat_long_points,
6145                                          time_step_count, time_step)
6146        for file in files:
6147            urs = Urs_points(file)
6148            assert time_step_count == urs.time_step_count
6149            assert time_step == urs.time_step
6150
6151            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
6152                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
6153                    assert allclose(n, dep[2])
6154                       
6155            count = 0
6156            for slice in urs:
6157                count += 1
6158                #print slice
6159                for lat_lon, quantity in map(None, lat_long_points, slice):
6160                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
6161                    #print "quantity", quantity
6162                    #print "e", e
6163                    #print "n", n
6164                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
6165                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
6166                        assert allclose(e, quantity)
6167                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
6168                        assert allclose(n, quantity)
6169            assert count == time_step_count
6170                     
6171        self.delete_mux(files)
6172
6173    def test_urs_ungridded2sww (self):
6174       
6175        #Zone:   50   
6176        #Easting:  240992.578  Northing: 7620442.472
6177        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6178        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6179        time_step_count = 2
6180        time_step = 400
6181        tide = 9000000
6182        base_name, files = self.write_mux(lat_long,
6183                                          time_step_count, time_step)
6184        urs_ungridded2sww(base_name, mean_stage=tide,
6185                          verbose=self.verbose)
6186       
6187        # now I want to check the sww file ...
6188        sww_file = base_name + '.sww'
6189       
6190        #Let's interigate the sww file
6191        # Note, the sww info is not gridded.  It is point data.
6192        fid = NetCDFFile(sww_file)
6193       
6194        # Make x and y absolute
6195        x = fid.variables['x'][:]
6196        y = fid.variables['y'][:]
6197        geo_reference = Geo_reference(NetCDFObject=fid)
6198        points = geo_reference.get_absolute(map(None, x, y))
6199        points = ensure_numeric(points)
6200        x = points[:,0]
6201        y = points[:,1]
6202       
6203        #Check that first coordinate is correctly represented       
6204        #Work out the UTM coordinates for first point
6205        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6206        assert allclose([x[0],y[0]], [e,n])
6207
6208        #Check the time vector
6209        times = fid.variables['time'][:]
6210       
6211        times_actual = []
6212        for i in range(time_step_count):
6213            times_actual.append(time_step * i)
6214       
6215        assert allclose(ensure_numeric(times),
6216                        ensure_numeric(times_actual))
6217       
6218        #Check first value
6219        stage = fid.variables['stage'][:]
6220        xmomentum = fid.variables['xmomentum'][:]
6221        ymomentum = fid.variables['ymomentum'][:]
6222        elevation = fid.variables['elevation'][:]
6223        assert allclose(stage[0,0], e +tide)  #Meters
6224
6225
6226        #Check the momentums - ua
6227        #momentum = velocity*(stage-elevation)
6228        # elevation = - depth
6229        #momentum = velocity_ua *(stage+depth)
6230        # = n*(e+tide+n) based on how I'm writing these files
6231        #
6232        answer_x = n*(e+tide+n)
6233        actual_x = xmomentum[0,0]
6234        #print "answer_x",answer_x
6235        #print "actual_x",actual_x
6236        assert allclose(answer_x, actual_x)  #Meters
6237       
6238        #Check the momentums - va
6239        #momentum = velocity*(stage-elevation)
6240        # elevation = - depth
6241        #momentum = velocity_va *(stage+depth)
6242        # = e*(e+tide+n) based on how I'm writing these files
6243        #
6244        answer_y = -1*e*(e+tide+n)
6245        actual_y = ymomentum[0,0]
6246        #print "answer_y",answer_y
6247        #print "actual_y",actual_y
6248        assert allclose(answer_y, actual_y)  #Meters
6249
6250        # check the stage values, first time step.
6251        # These arrays are equal since the Easting values were used as
6252        # the stage
6253        assert allclose(stage[0], x +tide)  #Meters
6254        # check the elevation values.
6255        # -ve since urs measures depth, sww meshers height,
6256        # these arrays are equal since the northing values were used as
6257        # the elevation
6258        assert allclose(-elevation, y)  #Meters
6259       
6260        fid.close()
6261        self.delete_mux(files)
6262        os.remove(sww_file)
6263 
6264    def test_urs_ungridded2swwII (self):
6265       
6266        #Zone:   50   
6267        #Easting:  240992.578  Northing: 7620442.472
6268        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6269        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6270        time_step_count = 2
6271        time_step = 400
6272        tide = 9000000
6273        geo_reference = Geo_reference(50, 3434543,34534543)
6274        base_name, files = self.write_mux(lat_long,
6275                                          time_step_count, time_step)
6276        urs_ungridded2sww(base_name, mean_stage=tide, origin = geo_reference,
6277                          verbose=self.verbose)
6278       
6279        # now I want to check the sww file ...
6280        sww_file = base_name + '.sww'
6281       
6282        #Let's interigate the sww file
6283        # Note, the sww info is not gridded.  It is point data.
6284        fid = NetCDFFile(sww_file)
6285       
6286        # Make x and y absolute
6287        x = fid.variables['x'][:]
6288        y = fid.variables['y'][:]
6289        geo_reference = Geo_reference(NetCDFObject=fid)
6290        points = geo_reference.get_absolute(map(None, x, y))
6291        points = ensure_numeric(points)
6292        x = points[:,0]
6293        y = points[:,1]
6294       
6295        #Check that first coordinate is correctly represented       
6296        #Work out the UTM coordinates for first point
6297        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6298        assert allclose([x[0],y[0]], [e,n])
6299
6300        #Check the time vector
6301        times = fid.variables['time'][:]
6302       
6303        times_actual = []
6304        for i in range(time_step_count):
6305            times_actual.append(time_step * i)
6306       
6307        assert allclose(ensure_numeric(times),
6308                        ensure_numeric(times_actual))
6309       
6310        #Check first value
6311        stage = fid.variables['stage'][:]
6312        xmomentum = fid.variables['xmomentum'][:]
6313        ymomentum = fid.variables['ymomentum'][:]
6314        elevation = fid.variables['elevation'][:]
6315        assert allclose(stage[0,0], e +tide)  #Meters
6316
6317        #Check the momentums - ua
6318        #momentum = velocity*(stage-elevation)
6319        # elevation = - depth
6320        #momentum = velocity_ua *(stage+depth)
6321        # = n*(e+tide+n) based on how I'm writing these files
6322        #
6323        answer_x = n*(e+tide+n)
6324        actual_x = xmomentum[0,0]
6325        #print "answer_x",answer_x
6326        #print "actual_x",actual_x
6327        assert allclose(answer_x, actual_x)  #Meters
6328       
6329        #Check the momentums - va
6330        #momentum = velocity*(stage-elevation)
6331        # elevation = - depth
6332        #momentum = velocity_va *(stage+depth)
6333        # = e*(e+tide+n) based on how I'm writing these files
6334        #
6335        answer_y = -1*e*(e+tide+n)
6336        actual_y = ymomentum[0,0]
6337        #print "answer_y",answer_y
6338        #print "actual_y",actual_y
6339        assert allclose(answer_y, actual_y)  #Meters
6340
6341        # check the stage values, first time step.
6342        # These arrays are equal since the Easting values were used as
6343        # the stage
6344        assert allclose(stage[0], x +tide)  #Meters
6345        # check the elevation values.
6346        # -ve since urs measures depth, sww meshers height,
6347        # these arrays are equal since the northing values were used as
6348        # the elevation
6349        assert allclose(-elevation, y)  #Meters
6350       
6351        fid.close()
6352        self.delete_mux(files)
6353        os.remove(sww_file)
6354 
6355    def test_urs_ungridded2swwIII (self):
6356       
6357        #Zone:   50   
6358        #Easting:  240992.578  Northing: 7620442.472
6359        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6360        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6361        time_step_count = 2
6362        time_step = 400
6363        tide = 9000000
6364        base_name, files = self.write_mux(lat_long,
6365                                          time_step_count, time_step)
6366        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6367                          verbose=self.verbose)
6368       
6369        # now I want to check the sww file ...
6370        sww_file = base_name + '.sww'
6371       
6372        #Let's interigate the sww file
6373        # Note, the sww info is not gridded.  It is point data.
6374        fid = NetCDFFile(sww_file)
6375       
6376        # Make x and y absolute
6377        x = fid.variables['x'][:]
6378        y = fid.variables['y'][:]
6379        geo_reference = Geo_reference(NetCDFObject=fid)
6380        points = geo_reference.get_absolute(map(None, x, y))
6381        points = ensure_numeric(points)
6382        x = points[:,0]
6383        y = points[:,1]
6384       
6385        #Check that first coordinate is correctly represented       
6386        #Work out the UTM coordinates for first point
6387        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6388        assert allclose([x[0],y[0]], [e,n])
6389
6390        #Check the time vector
6391        times = fid.variables['time'][:]
6392       
6393        times_actual = []
6394        for i in range(time_step_count):
6395            times_actual.append(time_step * i)
6396       
6397        assert allclose(ensure_numeric(times),
6398                        ensure_numeric(times_actual))
6399       
6400        #Check first value
6401        stage = fid.variables['stage'][:]
6402        xmomentum = fid.variables['xmomentum'][:]
6403        ymomentum = fid.variables['ymomentum'][:]
6404        elevation = fid.variables['elevation'][:]
6405        assert allclose(stage[0,0], e +tide)  #Meters
6406
6407        #Check the momentums - ua
6408        #momentum = velocity*(stage-elevation)
6409        # elevation = - depth
6410        #momentum = velocity_ua *(stage+depth)
6411        # = n*(e+tide+n) based on how I'm writing these files
6412        #
6413        answer_x = n*(e+tide+n)
6414        actual_x = xmomentum[0,0]
6415        #print "answer_x",answer_x
6416        #print "actual_x",actual_x
6417        assert allclose(answer_x, actual_x)  #Meters
6418       
6419        #Check the momentums - va
6420        #momentum = velocity*(stage-elevation)
6421        # elevation = - depth
6422        #momentum = velocity_va *(stage+depth)
6423        # = e*(e+tide+n) based on how I'm writing these files
6424        #
6425        answer_y = -1*e*(e+tide+n)
6426        actual_y = ymomentum[0,0]
6427        #print "answer_y",answer_y
6428        #print "actual_y",actual_y
6429        assert allclose(answer_y, actual_y)  #Meters
6430
6431        # check the stage values, first time step.
6432        # These arrays are equal since the Easting values were used as
6433        # the stage
6434        assert allclose(stage[0], x +tide)  #Meters
6435        # check the elevation values.
6436        # -ve since urs measures depth, sww meshers height,
6437        # these arrays are equal since the northing values were used as
6438        # the elevation
6439        assert allclose(-elevation, y)  #Meters
6440       
6441        fid.close()
6442        self.delete_mux(files)
6443        os.remove(sww_file)
6444
6445       
6446    def test_urs_ungridded_hole (self):
6447       
6448        #Zone:   50   
6449        #Easting:  240992.578  Northing: 7620442.472
6450        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6451        lat_long = [[-20.5, 114.5],
6452                    [-20.6, 114.6],
6453                    [-20.5, 115.],
6454                    [-20.6, 115.],
6455                    [-20.5, 115.5],
6456                    [-20.6, 115.4],
6457                   
6458                    [-21., 114.5],
6459                    [-21., 114.6],
6460                    [-21., 115.5],
6461                    [-21., 115.4],
6462                   
6463                    [-21.5, 114.5],
6464                    [-21.4, 114.6],
6465                    [-21.5, 115.],
6466                    [-21.4, 115.],
6467                    [-21.5, 115.5],
6468                    [-21.4, 115.4]
6469                    ]
6470        time_step_count = 6
6471        time_step = 100
6472        tide = 9000000
6473        base_name, files = self.write_mux(lat_long,
6474                                          time_step_count, time_step)
6475        #Easting:  292110.784  Northing: 7676551.710
6476        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
6477
6478        urs_ungridded2sww(base_name, mean_stage=-240992.0,
6479                          hole_points_UTM=[ 292110.784, 7676551.710 ],
6480                          verbose=self.verbose)
6481       
6482        # now I want to check the sww file ...
6483        sww_file = base_name + '.sww'
6484       
6485        #Let's interigate the sww file
6486        # Note, the sww info is not gridded.  It is point data.
6487        fid = NetCDFFile(sww_file)
6488       
6489        number_of_volumes = fid.variables['volumes']
6490        #print "number_of_volumes",len(number_of_volumes)
6491        assert allclose(16, len(number_of_volumes))
6492       
6493        fid.close()
6494        self.delete_mux(files)
6495        #print "sww_file", sww_file
6496        os.remove(sww_file)
6497       
6498    def test_urs_ungridded_holeII(self):
6499
6500        # Check that if using a hole that returns no triangles,
6501        # urs_ungridded2sww removes the hole label.
6502       
6503        lat_long = [[-20.5, 114.5],
6504                    [-20.6, 114.6],
6505                    [-20.5, 115.5],
6506                    [-20.6, 115.4],
6507                   
6508                   
6509                    [-21.5, 114.5],
6510                    [-21.4, 114.6],
6511                    [-21.5, 115.5],
6512                    [-21.4, 115.4]
6513                    ]
6514        time_step_count = 6
6515        time_step = 100
6516        tide = 9000000
6517        base_name, files = self.write_mux(lat_long,
6518                                          time_step_count, time_step)
6519        #Easting:  292110.784  Northing: 7676551.710
6520        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
6521
6522        urs_ungridded2sww(base_name, mean_stage=-240992.0,
6523                          hole_points_UTM=[ 292110.784, 7676551.710 ],
6524                          verbose=self.verbose)
6525       
6526        # now I want to check the sww file ...
6527        sww_file = base_name + '.sww'
6528        fid = NetCDFFile(sww_file)
6529       
6530        volumes = fid.variables['volumes']
6531        #print "number_of_volumes",len(volumes)
6532
6533        fid.close()
6534        os.remove(sww_file)
6535       
6536        urs_ungridded2sww(base_name, mean_stage=-240992.0)
6537       
6538        # now I want to check the sww file ...
6539        sww_file = base_name + '.sww'
6540        fid = NetCDFFile(sww_file)
6541       
6542        volumes_again = fid.variables['volumes']
6543        #print "number_of_volumes",len(volumes_again)
6544        assert allclose(len(volumes_again),
6545                        len(volumes))
6546        fid.close()
6547        os.remove(sww_file)
6548        self.delete_mux(files) 
6549       
6550    def test_urs_ungridded2sww_mint_maxt (self):
6551       
6552        #Zone:   50   
6553        #Easting:  240992.578  Northing: 7620442.472
6554        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6555        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6556        time_step_count = 6
6557        time_step = 100
6558        tide = 9000000
6559        base_name, files = self.write_mux(lat_long,
6560                                          time_step_count, time_step)
6561        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6562                          mint=101, maxt=500,
6563                          verbose=self.verbose)
6564       
6565        # now I want to check the sww file ...
6566        sww_file = base_name + '.sww'
6567       
6568        #Let's interigate the sww file
6569        # Note, the sww info is not gridded.  It is point data.
6570        fid = NetCDFFile(sww_file)
6571       
6572        # Make x and y absolute
6573        x = fid.variables['x'][:]
6574        y = fid.variables['y'][:]
6575        geo_reference = Geo_reference(NetCDFObject=fid)
6576        points = geo_reference.get_absolute(map(None, x, y))
6577        points = ensure_numeric(points)
6578        x = points[:,0]
6579        y = points[:,1]
6580       
6581        #Check that first coordinate is correctly represented       
6582        #Work out the UTM coordinates for first point
6583        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
6584        assert allclose([x[0],y[0]], [e,n])
6585
6586        #Check the time vector
6587        times = fid.variables['time'][:]
6588       
6589        times_actual = [0,100,200,300]
6590       
6591        assert allclose(ensure_numeric(times),
6592                        ensure_numeric(times_actual))
6593       
6594               #Check first value
6595        stage = fid.variables['stage'][:]
6596        xmomentum = fid.variables['xmomentum'][:]
6597        ymomentum = fid.variables['ymomentum'][:]
6598        elevation = fid.variables['elevation'][:]
6599        assert allclose(stage[0,0], e +tide)  #Meters
6600
6601        #Check the momentums - ua
6602        #momentum = velocity*(stage-elevation)
6603        # elevation = - depth
6604        #momentum = velocity_ua *(stage+depth)
6605        # = n*(e+tide+n) based on how I'm writing these files
6606        #
6607        answer_x = n*(e+tide+n)
6608        actual_x = xmomentum[0,0]
6609        #print "answer_x",answer_x
6610        #print "actual_x",actual_x
6611        assert allclose(answer_x, actual_x)  #Meters
6612       
6613        #Check the momentums - va
6614        #momentum = velocity*(stage-elevation)
6615        # elevation = - depth
6616        #momentum = velocity_va *(stage+depth)
6617        # = e*(e+tide+n) based on how I'm writing these files
6618        #
6619        answer_y = -1*e*(e+tide+n)
6620        actual_y = ymomentum[0,0]
6621        #print "answer_y",answer_y
6622        #print "actual_y",actual_y
6623        assert allclose(answer_y, actual_y)  #Meters
6624
6625        # check the stage values, first time step.
6626        # These arrays are equal since the Easting values were used as
6627        # the stage
6628        assert allclose(stage[0], x +tide)  #Meters
6629        # check the elevation values.
6630        # -ve since urs measures depth, sww meshers height,
6631        # these arrays are equal since the northing values were used as
6632        # the elevation
6633        assert allclose(-elevation, y)  #Meters
6634       
6635        fid.close()
6636        self.delete_mux(files)
6637        os.remove(sww_file)
6638       
6639    def test_urs_ungridded2sww_mint_maxtII (self):
6640       
6641        #Zone:   50   
6642        #Easting:  240992.578  Northing: 7620442.472
6643        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6644        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6645        time_step_count = 6
6646        time_step = 100
6647        tide = 9000000
6648        base_name, files = self.write_mux(lat_long,
6649                                          time_step_count, time_step)
6650        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
6651                          mint=0, maxt=100000,
6652                          verbose=self.verbose)
6653       
6654        # now I want to check the sww file ...
6655        sww_file = base_name + '.sww'
6656       
6657        #Let's interigate the sww file
6658        # Note, the sww info is not gridded.  It is point data.
6659        fid = NetCDFFile(sww_file)
6660       
6661        # Make x and y absolute
6662        geo_reference = Geo_reference(NetCDFObject=fid)
6663        points = geo_reference.get_absolute(map(None, fid.variables['x'][:],
6664                                                fid.variables['y'][:]))
6665        points = ensure_numeric(points)
6666        x = points[:,0]
6667       
6668        #Check the time vector
6669        times = fid.variables['time'][:]
6670       
6671        times_actual = [0,100,200,300,400,500]
6672        assert allclose(ensure_numeric(times),
6673                        ensure_numeric(times_actual))
6674       
6675        #Check first value
6676        stage = fid.variables['stage'][:]
6677        assert allclose(stage[0], x +tide)
6678       
6679        fid.close()
6680        self.delete_mux(files)
6681        os.remove(sww_file)
6682       
6683    def test_urs_ungridded2sww_mint_maxtIII (self):
6684       
6685        #Zone:   50   
6686        #Easting:  240992.578  Northing: 7620442.472
6687        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6688        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6689        time_step_count = 6
6690        time_step = 100
6691        tide = 9000000
6692        base_name, files = self.write_mux(lat_long,
6693                                          time_step_count, time_step)
6694        try:
6695            urs_ungridded2sww(base_name, mean_stage=tide,
6696                          origin =(50,23432,4343),
6697                          mint=301, maxt=399,
6698                              verbose=self.verbose)
6699        except: 
6700            pass
6701        else:
6702            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
6703
6704        self.delete_mux(files)
6705       
6706    def test_urs_ungridded2sww_mint_maxt_bad (self):       
6707        #Zone:   50   
6708        #Easting:  240992.578  Northing: 7620442.472
6709        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
6710        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
6711        time_step_count = 6
6712        time_step = 100
6713        tide = 9000000
6714        base_name, files = self.write_mux(lat_long,
6715                                          time_step_count, time_step)
6716        try:
6717            urs_ungridded2sww(base_name, mean_stage=tide,
6718                          origin =(50,23432,4343),
6719                          mint=301, maxt=301,
6720                              verbose=self.verbose)
6721        except: 
6722            pass
6723        else:
6724            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
6725
6726        self.delete_mux(files)
6727
6728       
6729    def test_URS_points_needed_and_urs_ungridded2sww(self):
6730        # This doesn't actually check anything
6731       
6732        ll_lat = -21.5
6733        ll_long = 114.5
6734        grid_spacing = 1./60.
6735        lat_amount = 30
6736        long_amount = 30
6737        time_step_count = 2
6738        time_step = 400
6739        tide = -200000
6740        zone = 50
6741
6742        boundary_polygon = [[250000,7660000],[280000,7660000],
6743                             [280000,7630000],[250000,7630000]]
6744        geo=URS_points_needed(boundary_polygon, zone,
6745                              ll_lat, ll_long, grid_spacing, 
6746                              lat_amount, long_amount,
6747                              verbose=self.verbose)
6748        lat_long = geo.get_data_points(as_lat_long=True)
6749        base_name, files = self.write_mux(lat_long,
6750                                          time_step_count, time_step)
6751        urs_ungridded2sww(base_name, mean_stage=tide,
6752                          verbose=self.verbose)
6753        self.delete_mux(files)
6754        os.remove( base_name + '.sww')
6755   
6756    def cache_test_URS_points_needed_and_urs_ungridded2sww(self):
6757       
6758        ll_lat = -21.5
6759        ll_long = 114.5
6760        grid_spacing = 1./60.
6761        lat_amount = 30
6762        long_amount = 30
6763        time_step_count = 2
6764        time_step = 400
6765        tide = -200000
6766        zone = 50
6767
6768        boundary_polygon = [[250000,7660000],[270000,7650000],
6769                             [280000,7630000],[250000,7630000]]
6770        geo=URS_points_needed(boundary_polygon, zone,
6771                              ll_lat, ll_long, grid_spacing, 
6772                              lat_amount, long_amount, use_cache=True,
6773                              verbose=True)
6774       
6775    def visual_test_URS_points_needed_and_urs_ungridded2sww(self):
6776       
6777        ll_lat = -21.5
6778        ll_long = 114.5
6779        grid_spacing = 1./60.
6780        lat_amount = 30
6781        long_amount = 30
6782        time_step_count = 2
6783        time_step = 400
6784        tide = -200000
6785        zone = 50
6786
6787        boundary_polygon = [[250000,7660000],[270000,7650000],
6788                             [280000,7630000],[250000,7630000]]
6789        geo=URS_points_needed(boundary_polygon, zone,
6790                              ll_lat, ll_long, grid_spacing, 
6791                              lat_amount, long_amount)
6792        lat_long = geo.get_data_points(as_lat_long=True)
6793        base_name, files = self.write_mux(lat_long,
6794                                          time_step_count, time_step)
6795        urs_ungridded2sww(base_name, mean_stage=tide)
6796        self.delete_mux(files)
6797        os.remove( base_name + '.sww')
6798        # extend this so it interpolates onto the boundary.
6799        # have it fail if there is NaN
6800
6801    def test_triangulation(self):
6802        #
6803       
6804       
6805        filename = tempfile.mktemp("_data_manager.sww")
6806        outfile = NetCDFFile(filename, "w")
6807        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6808        volumes = (0,1,2)
6809        elevation = [0,1,2]
6810        new_origin = None
6811        new_origin = Geo_reference(56, 0, 0)
6812        times = [0, 10]
6813        number_of_volumes = len(volumes)
6814        number_of_points = len(points_utm)
6815        sww = Write_sww()
6816        sww.store_header(outfile, times, number_of_volumes,
6817                         number_of_points, description='fully sick testing',
6818                         verbose=self.verbose,sww_precision=Float)
6819        sww.store_triangulation(outfile, points_utm, volumes,
6820                                elevation,  new_origin=new_origin,
6821                                verbose=self.verbose)       
6822        outfile.close()
6823        fid = NetCDFFile(filename)
6824
6825        x = fid.variables['x'][:]
6826        y = fid.variables['y'][:]
6827        fid.close()
6828
6829        assert allclose(array(map(None, x,y)), points_utm)
6830        os.remove(filename)
6831
6832       
6833    def test_triangulationII(self):
6834        #
6835       
6836       
6837        filename = tempfile.mktemp("_data_manager.sww")
6838        outfile = NetCDFFile(filename, "w")
6839        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6840        volumes = (0,1,2)
6841        elevation = [0,1,2]
6842        new_origin = None
6843        #new_origin = Geo_reference(56, 0, 0)
6844        times = [0, 10]
6845        number_of_volumes = len(volumes)
6846        number_of_points = len(points_utm)
6847        sww = Write_sww()
6848        sww.store_header(outfile, times, number_of_volumes,
6849                         number_of_points, description='fully sick testing',
6850                         verbose=self.verbose,sww_precision=Float)
6851        sww.store_triangulation(outfile, points_utm, volumes,
6852                                elevation,  new_origin=new_origin,
6853                                verbose=self.verbose)
6854        outfile.close()
6855        fid = NetCDFFile(filename)
6856
6857        x = fid.variables['x'][:]
6858        y = fid.variables['y'][:]
6859        results_georef = Geo_reference()
6860        results_georef.read_NetCDF(fid)
6861        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
6862        fid.close()
6863
6864        assert allclose(array(map(None, x,y)), points_utm)
6865        os.remove(filename)
6866
6867       
6868    def test_triangulation_new_origin(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, 1, 554354)
6879        points_utm = new_origin.change_points_geo_ref(points_utm)
6880        times = [0, 10]
6881        number_of_volumes = len(volumes)
6882        number_of_points = len(points_utm)
6883        sww = Write_sww()
6884        sww.store_header(outfile, times, number_of_volumes,
6885                         number_of_points, description='fully sick testing',
6886                         verbose=self.verbose,sww_precision=Float)
6887        sww.store_triangulation(outfile, points_utm, volumes,
6888                                elevation,  new_origin=new_origin,
6889                                verbose=self.verbose)
6890        outfile.close()
6891        fid = NetCDFFile(filename)
6892
6893        x = fid.variables['x'][:]
6894        y = fid.variables['y'][:]
6895        results_georef = Geo_reference()
6896        results_georef.read_NetCDF(fid)
6897        assert results_georef == new_origin
6898        fid.close()
6899
6900        absolute = Geo_reference(56, 0,0)
6901        assert allclose(array( \
6902            absolute.change_points_geo_ref(map(None, x,y),
6903                                           new_origin)),points_utm)
6904       
6905        os.remove(filename)
6906       
6907    def test_triangulation_points_georeference(self):
6908        #
6909       
6910       
6911        filename = tempfile.mktemp("_data_manager.sww")
6912        outfile = NetCDFFile(filename, "w")
6913        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6914        volumes = (0,1,2)
6915        elevation = [0,1,2]
6916        new_origin = None
6917        points_georeference = Geo_reference(56, 1, 554354)
6918        points_utm = points_georeference.change_points_geo_ref(points_utm)
6919        times = [0, 10]
6920        number_of_volumes = len(volumes)
6921        number_of_points = len(points_utm)
6922        sww = Write_sww()
6923        sww.store_header(outfile, times, number_of_volumes,
6924                         number_of_points, description='fully sick testing',
6925                         verbose=self.verbose,sww_precision=Float)
6926        sww.store_triangulation(outfile, points_utm, volumes,
6927                                elevation,  new_origin=new_origin,
6928                                points_georeference=points_georeference,
6929                                verbose=self.verbose)       
6930        outfile.close()
6931        fid = NetCDFFile(filename)
6932
6933        x = fid.variables['x'][:]
6934        y = fid.variables['y'][:]
6935        results_georef = Geo_reference()
6936        results_georef.read_NetCDF(fid)
6937        assert results_georef == points_georeference
6938        fid.close()
6939
6940        assert allclose(array(map(None, x,y)), points_utm)
6941        os.remove(filename)
6942       
6943    def test_triangulation_2_geo_refs(self):
6944        #
6945       
6946       
6947        filename = tempfile.mktemp("_data_manager.sww")
6948        outfile = NetCDFFile(filename, "w")
6949        points_utm = array([[0.,0.],[1.,1.], [0.,1.]])
6950        volumes = (0,1,2)
6951        elevation = [0,1,2]
6952        new_origin = Geo_reference(56, 1, 1)
6953        points_georeference = Geo_reference(56, 0, 0)
6954        points_utm = points_georeference.change_points_geo_ref(points_utm)
6955        times = [0, 10]
6956        number_of_volumes = len(volumes)
6957        number_of_points = len(points_utm)
6958        sww = Write_sww()
6959        sww.store_header(outfile, times, number_of_volumes,
6960                         number_of_points, description='fully sick testing',
6961                         verbose=self.verbose,sww_precision=Float)
6962        sww.store_triangulation(outfile, points_utm, volumes,
6963                                elevation,  new_origin=new_origin,
6964                                points_georeference=points_georeference,
6965                                verbose=self.verbose)       
6966        outfile.close()
6967        fid = NetCDFFile(filename)
6968
6969        x = fid.variables['x'][:]
6970        y = fid.variables['y'][:]
6971        results_georef = Geo_reference()
6972        results_georef.read_NetCDF(fid)
6973        assert results_georef == new_origin
6974        fid.close()
6975
6976
6977        absolute = Geo_reference(56, 0,0)
6978        assert allclose(array( \
6979            absolute.change_points_geo_ref(map(None, x,y),
6980                                           new_origin)),points_utm)
6981        os.remove(filename)
6982       
6983    def test_get_data_from_file(self):
6984#    from anuga.abstract_2d_finite_volumes.util import get_data_from_file
6985       
6986        import os
6987       
6988        fileName = tempfile.mktemp(".txt")
6989#        print"filename",fileName
6990        file = open(fileName,"w")
6991        file.write("elevation, stage\n\
69921.0, 3  \n\
69930.0, 4 \n\
69944.0, 3 \n\
69951.0, 6 \n")
6996        file.close()
6997       
6998        header,x = get_data_from_file(fileName)
6999#        print 'x',x
7000        os.remove(fileName)
7001       
7002        assert allclose(x[:,0], [1.0, 0.0,4.0, 1.0])
7003       
7004    def test_get_data_from_file1(self):
7005#    from anuga.abstract_2d_finite_volumes.util import get_data_from_file
7006       
7007        import os
7008       
7009        fileName = tempfile.mktemp(".txt")
7010#        print"filename",fileName
7011        file = open(fileName,"w")
7012        file.write("elevation stage\n\
70131.3 3  \n\
70140.0 4 \n\
70154.5 3.5 \n\
70161.0 6 \n")
7017        file.close()
7018       
7019        header, x = get_data_from_file(fileName,separator_value=' ')
7020        os.remove(fileName)
7021#        x = get_data_from_file(fileName)
7022#        print '1x',x[:,0]
7023       
7024        assert allclose(x[:,0], [1.3, 0.0,4.5, 1.0])
7025       
7026    def test_store_parameters(self):
7027        """tests store temporary file
7028        """
7029       
7030        from os import sep, getenv
7031       
7032        output_dir=''
7033        file_name='details.csv'
7034       
7035        kwargs = {'file_name':'new2.txt',
7036                  'output_dir':output_dir,
7037                  'file_name':file_name,
7038                  'who':'me',
7039                  'what':'detail',
7040                  'how':2,
7041                  'why':241,
7042#                  'completed':345
7043                  }
7044        store_parameters(verbose=False,**kwargs)
7045
7046        temp='detail_temp.csv'
7047        fid = open(temp)
7048        file_header = fid.readline()
7049        file_line = fid.readline()
7050        fid.close()
7051       
7052       
7053        keys = kwargs.keys()
7054        keys.sort()
7055        line=''
7056        header=''
7057        count=0
7058        #used the sorted keys to create the header and line data
7059        for k in keys:
7060#            print "%s = %s" %(k, kwargs[k])
7061            header = header+str(k)
7062            line = line+str(kwargs[k])
7063            count+=1
7064            if count <len(kwargs):
7065                header = header+','
7066                line = line+','
7067        header+='\n'
7068        line+='\n'
7069       
7070       
7071        #file exists
7072        assert access(temp,F_OK)
7073        assert header == file_header
7074        assert line == file_line
7075       
7076        os.remove(temp)
7077       
7078    def test_store_parameters1(self):
7079        """tests store in temporary file and other file
7080        """
7081       
7082        from os import sep, getenv
7083       
7084        output_dir=''
7085        file_name='details.csv'
7086       
7087        kwargs = {'file_name':'new2.txt',
7088                  'output_dir':output_dir,
7089                  'file_name':file_name,
7090                  'who':'me',
7091                  'what':'detail',
7092                  'how':2,
7093                  'why':241,
7094#                  'completed':345
7095                  }
7096        store_parameters(verbose=False,**kwargs)
7097       
7098        kwargs['how']=55
7099        kwargs['completed']=345
7100
7101        keys = kwargs.keys()
7102        keys.sort()
7103        line=''
7104        header=''
7105        count=0
7106        #used the sorted keys to create the header and line data
7107        for k in keys:
7108#            print "%s = %s" %(k, kwargs[k])
7109            header = header+str(k)
7110            line = line+str(kwargs[k])
7111            count+=1
7112            if count <len(kwargs):
7113                header = header+','
7114                line = line+','
7115        header+='\n'
7116        line+='\n'
7117       
7118        kwargs['how']=55
7119        kwargs['completed']=345
7120       
7121        store_parameters(verbose=False,**kwargs)
7122       
7123#        temp='detail_temp.csv'
7124        fid = open(file_name)
7125        file_header = fid.readline()
7126        file_line1 = fid.readline()
7127        file_line2 = fid.readline()
7128        fid.close()
7129       
7130       
7131        #file exists
7132#        print 'header',header,'line',line
7133#        print 'file_header',file_header,'file_line1',file_line1,'file_line2',file_line2
7134        assert access(file_name,F_OK)
7135        assert header == file_header
7136        assert line == file_line1
7137       
7138        temp='detail_temp.csv'
7139        os.remove(temp)
7140        os.remove(file_name)       
7141       
7142    def test_store_parameters2(self):
7143        """tests appending the data to the end of an existing file
7144        """
7145       
7146        from os import sep, getenv
7147       
7148        output_dir=''
7149        file_name='details.csv'
7150       
7151        kwargs = {'file_name':'new2.txt',
7152                  'output_dir':output_dir,
7153                  'file_name':file_name,
7154                  'who':'me',
7155                  'what':'detail',
7156                  'how':2,
7157                  'why':241,
7158                  'completed':345
7159                  }
7160        store_parameters(verbose=False,**kwargs)
7161       
7162        kwargs['how']=55
7163        kwargs['completed']=23.54532
7164       
7165        store_parameters(verbose=False,**kwargs)
7166       
7167        keys = kwargs.keys()
7168        keys.sort()
7169        line=''
7170        header=''
7171        count=0
7172        #used the sorted keys to create the header and line data
7173        for k in keys:
7174#            print "%s = %s" %(k, kwargs[k])
7175            header = header+str(k)
7176            line = line+str(kwargs[k])
7177            count+=1
7178            if count <len(kwargs):
7179                header = header+','
7180                line = line+','
7181        header+='\n'
7182        line+='\n'
7183       
7184        fid = open(file_name)
7185        file_header = fid.readline()
7186        file_line1 = fid.readline()
7187        file_line2 = fid.readline()
7188        fid.close()
7189       
7190        assert access(file_name,F_OK)
7191        assert header == file_header
7192        assert line == file_line2
7193       
7194        os.remove(file_name)       
7195       
7196
7197    def test_get_maximum_inundation(self):
7198        """Test that sww information can be converted correctly to maximum
7199        runup elevation and location (without and with georeferencing)
7200
7201        This test creates a slope and a runup which is maximal (~11m) at around 10s
7202        and levels out to the boundary condition (1m) at about 30s.
7203        """
7204
7205        import time, os
7206        from Numeric import array, zeros, allclose, Float, concatenate
7207        from Scientific.IO.NetCDF import NetCDFFile
7208
7209        #Setup
7210
7211        from mesh_factory import rectangular
7212
7213        # Create basic mesh (100m x 100m)
7214        points, vertices, boundary = rectangular(20, 5, 100, 50)
7215
7216        # Create shallow water domain
7217        domain = Domain(points, vertices, boundary)
7218        domain.default_order = 2
7219        domain.set_minimum_storable_height(0.01)
7220
7221        domain.set_name('runuptest')
7222        swwfile = domain.get_name() + '.sww'
7223
7224        domain.set_datadir('.')
7225        domain.format = 'sww'
7226        domain.smooth = True
7227
7228        # FIXME (Ole): Backwards compatibility
7229        # Look at sww file and see what happens when
7230        # domain.tight_slope_limiters = 1
7231        domain.tight_slope_limiters = 0 
7232       
7233        Br = Reflective_boundary(domain)
7234        Bd = Dirichlet_boundary([1.0,0,0])
7235
7236
7237        #---------- First run without geo referencing
7238       
7239        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
7240        domain.set_quantity('stage', -6)
7241        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
7242
7243        for t in domain.evolve(yieldstep=1, finaltime = 50):
7244            pass
7245
7246
7247        # Check maximal runup
7248        runup = get_maximum_inundation_elevation(swwfile)
7249        location = get_maximum_inundation_location(swwfile)
7250        #print 'Runup, location', runup, location
7251        assert allclose(runup, 11) or allclose(runup, 12) # old limiters
7252        assert allclose(location[0], 15) or allclose(location[0], 10)
7253
7254        # Check final runup
7255        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
7256        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
7257        # print 'Runup, location:',runup, location       
7258        assert allclose(runup, 1)
7259        assert allclose(location[0], 65)
7260
7261        # Check runup restricted to a polygon
7262        p = [[50,1], [99,1], [99,49], [50,49]]
7263        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
7264        location = get_maximum_inundation_location(swwfile, polygon=p)
7265        #print runup, location       
7266        assert allclose(runup, 4)
7267        assert allclose(location[0], 50)               
7268
7269        # Check that mimimum_storable_height works
7270        fid = NetCDFFile(swwfile, 'r') # Open existing file
7271       
7272        stage = fid.variables['stage'][:]
7273        z = fid.variables['elevation'][:]
7274        xmomentum = fid.variables['xmomentum'][:]
7275        ymomentum = fid.variables['ymomentum'][:]       
7276
7277        h = stage-z
7278        for i in range(len(stage)):
7279            if h[i] == 0.0:
7280                assert xmomentum[i] == 0.0
7281                assert ymomentum[i] == 0.0               
7282            else:
7283                assert h[i] >= domain.minimum_storable_height
7284       
7285        fid.close()
7286
7287        # Cleanup
7288        os.remove(swwfile)
7289       
7290
7291
7292        #------------- Now the same with georeferencing
7293
7294        domain.time=0.0
7295        E = 308500
7296        N = 6189000
7297        #E = N = 0
7298        domain.geo_reference = Geo_reference(56, E, N)
7299
7300        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
7301        domain.set_quantity('stage', -6)
7302        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
7303
7304        for t in domain.evolve(yieldstep=1, finaltime = 50):
7305            pass
7306
7307        # Check maximal runup
7308        runup = get_maximum_inundation_elevation(swwfile)
7309        location = get_maximum_inundation_location(swwfile)
7310        assert allclose(runup, 11) or allclose(runup, 12) # old limiters
7311        assert allclose(location[0], 15+E) or allclose(location[0], 10+E)
7312
7313        # Check final runup
7314        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
7315        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
7316        assert allclose(runup, 1)
7317        assert allclose(location[0], 65+E)
7318
7319        # Check runup restricted to a polygon
7320        p = array([[50,1], [99,1], [99,49], [50,49]]) + array([E, N])
7321
7322        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
7323        location = get_maximum_inundation_location(swwfile, polygon=p)
7324        assert allclose(runup, 4)
7325        assert allclose(location[0], 50+E)               
7326
7327
7328        # Cleanup
7329        os.remove(swwfile)
7330
7331
7332    def NOtest_get_flow_through_cross_section(self):
7333        """test_get_flow_through_cross_section(self):
7334
7335        Test that the total flow through a cross section can be
7336        correctly obtained from an sww file.
7337       
7338        This test creates a flat bed with a known flow through it and tests
7339        that the function correctly returns the expected flow.
7340
7341        The specifics are
7342        u = 2 m/s
7343        h = 1 m
7344        w = 5 m (width of channel)
7345
7346        q = u*h*w = 10 m^3/s
7347       
7348        """
7349
7350        import time, os
7351        from Numeric import array, zeros, allclose, Float, concatenate
7352        from Scientific.IO.NetCDF import NetCDFFile
7353
7354        # Setup
7355        from mesh_factory import rectangular
7356
7357        # Create basic mesh (100m x 5m)
7358        width = 5
7359        len = 100
7360        points, vertices, boundary = rectangular(len, width, 100, 5)
7361
7362        # Create shallow water domain
7363        domain = Domain(points, vertices, boundary)
7364        domain.default_order = 2
7365        domain.set_minimum_storable_height(0.01)
7366
7367        domain.set_name('flowtest')
7368        swwfile = domain.get_name() + '.sww'
7369
7370        domain.set_datadir('.')
7371        domain.format = 'sww'
7372        domain.smooth = True
7373
7374        h = 1.0
7375        u = 2.0
7376        uh = u*h
7377
7378        Br = Reflective_boundary(domain)     # Side walls
7379        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet:
7380
7381
7382        #---------- First run without geo referencing
7383       
7384        domain.set_quantity('elevation', 0.0)
7385        domain.set_quantity('stage', h)
7386        domain.set_quantity('xmomentum', uh)
7387        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
7388
7389        for t in domain.evolve(yieldstep=1, finaltime = 50):
7390            pass
7391
7392        # Check that momentum is as it should be in the interior
7393        f = file_function(swwfile,
7394                          quantities=['stage', 'xmomentum', 'ymomentum'],
7395                          interpolation_points=[[0,width/2],
7396                                                [len/2, width/2],
7397                                                [len, width/2]],
7398                          verbose=False)
7399
7400        for t in range(50):
7401            for i in range(3):
7402                assert allclose(f(t, i), [1, 2, 0])
7403           
7404
7405
7406        # Check flow through the middle
7407        cross_section = [[len/2,0], [len/2,width]]
7408        Q = get_flow_through_cross_section(swwfile,
7409                                           cross_section,
7410                                           verbose=True)
7411
7412        assert allclose(Q, uh*width) 
7413                                     
7414
7415
7416
7417
7418       
7419       
7420    def test_get_all_swwfiles(self):
7421        try:
7422            swwfiles = get_all_swwfiles('','test.txt')  #Invalid
7423        except IOError:
7424            pass
7425        else:
7426            raise 'Should have raised exception' 
7427       
7428    def test_get_all_swwfiles1(self):
7429       
7430        temp_dir = tempfile.mkdtemp('','sww_test')
7431        filename0 = tempfile.mktemp('.sww','test',temp_dir)
7432        filename1 = tempfile.mktemp('.sww','test',temp_dir)
7433        filename2 = tempfile.mktemp('.sww','test',temp_dir)
7434        filename3 = tempfile.mktemp('.sww','test',temp_dir)
7435       
7436        #print'filename', filename0,filename1,filename2,filename3
7437       
7438        fid0 = open(filename0, 'w')
7439        fid1 = open(filename1, 'w')
7440        fid2 = open(filename2, 'w')
7441        fid3 = open(filename3, 'w')
7442
7443        fid0.write('hello')
7444        fid1.write('hello')
7445        fid2.write('hello')
7446        fid3.write('hello')
7447       
7448        fid0.close()
7449        fid1.close()
7450        fid2.close()
7451        fid3.close()
7452       
7453       
7454        dir, name0 = os.path.split(filename0)
7455        #print 'dir',dir,name0
7456       
7457        iterate=get_all_swwfiles(dir,'test')
7458       
7459        del_dir(temp_dir)
7460#        removeall(temp_dir)
7461
7462        _, name0 = os.path.split(filename0) 
7463        #print'name0',name0[:-4],iterate[0]   
7464        _, name1 = os.path.split(filename1)       
7465        _, name2 = os.path.split(filename2)       
7466        _, name3 = os.path.split(filename3)       
7467
7468        assert name0[:-4] in iterate
7469        assert name1[:-4] in iterate
7470        assert name2[:-4] in iterate
7471        assert name3[:-4] in iterate
7472       
7473        assert len(iterate)==4
7474
7475    def test_screen_catcher(self):
7476   
7477        stdout_orginal = sys.stdout
7478        stderr_orginal = sys.stderr
7479       
7480        output_dir = tempfile.mkdtemp('','output_test')
7481        #print output_dir
7482        start_screen_catcher(output_dir+sep,verbose=False)
7483       
7484        print 'hello screen catcher'
7485        print 'goodbye screen catcher'
7486       
7487        sys.stdout = stdout_orginal
7488        sys.stderr = stderr_orginal
7489       
7490        output_file = output_dir+sep+'screen_output.txt'
7491        assert access(output_file,F_OK)
7492       
7493        fid = open(output_file)
7494        file_line1 = fid.readline()
7495        file_line2 = fid.readline()
7496       
7497        #print 'file contents',file_line1, file_line2
7498        assert (file_line1 == 'hello screen catcher\n')
7499        assert (file_line2 == 'goodbye screen catcher\n')
7500       
7501 
7502    def test_points2polygon(self): 
7503        att_dict = {}
7504        pointlist = array([[1.0, 0.0],[0.0, 1.0],[0.0, 0.0]])
7505        att_dict['elevation'] = array([10.1, 0.0, 10.4])
7506        att_dict['brightness'] = array([10.0, 1.0, 10.4])
7507       
7508        fileName = tempfile.mktemp(".csv")
7509       
7510        G = Geospatial_data(pointlist, att_dict)
7511       
7512        G.export_points_file(fileName)
7513       
7514        polygon = points2polygon(fileName)
7515       
7516        # This test may fail if the order changes
7517        assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])
7518       
7519
7520
7521
7522#-------------------------------------------------------------
7523if __name__ == "__main__":
7524
7525    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
7526#    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
7527    suite = unittest.makeSuite(Test_Data_Manager,'test')
7528    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
7529    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_holeII')
7530
7531   
7532    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
7533        Test_Data_Manager.verbose=True
7534        saveout = sys.stdout   
7535        filename = ".temp_verbose"
7536        fid = open(filename, 'w')
7537        sys.stdout = fid
7538    else:
7539        pass
7540    runner = unittest.TextTestRunner() # verbosity=2)
7541    runner.run(suite)
7542
7543    # Cleaning up
7544    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
7545        sys.stdout = saveout
7546        fid.close() 
7547        os.remove(filename)
7548
7549
Note: See TracBrowser for help on using the repository browser.