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

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

Cosmetics while pondering over storage of extrema

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