source: inundation/pyvolution/test_data_manager.py @ 2022

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

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