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

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

Addition to URS_points_needed_to_file so it can do the Northern hemisphere.

File size: 249.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