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

Last change on this file since 4414 was 4414, checked in by duncan, 17 years ago

checking in a set up of testing verbose sections of the code. May not be the final version.

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