source: inundation/pyvolution/test_data_manager.py @ 2131

Last change on this file since 2131 was 2074, checked in by duncan, 19 years ago

Removing file leaks in tests, moving tsh2sww to data_manager
File size: 113.5 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import copy
6from Numeric import zeros, array, allclose, Float
7from util import mean
8import tempfile
9import os
10from Scientific.IO.NetCDF import NetCDFFile
11
12from data_manager import *
13from shallow_water import *
14from config import epsilon
15import data_manager
16
17from coordinate_transforms.geo_reference import Geo_reference
18
19class Test_Data_Manager(unittest.TestCase):
20    def setUp(self):
21        import time
22        from mesh_factory import rectangular
23
24        #Create basic mesh
25        points, vertices, boundary = rectangular(2, 2)
26
27        #Create shallow water domain
28        domain = Domain(points, vertices, boundary)
29        domain.default_order=2
30
31
32        #Set some field values
33        domain.set_quantity('elevation', lambda x,y: -x)
34        domain.set_quantity('friction', 0.03)
35
36
37        ######################
38        # Boundary conditions
39        B = Transmissive_boundary(domain)
40        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
41
42
43        ######################
44        #Initial condition - with jumps
45
46
47        bed = domain.quantities['elevation'].vertex_values
48        stage = zeros(bed.shape, Float)
49
50        h = 0.3
51        for i in range(stage.shape[0]):
52            if i % 2 == 0:
53                stage[i,:] = bed[i,:] + h
54            else:
55                stage[i,:] = bed[i,:]
56
57        domain.set_quantity('stage', stage)
58        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
59
60        domain.distribute_to_vertices_and_edges()
61
62
63        self.domain = domain
64
65        C = domain.get_vertex_coordinates()
66        self.X = C[:,0:6:2].copy()
67        self.Y = C[:,1:6:2].copy()
68
69        self.F = bed
70
71
72    def tearDown(self):
73        pass
74
75    def test_sww_constant(self):
76        """Test that constant sww information can be written correctly
77        (non smooth)
78        """
79
80        import time, os
81        from Numeric import array, zeros, allclose, Float, concatenate
82        from Scientific.IO.NetCDF import NetCDFFile
83
84        self.domain.filename = 'datatest' + str(id(self))
85        self.domain.format = 'sww'
86        self.domain.smooth = False
87
88        sww = get_dataobject(self.domain)
89        sww.store_connectivity()
90
91        #Check contents
92        #Get NetCDF
93        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
94
95        # Get the variables
96        x = fid.variables['x']
97        y = fid.variables['y']
98        z = fid.variables['elevation']
99
100        volumes = fid.variables['volumes']
101
102
103        assert allclose (x[:], self.X.flat)
104        assert allclose (y[:], self.Y.flat)
105        assert allclose (z[:], self.F.flat)
106
107        V = volumes
108
109        P = len(self.domain)
110        for k in range(P):
111            assert V[k, 0] == 3*k
112            assert V[k, 1] == 3*k+1
113            assert V[k, 2] == 3*k+2
114
115
116        fid.close()
117
118        #Cleanup
119        os.remove(sww.filename)
120
121
122    def test_sww_constant_smooth(self):
123        """Test that constant sww information can be written correctly
124        (non smooth)
125        """
126
127        import time, os
128        from Numeric import array, zeros, allclose, Float, concatenate
129        from Scientific.IO.NetCDF import NetCDFFile
130
131        self.domain.filename = 'datatest' + str(id(self))
132        self.domain.format = 'sww'
133        self.domain.smooth = True
134
135        sww = get_dataobject(self.domain)
136        sww.store_connectivity()
137
138        #Check contents
139        #Get NetCDF
140        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
141
142        # Get the variables
143        x = fid.variables['x']
144        y = fid.variables['y']
145        z = fid.variables['elevation']
146
147        volumes = fid.variables['volumes']
148
149        X = x[:]
150        Y = y[:]
151
152        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
153        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
154        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
155
156        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
157
158        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
159
160        Z = z[:]
161        assert Z[4] == -0.5
162
163        V = volumes
164        assert V[2,0] == 4
165        assert V[2,1] == 5
166        assert V[2,2] == 1
167
168        assert V[4,0] == 6
169        assert V[4,1] == 7
170        assert V[4,2] == 3
171
172
173        fid.close()
174
175        #Cleanup
176        os.remove(sww.filename)
177
178
179
180    def test_sww_variable(self):
181        """Test that sww information can be written correctly
182        """
183
184        import time, os
185        from Numeric import array, zeros, allclose, Float, concatenate
186        from Scientific.IO.NetCDF import NetCDFFile
187
188        self.domain.filename = 'datatest' + str(id(self))
189        self.domain.format = 'sww'
190        self.domain.smooth = True
191        self.domain.reduction = mean
192
193        sww = get_dataobject(self.domain)
194        sww.store_connectivity()
195        sww.store_timestep('stage')
196
197        #Check contents
198        #Get NetCDF
199        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
200
201
202        # Get the variables
203        x = fid.variables['x']
204        y = fid.variables['y']
205        z = fid.variables['elevation']
206        time = fid.variables['time']
207        stage = fid.variables['stage']
208
209
210        Q = self.domain.quantities['stage']
211        Q0 = Q.vertex_values[:,0]
212        Q1 = Q.vertex_values[:,1]
213        Q2 = Q.vertex_values[:,2]
214
215        A = stage[0,:]
216        #print A[0], (Q2[0,0] + Q1[1,0])/2
217        assert allclose(A[0], (Q2[0] + Q1[1])/2)
218        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
219        assert allclose(A[2], Q0[3])
220        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
221
222        #Center point
223        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
224                                 Q0[5] + Q2[6] + Q1[7])/6)
225
226
227
228        fid.close()
229
230        #Cleanup
231        os.remove(sww.filename)
232
233
234    def test_sww_variable2(self):
235        """Test that sww information can be written correctly
236        multiple timesteps. Use average as reduction operator
237        """
238
239        import time, os
240        from Numeric import array, zeros, allclose, Float, concatenate
241        from Scientific.IO.NetCDF import NetCDFFile
242
243        self.domain.filename = 'datatest' + str(id(self))
244        self.domain.format = 'sww'
245        self.domain.smooth = True
246
247        self.domain.reduction = mean
248
249        sww = get_dataobject(self.domain)
250        sww.store_connectivity()
251        sww.store_timestep('stage')
252        self.domain.evolve_to_end(finaltime = 0.01)
253        sww.store_timestep('stage')
254
255
256        #Check contents
257        #Get NetCDF
258        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
259
260        # Get the variables
261        x = fid.variables['x']
262        y = fid.variables['y']
263        z = fid.variables['elevation']
264        time = fid.variables['time']
265        stage = fid.variables['stage']
266
267        #Check values
268        Q = self.domain.quantities['stage']
269        Q0 = Q.vertex_values[:,0]
270        Q1 = Q.vertex_values[:,1]
271        Q2 = Q.vertex_values[:,2]
272
273        A = stage[1,:]
274        assert allclose(A[0], (Q2[0] + Q1[1])/2)
275        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
276        assert allclose(A[2], Q0[3])
277        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
278
279        #Center point
280        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
281                                 Q0[5] + Q2[6] + Q1[7])/6)
282
283
284        fid.close()
285
286        #Cleanup
287        os.remove(sww.filename)
288
289    def test_sww_variable3(self):
290        """Test that sww information can be written correctly
291        multiple timesteps using a different reduction operator (min)
292        """
293
294        import time, os
295        from Numeric import array, zeros, allclose, Float, concatenate
296        from Scientific.IO.NetCDF import NetCDFFile
297
298        self.domain.filename = 'datatest' + str(id(self))
299        self.domain.format = 'sww'
300        self.domain.smooth = True
301        self.domain.reduction = min
302
303        sww = get_dataobject(self.domain)
304        sww.store_connectivity()
305        sww.store_timestep('stage')
306
307        self.domain.evolve_to_end(finaltime = 0.01)
308        sww.store_timestep('stage')
309
310
311        #Check contents
312        #Get NetCDF
313        fid = NetCDFFile(sww.filename, 'r')
314
315
316        # Get the variables
317        x = fid.variables['x']
318        y = fid.variables['y']
319        z = fid.variables['elevation']
320        time = fid.variables['time']
321        stage = fid.variables['stage']
322
323        #Check values
324        Q = self.domain.quantities['stage']
325        Q0 = Q.vertex_values[:,0]
326        Q1 = Q.vertex_values[:,1]
327        Q2 = Q.vertex_values[:,2]
328
329        A = stage[1,:]
330        assert allclose(A[0], min(Q2[0], Q1[1]))
331        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
332        assert allclose(A[2], Q0[3])
333        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
334
335        #Center point
336        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
337                                  Q0[5], Q2[6], Q1[7]))
338
339
340        fid.close()
341
342        #Cleanup
343        os.remove(sww.filename)
344
345
346    def test_sync(self):
347        """Test info stored at each timestep is as expected (incl initial condition)
348        """
349
350        import time, os, config
351        from Numeric import array, zeros, allclose, Float, concatenate
352        from Scientific.IO.NetCDF import NetCDFFile
353
354        self.domain.filename = 'synctest'
355        self.domain.format = 'sww'
356        self.domain.smooth = False
357        self.domain.store = True
358        self.domain.beta_h = 0
359
360        #Evolution
361        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
362            stage = self.domain.quantities['stage'].vertex_values
363
364            #Get NetCDF
365            fid = NetCDFFile(self.domain.writer.filename, 'r')
366            stage_file = fid.variables['stage']
367
368            if t == 0.0:
369                assert allclose(stage, self.initial_stage)
370                assert allclose(stage_file[:], stage.flat)
371            else:
372                assert not allclose(stage, self.initial_stage)
373                assert not allclose(stage_file[:], stage.flat)
374
375            fid.close()
376
377        os.remove(self.domain.writer.filename)
378
379
380
381    def test_sww_DSG(self):
382        """Not a test, rather a look at the sww format
383        """
384
385        import time, os
386        from Numeric import array, zeros, allclose, Float, concatenate
387        from Scientific.IO.NetCDF import NetCDFFile
388
389        self.domain.filename = 'datatest' + str(id(self))
390        self.domain.format = 'sww'
391        self.domain.smooth = True
392        self.domain.reduction = mean
393
394        sww = get_dataobject(self.domain)
395        sww.store_connectivity()
396        sww.store_timestep('stage')
397
398        #Check contents
399        #Get NetCDF
400        fid = NetCDFFile(sww.filename, 'r')
401
402        # Get the variables
403        x = fid.variables['x']
404        y = fid.variables['y']
405        z = fid.variables['elevation']
406
407        volumes = fid.variables['volumes']
408        time = fid.variables['time']
409
410        # 2D
411        stage = fid.variables['stage']
412
413        X = x[:]
414        Y = y[:]
415        Z = z[:]
416        V = volumes[:]
417        T = time[:]
418        S = stage[:,:]
419
420#         print "****************************"
421#         print "X ",X
422#         print "****************************"
423#         print "Y ",Y
424#         print "****************************"
425#         print "Z ",Z
426#         print "****************************"
427#         print "V ",V
428#         print "****************************"
429#         print "Time ",T
430#         print "****************************"
431#         print "Stage ",S
432#         print "****************************"
433
434
435        fid.close()
436
437        #Cleanup
438        os.remove(sww.filename)
439
440
441    def test_write_pts(self):
442        #Get (enough) datapoints
443
444        from Numeric import array
445        points = array([[ 0.66666667, 0.66666667],
446                        [ 1.33333333, 1.33333333],
447                        [ 2.66666667, 0.66666667],
448                        [ 0.66666667, 2.66666667],
449                        [ 0.0, 1.0],
450                        [ 0.0, 3.0],
451                        [ 1.0, 0.0],
452                        [ 1.0, 1.0],
453                        [ 1.0, 2.0],
454                        [ 1.0, 3.0],
455                        [ 2.0, 1.0],
456                        [ 3.0, 0.0],
457                        [ 3.0, 1.0]])
458
459        z = points[:,0] + 2*points[:,1]
460
461        ptsfile = 'testptsfile.pts'
462        write_ptsfile(ptsfile, points, z,
463                      attribute_name = 'linear_combination')
464
465        #Check contents
466        #Get NetCDF
467        from Scientific.IO.NetCDF import NetCDFFile       
468        fid = NetCDFFile(ptsfile, 'r')
469
470        # Get the variables
471        #print fid.variables.keys()
472        points1 = fid.variables['points']
473        z1 = fid.variables['linear_combination']
474
475        #Check values
476
477        #print points[:]
478        #print ref_points
479        assert allclose(points, points1)
480
481        #print attributes[:]
482        #print ref_elevation
483        assert allclose(z, z1)
484
485        #Cleanup
486        fid.close()
487
488        import os
489        os.remove(ptsfile)
490       
491
492
493
494    def test_dem2pts(self):
495        """Test conversion from dem in ascii format to native NetCDF xya format
496        """
497
498        import time, os
499        from Numeric import array, zeros, allclose, Float, concatenate
500        from Scientific.IO.NetCDF import NetCDFFile
501
502        #Write test asc file
503        root = 'demtest'
504
505        filename = root+'.asc'
506        fid = open(filename, 'w')
507        fid.write("""ncols         5
508nrows         6
509xllcorner     2000.5
510yllcorner     3000.5
511cellsize      25
512NODATA_value  -9999
513""")
514        #Create linear function
515
516        ref_points = []
517        ref_elevation = []
518        for i in range(6):
519            y = (6-i)*25.0
520            for j in range(5):
521                x = j*25.0
522                z = x+2*y
523
524                ref_points.append( [x,y] )
525                ref_elevation.append(z)
526                fid.write('%f ' %z)
527            fid.write('\n')
528
529        fid.close()
530
531        #Write prj file with metadata
532        metafilename = root+'.prj'
533        fid = open(metafilename, 'w')
534
535
536        fid.write("""Projection UTM
537Zone 56
538Datum WGS84
539Zunits NO
540Units METERS
541Spheroid WGS84
542Xshift 0.0000000000
543Yshift 10000000.0000000000
544Parameters
545""")
546        fid.close()
547
548        #Convert to NetCDF pts
549        convert_dem_from_ascii2netcdf(root)
550        dem2pts(root)
551
552        #Check contents
553        #Get NetCDF
554        fid = NetCDFFile(root+'.pts', 'r')
555
556        # Get the variables
557        #print fid.variables.keys()
558        points = fid.variables['points']
559        elevation = fid.variables['elevation']
560
561        #Check values
562
563        #print points[:]
564        #print ref_points
565        assert allclose(points, ref_points)
566
567        #print attributes[:]
568        #print ref_elevation
569        assert allclose(elevation, ref_elevation)
570
571        #Cleanup
572        fid.close()
573
574
575        os.remove(root + '.pts')
576        os.remove(root + '.dem')
577        os.remove(root + '.asc')
578        os.remove(root + '.prj')
579
580
581
582    def test_dem2pts_bounding_box(self):
583        """Test conversion from dem in ascii format to native NetCDF xya format
584        """
585
586        import time, os
587        from Numeric import array, zeros, allclose, Float, concatenate
588        from Scientific.IO.NetCDF import NetCDFFile
589
590        #Write test asc file
591        root = 'demtest'
592
593        filename = root+'.asc'
594        fid = open(filename, 'w')
595        fid.write("""ncols         5
596nrows         6
597xllcorner     2000.5
598yllcorner     3000.5
599cellsize      25
600NODATA_value  -9999
601""")
602        #Create linear function
603
604        ref_points = []
605        ref_elevation = []
606        for i in range(6):
607            y = (6-i)*25.0
608            for j in range(5):
609                x = j*25.0
610                z = x+2*y
611
612                ref_points.append( [x,y] )
613                ref_elevation.append(z)
614                fid.write('%f ' %z)
615            fid.write('\n')
616
617        fid.close()
618
619        #Write prj file with metadata
620        metafilename = root+'.prj'
621        fid = open(metafilename, 'w')
622
623
624        fid.write("""Projection UTM
625Zone 56
626Datum WGS84
627Zunits NO
628Units METERS
629Spheroid WGS84
630Xshift 0.0000000000
631Yshift 10000000.0000000000
632Parameters
633""")
634        fid.close()
635
636        #Convert to NetCDF pts
637        convert_dem_from_ascii2netcdf(root)
638        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
639                northing_min=3035.0, northing_max=3125.5)
640
641        #Check contents
642        #Get NetCDF
643        fid = NetCDFFile(root+'.pts', 'r')
644
645        # Get the variables
646        #print fid.variables.keys()
647        points = fid.variables['points']
648        elevation = fid.variables['elevation']
649
650        #Check values
651        assert fid.xllcorner[0] == 2010.0
652        assert fid.yllcorner[0] == 3035.0
653
654        #create new reference points
655        ref_points = []
656        ref_elevation = []
657        for i in range(4):
658            y = (4-i)*25.0 + 25.0
659            y_new = y + 3000.5 - 3035.0
660            for j in range(4):
661                x = j*25.0 + 25.0
662                x_new = x + 2000.5 - 2010.0
663                z = x+2*y
664
665                ref_points.append( [x_new,y_new] )
666                ref_elevation.append(z)
667
668        #print points[:]
669        #print ref_points
670        assert allclose(points, ref_points)
671
672        #print attributes[:]
673        #print ref_elevation
674        assert allclose(elevation, ref_elevation)
675
676        #Cleanup
677        fid.close()
678
679
680        os.remove(root + '.pts')
681        os.remove(root + '.dem')
682        os.remove(root + '.asc')
683        os.remove(root + '.prj')
684
685
686
687    def test_hecras_cross_sections2pts(self):
688        """Test conversion from HECRAS cross sections in ascii format
689        to native NetCDF pts format
690        """
691
692        import time, os
693        from Numeric import array, zeros, allclose, Float, concatenate
694        from Scientific.IO.NetCDF import NetCDFFile
695
696        #Write test asc file
697        root = 'hecrastest'
698
699        filename = root+'.sdf'
700        fid = open(filename, 'w')
701        fid.write("""
702# RAS export file created on Mon 15Aug2005 11:42
703# by HEC-RAS Version 3.1.1
704
705BEGIN HEADER:
706  UNITS: METRIC
707  DTM TYPE: TIN
708  DTM: v:\1\cit\perth_topo\river_tin
709  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
710  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
711  MAP PROJECTION: UTM
712  PROJECTION ZONE: 50
713  DATUM: AGD66
714  VERTICAL DATUM:
715  NUMBER OF REACHES:  19
716  NUMBER OF CROSS-SECTIONS:  2
717END HEADER:
718
719
720BEGIN CROSS-SECTIONS:
721
722  CROSS-SECTION:
723    STREAM ID:Southern-Wungong
724    REACH ID:Southern-Wungong
725    STATION:21410   
726    CUT LINE:
727      407546.08 , 6437277.542
728      407329.32 , 6437489.482
729      407283.11 , 6437541.232
730    SURFACE LINE:
731     407546.08,   6437277.54,   52.14
732     407538.88,   6437284.58,   51.07
733     407531.68,   6437291.62,   50.56
734     407524.48,   6437298.66,   49.58
735     407517.28,   6437305.70,   49.09
736     407510.08,   6437312.74,   48.76
737  END:
738
739  CROSS-SECTION:
740    STREAM ID:Swan River     
741    REACH ID:Swan Mouth     
742    STATION:840.*   
743    CUT LINE:
744      381178.0855 , 6452559.0685
745      380485.4755 , 6453169.272
746    SURFACE LINE:
747     381178.09,   6452559.07,   4.17
748     381169.49,   6452566.64,   4.26
749     381157.78,   6452576.96,   4.34
750     381155.97,   6452578.56,   4.35
751     381143.72,   6452589.35,   4.43
752     381136.69,   6452595.54,   4.58
753     381114.74,   6452614.88,   4.41
754     381075.53,   6452649.43,   4.17
755     381071.47,   6452653.00,   3.99
756     381063.46,   6452660.06,   3.67
757     381054.41,   6452668.03,   3.67
758  END: 
759END CROSS-SECTIONS:
760""")
761
762        fid.close()
763
764
765        #Convert to NetCDF pts
766        hecras_cross_sections2pts(root)
767
768        #Check contents
769        #Get NetCDF
770        fid = NetCDFFile(root+'.pts', 'r')
771
772        # Get the variables
773        #print fid.variables.keys()
774        points = fid.variables['points']
775        elevation = fid.variables['elevation']
776
777        #Check values
778        ref_points = [[407546.08, 6437277.54],
779                      [407538.88, 6437284.58],
780                      [407531.68, 6437291.62],
781                      [407524.48, 6437298.66],
782                      [407517.28, 6437305.70],
783                      [407510.08, 6437312.74]]
784
785        ref_points += [[381178.09, 6452559.07], 
786                       [381169.49, 6452566.64], 
787                       [381157.78, 6452576.96], 
788                       [381155.97, 6452578.56], 
789                       [381143.72, 6452589.35], 
790                       [381136.69, 6452595.54], 
791                       [381114.74, 6452614.88], 
792                       [381075.53, 6452649.43], 
793                       [381071.47, 6452653.00], 
794                       [381063.46, 6452660.06], 
795                       [381054.41, 6452668.03]] 
796       
797                     
798        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
799        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
800
801        #print points[:]
802        #print ref_points
803        assert allclose(points, ref_points)
804
805        #print attributes[:]
806        #print ref_elevation
807        assert allclose(elevation, ref_elevation)
808
809        #Cleanup
810        fid.close()
811
812
813        os.remove(root + '.sdf')
814        os.remove(root + '.pts')
815
816
817
818    def test_sww2dem_asc_elevation(self):
819        """Test that sww information can be converted correctly to asc/prj
820        format readable by e.g. ArcView
821        """
822
823        import time, os
824        from Numeric import array, zeros, allclose, Float, concatenate
825        from Scientific.IO.NetCDF import NetCDFFile
826
827        #Setup
828        self.domain.filename = 'datatest'
829
830        prjfile = self.domain.filename + '_elevation.prj'
831        ascfile = self.domain.filename + '_elevation.asc'
832        swwfile = self.domain.filename + '.sww'
833
834        self.domain.set_datadir('.')
835        self.domain.format = 'sww'
836        self.domain.smooth = True
837        self.domain.set_quantity('elevation', lambda x,y: -x-y)
838
839        self.domain.geo_reference = Geo_reference(56,308500,6189000)
840
841        sww = get_dataobject(self.domain)
842        sww.store_connectivity()
843        sww.store_timestep('stage')
844
845        self.domain.evolve_to_end(finaltime = 0.01)
846        sww.store_timestep('stage')
847
848        cellsize = 0.25
849        #Check contents
850        #Get NetCDF
851
852        fid = NetCDFFile(sww.filename, 'r')
853
854        # Get the variables
855        x = fid.variables['x'][:]
856        y = fid.variables['y'][:]
857        z = fid.variables['elevation'][:]
858        time = fid.variables['time'][:]
859        stage = fid.variables['stage'][:]
860
861
862        #Export to ascii/prj files
863        sww2dem(self.domain.filename,
864                quantity = 'elevation',
865                cellsize = cellsize,
866                verbose = False,
867                format = 'asc')
868               
869        #Check prj (meta data)
870        prjid = open(prjfile)
871        lines = prjid.readlines()
872        prjid.close()
873
874        L = lines[0].strip().split()
875        assert L[0].strip().lower() == 'projection'
876        assert L[1].strip().lower() == 'utm'
877
878        L = lines[1].strip().split()
879        assert L[0].strip().lower() == 'zone'
880        assert L[1].strip().lower() == '56'
881
882        L = lines[2].strip().split()
883        assert L[0].strip().lower() == 'datum'
884        assert L[1].strip().lower() == 'wgs84'
885
886        L = lines[3].strip().split()
887        assert L[0].strip().lower() == 'zunits'
888        assert L[1].strip().lower() == 'no'
889
890        L = lines[4].strip().split()
891        assert L[0].strip().lower() == 'units'
892        assert L[1].strip().lower() == 'meters'
893
894        L = lines[5].strip().split()
895        assert L[0].strip().lower() == 'spheroid'
896        assert L[1].strip().lower() == 'wgs84'
897
898        L = lines[6].strip().split()
899        assert L[0].strip().lower() == 'xshift'
900        assert L[1].strip().lower() == '500000'
901
902        L = lines[7].strip().split()
903        assert L[0].strip().lower() == 'yshift'
904        assert L[1].strip().lower() == '10000000'
905
906        L = lines[8].strip().split()
907        assert L[0].strip().lower() == 'parameters'
908
909
910        #Check asc file
911        ascid = open(ascfile)
912        lines = ascid.readlines()
913        ascid.close()
914
915        L = lines[0].strip().split()
916        assert L[0].strip().lower() == 'ncols'
917        assert L[1].strip().lower() == '5'
918
919        L = lines[1].strip().split()
920        assert L[0].strip().lower() == 'nrows'
921        assert L[1].strip().lower() == '5'
922
923        L = lines[2].strip().split()
924        assert L[0].strip().lower() == 'xllcorner'
925        assert allclose(float(L[1].strip().lower()), 308500)
926
927        L = lines[3].strip().split()
928        assert L[0].strip().lower() == 'yllcorner'
929        assert allclose(float(L[1].strip().lower()), 6189000)
930
931        L = lines[4].strip().split()
932        assert L[0].strip().lower() == 'cellsize'
933        assert allclose(float(L[1].strip().lower()), cellsize)
934
935        L = lines[5].strip().split()
936        assert L[0].strip() == 'NODATA_value'
937        assert L[1].strip().lower() == '-9999'
938
939        #Check grid values
940        for j in range(5):
941            L = lines[6+j].strip().split()
942            y = (4-j) * cellsize
943            for i in range(5):
944                assert allclose(float(L[i]), -i*cellsize - y)
945
946
947        fid.close()
948
949        #Cleanup
950        os.remove(prjfile)
951        os.remove(ascfile)
952        os.remove(swwfile)
953
954
955
956    def test_sww2dem_larger(self):
957        """Test that sww information can be converted correctly to asc/prj
958        format readable by e.g. ArcView. Here:
959
960        ncols         11
961        nrows         11
962        xllcorner     308500
963        yllcorner     6189000
964        cellsize      10.000000
965        NODATA_value  -9999
966        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
967         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
968         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
969         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
970         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
971         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
972         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
973         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
974         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
975         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
976           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
977       
978        """
979
980        import time, os
981        from Numeric import array, zeros, allclose, Float, concatenate
982        from Scientific.IO.NetCDF import NetCDFFile
983
984        #Setup
985
986        from mesh_factory import rectangular
987
988        #Create basic mesh (100m x 100m)
989        points, vertices, boundary = rectangular(2, 2, 100, 100)
990
991        #Create shallow water domain
992        domain = Domain(points, vertices, boundary)
993        domain.default_order = 2
994
995        domain.filename = 'datatest'
996
997        prjfile = domain.filename + '_elevation.prj'
998        ascfile = domain.filename + '_elevation.asc'
999        swwfile = domain.filename + '.sww'
1000
1001        domain.set_datadir('.')
1002        domain.format = 'sww'
1003        domain.smooth = True
1004        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1005
1006        #
1007        domain.set_quantity('elevation', lambda x,y: -x-y)
1008        domain.set_quantity('stage', 0)       
1009
1010        B = Transmissive_boundary(domain)
1011        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1012
1013
1014        #
1015        sww = get_dataobject(domain)
1016        sww.store_connectivity()
1017        sww.store_timestep('stage')
1018
1019        domain.evolve_to_end(finaltime = 0.01)
1020        sww.store_timestep('stage')
1021
1022        cellsize = 10  #10m grid
1023
1024       
1025        #Check contents
1026        #Get NetCDF
1027
1028        fid = NetCDFFile(sww.filename, 'r')
1029
1030        # Get the variables
1031        x = fid.variables['x'][:]
1032        y = fid.variables['y'][:]
1033        z = fid.variables['elevation'][:]
1034        time = fid.variables['time'][:]
1035        stage = fid.variables['stage'][:]
1036
1037
1038        #Export to ascii/prj files
1039        sww2dem(domain.filename,
1040                quantity = 'elevation',
1041                cellsize = cellsize,
1042                verbose = False,
1043                format = 'asc')
1044       
1045               
1046        #Check prj (meta data)
1047        prjid = open(prjfile)
1048        lines = prjid.readlines()
1049        prjid.close()
1050
1051        L = lines[0].strip().split()
1052        assert L[0].strip().lower() == 'projection'
1053        assert L[1].strip().lower() == 'utm'
1054
1055        L = lines[1].strip().split()
1056        assert L[0].strip().lower() == 'zone'
1057        assert L[1].strip().lower() == '56'
1058
1059        L = lines[2].strip().split()
1060        assert L[0].strip().lower() == 'datum'
1061        assert L[1].strip().lower() == 'wgs84'
1062
1063        L = lines[3].strip().split()
1064        assert L[0].strip().lower() == 'zunits'
1065        assert L[1].strip().lower() == 'no'
1066
1067        L = lines[4].strip().split()
1068        assert L[0].strip().lower() == 'units'
1069        assert L[1].strip().lower() == 'meters'
1070
1071        L = lines[5].strip().split()
1072        assert L[0].strip().lower() == 'spheroid'
1073        assert L[1].strip().lower() == 'wgs84'
1074
1075        L = lines[6].strip().split()
1076        assert L[0].strip().lower() == 'xshift'
1077        assert L[1].strip().lower() == '500000'
1078
1079        L = lines[7].strip().split()
1080        assert L[0].strip().lower() == 'yshift'
1081        assert L[1].strip().lower() == '10000000'
1082
1083        L = lines[8].strip().split()
1084        assert L[0].strip().lower() == 'parameters'
1085
1086
1087        #Check asc file
1088        ascid = open(ascfile)
1089        lines = ascid.readlines()
1090        ascid.close()
1091
1092        L = lines[0].strip().split()
1093        assert L[0].strip().lower() == 'ncols'
1094        assert L[1].strip().lower() == '11'
1095
1096        L = lines[1].strip().split()
1097        assert L[0].strip().lower() == 'nrows'
1098        assert L[1].strip().lower() == '11'
1099
1100        L = lines[2].strip().split()
1101        assert L[0].strip().lower() == 'xllcorner'
1102        assert allclose(float(L[1].strip().lower()), 308500)
1103
1104        L = lines[3].strip().split()
1105        assert L[0].strip().lower() == 'yllcorner'
1106        assert allclose(float(L[1].strip().lower()), 6189000)
1107
1108        L = lines[4].strip().split()
1109        assert L[0].strip().lower() == 'cellsize'
1110        assert allclose(float(L[1].strip().lower()), cellsize)
1111
1112        L = lines[5].strip().split()
1113        assert L[0].strip() == 'NODATA_value'
1114        assert L[1].strip().lower() == '-9999'
1115
1116        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1117        for i, line in enumerate(lines[6:]):
1118            for j, value in enumerate( line.split() ):
1119                #assert allclose(float(value), -(10-i+j)*cellsize)
1120                assert float(value) == -(10-i+j)*cellsize
1121       
1122
1123        fid.close()
1124
1125        #Cleanup
1126        os.remove(prjfile)
1127        os.remove(ascfile)
1128        os.remove(swwfile)
1129
1130
1131
1132    def test_sww2dem_boundingbox(self):
1133        """Test that sww information can be converted correctly to asc/prj
1134        format readable by e.g. ArcView.
1135        This will test that mesh can be restricted by bounding box
1136
1137        Original extent is 100m x 100m:
1138
1139        Eastings:   308500 -  308600
1140        Northings: 6189000 - 6189100
1141
1142        Bounding box changes this to the 50m x 50m square defined by
1143
1144        Eastings:   308530 -  308570
1145        Northings: 6189050 - 6189100
1146
1147        The cropped values should be
1148
1149         -130 -140 -150 -160 -170 
1150         -120 -130 -140 -150 -160 
1151         -110 -120 -130 -140 -150 
1152         -100 -110 -120 -130 -140 
1153          -90 -100 -110 -120 -130 
1154          -80  -90 -100 -110 -120
1155
1156        and the new lower reference point should be   
1157        Eastings:   308530
1158        Northings: 6189050
1159       
1160        Original dataset is the same as in test_sww2dem_larger()
1161       
1162        """
1163
1164        import time, os
1165        from Numeric import array, zeros, allclose, Float, concatenate
1166        from Scientific.IO.NetCDF import NetCDFFile
1167
1168        #Setup
1169
1170        from mesh_factory import rectangular
1171
1172        #Create basic mesh (100m x 100m)
1173        points, vertices, boundary = rectangular(2, 2, 100, 100)
1174
1175        #Create shallow water domain
1176        domain = Domain(points, vertices, boundary)
1177        domain.default_order = 2
1178
1179        domain.filename = 'datatest'
1180
1181        prjfile = domain.filename + '_elevation.prj'
1182        ascfile = domain.filename + '_elevation.asc'
1183        swwfile = domain.filename + '.sww'
1184
1185        domain.set_datadir('.')
1186        domain.format = 'sww'
1187        domain.smooth = True
1188        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1189
1190        #
1191        domain.set_quantity('elevation', lambda x,y: -x-y)
1192        domain.set_quantity('stage', 0)       
1193
1194        B = Transmissive_boundary(domain)
1195        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1196
1197
1198        #
1199        sww = get_dataobject(domain)
1200        sww.store_connectivity()
1201        sww.store_timestep('stage')
1202
1203        domain.evolve_to_end(finaltime = 0.01)
1204        sww.store_timestep('stage')
1205
1206        cellsize = 10  #10m grid
1207
1208       
1209        #Check contents
1210        #Get NetCDF
1211
1212        fid = NetCDFFile(sww.filename, 'r')
1213
1214        # Get the variables
1215        x = fid.variables['x'][:]
1216        y = fid.variables['y'][:]
1217        z = fid.variables['elevation'][:]
1218        time = fid.variables['time'][:]
1219        stage = fid.variables['stage'][:]
1220
1221
1222        #Export to ascii/prj files
1223        sww2dem(domain.filename,
1224                quantity = 'elevation',
1225                cellsize = cellsize,
1226                easting_min = 308530,
1227                easting_max = 308570,
1228                northing_min = 6189050, 
1229                northing_max = 6189100,                               
1230                verbose = False,
1231                format = 'asc')
1232       
1233        fid.close()
1234
1235       
1236        #Check prj (meta data)
1237        prjid = open(prjfile)
1238        lines = prjid.readlines()
1239        prjid.close()
1240
1241        L = lines[0].strip().split()
1242        assert L[0].strip().lower() == 'projection'
1243        assert L[1].strip().lower() == 'utm'
1244
1245        L = lines[1].strip().split()
1246        assert L[0].strip().lower() == 'zone'
1247        assert L[1].strip().lower() == '56'
1248
1249        L = lines[2].strip().split()
1250        assert L[0].strip().lower() == 'datum'
1251        assert L[1].strip().lower() == 'wgs84'
1252
1253        L = lines[3].strip().split()
1254        assert L[0].strip().lower() == 'zunits'
1255        assert L[1].strip().lower() == 'no'
1256
1257        L = lines[4].strip().split()
1258        assert L[0].strip().lower() == 'units'
1259        assert L[1].strip().lower() == 'meters'
1260
1261        L = lines[5].strip().split()
1262        assert L[0].strip().lower() == 'spheroid'
1263        assert L[1].strip().lower() == 'wgs84'
1264
1265        L = lines[6].strip().split()
1266        assert L[0].strip().lower() == 'xshift'
1267        assert L[1].strip().lower() == '500000'
1268
1269        L = lines[7].strip().split()
1270        assert L[0].strip().lower() == 'yshift'
1271        assert L[1].strip().lower() == '10000000'
1272
1273        L = lines[8].strip().split()
1274        assert L[0].strip().lower() == 'parameters'
1275
1276
1277        #Check asc file
1278        ascid = open(ascfile)
1279        lines = ascid.readlines()
1280        ascid.close()
1281
1282        L = lines[0].strip().split()
1283        assert L[0].strip().lower() == 'ncols'
1284        assert L[1].strip().lower() == '5'
1285
1286        L = lines[1].strip().split()
1287        assert L[0].strip().lower() == 'nrows'
1288        assert L[1].strip().lower() == '6'
1289
1290        L = lines[2].strip().split()
1291        assert L[0].strip().lower() == 'xllcorner'
1292        assert allclose(float(L[1].strip().lower()), 308530)
1293
1294        L = lines[3].strip().split()
1295        assert L[0].strip().lower() == 'yllcorner'
1296        assert allclose(float(L[1].strip().lower()), 6189050)
1297
1298        L = lines[4].strip().split()
1299        assert L[0].strip().lower() == 'cellsize'
1300        assert allclose(float(L[1].strip().lower()), cellsize)
1301
1302        L = lines[5].strip().split()
1303        assert L[0].strip() == 'NODATA_value'
1304        assert L[1].strip().lower() == '-9999'
1305
1306        #Check grid values
1307        for i, line in enumerate(lines[6:]):
1308            for j, value in enumerate( line.split() ):
1309                #assert float(value) == -(10-i+j)*cellsize               
1310                assert float(value) == -(10-i+j+3)*cellsize
1311       
1312
1313
1314        #Cleanup
1315        os.remove(prjfile)
1316        os.remove(ascfile)
1317        os.remove(swwfile)
1318
1319
1320
1321    def test_sww2dem_asc_stage_reduction(self):
1322        """Test that sww information can be converted correctly to asc/prj
1323        format readable by e.g. ArcView
1324
1325        This tests the reduction of quantity stage using min
1326        """
1327
1328        import time, os
1329        from Numeric import array, zeros, allclose, Float, concatenate
1330        from Scientific.IO.NetCDF import NetCDFFile
1331
1332        #Setup
1333        self.domain.filename = 'datatest'
1334
1335        prjfile = self.domain.filename + '_stage.prj'
1336        ascfile = self.domain.filename + '_stage.asc'
1337        swwfile = self.domain.filename + '.sww'
1338
1339        self.domain.set_datadir('.')
1340        self.domain.format = 'sww'
1341        self.domain.smooth = True
1342        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1343
1344        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1345
1346
1347        sww = get_dataobject(self.domain)
1348        sww.store_connectivity()
1349        sww.store_timestep('stage')
1350
1351        self.domain.evolve_to_end(finaltime = 0.01)
1352        sww.store_timestep('stage')
1353
1354        cellsize = 0.25
1355        #Check contents
1356        #Get NetCDF
1357
1358        fid = NetCDFFile(sww.filename, 'r')
1359
1360        # Get the variables
1361        x = fid.variables['x'][:]
1362        y = fid.variables['y'][:]
1363        z = fid.variables['elevation'][:]
1364        time = fid.variables['time'][:]
1365        stage = fid.variables['stage'][:]
1366
1367
1368        #Export to ascii/prj files
1369        sww2dem(self.domain.filename,
1370                quantity = 'stage',
1371                cellsize = cellsize,
1372                reduction = min,
1373                format = 'asc')
1374
1375
1376        #Check asc file
1377        ascid = open(ascfile)
1378        lines = ascid.readlines()
1379        ascid.close()
1380
1381        L = lines[0].strip().split()
1382        assert L[0].strip().lower() == 'ncols'
1383        assert L[1].strip().lower() == '5'
1384
1385        L = lines[1].strip().split()
1386        assert L[0].strip().lower() == 'nrows'
1387        assert L[1].strip().lower() == '5'
1388
1389        L = lines[2].strip().split()
1390        assert L[0].strip().lower() == 'xllcorner'
1391        assert allclose(float(L[1].strip().lower()), 308500)
1392
1393        L = lines[3].strip().split()
1394        assert L[0].strip().lower() == 'yllcorner'
1395        assert allclose(float(L[1].strip().lower()), 6189000)
1396
1397        L = lines[4].strip().split()
1398        assert L[0].strip().lower() == 'cellsize'
1399        assert allclose(float(L[1].strip().lower()), cellsize)
1400
1401        L = lines[5].strip().split()
1402        assert L[0].strip() == 'NODATA_value'
1403        assert L[1].strip().lower() == '-9999'
1404
1405
1406        #Check grid values (where applicable)
1407        for j in range(5):
1408            if j%2 == 0:
1409                L = lines[6+j].strip().split()
1410                jj = 4-j
1411                for i in range(5):
1412                    if i%2 == 0:
1413                        index = jj/2 + i/2*3
1414                        val0 = stage[0,index]
1415                        val1 = stage[1,index]
1416
1417                        #print i, j, index, ':', L[i], val0, val1
1418                        assert allclose(float(L[i]), min(val0, val1))
1419
1420
1421        fid.close()
1422
1423        #Cleanup
1424        os.remove(prjfile)
1425        os.remove(ascfile)
1426        #os.remove(swwfile)
1427
1428
1429
1430    def test_sww2dem_asc_derived_quantity(self):
1431        """Test that sww information can be converted correctly to asc/prj
1432        format readable by e.g. ArcView
1433
1434        This tests the use of derived quantities
1435        """
1436
1437        import time, os
1438        from Numeric import array, zeros, allclose, Float, concatenate
1439        from Scientific.IO.NetCDF import NetCDFFile
1440
1441        #Setup
1442        self.domain.filename = 'datatest'
1443
1444        prjfile = self.domain.filename + '_depth.prj'
1445        ascfile = self.domain.filename + '_depth.asc'
1446        swwfile = self.domain.filename + '.sww'
1447
1448        self.domain.set_datadir('.')
1449        self.domain.format = 'sww'
1450        self.domain.smooth = True
1451        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1452        self.domain.set_quantity('stage', 0.0)
1453
1454        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1455
1456
1457        sww = get_dataobject(self.domain)
1458        sww.store_connectivity()
1459        sww.store_timestep('stage')
1460
1461        self.domain.evolve_to_end(finaltime = 0.01)
1462        sww.store_timestep('stage')
1463
1464        cellsize = 0.25
1465        #Check contents
1466        #Get NetCDF
1467
1468        fid = NetCDFFile(sww.filename, 'r')
1469
1470        # Get the variables
1471        x = fid.variables['x'][:]
1472        y = fid.variables['y'][:]
1473        z = fid.variables['elevation'][:]
1474        time = fid.variables['time'][:]
1475        stage = fid.variables['stage'][:]
1476
1477
1478        #Export to ascii/prj files
1479        sww2dem(self.domain.filename,
1480                basename_out = 'datatest_depth',
1481                quantity = 'stage - elevation',
1482                cellsize = cellsize,
1483                reduction = min,
1484                format = 'asc',
1485                verbose = False)
1486
1487
1488        #Check asc file
1489        ascid = open(ascfile)
1490        lines = ascid.readlines()
1491        ascid.close()
1492
1493        L = lines[0].strip().split()
1494        assert L[0].strip().lower() == 'ncols'
1495        assert L[1].strip().lower() == '5'
1496
1497        L = lines[1].strip().split()
1498        assert L[0].strip().lower() == 'nrows'
1499        assert L[1].strip().lower() == '5'
1500
1501        L = lines[2].strip().split()
1502        assert L[0].strip().lower() == 'xllcorner'
1503        assert allclose(float(L[1].strip().lower()), 308500)
1504
1505        L = lines[3].strip().split()
1506        assert L[0].strip().lower() == 'yllcorner'
1507        assert allclose(float(L[1].strip().lower()), 6189000)
1508
1509        L = lines[4].strip().split()
1510        assert L[0].strip().lower() == 'cellsize'
1511        assert allclose(float(L[1].strip().lower()), cellsize)
1512
1513        L = lines[5].strip().split()
1514        assert L[0].strip() == 'NODATA_value'
1515        assert L[1].strip().lower() == '-9999'
1516
1517
1518        #Check grid values (where applicable)
1519        for j in range(5):
1520            if j%2 == 0:
1521                L = lines[6+j].strip().split()
1522                jj = 4-j
1523                for i in range(5):
1524                    if i%2 == 0:
1525                        index = jj/2 + i/2*3
1526                        val0 = stage[0,index] - z[index]
1527                        val1 = stage[1,index] - z[index]
1528
1529                        #print i, j, index, ':', L[i], val0, val1
1530                        assert allclose(float(L[i]), min(val0, val1))
1531
1532
1533        fid.close()
1534
1535        #Cleanup
1536        os.remove(prjfile)
1537        os.remove(ascfile)
1538        #os.remove(swwfile)
1539
1540
1541
1542
1543
1544    def test_sww2dem_asc_missing_points(self):
1545        """Test that sww information can be converted correctly to asc/prj
1546        format readable by e.g. ArcView
1547
1548        This test includes the writing of missing values
1549        """
1550
1551        import time, os
1552        from Numeric import array, zeros, allclose, Float, concatenate
1553        from Scientific.IO.NetCDF import NetCDFFile
1554
1555        #Setup mesh not coinciding with rectangle.
1556        #This will cause missing values to occur in gridded data
1557
1558
1559        points = [                        [1.0, 1.0],
1560                              [0.5, 0.5], [1.0, 0.5],
1561                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1562
1563        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1564
1565        #Create shallow water domain
1566        domain = Domain(points, vertices)
1567        domain.default_order=2
1568
1569
1570        #Set some field values
1571        domain.set_quantity('elevation', lambda x,y: -x-y)
1572        domain.set_quantity('friction', 0.03)
1573
1574
1575        ######################
1576        # Boundary conditions
1577        B = Transmissive_boundary(domain)
1578        domain.set_boundary( {'exterior': B} )
1579
1580
1581        ######################
1582        #Initial condition - with jumps
1583
1584        bed = domain.quantities['elevation'].vertex_values
1585        stage = zeros(bed.shape, Float)
1586
1587        h = 0.3
1588        for i in range(stage.shape[0]):
1589            if i % 2 == 0:
1590                stage[i,:] = bed[i,:] + h
1591            else:
1592                stage[i,:] = bed[i,:]
1593
1594        domain.set_quantity('stage', stage)
1595        domain.distribute_to_vertices_and_edges()
1596
1597        domain.filename = 'datatest'
1598
1599        prjfile = domain.filename + '_elevation.prj'
1600        ascfile = domain.filename + '_elevation.asc'
1601        swwfile = domain.filename + '.sww'
1602
1603        domain.set_datadir('.')
1604        domain.format = 'sww'
1605        domain.smooth = True
1606
1607        domain.geo_reference = Geo_reference(56,308500,6189000)
1608
1609        sww = get_dataobject(domain)
1610        sww.store_connectivity()
1611        sww.store_timestep('stage')
1612
1613        cellsize = 0.25
1614        #Check contents
1615        #Get NetCDF
1616
1617        fid = NetCDFFile(swwfile, 'r')
1618
1619        # Get the variables
1620        x = fid.variables['x'][:]
1621        y = fid.variables['y'][:]
1622        z = fid.variables['elevation'][:]
1623        time = fid.variables['time'][:]
1624
1625        try:
1626            geo_reference = Geo_reference(NetCDFObject=fid)
1627        except AttributeError, e:
1628            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1629
1630        #Export to ascii/prj files
1631        sww2dem(domain.filename,
1632                quantity = 'elevation',
1633                cellsize = cellsize,
1634                verbose = False,
1635                format = 'asc')
1636
1637
1638        #Check asc file
1639        ascid = open(ascfile)
1640        lines = ascid.readlines()
1641        ascid.close()
1642
1643        L = lines[0].strip().split()
1644        assert L[0].strip().lower() == 'ncols'
1645        assert L[1].strip().lower() == '5'
1646
1647        L = lines[1].strip().split()
1648        assert L[0].strip().lower() == 'nrows'
1649        assert L[1].strip().lower() == '5'
1650
1651        L = lines[2].strip().split()
1652        assert L[0].strip().lower() == 'xllcorner'
1653        assert allclose(float(L[1].strip().lower()), 308500)
1654
1655        L = lines[3].strip().split()
1656        assert L[0].strip().lower() == 'yllcorner'
1657        assert allclose(float(L[1].strip().lower()), 6189000)
1658
1659        L = lines[4].strip().split()
1660        assert L[0].strip().lower() == 'cellsize'
1661        assert allclose(float(L[1].strip().lower()), cellsize)
1662
1663        L = lines[5].strip().split()
1664        assert L[0].strip() == 'NODATA_value'
1665        assert L[1].strip().lower() == '-9999'
1666
1667        #Check grid values
1668        for j in range(5):
1669            L = lines[6+j].strip().split()
1670            assert len(L) == 5
1671            y = (4-j) * cellsize
1672
1673            for i in range(5):
1674                #print i
1675                if i+j >= 4:
1676                    assert allclose(float(L[i]), -i*cellsize - y)
1677                else:
1678                    #Missing values
1679                    assert allclose(float(L[i]), -9999)
1680
1681
1682
1683        fid.close()
1684
1685        #Cleanup
1686        os.remove(prjfile)
1687        os.remove(ascfile)
1688        os.remove(swwfile)
1689
1690    def test_sww2ers_simple(self):
1691        """Test that sww information can be converted correctly to asc/prj
1692        format readable by e.g. ArcView
1693        """
1694
1695        import time, os
1696        from Numeric import array, zeros, allclose, Float, concatenate
1697        from Scientific.IO.NetCDF import NetCDFFile
1698
1699
1700        NODATA_value = 1758323
1701
1702        #Setup
1703        self.domain.filename = 'datatest'
1704
1705        headerfile = self.domain.filename + '.ers'
1706        swwfile = self.domain.filename + '.sww'
1707
1708        self.domain.set_datadir('.')
1709        self.domain.format = 'sww'
1710        self.domain.smooth = True
1711        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1712
1713        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1714
1715        sww = get_dataobject(self.domain)
1716        sww.store_connectivity()
1717        sww.store_timestep('stage')
1718
1719        self.domain.evolve_to_end(finaltime = 0.01)
1720        sww.store_timestep('stage')
1721
1722        cellsize = 0.25
1723        #Check contents
1724        #Get NetCDF
1725
1726        fid = NetCDFFile(sww.filename, 'r')
1727
1728        # Get the variables
1729        x = fid.variables['x'][:]
1730        y = fid.variables['y'][:]
1731        z = fid.variables['elevation'][:]
1732        time = fid.variables['time'][:]
1733        stage = fid.variables['stage'][:]
1734
1735
1736        #Export to ers files
1737        #sww2ers(self.domain.filename,
1738        #        quantity = 'elevation',
1739        #        cellsize = cellsize,
1740        #        verbose = False)
1741               
1742        sww2dem(self.domain.filename,
1743                quantity = 'elevation',
1744                cellsize = cellsize,
1745                NODATA_value = NODATA_value,
1746                verbose = False,
1747                format = 'ers')
1748
1749        #Check header data
1750        from ermapper_grids import read_ermapper_header, read_ermapper_data
1751       
1752        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
1753        #print header
1754        assert header['projection'].lower() == '"utm-56"'
1755        assert header['datum'].lower() == '"wgs84"'
1756        assert header['units'].lower() == '"meters"'   
1757        assert header['value'].lower() == '"elevation"'         
1758        assert header['xdimension'] == '0.25'
1759        assert header['ydimension'] == '0.25'   
1760        assert float(header['eastings']) == 308500.0   #xllcorner
1761        assert float(header['northings']) == 6189000.0 #yllcorner       
1762        assert int(header['nroflines']) == 5
1763        assert int(header['nrofcellsperline']) == 5     
1764        assert int(header['nullcellvalue']) == NODATA_value
1765        #FIXME - there is more in the header                   
1766
1767               
1768        #Check grid data               
1769        grid = read_ermapper_data(self.domain.filename + '_elevation') 
1770       
1771        #FIXME (Ole): Why is this the desired reference grid for -x-y?
1772        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
1773                    -1,    -1.25, -1.5,  -1.75, -2.0,
1774                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
1775                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1776                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
1777                                         
1778
1779        #print grid
1780        assert allclose(grid, ref_grid)
1781
1782        fid.close()
1783       
1784        #Cleanup
1785        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1786        #Done (Ole) - it was because sww2ers didn't close it's sww file
1787        os.remove(sww.filename)
1788
1789
1790
1791    def test_ferret2sww(self):
1792        """Test that georeferencing etc works when converting from
1793        ferret format (lat/lon) to sww format (UTM)
1794        """
1795        from Scientific.IO.NetCDF import NetCDFFile
1796
1797        #The test file has
1798        # LON = 150.66667, 150.83334, 151, 151.16667
1799        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1800        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1801        #
1802        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1803        # Fourth value (index==3) is -6.50198 cm
1804
1805
1806        from coordinate_transforms.redfearn import redfearn
1807
1808        fid = NetCDFFile('small_ha.nc')
1809        first_value = fid.variables['HA'][:][0,0,0]
1810        fourth_value = fid.variables['HA'][:][0,0,3]
1811
1812
1813        #Call conversion (with zero origin)
1814        ferret2sww('small', verbose=False,
1815                   origin = (56, 0, 0))
1816
1817
1818        #Work out the UTM coordinates for first point
1819        zone, e, n = redfearn(-34.5, 150.66667)
1820        #print zone, e, n
1821
1822        #Read output file 'small.sww'
1823        fid = NetCDFFile('small.sww')
1824
1825        x = fid.variables['x'][:]
1826        y = fid.variables['y'][:]
1827
1828        #Check that first coordinate is correctly represented
1829        assert allclose(x[0], e)
1830        assert allclose(y[0], n)
1831
1832        #Check first value
1833        stage = fid.variables['stage'][:]
1834        xmomentum = fid.variables['xmomentum'][:]
1835        ymomentum = fid.variables['ymomentum'][:]
1836
1837        #print ymomentum
1838
1839        assert allclose(stage[0,0], first_value/100)  #Meters
1840
1841        #Check fourth value
1842        assert allclose(stage[0,3], fourth_value/100)  #Meters
1843
1844        fid.close()
1845
1846        #Cleanup
1847        import os
1848        os.remove('small.sww')
1849
1850
1851
1852    def test_ferret2sww_2(self):
1853        """Test that georeferencing etc works when converting from
1854        ferret format (lat/lon) to sww format (UTM)
1855        """
1856        from Scientific.IO.NetCDF import NetCDFFile
1857
1858        #The test file has
1859        # LON = 150.66667, 150.83334, 151, 151.16667
1860        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1861        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1862        #
1863        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1864        # Fourth value (index==3) is -6.50198 cm
1865
1866
1867        from coordinate_transforms.redfearn import redfearn
1868
1869        fid = NetCDFFile('small_ha.nc')
1870
1871        #Pick a coordinate and a value
1872
1873        time_index = 1
1874        lat_index = 0
1875        lon_index = 2
1876
1877        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
1878        test_time = fid.variables['TIME'][:][time_index]
1879        test_lat = fid.variables['LAT'][:][lat_index]
1880        test_lon = fid.variables['LON'][:][lon_index]
1881
1882        linear_point_index = lat_index*4 + lon_index
1883        fid.close()
1884
1885        #Call conversion (with zero origin)
1886        ferret2sww('small', verbose=False,
1887                   origin = (56, 0, 0))
1888
1889
1890        #Work out the UTM coordinates for test point
1891        zone, e, n = redfearn(test_lat, test_lon)
1892
1893        #Read output file 'small.sww'
1894        fid = NetCDFFile('small.sww')
1895
1896        x = fid.variables['x'][:]
1897        y = fid.variables['y'][:]
1898
1899        #Check that test coordinate is correctly represented
1900        assert allclose(x[linear_point_index], e)
1901        assert allclose(y[linear_point_index], n)
1902
1903        #Check test value
1904        stage = fid.variables['stage'][:]
1905
1906        assert allclose(stage[time_index, linear_point_index], test_value/100)
1907
1908        fid.close()
1909
1910        #Cleanup
1911        import os
1912        os.remove('small.sww')
1913
1914
1915
1916    def test_ferret2sww3(self):
1917        """Elevation included
1918        """
1919        from Scientific.IO.NetCDF import NetCDFFile
1920
1921        #The test file has
1922        # LON = 150.66667, 150.83334, 151, 151.16667
1923        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1924        # ELEVATION = [-1 -2 -3 -4
1925        #             -5 -6 -7 -8
1926        #              ...
1927        #              ...      -16]
1928        # where the top left corner is -1m,
1929        # and the ll corner is -13.0m
1930        #
1931        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1932        # Fourth value (index==3) is -6.50198 cm
1933
1934        from coordinate_transforms.redfearn import redfearn
1935        import os
1936        fid1 = NetCDFFile('test_ha.nc','w')
1937        fid2 = NetCDFFile('test_ua.nc','w')
1938        fid3 = NetCDFFile('test_va.nc','w')
1939        fid4 = NetCDFFile('test_e.nc','w')
1940
1941        h1_list = [150.66667,150.83334,151.]
1942        h2_list = [-34.5,-34.33333]
1943
1944        long_name = 'LON'
1945        lat_name = 'LAT'
1946
1947        nx = 3
1948        ny = 2
1949
1950        for fid in [fid1,fid2,fid3]:
1951            fid.createDimension(long_name,nx)
1952            fid.createVariable(long_name,'d',(long_name,))
1953            fid.variables[long_name].point_spacing='uneven'
1954            fid.variables[long_name].units='degrees_east'
1955            fid.variables[long_name].assignValue(h1_list)
1956
1957            fid.createDimension(lat_name,ny)
1958            fid.createVariable(lat_name,'d',(lat_name,))
1959            fid.variables[lat_name].point_spacing='uneven'
1960            fid.variables[lat_name].units='degrees_north'
1961            fid.variables[lat_name].assignValue(h2_list)
1962
1963            fid.createDimension('TIME',2)
1964            fid.createVariable('TIME','d',('TIME',))
1965            fid.variables['TIME'].point_spacing='uneven'
1966            fid.variables['TIME'].units='seconds'
1967            fid.variables['TIME'].assignValue([0.,1.])
1968            if fid == fid3: break
1969
1970
1971        for fid in [fid4]:
1972            fid.createDimension(long_name,nx)
1973            fid.createVariable(long_name,'d',(long_name,))
1974            fid.variables[long_name].point_spacing='uneven'
1975            fid.variables[long_name].units='degrees_east'
1976            fid.variables[long_name].assignValue(h1_list)
1977
1978            fid.createDimension(lat_name,ny)
1979            fid.createVariable(lat_name,'d',(lat_name,))
1980            fid.variables[lat_name].point_spacing='uneven'
1981            fid.variables[lat_name].units='degrees_north'
1982            fid.variables[lat_name].assignValue(h2_list)
1983
1984        name = {}
1985        name[fid1]='HA'
1986        name[fid2]='UA'
1987        name[fid3]='VA'
1988        name[fid4]='ELEVATION'
1989
1990        units = {}
1991        units[fid1]='cm'
1992        units[fid2]='cm/s'
1993        units[fid3]='cm/s'
1994        units[fid4]='m'
1995
1996        values = {}
1997        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
1998        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
1999        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2000        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2001
2002        for fid in [fid1,fid2,fid3]:
2003          fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
2004          fid.variables[name[fid]].point_spacing='uneven'
2005          fid.variables[name[fid]].units=units[fid]
2006          fid.variables[name[fid]].assignValue(values[fid])
2007          fid.variables[name[fid]].missing_value = -99999999.
2008          if fid == fid3: break
2009
2010        for fid in [fid4]:
2011            fid.createVariable(name[fid],'d',(lat_name,long_name))
2012            fid.variables[name[fid]].point_spacing='uneven'
2013            fid.variables[name[fid]].units=units[fid]
2014            fid.variables[name[fid]].assignValue(values[fid])
2015            fid.variables[name[fid]].missing_value = -99999999.
2016
2017
2018        fid1.sync(); fid1.close()
2019        fid2.sync(); fid2.close()
2020        fid3.sync(); fid3.close()
2021        fid4.sync(); fid4.close()
2022
2023        fid1 = NetCDFFile('test_ha.nc','r')
2024        fid2 = NetCDFFile('test_e.nc','r')
2025        fid3 = NetCDFFile('test_va.nc','r')
2026
2027
2028        first_amp = fid1.variables['HA'][:][0,0,0]
2029        third_amp = fid1.variables['HA'][:][0,0,2]
2030        first_elevation = fid2.variables['ELEVATION'][0,0]
2031        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2032        first_speed = fid3.variables['VA'][0,0,0]
2033        third_speed = fid3.variables['VA'][:][0,0,2]
2034
2035        fid1.close()
2036        fid2.close()
2037        fid3.close()
2038
2039        #Call conversion (with zero origin)
2040        ferret2sww('test', verbose=False,
2041                   origin = (56, 0, 0))
2042
2043        os.remove('test_va.nc')
2044        os.remove('test_ua.nc')
2045        os.remove('test_ha.nc')
2046        os.remove('test_e.nc')
2047
2048        #Read output file 'test.sww'
2049        fid = NetCDFFile('test.sww')
2050
2051
2052        #Check first value
2053        elevation = fid.variables['elevation'][:]
2054        stage = fid.variables['stage'][:]
2055        xmomentum = fid.variables['xmomentum'][:]
2056        ymomentum = fid.variables['ymomentum'][:]
2057
2058        #print ymomentum
2059        first_height = first_amp/100 - first_elevation
2060        third_height = third_amp/100 - third_elevation
2061        first_momentum=first_speed*first_height/100
2062        third_momentum=third_speed*third_height/100
2063
2064        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2065        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2066
2067        fid.close()
2068
2069        #Cleanup
2070        os.remove('test.sww')
2071
2072
2073
2074
2075    def test_ferret2sww_nz_origin(self):
2076        from Scientific.IO.NetCDF import NetCDFFile
2077        from coordinate_transforms.redfearn import redfearn
2078
2079        #Call conversion (with nonzero origin)
2080        ferret2sww('small', verbose=False,
2081                   origin = (56, 100000, 200000))
2082
2083
2084        #Work out the UTM coordinates for first point
2085        zone, e, n = redfearn(-34.5, 150.66667)
2086
2087        #Read output file 'small.sww'
2088        fid = NetCDFFile('small.sww', 'r')
2089
2090        x = fid.variables['x'][:]
2091        y = fid.variables['y'][:]
2092
2093        #Check that first coordinate is correctly represented
2094        assert allclose(x[0], e-100000)
2095        assert allclose(y[0], n-200000)
2096
2097        fid.close()
2098
2099        #Cleanup
2100        os.remove('small.sww')
2101
2102
2103
2104    def test_sww_extent(self):
2105        """Not a test, rather a look at the sww format
2106        """
2107
2108        import time, os
2109        from Numeric import array, zeros, allclose, Float, concatenate
2110        from Scientific.IO.NetCDF import NetCDFFile
2111
2112        self.domain.filename = 'datatest' + str(id(self))
2113        self.domain.format = 'sww'
2114        self.domain.smooth = True
2115        self.domain.reduction = mean
2116        self.domain.set_datadir('.')
2117
2118
2119        sww = get_dataobject(self.domain)
2120        sww.store_connectivity()
2121        sww.store_timestep('stage')
2122        self.domain.time = 2.
2123
2124        #Modify stage at second timestep
2125        stage = self.domain.quantities['stage'].vertex_values
2126        self.domain.set_quantity('stage', stage/2)
2127
2128        sww.store_timestep('stage')
2129
2130        file_and_extension_name = self.domain.filename + ".sww"
2131        #print "file_and_extension_name",file_and_extension_name
2132        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2133               extent_sww(file_and_extension_name )
2134
2135        assert allclose(xmin, 0.0)
2136        assert allclose(xmax, 1.0)
2137        assert allclose(ymin, 0.0)
2138        assert allclose(ymax, 1.0)
2139        assert allclose(stagemin, -0.85)
2140        assert allclose(stagemax, 0.15)
2141
2142
2143        #Cleanup
2144        os.remove(sww.filename)
2145
2146
2147
2148    def test_sww2domain(self):
2149        ################################################
2150        #Create a test domain, and evolve and save it.
2151        ################################################
2152        from mesh_factory import rectangular
2153        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2154             Constant_height, Time_boundary, Transmissive_boundary
2155        from Numeric import array
2156
2157        #Create basic mesh
2158
2159        yiel=0.01
2160        points, vertices, boundary = rectangular(10,10)
2161
2162        #Create shallow water domain
2163        domain = Domain(points, vertices, boundary)
2164        domain.geo_reference = Geo_reference(56,11,11)
2165        domain.smooth = False
2166        domain.visualise = False
2167        domain.store = True
2168        domain.filename = 'bedslope'
2169        domain.default_order=2
2170        #Bed-slope and friction
2171        domain.set_quantity('elevation', lambda x,y: -x/3)
2172        domain.set_quantity('friction', 0.1)
2173        # Boundary conditions
2174        from math import sin, pi
2175        Br = Reflective_boundary(domain)
2176        Bt = Transmissive_boundary(domain)
2177        Bd = Dirichlet_boundary([0.2,0.,0.])
2178        Bw = Time_boundary(domain=domain,
2179                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2180
2181        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2182        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2183
2184        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2185        #Initial condition
2186        h = 0.05
2187        elevation = domain.quantities['elevation'].vertex_values
2188        domain.set_quantity('stage', elevation + h)
2189        #elevation = domain.get_quantity('elevation')
2190        #domain.set_quantity('stage', elevation + h)
2191
2192        domain.check_integrity()
2193        #Evolution
2194        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2195        #    domain.write_time()
2196            pass
2197
2198
2199        ##########################################
2200        #Import the example's file as a new domain
2201        ##########################################
2202        from data_manager import sww2domain
2203        from Numeric import allclose
2204        import os
2205
2206        filename = domain.datadir+os.sep+domain.filename+'.sww'
2207        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2208        #points, vertices, boundary = rectangular(15,15)
2209        #domain2.boundary = boundary
2210        ###################
2211        ##NOW TEST IT!!!
2212        ###################
2213
2214        bits = ['vertex_coordinates']
2215        for quantity in ['elevation']+domain.quantities_to_be_stored:
2216            bits.append('quantities["%s"].get_integral()'%quantity)
2217            bits.append('get_quantity("%s")'%quantity)
2218
2219        for bit in bits:
2220            #print 'testing that domain.'+bit+' has been restored'
2221            #print bit
2222        #print 'done'
2223            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2224
2225        ######################################
2226        #Now evolve them both, just to be sure
2227        ######################################x
2228        visualise = False
2229        #visualise = True
2230        domain.visualise = visualise
2231        domain.time = 0.
2232        from time import sleep
2233
2234        final = .1
2235        domain.set_quantity('friction', 0.1)
2236        domain.store = False
2237        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2238
2239
2240        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2241            if visualise: sleep(1.)
2242            #domain.write_time()
2243            pass
2244
2245        final = final - (domain2.starttime-domain.starttime)
2246        #BUT since domain1 gets time hacked back to 0:
2247        final = final + (domain2.starttime-domain.starttime)
2248
2249        domain2.smooth = False
2250        domain2.visualise = visualise
2251        domain2.store = False
2252        domain2.default_order=2
2253        domain2.set_quantity('friction', 0.1)
2254        #Bed-slope and friction
2255        # Boundary conditions
2256        Bd2=Dirichlet_boundary([0.2,0.,0.])
2257        domain2.boundary = domain.boundary
2258        #print 'domain2.boundary'
2259        #print domain2.boundary
2260        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2261        #domain2.set_boundary({'exterior': Bd})
2262
2263        domain2.check_integrity()
2264
2265        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2266            if visualise: sleep(1.)
2267            #domain2.write_time()
2268            pass
2269
2270        ###################
2271        ##NOW TEST IT!!!
2272        ##################
2273
2274        bits = [ 'vertex_coordinates']
2275
2276        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2277            bits.append('quantities["%s"].get_integral()'%quantity)
2278            bits.append('get_quantity("%s")'%quantity)
2279
2280        for bit in bits:
2281            #print bit
2282            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2283
2284
2285    def test_sww2domain2(self):
2286        ##################################################################
2287        #Same as previous test, but this checks how NaNs are handled.
2288        ##################################################################
2289
2290
2291        from mesh_factory import rectangular
2292        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2293             Constant_height, Time_boundary, Transmissive_boundary
2294        from Numeric import array
2295
2296        #Create basic mesh
2297        points, vertices, boundary = rectangular(2,2)
2298
2299        #Create shallow water domain
2300        domain = Domain(points, vertices, boundary)
2301        domain.smooth = False
2302        domain.visualise = False
2303        domain.store = True
2304        domain.filename = 'bedslope'
2305        domain.default_order=2
2306        domain.quantities_to_be_stored=['stage']
2307
2308        domain.set_quantity('elevation', lambda x,y: -x/3)
2309        domain.set_quantity('friction', 0.1)
2310
2311        from math import sin, pi
2312        Br = Reflective_boundary(domain)
2313        Bt = Transmissive_boundary(domain)
2314        Bd = Dirichlet_boundary([0.2,0.,0.])
2315        Bw = Time_boundary(domain=domain,
2316                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2317
2318        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2319
2320        h = 0.05
2321        elevation = domain.quantities['elevation'].vertex_values
2322        domain.set_quantity('stage', elevation + h)
2323
2324        domain.check_integrity()
2325
2326        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2327            pass
2328            #domain.write_time()
2329
2330
2331
2332        ##################################
2333        #Import the file as a new domain
2334        ##################################
2335        from data_manager import sww2domain
2336        from Numeric import allclose
2337        import os
2338
2339        filename = domain.datadir+os.sep+domain.filename+'.sww'
2340
2341        #Fail because NaNs are present
2342        try:
2343            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2344            assert True == False
2345        except:
2346            #Now import it, filling NaNs to be 0
2347            filler = 0
2348            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2349        bits = [ 'geo_reference.get_xllcorner()',
2350                'geo_reference.get_yllcorner()',
2351                'vertex_coordinates']
2352
2353        for quantity in ['elevation']+domain.quantities_to_be_stored:
2354            bits.append('quantities["%s"].get_integral()'%quantity)
2355            bits.append('get_quantity("%s")'%quantity)
2356
2357        for bit in bits:
2358        #    print 'testing that domain.'+bit+' has been restored'
2359            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2360
2361        assert max(max(domain2.get_quantity('xmomentum')))==filler
2362        assert min(min(domain2.get_quantity('xmomentum')))==filler
2363        assert max(max(domain2.get_quantity('ymomentum')))==filler
2364        assert min(min(domain2.get_quantity('ymomentum')))==filler
2365
2366        #print 'passed'
2367
2368        #cleanup
2369        #import os
2370        #os.remove(domain.datadir+'/'+domain.filename+'.sww')
2371
2372
2373    #def test_weed(self):
2374        from data_manager import weed
2375
2376        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2377        volumes1 = [[0,1,2],[3,4,5]]
2378        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2379        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2380
2381        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2382
2383        assert len(points2)==len(coordinates2)
2384        for i in range(len(coordinates2)):
2385            coordinate = tuple(coordinates2[i])
2386            assert points2.has_key(coordinate)
2387            points2[coordinate]=i
2388
2389        for triangle in volumes1:
2390            for coordinate in triangle:
2391                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2392                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2393
2394
2395    #FIXME This fails - smooth makes the comparism too hard for allclose
2396    def ztest_sww2domain3(self):
2397        ################################################
2398        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2399        ################################################
2400        from mesh_factory import rectangular
2401        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2402             Constant_height, Time_boundary, Transmissive_boundary
2403        from Numeric import array
2404        #Create basic mesh
2405
2406        yiel=0.01
2407        points, vertices, boundary = rectangular(10,10)
2408
2409        #Create shallow water domain
2410        domain = Domain(points, vertices, boundary)
2411        domain.geo_reference = Geo_reference(56,11,11)
2412        domain.smooth = True
2413        domain.visualise = False
2414        domain.store = True
2415        domain.filename = 'bedslope'
2416        domain.default_order=2
2417        #Bed-slope and friction
2418        domain.set_quantity('elevation', lambda x,y: -x/3)
2419        domain.set_quantity('friction', 0.1)
2420        # Boundary conditions
2421        from math import sin, pi
2422        Br = Reflective_boundary(domain)
2423        Bt = Transmissive_boundary(domain)
2424        Bd = Dirichlet_boundary([0.2,0.,0.])
2425        Bw = Time_boundary(domain=domain,
2426                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2427
2428        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2429
2430        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2431        #Initial condition
2432        h = 0.05
2433        elevation = domain.quantities['elevation'].vertex_values
2434        domain.set_quantity('stage', elevation + h)
2435        #elevation = domain.get_quantity('elevation')
2436        #domain.set_quantity('stage', elevation + h)
2437
2438        domain.check_integrity()
2439        #Evolution
2440        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2441        #    domain.write_time()
2442            pass
2443
2444
2445        ##########################################
2446        #Import the example's file as a new domain
2447        ##########################################
2448        from data_manager import sww2domain
2449        from Numeric import allclose
2450        import os
2451
2452        filename = domain.datadir+os.sep+domain.filename+'.sww'
2453        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2454        #points, vertices, boundary = rectangular(15,15)
2455        #domain2.boundary = boundary
2456        ###################
2457        ##NOW TEST IT!!!
2458        ###################
2459
2460        #FIXME smooth domain so that they can be compared
2461
2462
2463        bits = []#'vertex_coordinates']
2464        for quantity in ['elevation']+domain.quantities_to_be_stored:
2465            bits.append('quantities["%s"].get_integral()'%quantity)
2466            #bits.append('get_quantity("%s")'%quantity)
2467
2468        for bit in bits:
2469            #print 'testing that domain.'+bit+' has been restored'
2470            #print bit
2471            #print 'done'
2472            #print ('domain.'+bit), eval('domain.'+bit)
2473            #print ('domain2.'+bit), eval('domain2.'+bit)
2474            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2475            pass
2476
2477        ######################################
2478        #Now evolve them both, just to be sure
2479        ######################################x
2480        visualise = False
2481        visualise = True
2482        domain.visualise = visualise
2483        domain.time = 0.
2484        from time import sleep
2485
2486        final = .5
2487        domain.set_quantity('friction', 0.1)
2488        domain.store = False
2489        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2490
2491        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2492            if visualise: sleep(.03)
2493            #domain.write_time()
2494            pass
2495
2496        domain2.smooth = True
2497        domain2.visualise = visualise
2498        domain2.store = False
2499        domain2.default_order=2
2500        domain2.set_quantity('friction', 0.1)
2501        #Bed-slope and friction
2502        # Boundary conditions
2503        Bd2=Dirichlet_boundary([0.2,0.,0.])
2504        Br2 = Reflective_boundary(domain2)
2505        domain2.boundary = domain.boundary
2506        #print 'domain2.boundary'
2507        #print domain2.boundary
2508        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2509        #domain2.boundary = domain.boundary
2510        #domain2.set_boundary({'exterior': Bd})
2511
2512        domain2.check_integrity()
2513
2514        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2515            if visualise: sleep(.03)
2516            #domain2.write_time()
2517            pass
2518
2519        ###################
2520        ##NOW TEST IT!!!
2521        ##################
2522
2523        bits = [ 'vertex_coordinates']
2524
2525        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2526            #bits.append('quantities["%s"].get_integral()'%quantity)
2527            bits.append('get_quantity("%s")'%quantity)
2528
2529        for bit in bits:
2530            print bit
2531            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2532
2533
2534    def test_decimate_dem(self):
2535        """Test decimation of dem file
2536        """
2537
2538        import os
2539        from Numeric import ones, allclose, Float, arange
2540        from Scientific.IO.NetCDF import NetCDFFile
2541
2542        #Write test dem file
2543        root = 'decdemtest'
2544
2545        filename = root + '.dem'
2546        fid = NetCDFFile(filename, 'w')
2547
2548        fid.institution = 'Geoscience Australia'
2549        fid.description = 'NetCDF DEM format for compact and portable ' +\
2550                          'storage of spatial point data'
2551
2552        nrows = 15
2553        ncols = 18
2554
2555        fid.ncols = ncols
2556        fid.nrows = nrows
2557        fid.xllcorner = 2000.5
2558        fid.yllcorner = 3000.5
2559        fid.cellsize = 25
2560        fid.NODATA_value = -9999
2561
2562        fid.zone = 56
2563        fid.false_easting = 0.0
2564        fid.false_northing = 0.0
2565        fid.projection = 'UTM'
2566        fid.datum = 'WGS84'
2567        fid.units = 'METERS'
2568
2569        fid.createDimension('number_of_points', nrows*ncols)
2570
2571        fid.createVariable('elevation', Float, ('number_of_points',))
2572
2573        elevation = fid.variables['elevation']
2574
2575        elevation[:] = (arange(nrows*ncols))
2576
2577        fid.close()
2578
2579        #generate the elevation values expected in the decimated file
2580        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2581                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2582                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2583                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2584                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2585                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2586                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2587                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2588                         (144+145+146+162+163+164+180+181+182) / 9.0,
2589                         (148+149+150+166+167+168+184+185+186) / 9.0,
2590                         (152+153+154+170+171+172+188+189+190) / 9.0,
2591                         (156+157+158+174+175+176+192+193+194) / 9.0,
2592                         (216+217+218+234+235+236+252+253+254) / 9.0,
2593                         (220+221+222+238+239+240+256+257+258) / 9.0,
2594                         (224+225+226+242+243+244+260+261+262) / 9.0,
2595                         (228+229+230+246+247+248+264+265+266) / 9.0]
2596
2597        #generate a stencil for computing the decimated values
2598        stencil = ones((3,3), Float) / 9.0
2599
2600        decimate_dem(root, stencil=stencil, cellsize_new=100)
2601
2602        #Open decimated NetCDF file
2603        fid = NetCDFFile(root + '_100.dem', 'r')
2604
2605        # Get decimated elevation
2606        elevation = fid.variables['elevation']
2607
2608        #Check values
2609        assert allclose(elevation, ref_elevation)
2610
2611        #Cleanup
2612        fid.close()
2613
2614        os.remove(root + '.dem')
2615        os.remove(root + '_100.dem')
2616
2617    def test_decimate_dem_NODATA(self):
2618        """Test decimation of dem file that includes NODATA values
2619        """
2620
2621        import os
2622        from Numeric import ones, allclose, Float, arange, reshape
2623        from Scientific.IO.NetCDF import NetCDFFile
2624
2625        #Write test dem file
2626        root = 'decdemtest'
2627
2628        filename = root + '.dem'
2629        fid = NetCDFFile(filename, 'w')
2630
2631        fid.institution = 'Geoscience Australia'
2632        fid.description = 'NetCDF DEM format for compact and portable ' +\
2633                          'storage of spatial point data'
2634
2635        nrows = 15
2636        ncols = 18
2637        NODATA_value = -9999
2638
2639        fid.ncols = ncols
2640        fid.nrows = nrows
2641        fid.xllcorner = 2000.5
2642        fid.yllcorner = 3000.5
2643        fid.cellsize = 25
2644        fid.NODATA_value = NODATA_value
2645
2646        fid.zone = 56
2647        fid.false_easting = 0.0
2648        fid.false_northing = 0.0
2649        fid.projection = 'UTM'
2650        fid.datum = 'WGS84'
2651        fid.units = 'METERS'
2652
2653        fid.createDimension('number_of_points', nrows*ncols)
2654
2655        fid.createVariable('elevation', Float, ('number_of_points',))
2656
2657        elevation = fid.variables['elevation']
2658
2659        #generate initial elevation values
2660        elevation_tmp = (arange(nrows*ncols))
2661        #add some NODATA values
2662        elevation_tmp[0]   = NODATA_value
2663        elevation_tmp[95]  = NODATA_value
2664        elevation_tmp[188] = NODATA_value
2665        elevation_tmp[189] = NODATA_value
2666        elevation_tmp[190] = NODATA_value
2667        elevation_tmp[209] = NODATA_value
2668        elevation_tmp[252] = NODATA_value
2669
2670        elevation[:] = elevation_tmp
2671
2672        fid.close()
2673
2674        #generate the elevation values expected in the decimated file
2675        ref_elevation = [NODATA_value,
2676                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2677                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2678                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2679                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2680                         NODATA_value,
2681                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2682                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2683                         (144+145+146+162+163+164+180+181+182) / 9.0,
2684                         (148+149+150+166+167+168+184+185+186) / 9.0,
2685                         NODATA_value,
2686                         (156+157+158+174+175+176+192+193+194) / 9.0,
2687                         NODATA_value,
2688                         (220+221+222+238+239+240+256+257+258) / 9.0,
2689                         (224+225+226+242+243+244+260+261+262) / 9.0,
2690                         (228+229+230+246+247+248+264+265+266) / 9.0]
2691
2692        #generate a stencil for computing the decimated values
2693        stencil = ones((3,3), Float) / 9.0
2694
2695        decimate_dem(root, stencil=stencil, cellsize_new=100)
2696
2697        #Open decimated NetCDF file
2698        fid = NetCDFFile(root + '_100.dem', 'r')
2699
2700        # Get decimated elevation
2701        elevation = fid.variables['elevation']
2702
2703        #Check values
2704        assert allclose(elevation, ref_elevation)
2705
2706        #Cleanup
2707        fid.close()
2708
2709        os.remove(root + '.dem')
2710        os.remove(root + '_100.dem')
2711
2712    def xxxtestz_sww2ers_real(self):
2713        """Test that sww information can be converted correctly to asc/prj
2714        format readable by e.g. ArcView
2715        """
2716
2717        import time, os
2718        from Numeric import array, zeros, allclose, Float, concatenate
2719        from Scientific.IO.NetCDF import NetCDFFile
2720
2721        # the memory optimised least squares
2722        #  cellsize = 20,   # this one seems to hang
2723        #  cellsize = 200000, # Ran 1 test in 269.703s
2724                                #Ran 1 test in 267.344s
2725        #  cellsize = 20000,  # Ran 1 test in 460.922s
2726        #  cellsize = 2000   #Ran 1 test in 5340.250s
2727        #  cellsize = 200   #this one seems to hang, building matirx A
2728
2729        # not optimised
2730        # seems to hang
2731        #  cellsize = 2000   # Ran 1 test in 5334.563s
2732        #Export to ascii/prj files
2733        sww2dem('karratha_100m',
2734                quantity = 'depth',
2735                cellsize = 200000,
2736                verbose = True)
2737
2738    def test_read_asc(self):
2739        """Test conversion from dem in ascii format to native NetCDF xya 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        import data_manager
2747        #Write test asc file
2748        filename = tempfile.mktemp(".000")
2749        fid = open(filename, 'w')
2750        fid.write("""ncols         7
2751nrows         4
2752xllcorner     2000.5
2753yllcorner     3000.5
2754cellsize      25
2755NODATA_value  -9999
2756    97.921    99.285   125.588   180.830   258.645   342.872   415.836
2757   473.157   514.391   553.893   607.120   678.125   777.283   883.038
2758   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
2759   502.645   516.230   504.739   450.604   388.500   338.097   514.980
2760""")
2761        fid.close()
2762        bath_metadata, grid = \
2763                   data_manager._read_asc(filename, verbose=False)
2764        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
2765        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
2766        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
2767        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
2768        self.failUnless(grid[0][0]  == 97.921,  'Failed')
2769        self.failUnless(grid[3][6]  == 514.980,  'Failed')
2770
2771        os.remove(filename)
2772       
2773    def test_asc_csiro2sww(self):
2774        import tempfile
2775
2776        bath_dir = tempfile.mkdtemp()
2777        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
2778        #bath_dir = 'bath_data_manager_test'
2779        #print "os.getcwd( )",os.getcwd( )
2780        elevation_dir =  tempfile.mkdtemp()
2781        #elevation_dir = 'elev_expanded'
2782        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
2783        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
2784       
2785        fid = open(bath_dir_filename, 'w')
2786        fid.write(""" ncols             3
2787 nrows             2
2788 xllcorner    148.00000
2789 yllcorner    -38.00000
2790 cellsize       0.25
2791 nodata_value   -9999.0
2792    9000.000    -1000.000    3000.0
2793   -1000.000    9000.000  -1000.000
2794""")
2795        fid.close()
2796       
2797        fid = open(elevation_dir_filename1, 'w')
2798        fid.write(""" ncols             3
2799 nrows             2
2800 xllcorner    148.00000
2801 yllcorner    -38.00000
2802 cellsize       0.25
2803 nodata_value   -9999.0
2804    9000.000    0.000    3000.0
2805     0.000     9000.000     0.000
2806""")
2807        fid.close()
2808
2809        fid = open(elevation_dir_filename2, 'w')
2810        fid.write(""" ncols             3
2811 nrows             2
2812 xllcorner    148.00000
2813 yllcorner    -38.00000
2814 cellsize       0.25
2815 nodata_value   -9999.0
2816    9000.000    4000.000    4000.0
2817    4000.000    9000.000    4000.000
2818""")
2819        fid.close()
2820         
2821        ucur_dir =  tempfile.mkdtemp()
2822        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
2823        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
2824       
2825        fid = open(ucur_dir_filename1, 'w')
2826        fid.write(""" ncols             3
2827 nrows             2
2828 xllcorner    148.00000
2829 yllcorner    -38.00000
2830 cellsize       0.25
2831 nodata_value   -9999.0
2832    90.000    60.000    30.0
2833    10.000    10.000    10.000
2834""")
2835        fid.close()
2836        fid = open(ucur_dir_filename2, 'w')
2837        fid.write(""" ncols             3
2838 nrows             2
2839 xllcorner    148.00000
2840 yllcorner    -38.00000
2841 cellsize       0.25
2842 nodata_value   -9999.0
2843    90.000    60.000    30.0
2844    10.000    10.000    10.000
2845""")
2846        fid.close()
2847       
2848        vcur_dir =  tempfile.mkdtemp()
2849        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
2850        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
2851       
2852        fid = open(vcur_dir_filename1, 'w')
2853        fid.write(""" ncols             3
2854 nrows             2
2855 xllcorner    148.00000
2856 yllcorner    -38.00000
2857 cellsize       0.25
2858 nodata_value   -9999.0
2859    90.000    60.000    30.0
2860    10.000    10.000    10.000
2861""")
2862        fid.close()
2863        fid = open(vcur_dir_filename2, 'w')
2864        fid.write(""" ncols             3
2865 nrows             2
2866 xllcorner    148.00000
2867 yllcorner    -38.00000
2868 cellsize       0.25
2869 nodata_value   -9999.0
2870    90.000    60.000    30.0
2871    10.000    10.000    10.000
2872""")
2873        fid.close()
2874       
2875        sww_file = 'a_test.sww'
2876        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
2877
2878        # check the sww file
2879       
2880        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
2881        x = fid.variables['x'][:]
2882        y = fid.variables['y'][:]
2883        z = fid.variables['z'][:]
2884        stage = fid.variables['stage'][:]
2885        xmomentum = fid.variables['xmomentum'][:]
2886        geo_ref = Geo_reference(NetCDFObject=fid)
2887        #print "geo_ref",geo_ref
2888        x_ref = geo_ref.get_xllcorner()
2889        y_ref = geo_ref.get_yllcorner()
2890        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
2891        assert allclose(x_ref, 587798.418) # (-38, 148)
2892        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
2893       
2894        #Zone:   55   
2895        #Easting:  588095.674  Northing: 5821451.722
2896        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
2897        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
2898
2899        #Zone:   55   
2900        #Easting:  632145.632  Northing: 5820863.269
2901        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
2902        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
2903
2904        #Zone:   55   
2905        #Easting:  609748.788  Northing: 5793447.860
2906        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
2907        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
2908
2909        assert allclose(z[0],9000.0 )
2910        assert allclose(stage[0][1],0.0 )
2911
2912        #(4000+1000)*60
2913        assert allclose(xmomentum[1][1],300000.0 )
2914       
2915       
2916        fid.close()
2917   
2918        #tidy up
2919        os.remove(bath_dir_filename)
2920        os.rmdir(bath_dir)
2921
2922        os.remove(elevation_dir_filename1)
2923        os.remove(elevation_dir_filename2)
2924        os.rmdir(elevation_dir)
2925       
2926        os.remove(ucur_dir_filename1)
2927        os.remove(ucur_dir_filename2)
2928        os.rmdir(ucur_dir)
2929       
2930        os.remove(vcur_dir_filename1)
2931        os.remove(vcur_dir_filename2)
2932        os.rmdir(vcur_dir)
2933
2934
2935        # remove sww file
2936        os.remove(sww_file)
2937       
2938    def test_asc_csiro2sww2(self):
2939        import tempfile
2940
2941        bath_dir = tempfile.mkdtemp()
2942        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
2943        #bath_dir = 'bath_data_manager_test'
2944        #print "os.getcwd( )",os.getcwd( )
2945        elevation_dir =  tempfile.mkdtemp()
2946        #elevation_dir = 'elev_expanded'
2947        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
2948        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
2949       
2950        fid = open(bath_dir_filename, 'w')
2951        fid.write(""" ncols             3
2952 nrows             2
2953 xllcorner    148.00000
2954 yllcorner    -38.00000
2955 cellsize       0.25
2956 nodata_value   -9999.0
2957    9000.000    -1000.000    3000.0
2958   -1000.000    9000.000  -1000.000
2959""")
2960        fid.close()
2961       
2962        fid = open(elevation_dir_filename1, 'w')
2963        fid.write(""" ncols             3
2964 nrows             2
2965 xllcorner    148.00000
2966 yllcorner    -38.00000
2967 cellsize       0.25
2968 nodata_value   -9999.0
2969    9000.000    0.000    3000.0
2970     0.000     -9999.000     -9999.000
2971""")
2972        fid.close()
2973
2974        fid = open(elevation_dir_filename2, 'w')
2975        fid.write(""" ncols             3
2976 nrows             2
2977 xllcorner    148.00000
2978 yllcorner    -38.00000
2979 cellsize       0.25
2980 nodata_value   -9999.0
2981    9000.000    4000.000    4000.0
2982    4000.000    9000.000    4000.000
2983""")
2984        fid.close()
2985         
2986        ucur_dir =  tempfile.mkdtemp()
2987        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
2988        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
2989       
2990        fid = open(ucur_dir_filename1, 'w')
2991        fid.write(""" ncols             3
2992 nrows             2
2993 xllcorner    148.00000
2994 yllcorner    -38.00000
2995 cellsize       0.25
2996 nodata_value   -9999.0
2997    90.000    60.000    30.0
2998    10.000    10.000    10.000
2999""")
3000        fid.close()
3001        fid = open(ucur_dir_filename2, 'w')
3002        fid.write(""" ncols             3
3003 nrows             2
3004 xllcorner    148.00000
3005 yllcorner    -38.00000
3006 cellsize       0.25
3007 nodata_value   -9999.0
3008    90.000    60.000    30.0
3009    10.000    10.000    10.000
3010""")
3011        fid.close()
3012       
3013        vcur_dir =  tempfile.mkdtemp()
3014        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3015        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3016       
3017        fid = open(vcur_dir_filename1, 'w')
3018        fid.write(""" ncols             3
3019 nrows             2
3020 xllcorner    148.00000
3021 yllcorner    -38.00000
3022 cellsize       0.25
3023 nodata_value   -9999.0
3024    90.000    60.000    30.0
3025    10.000    10.000    10.000
3026""")
3027        fid.close()
3028        fid = open(vcur_dir_filename2, 'w')
3029        fid.write(""" ncols             3
3030 nrows             2
3031 xllcorner    148.00000
3032 yllcorner    -38.00000
3033 cellsize       0.25
3034 nodata_value   -9999.0
3035    90.000    60.000    30.0
3036    10.000    10.000    10.000
3037""")
3038        fid.close()
3039       
3040        try:
3041            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3042                          vcur_dir, sww_file)
3043        except:
3044            #tidy up
3045            os.remove(bath_dir_filename)
3046            os.rmdir(bath_dir)
3047
3048            os.remove(elevation_dir_filename1)
3049            os.remove(elevation_dir_filename2)
3050            os.rmdir(elevation_dir)
3051       
3052            os.remove(ucur_dir_filename1)
3053            os.remove(ucur_dir_filename2)
3054            os.rmdir(ucur_dir)
3055       
3056            os.remove(vcur_dir_filename1)
3057            os.remove(vcur_dir_filename2)
3058            os.rmdir(vcur_dir)
3059        else:
3060            #tidy up
3061            os.remove(bath_dir_filename)
3062            os.rmdir(bath_dir)
3063
3064            os.remove(elevation_dir_filename1)
3065            os.remove(elevation_dir_filename2)
3066            os.rmdir(elevation_dir)
3067            raise 'Should raise exception'
3068       
3069            os.remove(ucur_dir_filename1)
3070            os.remove(ucur_dir_filename2)
3071            os.rmdir(ucur_dir)
3072       
3073            os.remove(vcur_dir_filename1)
3074            os.remove(vcur_dir_filename2)
3075            os.rmdir(vcur_dir)
3076
3077       
3078     
3079    def test_asc_csiro2sww3(self):
3080        import tempfile
3081
3082        bath_dir = tempfile.mkdtemp()
3083        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3084        #bath_dir = 'bath_data_manager_test'
3085        #print "os.getcwd( )",os.getcwd( )
3086        elevation_dir =  tempfile.mkdtemp()
3087        #elevation_dir = 'elev_expanded'
3088        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3089        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3090       
3091        fid = open(bath_dir_filename, 'w')
3092        fid.write(""" ncols             3
3093 nrows             2
3094 xllcorner    148.00000
3095 yllcorner    -38.00000
3096 cellsize       0.25
3097 nodata_value   -9999.0
3098    9000.000    -1000.000    3000.0
3099   -1000.000    9000.000  -1000.000
3100""")
3101        fid.close()
3102       
3103        fid = open(elevation_dir_filename1, 'w')
3104        fid.write(""" ncols             3
3105 nrows             2
3106 xllcorner    148.00000
3107 yllcorner    -38.00000
3108 cellsize       0.25
3109 nodata_value   -9999.0
3110    9000.000    0.000    3000.0
3111     0.000     -9999.000     -9999.000
3112""")
3113        fid.close()
3114
3115        fid = open(elevation_dir_filename2, 'w')
3116        fid.write(""" ncols             3
3117 nrows             2
3118 xllcorner    148.00000
3119 yllcorner    -38.00000
3120 cellsize       0.25
3121 nodata_value   -9999.0
3122    9000.000    4000.000    4000.0
3123    4000.000    9000.000    4000.000
3124""")
3125        fid.close()
3126         
3127        ucur_dir =  tempfile.mkdtemp()
3128        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3129        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3130       
3131        fid = open(ucur_dir_filename1, 'w')
3132        fid.write(""" ncols             3
3133 nrows             2
3134 xllcorner    148.00000
3135 yllcorner    -38.00000
3136 cellsize       0.25
3137 nodata_value   -9999.0
3138    90.000    60.000    30.0
3139    10.000    10.000    10.000
3140""")
3141        fid.close()
3142        fid = open(ucur_dir_filename2, 'w')
3143        fid.write(""" ncols             3
3144 nrows             2
3145 xllcorner    148.00000
3146 yllcorner    -38.00000
3147 cellsize       0.25
3148 nodata_value   -9999.0
3149    90.000    60.000    30.0
3150    10.000    10.000    10.000
3151""")
3152        fid.close()
3153       
3154        vcur_dir =  tempfile.mkdtemp()
3155        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3156        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3157       
3158        fid = open(vcur_dir_filename1, 'w')
3159        fid.write(""" ncols             3
3160 nrows             2
3161 xllcorner    148.00000
3162 yllcorner    -38.00000
3163 cellsize       0.25
3164 nodata_value   -9999.0
3165    90.000    60.000    30.0
3166    10.000    10.000    10.000
3167""")
3168        fid.close()
3169        fid = open(vcur_dir_filename2, 'w')
3170        fid.write(""" ncols             3
3171 nrows             2
3172 xllcorner    148.00000
3173 yllcorner    -38.00000
3174 cellsize       0.25
3175 nodata_value   -9999.0
3176    90.000    60.000    30.0
3177    10.000    10.000    10.000
3178""")
3179        fid.close()
3180       
3181        sww_file = 'a_test.sww'
3182        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3183                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3184                      mean_stage = 100)
3185
3186        # check the sww file
3187       
3188        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3189        x = fid.variables['x'][:]
3190        y = fid.variables['y'][:]
3191        z = fid.variables['z'][:]
3192        stage = fid.variables['stage'][:]
3193        xmomentum = fid.variables['xmomentum'][:]
3194        geo_ref = Geo_reference(NetCDFObject=fid)
3195        #print "geo_ref",geo_ref
3196        x_ref = geo_ref.get_xllcorner()
3197        y_ref = geo_ref.get_yllcorner()
3198        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3199        assert allclose(x_ref, 587798.418) # (-38, 148)
3200        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3201       
3202        #Zone:   55   
3203        #Easting:  588095.674  Northing: 5821451.722
3204        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3205        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3206
3207        #Zone:   55   
3208        #Easting:  632145.632  Northing: 5820863.269
3209        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3210        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3211
3212        #Zone:   55   
3213        #Easting:  609748.788  Northing: 5793447.860
3214        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3215        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3216
3217        assert allclose(z[0],9000.0 )
3218        assert allclose(stage[0][4],100.0 )
3219        assert allclose(stage[0][5],100.0 )
3220
3221        #(100.0 - 9000)*10
3222        assert allclose(xmomentum[0][4], -89000.0 )
3223       
3224        #(100.0 - -1000.000)*10
3225        assert allclose(xmomentum[0][5], 11000.0 )
3226       
3227        fid.close()
3228   
3229        #tidy up
3230        os.remove(bath_dir_filename)
3231        os.rmdir(bath_dir)
3232
3233        os.remove(elevation_dir_filename1)
3234        os.remove(elevation_dir_filename2)
3235        os.rmdir(elevation_dir)
3236       
3237        os.remove(ucur_dir_filename1)
3238        os.remove(ucur_dir_filename2)
3239        os.rmdir(ucur_dir)
3240       
3241        os.remove(vcur_dir_filename1)
3242        os.remove(vcur_dir_filename2)
3243        os.rmdir(vcur_dir)
3244
3245        # remove sww file
3246        os.remove(sww_file)
3247   
3248     
3249    def test_asc_csiro2sww4(self):
3250        """
3251        Test specifying the extent
3252        """
3253       
3254        import tempfile
3255
3256        bath_dir = tempfile.mkdtemp()
3257        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3258        #bath_dir = 'bath_data_manager_test'
3259        #print "os.getcwd( )",os.getcwd( )
3260        elevation_dir =  tempfile.mkdtemp()
3261        #elevation_dir = 'elev_expanded'
3262        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3263        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3264       
3265        fid = open(bath_dir_filename, 'w')
3266        fid.write(""" ncols             4
3267 nrows             4
3268 xllcorner    148.00000
3269 yllcorner    -38.00000
3270 cellsize       0.25
3271 nodata_value   -9999.0
3272   -9000.000    -1000.000   -3000.0 -2000.000
3273   -1000.000    9000.000  -1000.000 -3000.000
3274   -4000.000    6000.000   2000.000 -5000.000
3275   -9000.000    -1000.000   -3000.0 -2000.000
3276""")
3277        fid.close()
3278       
3279        fid = open(elevation_dir_filename1, 'w')
3280        fid.write(""" ncols             4
3281 nrows             4
3282 xllcorner    148.00000
3283 yllcorner    -38.00000
3284 cellsize       0.25
3285 nodata_value   -9999.0
3286   -900.000    -100.000   -300.0 -200.000
3287   -100.000    900.000  -100.000 -300.000
3288   -400.000    600.000   200.000 -500.000
3289   -900.000    -100.000   -300.0 -200.000
3290""")
3291        fid.close()
3292
3293        fid = open(elevation_dir_filename2, 'w')
3294        fid.write(""" ncols             4
3295 nrows             4
3296 xllcorner    148.00000
3297 yllcorner    -38.00000
3298 cellsize       0.25
3299 nodata_value   -9999.0
3300   -990.000    -110.000   -330.0 -220.000
3301   -110.000    990.000  -110.000 -330.000
3302   -440.000    660.000   220.000 -550.000
3303   -990.000    -110.000   -330.0 -220.000
3304""")
3305        fid.close()
3306         
3307        ucur_dir =  tempfile.mkdtemp()
3308        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3309        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3310       
3311        fid = open(ucur_dir_filename1, 'w')
3312        fid.write(""" ncols             4
3313 nrows             4
3314 xllcorner    148.00000
3315 yllcorner    -38.00000
3316 cellsize       0.25
3317 nodata_value   -9999.0
3318   -90.000    -10.000   -30.0 -20.000
3319   -10.000    90.000  -10.000 -30.000
3320   -40.000    60.000   20.000 -50.000
3321   -90.000    -10.000   -30.0 -20.000
3322""")
3323        fid.close()
3324        fid = open(ucur_dir_filename2, 'w')
3325        fid.write(""" ncols             4
3326 nrows             4
3327 xllcorner    148.00000
3328 yllcorner    -38.00000
3329 cellsize       0.25
3330 nodata_value   -9999.0
3331   -90.000    -10.000   -30.0 -20.000
3332   -10.000    99.000  -11.000 -30.000
3333   -40.000    66.000   22.000 -50.000
3334   -90.000    -10.000   -30.0 -20.000
3335""")
3336        fid.close()
3337       
3338        vcur_dir =  tempfile.mkdtemp()
3339        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3340        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3341       
3342        fid = open(vcur_dir_filename1, 'w')
3343        fid.write(""" ncols             4
3344 nrows             4
3345 xllcorner    148.00000
3346 yllcorner    -38.00000
3347 cellsize       0.25
3348 nodata_value   -9999.0
3349   -90.000    -10.000   -30.0 -20.000
3350   -10.000    80.000  -20.000 -30.000
3351   -40.000    50.000   10.000 -50.000
3352   -90.000    -10.000   -30.0 -20.000
3353""")
3354        fid.close()
3355        fid = open(vcur_dir_filename2, 'w')
3356        fid.write(""" ncols             4
3357 nrows             4
3358 xllcorner    148.00000
3359 yllcorner    -38.00000
3360 cellsize       0.25
3361 nodata_value   -9999.0
3362   -90.000    -10.000   -30.0 -20.000
3363   -10.000    88.000  -22.000 -30.000
3364   -40.000    55.000   11.000 -50.000
3365   -90.000    -10.000   -30.0 -20.000
3366""")
3367        fid.close()
3368       
3369        sww_file = tempfile.mktemp(".sww")
3370        #sww_file = 'a_test.sww'
3371        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3372                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3373                      mean_stage = 100,
3374                       minlat = -37.6, maxlat = -37.6,
3375                  minlon = 148.3, maxlon = 148.3
3376                      #,verbose = True
3377                      )
3378
3379        # check the sww file
3380       
3381        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3382        x = fid.variables['x'][:]
3383        y = fid.variables['y'][:]
3384        z = fid.variables['z'][:]
3385        stage = fid.variables['stage'][:]
3386        xmomentum = fid.variables['xmomentum'][:]
3387        ymomentum = fid.variables['ymomentum'][:]
3388        geo_ref = Geo_reference(NetCDFObject=fid)
3389        #print "geo_ref",geo_ref
3390        x_ref = geo_ref.get_xllcorner()
3391        y_ref = geo_ref.get_yllcorner()
3392        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3393
3394        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3395        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3396        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3397
3398        #Easting:  632145.632  Northing: 5820863.269
3399        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3400       
3401        #print "x",x
3402        #print "y",y
3403        self.failUnless(len(x) == 4,'failed') # 2*2
3404        self.failUnless(len(x) == 4,'failed') # 2*2
3405       
3406        #Zone:   55
3407        #Easting:  632145.632  Northing: 5820863.269
3408        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3409        # magic number - y is close enough for me.
3410        assert allclose(x[3], 632145.63 - x_ref)
3411        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3412
3413        assert allclose(z[0],9000.0 ) #z is elevation info
3414        #print "z",z
3415        # 2 time steps, 4 points
3416        self.failUnless(xmomentum.shape == (2,4), 'failed') 
3417        self.failUnless(ymomentum.shape == (2,4), 'failed') 
3418       
3419        #(100.0 - -1000.000)*10
3420        #assert allclose(xmomentum[0][5], 11000.0 )
3421       
3422        fid.close()
3423       
3424        # is the sww file readable?
3425        #Lets see if we can convert it to a dem!
3426        #print "sww_file",sww_file
3427        #dem_file = tempfile.mktemp(".dem")
3428        domain = sww2domain(sww_file) ###, dem_file)
3429        domain.check_integrity()
3430
3431        #tidy up
3432        os.remove(bath_dir_filename)
3433        os.rmdir(bath_dir)
3434
3435        os.remove(elevation_dir_filename1)
3436        os.remove(elevation_dir_filename2)
3437        os.rmdir(elevation_dir)
3438       
3439        os.remove(ucur_dir_filename1)
3440        os.remove(ucur_dir_filename2)
3441        os.rmdir(ucur_dir)
3442       
3443        os.remove(vcur_dir_filename1)
3444        os.remove(vcur_dir_filename2)
3445        os.rmdir(vcur_dir)
3446
3447
3448     
3449       
3450        # remove sww file
3451        os.remove(sww_file)
3452       
3453        # remove dem file
3454        #os.remove(dem_file)
3455
3456    def test_get_min_max_indexes(self):
3457        latitudes = [3,2,1,0]
3458        longitudes = [0,10,20,30]
3459
3460        # k - lat
3461        # l - lon
3462        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3463            latitudes,longitudes,
3464            -10,4,-10,31)
3465
3466        #print "kmin",kmin;print "kmax",kmax
3467        #print "lmin",lmin;print "lmax",lmax
3468        latitudes_new = latitudes[kmin:kmax]
3469        longitudes_news = longitudes[lmin:lmax]
3470        #print "latitudes_new", latitudes_new
3471        #print "longitudes_news",longitudes_news
3472        self.failUnless(latitudes == latitudes_new and \
3473                        longitudes == longitudes_news,
3474                         'failed')
3475       
3476        ## 2nd test
3477        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3478            latitudes,longitudes,
3479            0.5,2.5,5,25)
3480        #print "kmin",kmin;print "kmax",kmax
3481        #print "lmin",lmin;print "lmax",lmax
3482        latitudes_new = latitudes[kmin:kmax]
3483        longitudes_news = longitudes[lmin:lmax]
3484        #print "latitudes_new", latitudes_new
3485        #print "longitudes_news",longitudes_news
3486       
3487        self.failUnless(latitudes == latitudes_new and \
3488                        longitudes == longitudes_news,
3489                         'failed')
3490       
3491        ## 3rd test
3492        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3493                                                                    longitudes,
3494                                                      1.1,1.9,12,17)
3495        #print "kmin",kmin;print "kmax",kmax
3496        #print "lmin",lmin;print "lmax",lmax
3497        latitudes_new = latitudes[kmin:kmax]
3498        longitudes_news = longitudes[lmin:lmax]
3499        #print "latitudes_new", latitudes_new
3500        #print "longitudes_news",longitudes_news
3501       
3502        self.failUnless(latitudes_new == [2, 1] and \
3503                        longitudes_news == [10, 20],
3504                         'failed')
3505
3506       
3507        ## 4th test
3508        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3509            latitudes,longitudes,
3510                                                      -0.1,1.9,-2,17)
3511        #print "kmin",kmin;print "kmax",kmax
3512        #print "lmin",lmin;print "lmax",lmax
3513        latitudes_new = latitudes[kmin:kmax]
3514        longitudes_news = longitudes[lmin:lmax]
3515        #print "latitudes_new", latitudes_new
3516        #print "longitudes_news",longitudes_news
3517       
3518        self.failUnless(latitudes_new == [2, 1, 0] and \
3519                        longitudes_news == [0, 10, 20],
3520                         'failed')
3521        ## 5th test
3522        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3523            latitudes,longitudes,
3524            0.1,1.9,2,17)
3525        #print "kmin",kmin;print "kmax",kmax
3526        #print "lmin",lmin;print "lmax",lmax
3527        latitudes_new = latitudes[kmin:kmax]
3528        longitudes_news = longitudes[lmin:lmax]
3529        #print "latitudes_new", latitudes_new
3530        #print "longitudes_news",longitudes_news
3531       
3532        self.failUnless(latitudes_new == [2, 1, 0] and \
3533                        longitudes_news == [0, 10, 20],
3534                         'failed')
3535       
3536        ## 6th test
3537       
3538        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3539            latitudes,longitudes,
3540            1.5,4,18,32)
3541        #print "kmin",kmin;print "kmax",kmax
3542        #print "lmin",lmin;print "lmax",lmax
3543        latitudes_new = latitudes[kmin:kmax]
3544        longitudes_news = longitudes[lmin:lmax]
3545        #print "latitudes_new", latitudes_new
3546        #print "longitudes_news",longitudes_news
3547       
3548        self.failUnless(latitudes_new == [3, 2, 1] and \
3549                        longitudes_news == [10, 20, 30],
3550                         'failed')
3551       
3552       
3553        ## 7th test
3554        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3555        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3556            latitudes,longitudes,
3557            1.5,1.5,15,15)
3558        #print "kmin",kmin;print "kmax",kmax
3559        #print "lmin",lmin;print "lmax",lmax
3560        latitudes_new = latitudes[kmin:kmax]
3561        longitudes_news = longitudes[lmin:lmax]
3562        m2d = m2d[kmin:kmax,lmin:lmax]
3563        #print "m2d", m2d
3564        #print "latitudes_new", latitudes_new
3565        #print "longitudes_news",longitudes_news
3566       
3567        self.failUnless(latitudes_new == [2, 1] and \
3568                        longitudes_news == [10, 20],
3569                         'failed')
3570       
3571        self.failUnless(m2d == [[5,6],[9,10]],
3572                         'failed')
3573       
3574    def test_get_min_max_indexes2(self):
3575        latitudes = [-30,-35,-40,-45]
3576        longitudes = [148,149,150,151]
3577
3578        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3579       
3580        # k - lat
3581        # l - lon
3582        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3583            latitudes,longitudes,
3584            -37,-27,147,149.5)
3585
3586        #print "kmin",kmin;print "kmax",kmax
3587        #print "lmin",lmin;print "lmax",lmax
3588        #print "m2d", m2d
3589        #print "latitudes", latitudes
3590        #print "longitudes",longitudes
3591        #print "latitudes[kmax]", latitudes[kmax]
3592        latitudes_new = latitudes[kmin:kmax]
3593        longitudes_new = longitudes[lmin:lmax]
3594        m2d = m2d[kmin:kmax,lmin:lmax]
3595        #print "m2d", m2d
3596        #print "latitudes_new", latitudes_new
3597        #print "longitudes_new",longitudes_new
3598       
3599        self.failUnless(latitudes_new == [-30, -35, -40] and \
3600                        longitudes_new == [148, 149,150],
3601                         'failed')
3602        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
3603                         'failed')
3604       
3605    def test_get_min_max_indexes3(self):
3606        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3607        longitudes = [148,149,150,151]
3608
3609        # k - lat
3610        # l - lon
3611        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3612            latitudes,longitudes,
3613            -43,-37,148.5,149.5)
3614
3615       
3616        #print "kmin",kmin;print "kmax",kmax
3617        #print "lmin",lmin;print "lmax",lmax
3618        #print "latitudes", latitudes
3619        #print "longitudes",longitudes
3620        latitudes_new = latitudes[kmin:kmax]
3621        longitudes_news = longitudes[lmin:lmax]
3622        #print "latitudes_new", latitudes_new
3623        #print "longitudes_news",longitudes_news
3624       
3625        self.failUnless(latitudes_new == [-35, -40, -45] and \
3626                        longitudes_news == [148, 149,150],
3627                         'failed')
3628       
3629    def test_get_min_max_indexes4(self):
3630        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3631        longitudes = [148,149,150,151]
3632
3633        # k - lat
3634        # l - lon
3635        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3636            latitudes,longitudes)
3637       
3638        #print "kmin",kmin;print "kmax",kmax
3639        #print "lmin",lmin;print "lmax",lmax
3640        #print "latitudes", latitudes
3641        #print "longitudes",longitudes
3642        latitudes_new = latitudes[kmin:kmax]
3643        longitudes_news = longitudes[lmin:lmax]
3644        #print "latitudes_new", latitudes_new
3645        #print "longitudes_news",longitudes_news
3646       
3647        self.failUnless(latitudes_new == latitudes  and \
3648                        longitudes_news == longitudes,
3649                         'failed')
3650       
3651    def test_tsh2sww(self):
3652        import os
3653        import tempfile
3654       
3655        tsh_file = tempfile.mktemp(".tsh")
3656        file = open(tsh_file,"w")
3657        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
36580 0.0 0.0 0.0 0.0 0.01 \n \
36591 1.0 0.0 10.0 10.0 0.02  \n \
36602 0.0 1.0 0.0 10.0 0.03  \n \
36613 0.5 0.25 8.0 12.0 0.04  \n \
3662# Vert att title  \n \
3663elevation  \n \
3664stage  \n \
3665friction  \n \
36662 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
36670 0 3 2 -1  -1  1 dsg\n\
36681 0 1 3 -1  0 -1   ole nielsen\n\
36694 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
36700 1 0 2 \n\
36711 0 2 3 \n\
36722 2 3 \n\
36733 3 1 1 \n\
36743 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
36750 216.0 -86.0 \n \
36761 160.0 -167.0 \n \
36772 114.0 -91.0 \n \
36783 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
36790 0 1 0 \n \
36801 1 2 0 \n \
36812 2 0 0 \n \
36820 # <x> <y> ...Mesh Holes... \n \
36830 # <x> <y> <attribute>...Mesh Regions... \n \
36840 # <x> <y> <attribute>...Mesh Regions, area... \n\
3685#Geo reference \n \
368656 \n \
3687140 \n \
3688120 \n")
3689        file.close()
3690
3691        #sww_file = tempfile.mktemp(".sww")
3692        #print "sww_file",sww_file
3693        #print "sww_file",tsh_file
3694        tsh2sww(tsh_file)
3695
3696        os.remove(tsh_file)
3697        os.remove(tsh_file[:-4] + '.sww')
3698#-------------------------------------------------------------
3699if __name__ == "__main__":
3700    #suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
3701    suite = unittest.makeSuite(Test_Data_Manager,'test')
3702    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
3703    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')   
3704    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
3705    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
3706    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
3707    runner = unittest.TextTestRunner()
3708    runner.run(suite)
3709 
Note: See TracBrowser for help on using the repository browser.