source: inundation/pyvolution/test_data_manager.py @ 2060

Last change on this file since 2060 was 2060, checked in by ole, 18 years ago

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