source: inundation/pyvolution/test_data_manager.py @ 1919

Last change on this file since 1919 was 1919, checked in by ole, 19 years ago

Implemented arbitrary expressions for sww2dem using code from changeset:1916
File size: 79.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
8
9from data_manager import *
10from shallow_water import *
11from config import epsilon
12
13from coordinate_transforms.geo_reference import Geo_reference
14
15class Test_Data_Manager(unittest.TestCase):
16    def setUp(self):
17        import time
18        from mesh_factory import rectangular
19
20        #Create basic mesh
21        points, vertices, boundary = rectangular(2, 2)
22
23        #Create shallow water domain
24        domain = Domain(points, vertices, boundary)
25        domain.default_order=2
26
27
28        #Set some field values
29        domain.set_quantity('elevation', lambda x,y: -x)
30        domain.set_quantity('friction', 0.03)
31
32
33        ######################
34        # Boundary conditions
35        B = Transmissive_boundary(domain)
36        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
37
38
39        ######################
40        #Initial condition - with jumps
41
42
43        bed = domain.quantities['elevation'].vertex_values
44        stage = zeros(bed.shape, Float)
45
46        h = 0.3
47        for i in range(stage.shape[0]):
48            if i % 2 == 0:
49                stage[i,:] = bed[i,:] + h
50            else:
51                stage[i,:] = bed[i,:]
52
53        domain.set_quantity('stage', stage)
54        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
55
56        domain.distribute_to_vertices_and_edges()
57
58
59        self.domain = domain
60
61        C = domain.get_vertex_coordinates()
62        self.X = C[:,0:6:2].copy()
63        self.Y = C[:,1:6:2].copy()
64
65        self.F = bed
66
67
68    def tearDown(self):
69        pass
70
71
72
73
74#     def test_xya(self):
75#         import os
76#         from Numeric import concatenate
77
78#         import time, os
79#         from Numeric import array, zeros, allclose, Float, concatenate
80
81#         domain = self.domain
82
83#         domain.filename = 'datatest' + str(time.time())
84#         domain.format = 'xya'
85#         domain.smooth = True
86
87#         xya = get_dataobject(self.domain)
88#         xya.store_all()
89
90
91#         #Read back
92#         file = open(xya.filename)
93#         lFile = file.read().split('\n')
94#         lFile = lFile[:-1]
95
96#         file.close()
97#         os.remove(xya.filename)
98
99#         #Check contents
100#         if domain.smooth:
101#             self.failUnless(lFile[0] == '9 3 # <vertex #> <x> <y> [attributes]')
102#         else:
103#             self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]')
104
105#         #Get smoothed field values with X and Y
106#         X,Y,F,V = domain.get_vertex_values(xy=True, value_array='field_values',
107#                                            indices = (0,1), precision = Float)
108
109
110#         Q,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
111#                                            indices = (0,), precision = Float)
112
113
114
115#         for i, line in enumerate(lFile[1:]):
116#             fields = line.split()
117
118#             assert len(fields) == 5
119
120#             assert allclose(float(fields[0]), X[i])
121#             assert allclose(float(fields[1]), Y[i])
122#             assert allclose(float(fields[2]), F[i,0])
123#             assert allclose(float(fields[3]), Q[i,0])
124#             assert allclose(float(fields[4]), F[i,1])
125
126
127
128
129    def test_sww_constant(self):
130        """Test that constant sww information can be written correctly
131        (non smooth)
132        """
133
134        import time, os
135        from Numeric import array, zeros, allclose, Float, concatenate
136        from Scientific.IO.NetCDF import NetCDFFile
137
138        self.domain.filename = 'datatest' + str(id(self))
139        self.domain.format = 'sww'
140        self.domain.smooth = False
141
142        sww = get_dataobject(self.domain)
143        sww.store_connectivity()
144
145        #Check contents
146        #Get NetCDF
147        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
148
149        # Get the variables
150        x = fid.variables['x']
151        y = fid.variables['y']
152        z = fid.variables['elevation']
153
154        volumes = fid.variables['volumes']
155
156
157        assert allclose (x[:], self.X.flat)
158        assert allclose (y[:], self.Y.flat)
159        assert allclose (z[:], self.F.flat)
160
161        V = volumes
162
163        P = len(self.domain)
164        for k in range(P):
165            assert V[k, 0] == 3*k
166            assert V[k, 1] == 3*k+1
167            assert V[k, 2] == 3*k+2
168
169
170        fid.close()
171
172        #Cleanup
173        os.remove(sww.filename)
174
175
176    def test_sww_constant_smooth(self):
177        """Test that constant sww information can be written correctly
178        (non smooth)
179        """
180
181        import time, os
182        from Numeric import array, zeros, allclose, Float, concatenate
183        from Scientific.IO.NetCDF import NetCDFFile
184
185        self.domain.filename = 'datatest' + str(id(self))
186        self.domain.format = 'sww'
187        self.domain.smooth = True
188
189        sww = get_dataobject(self.domain)
190        sww.store_connectivity()
191
192        #Check contents
193        #Get NetCDF
194        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
195
196        # Get the variables
197        x = fid.variables['x']
198        y = fid.variables['y']
199        z = fid.variables['elevation']
200
201        volumes = fid.variables['volumes']
202
203        X = x[:]
204        Y = y[:]
205
206        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
207        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
208        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
209
210        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
211
212        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
213
214        Z = z[:]
215        assert Z[4] == -0.5
216
217        V = volumes
218        assert V[2,0] == 4
219        assert V[2,1] == 5
220        assert V[2,2] == 1
221
222        assert V[4,0] == 6
223        assert V[4,1] == 7
224        assert V[4,2] == 3
225
226
227        fid.close()
228
229        #Cleanup
230        os.remove(sww.filename)
231
232
233
234    def test_sww_variable(self):
235        """Test that sww information can be written correctly
236        """
237
238        import time, os
239        from Numeric import array, zeros, allclose, Float, concatenate
240        from Scientific.IO.NetCDF import NetCDFFile
241
242        self.domain.filename = 'datatest' + str(id(self))
243        self.domain.format = 'sww'
244        self.domain.smooth = True
245        self.domain.reduction = mean
246
247        sww = get_dataobject(self.domain)
248        sww.store_connectivity()
249        sww.store_timestep('stage')
250
251        #Check contents
252        #Get NetCDF
253        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
254
255
256        # Get the variables
257        x = fid.variables['x']
258        y = fid.variables['y']
259        z = fid.variables['elevation']
260        time = fid.variables['time']
261        stage = fid.variables['stage']
262
263
264        Q = self.domain.quantities['stage']
265        Q0 = Q.vertex_values[:,0]
266        Q1 = Q.vertex_values[:,1]
267        Q2 = Q.vertex_values[:,2]
268
269        A = stage[0,:]
270        #print A[0], (Q2[0,0] + Q1[1,0])/2
271        assert allclose(A[0], (Q2[0] + Q1[1])/2)
272        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
273        assert allclose(A[2], Q0[3])
274        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
275
276        #Center point
277        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
278                                 Q0[5] + Q2[6] + Q1[7])/6)
279
280
281
282        fid.close()
283
284        #Cleanup
285        os.remove(sww.filename)
286
287
288    def test_sww_variable2(self):
289        """Test that sww information can be written correctly
290        multiple timesteps. Use average as reduction operator
291        """
292
293        import time, os
294        from Numeric import array, zeros, allclose, Float, concatenate
295        from Scientific.IO.NetCDF import NetCDFFile
296
297        self.domain.filename = 'datatest' + str(id(self))
298        self.domain.format = 'sww'
299        self.domain.smooth = True
300
301        self.domain.reduction = mean
302
303        sww = get_dataobject(self.domain)
304        sww.store_connectivity()
305        sww.store_timestep('stage')
306        self.domain.evolve_to_end(finaltime = 0.01)
307        sww.store_timestep('stage')
308
309
310        #Check contents
311        #Get NetCDF
312        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
313
314        # Get the variables
315        x = fid.variables['x']
316        y = fid.variables['y']
317        z = fid.variables['elevation']
318        time = fid.variables['time']
319        stage = fid.variables['stage']
320
321        #Check values
322        Q = self.domain.quantities['stage']
323        Q0 = Q.vertex_values[:,0]
324        Q1 = Q.vertex_values[:,1]
325        Q2 = Q.vertex_values[:,2]
326
327        A = stage[1,:]
328        assert allclose(A[0], (Q2[0] + Q1[1])/2)
329        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
330        assert allclose(A[2], Q0[3])
331        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
332
333        #Center point
334        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
335                                 Q0[5] + Q2[6] + Q1[7])/6)
336
337
338        fid.close()
339
340        #Cleanup
341        os.remove(sww.filename)
342
343    def test_sww_variable3(self):
344        """Test that sww information can be written correctly
345        multiple timesteps using a different reduction operator (min)
346        """
347
348        import time, os
349        from Numeric import array, zeros, allclose, Float, concatenate
350        from Scientific.IO.NetCDF import NetCDFFile
351
352        self.domain.filename = 'datatest' + str(id(self))
353        self.domain.format = 'sww'
354        self.domain.smooth = True
355        self.domain.reduction = min
356
357        sww = get_dataobject(self.domain)
358        sww.store_connectivity()
359        sww.store_timestep('stage')
360
361        self.domain.evolve_to_end(finaltime = 0.01)
362        sww.store_timestep('stage')
363
364
365        #Check contents
366        #Get NetCDF
367        fid = NetCDFFile(sww.filename, 'r')
368
369
370        # Get the variables
371        x = fid.variables['x']
372        y = fid.variables['y']
373        z = fid.variables['elevation']
374        time = fid.variables['time']
375        stage = fid.variables['stage']
376
377        #Check values
378        Q = self.domain.quantities['stage']
379        Q0 = Q.vertex_values[:,0]
380        Q1 = Q.vertex_values[:,1]
381        Q2 = Q.vertex_values[:,2]
382
383        A = stage[1,:]
384        assert allclose(A[0], min(Q2[0], Q1[1]))
385        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
386        assert allclose(A[2], Q0[3])
387        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
388
389        #Center point
390        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
391                                  Q0[5], Q2[6], Q1[7]))
392
393
394        fid.close()
395
396        #Cleanup
397        os.remove(sww.filename)
398
399
400    def test_sync(self):
401        """Test info stored at each timestep is as expected (incl initial condition)
402        """
403
404        import time, os, config
405        from Numeric import array, zeros, allclose, Float, concatenate
406        from Scientific.IO.NetCDF import NetCDFFile
407
408        self.domain.filename = 'synctest'
409        self.domain.format = 'sww'
410        self.domain.smooth = False
411        self.domain.store = True
412        self.domain.beta_h = 0
413
414        #Evolution
415        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
416            stage = self.domain.quantities['stage'].vertex_values
417
418            #Get NetCDF
419            fid = NetCDFFile(self.domain.writer.filename, 'r')
420            stage_file = fid.variables['stage']
421
422            if t == 0.0:
423                assert allclose(stage, self.initial_stage)
424                assert allclose(stage_file[:], stage.flat)
425            else:
426                assert not allclose(stage, self.initial_stage)
427                assert not allclose(stage_file[:], stage.flat)
428
429            fid.close()
430
431        os.remove(self.domain.writer.filename)
432
433
434
435    def test_sww_DSG(self):
436        """Not a test, rather a look at the sww format
437        """
438
439        import time, os
440        from Numeric import array, zeros, allclose, Float, concatenate
441        from Scientific.IO.NetCDF import NetCDFFile
442
443        self.domain.filename = 'datatest' + str(id(self))
444        self.domain.format = 'sww'
445        self.domain.smooth = True
446        self.domain.reduction = mean
447
448        sww = get_dataobject(self.domain)
449        sww.store_connectivity()
450        sww.store_timestep('stage')
451
452        #Check contents
453        #Get NetCDF
454        fid = NetCDFFile(sww.filename, 'r')
455
456        # Get the variables
457        x = fid.variables['x']
458        y = fid.variables['y']
459        z = fid.variables['elevation']
460
461        volumes = fid.variables['volumes']
462        time = fid.variables['time']
463
464        # 2D
465        stage = fid.variables['stage']
466
467        X = x[:]
468        Y = y[:]
469        Z = z[:]
470        V = volumes[:]
471        T = time[:]
472        S = stage[:,:]
473
474#         print "****************************"
475#         print "X ",X
476#         print "****************************"
477#         print "Y ",Y
478#         print "****************************"
479#         print "Z ",Z
480#         print "****************************"
481#         print "V ",V
482#         print "****************************"
483#         print "Time ",T
484#         print "****************************"
485#         print "Stage ",S
486#         print "****************************"
487
488
489        fid.close()
490
491        #Cleanup
492        os.remove(sww.filename)
493
494
495    def test_write_pts(self):
496        #Get (enough) datapoints
497
498        from Numeric import array
499        points = array([[ 0.66666667, 0.66666667],
500                        [ 1.33333333, 1.33333333],
501                        [ 2.66666667, 0.66666667],
502                        [ 0.66666667, 2.66666667],
503                        [ 0.0, 1.0],
504                        [ 0.0, 3.0],
505                        [ 1.0, 0.0],
506                        [ 1.0, 1.0],
507                        [ 1.0, 2.0],
508                        [ 1.0, 3.0],
509                        [ 2.0, 1.0],
510                        [ 3.0, 0.0],
511                        [ 3.0, 1.0]])
512
513        z = points[:,0] + 2*points[:,1]
514
515        ptsfile = 'testptsfile.pts'
516        write_ptsfile(ptsfile, points, z,
517                      attribute_name = 'linear_combination')
518
519        #Check contents
520        #Get NetCDF
521        from Scientific.IO.NetCDF import NetCDFFile       
522        fid = NetCDFFile(ptsfile, 'r')
523
524        # Get the variables
525        #print fid.variables.keys()
526        points1 = fid.variables['points']
527        z1 = fid.variables['linear_combination']
528
529        #Check values
530
531        #print points[:]
532        #print ref_points
533        assert allclose(points, points1)
534
535        #print attributes[:]
536        #print ref_elevation
537        assert allclose(z, z1)
538
539        #Cleanup
540        fid.close()
541
542        import os
543        os.remove(ptsfile)
544       
545
546
547
548    def test_dem2pts(self):
549        """Test conversion from dem in ascii format to native NetCDF xya format
550        """
551
552        import time, os
553        from Numeric import array, zeros, allclose, Float, concatenate
554        from Scientific.IO.NetCDF import NetCDFFile
555
556        #Write test asc file
557        root = 'demtest'
558
559        filename = root+'.asc'
560        fid = open(filename, 'w')
561        fid.write("""ncols         5
562nrows         6
563xllcorner     2000.5
564yllcorner     3000.5
565cellsize      25
566NODATA_value  -9999
567""")
568        #Create linear function
569
570        ref_points = []
571        ref_elevation = []
572        for i in range(6):
573            y = (6-i)*25.0
574            for j in range(5):
575                x = j*25.0
576                z = x+2*y
577
578                ref_points.append( [x,y] )
579                ref_elevation.append(z)
580                fid.write('%f ' %z)
581            fid.write('\n')
582
583        fid.close()
584
585        #Write prj file with metadata
586        metafilename = root+'.prj'
587        fid = open(metafilename, 'w')
588
589
590        fid.write("""Projection UTM
591Zone 56
592Datum WGS84
593Zunits NO
594Units METERS
595Spheroid WGS84
596Xshift 0.0000000000
597Yshift 10000000.0000000000
598Parameters
599""")
600        fid.close()
601
602        #Convert to NetCDF pts
603        convert_dem_from_ascii2netcdf(root)
604        dem2pts(root)
605
606        #Check contents
607        #Get NetCDF
608        fid = NetCDFFile(root+'.pts', 'r')
609
610        # Get the variables
611        #print fid.variables.keys()
612        points = fid.variables['points']
613        elevation = fid.variables['elevation']
614
615        #Check values
616
617        #print points[:]
618        #print ref_points
619        assert allclose(points, ref_points)
620
621        #print attributes[:]
622        #print ref_elevation
623        assert allclose(elevation, ref_elevation)
624
625        #Cleanup
626        fid.close()
627
628
629        os.remove(root + '.pts')
630        os.remove(root + '.dem')
631        os.remove(root + '.asc')
632        os.remove(root + '.prj')
633
634
635
636    def test_dem2pts_bounding_box(self):
637        """Test conversion from dem in ascii format to native NetCDF xya format
638        """
639
640        import time, os
641        from Numeric import array, zeros, allclose, Float, concatenate
642        from Scientific.IO.NetCDF import NetCDFFile
643
644        #Write test asc file
645        root = 'demtest'
646
647        filename = root+'.asc'
648        fid = open(filename, 'w')
649        fid.write("""ncols         5
650nrows         6
651xllcorner     2000.5
652yllcorner     3000.5
653cellsize      25
654NODATA_value  -9999
655""")
656        #Create linear function
657
658        ref_points = []
659        ref_elevation = []
660        for i in range(6):
661            y = (6-i)*25.0
662            for j in range(5):
663                x = j*25.0
664                z = x+2*y
665
666                ref_points.append( [x,y] )
667                ref_elevation.append(z)
668                fid.write('%f ' %z)
669            fid.write('\n')
670
671        fid.close()
672
673        #Write prj file with metadata
674        metafilename = root+'.prj'
675        fid = open(metafilename, 'w')
676
677
678        fid.write("""Projection UTM
679Zone 56
680Datum WGS84
681Zunits NO
682Units METERS
683Spheroid WGS84
684Xshift 0.0000000000
685Yshift 10000000.0000000000
686Parameters
687""")
688        fid.close()
689
690        #Convert to NetCDF pts
691        convert_dem_from_ascii2netcdf(root)
692        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
693                northing_min=3035.0, northing_max=3125.5)
694
695        #Check contents
696        #Get NetCDF
697        fid = NetCDFFile(root+'.pts', 'r')
698
699        # Get the variables
700        #print fid.variables.keys()
701        points = fid.variables['points']
702        elevation = fid.variables['elevation']
703
704        #Check values
705        assert fid.xllcorner[0] == 2010.0
706        assert fid.yllcorner[0] == 3035.0
707
708        #create new reference points
709        ref_points = []
710        ref_elevation = []
711        for i in range(4):
712            y = (4-i)*25.0 + 25.0
713            y_new = y + 3000.5 - 3035.0
714            for j in range(4):
715                x = j*25.0 + 25.0
716                x_new = x + 2000.5 - 2010.0
717                z = x+2*y
718
719                ref_points.append( [x_new,y_new] )
720                ref_elevation.append(z)
721
722        #print points[:]
723        #print ref_points
724        assert allclose(points, ref_points)
725
726        #print attributes[:]
727        #print ref_elevation
728        assert allclose(elevation, ref_elevation)
729
730        #Cleanup
731        fid.close()
732
733
734        os.remove(root + '.pts')
735        os.remove(root + '.dem')
736        os.remove(root + '.asc')
737        os.remove(root + '.prj')
738
739
740
741    def test_sww2dem_asc_elevation(self):
742        """Test that sww information can be converted correctly to asc/prj
743        format readable by e.g. ArcView
744        """
745
746        import time, os
747        from Numeric import array, zeros, allclose, Float, concatenate
748        from Scientific.IO.NetCDF import NetCDFFile
749
750        #Setup
751        self.domain.filename = 'datatest'
752
753        prjfile = self.domain.filename + '_elevation.prj'
754        ascfile = self.domain.filename + '_elevation.asc'
755        swwfile = self.domain.filename + '.sww'
756
757        self.domain.set_datadir('.')
758        self.domain.format = 'sww'
759        self.domain.smooth = True
760        self.domain.set_quantity('elevation', lambda x,y: -x-y)
761
762        self.domain.geo_reference = Geo_reference(56,308500,6189000)
763
764        sww = get_dataobject(self.domain)
765        sww.store_connectivity()
766        sww.store_timestep('stage')
767
768        self.domain.evolve_to_end(finaltime = 0.01)
769        sww.store_timestep('stage')
770
771        cellsize = 0.25
772        #Check contents
773        #Get NetCDF
774
775        fid = NetCDFFile(sww.filename, 'r')
776
777        # Get the variables
778        x = fid.variables['x'][:]
779        y = fid.variables['y'][:]
780        z = fid.variables['elevation'][:]
781        time = fid.variables['time'][:]
782        stage = fid.variables['stage'][:]
783
784
785        #Export to ascii/prj files
786        sww2dem(self.domain.filename,
787                quantity = 'elevation',
788                cellsize = cellsize,
789                verbose = False,
790                format = 'asc')
791               
792        #Check prj (meta data)
793        prjid = open(prjfile)
794        lines = prjid.readlines()
795        prjid.close()
796
797        L = lines[0].strip().split()
798        assert L[0].strip().lower() == 'projection'
799        assert L[1].strip().lower() == 'utm'
800
801        L = lines[1].strip().split()
802        assert L[0].strip().lower() == 'zone'
803        assert L[1].strip().lower() == '56'
804
805        L = lines[2].strip().split()
806        assert L[0].strip().lower() == 'datum'
807        assert L[1].strip().lower() == 'wgs84'
808
809        L = lines[3].strip().split()
810        assert L[0].strip().lower() == 'zunits'
811        assert L[1].strip().lower() == 'no'
812
813        L = lines[4].strip().split()
814        assert L[0].strip().lower() == 'units'
815        assert L[1].strip().lower() == 'meters'
816
817        L = lines[5].strip().split()
818        assert L[0].strip().lower() == 'spheroid'
819        assert L[1].strip().lower() == 'wgs84'
820
821        L = lines[6].strip().split()
822        assert L[0].strip().lower() == 'xshift'
823        assert L[1].strip().lower() == '500000'
824
825        L = lines[7].strip().split()
826        assert L[0].strip().lower() == 'yshift'
827        assert L[1].strip().lower() == '10000000'
828
829        L = lines[8].strip().split()
830        assert L[0].strip().lower() == 'parameters'
831
832
833        #Check asc file
834        ascid = open(ascfile)
835        lines = ascid.readlines()
836        ascid.close()
837
838        L = lines[0].strip().split()
839        assert L[0].strip().lower() == 'ncols'
840        assert L[1].strip().lower() == '5'
841
842        L = lines[1].strip().split()
843        assert L[0].strip().lower() == 'nrows'
844        assert L[1].strip().lower() == '5'
845
846        L = lines[2].strip().split()
847        assert L[0].strip().lower() == 'xllcorner'
848        assert allclose(float(L[1].strip().lower()), 308500)
849
850        L = lines[3].strip().split()
851        assert L[0].strip().lower() == 'yllcorner'
852        assert allclose(float(L[1].strip().lower()), 6189000)
853
854        L = lines[4].strip().split()
855        assert L[0].strip().lower() == 'cellsize'
856        assert allclose(float(L[1].strip().lower()), cellsize)
857
858        L = lines[5].strip().split()
859        assert L[0].strip() == 'NODATA_value'
860        assert L[1].strip().lower() == '-9999'
861
862        #Check grid values
863        for j in range(5):
864            L = lines[6+j].strip().split()
865            y = (4-j) * cellsize
866            for i in range(5):
867                assert allclose(float(L[i]), -i*cellsize - y)
868
869
870        fid.close()
871
872        #Cleanup
873        os.remove(prjfile)
874        os.remove(ascfile)
875        os.remove(swwfile)
876
877
878
879    def test_sww2dem_larger(self):
880        """Test that sww information can be converted correctly to asc/prj
881        format readable by e.g. ArcView. Here:
882
883        ncols         11
884        nrows         11
885        xllcorner     308500
886        yllcorner     6189000
887        cellsize      10.000000
888        NODATA_value  -9999
889        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
890         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
891         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
892         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
893         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
894         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
895         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
896         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
897         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
898         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
899           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
900       
901        """
902
903        import time, os
904        from Numeric import array, zeros, allclose, Float, concatenate
905        from Scientific.IO.NetCDF import NetCDFFile
906
907        #Setup
908
909        from mesh_factory import rectangular
910
911        #Create basic mesh (100m x 100m)
912        points, vertices, boundary = rectangular(2, 2, 100, 100)
913
914        #Create shallow water domain
915        domain = Domain(points, vertices, boundary)
916        domain.default_order = 2
917
918        domain.filename = 'datatest'
919
920        prjfile = domain.filename + '_elevation.prj'
921        ascfile = domain.filename + '_elevation.asc'
922        swwfile = domain.filename + '.sww'
923
924        domain.set_datadir('.')
925        domain.format = 'sww'
926        domain.smooth = True
927        domain.geo_reference = Geo_reference(56, 308500, 6189000)
928
929        #
930        domain.set_quantity('elevation', lambda x,y: -x-y)
931        domain.set_quantity('stage', 0)       
932
933        B = Transmissive_boundary(domain)
934        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
935
936
937        #
938        sww = get_dataobject(domain)
939        sww.store_connectivity()
940        sww.store_timestep('stage')
941
942        domain.evolve_to_end(finaltime = 0.01)
943        sww.store_timestep('stage')
944
945        cellsize = 10  #10m grid
946
947       
948        #Check contents
949        #Get NetCDF
950
951        fid = NetCDFFile(sww.filename, 'r')
952
953        # Get the variables
954        x = fid.variables['x'][:]
955        y = fid.variables['y'][:]
956        z = fid.variables['elevation'][:]
957        time = fid.variables['time'][:]
958        stage = fid.variables['stage'][:]
959
960
961        #Export to ascii/prj files
962        sww2dem(domain.filename,
963                quantity = 'elevation',
964                cellsize = cellsize,
965                verbose = False,
966                format = 'asc')
967       
968               
969        #Check prj (meta data)
970        prjid = open(prjfile)
971        lines = prjid.readlines()
972        prjid.close()
973
974        L = lines[0].strip().split()
975        assert L[0].strip().lower() == 'projection'
976        assert L[1].strip().lower() == 'utm'
977
978        L = lines[1].strip().split()
979        assert L[0].strip().lower() == 'zone'
980        assert L[1].strip().lower() == '56'
981
982        L = lines[2].strip().split()
983        assert L[0].strip().lower() == 'datum'
984        assert L[1].strip().lower() == 'wgs84'
985
986        L = lines[3].strip().split()
987        assert L[0].strip().lower() == 'zunits'
988        assert L[1].strip().lower() == 'no'
989
990        L = lines[4].strip().split()
991        assert L[0].strip().lower() == 'units'
992        assert L[1].strip().lower() == 'meters'
993
994        L = lines[5].strip().split()
995        assert L[0].strip().lower() == 'spheroid'
996        assert L[1].strip().lower() == 'wgs84'
997
998        L = lines[6].strip().split()
999        assert L[0].strip().lower() == 'xshift'
1000        assert L[1].strip().lower() == '500000'
1001
1002        L = lines[7].strip().split()
1003        assert L[0].strip().lower() == 'yshift'
1004        assert L[1].strip().lower() == '10000000'
1005
1006        L = lines[8].strip().split()
1007        assert L[0].strip().lower() == 'parameters'
1008
1009
1010        #Check asc file
1011        ascid = open(ascfile)
1012        lines = ascid.readlines()
1013        ascid.close()
1014
1015        L = lines[0].strip().split()
1016        assert L[0].strip().lower() == 'ncols'
1017        assert L[1].strip().lower() == '11'
1018
1019        L = lines[1].strip().split()
1020        assert L[0].strip().lower() == 'nrows'
1021        assert L[1].strip().lower() == '11'
1022
1023        L = lines[2].strip().split()
1024        assert L[0].strip().lower() == 'xllcorner'
1025        assert allclose(float(L[1].strip().lower()), 308500)
1026
1027        L = lines[3].strip().split()
1028        assert L[0].strip().lower() == 'yllcorner'
1029        assert allclose(float(L[1].strip().lower()), 6189000)
1030
1031        L = lines[4].strip().split()
1032        assert L[0].strip().lower() == 'cellsize'
1033        assert allclose(float(L[1].strip().lower()), cellsize)
1034
1035        L = lines[5].strip().split()
1036        assert L[0].strip() == 'NODATA_value'
1037        assert L[1].strip().lower() == '-9999'
1038
1039        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1040        for i, line in enumerate(lines[6:]):
1041            for j, value in enumerate( line.split() ):
1042                #assert allclose(float(value), -(10-i+j)*cellsize)
1043                assert float(value) == -(10-i+j)*cellsize
1044       
1045
1046        fid.close()
1047
1048        #Cleanup
1049        os.remove(prjfile)
1050        os.remove(ascfile)
1051        os.remove(swwfile)
1052
1053
1054
1055    def test_sww2dem_boundingbox(self):
1056        """Test that sww information can be converted correctly to asc/prj
1057        format readable by e.g. ArcView.
1058        This will test that mesh can be restricted by bounding box
1059
1060        Original extent is 100m x 100m:
1061
1062        Eastings:   308500 -  308600
1063        Northings: 6189000 - 6189100
1064
1065        Bounding box changes this to the 50m x 50m square defined by
1066
1067        Eastings:   308530 -  308570
1068        Northings: 6189050 - 6189100
1069
1070        The cropped values should be
1071
1072         -130 -140 -150 -160 -170 
1073         -120 -130 -140 -150 -160 
1074         -110 -120 -130 -140 -150 
1075         -100 -110 -120 -130 -140 
1076          -90 -100 -110 -120 -130 
1077          -80  -90 -100 -110 -120
1078
1079        and the new lower reference point should be   
1080        Eastings:   308530
1081        Northings: 6189050
1082       
1083        Original dataset is the same as in test_sww2dem_larger()
1084       
1085        """
1086
1087        import time, os
1088        from Numeric import array, zeros, allclose, Float, concatenate
1089        from Scientific.IO.NetCDF import NetCDFFile
1090
1091        #Setup
1092
1093        from mesh_factory import rectangular
1094
1095        #Create basic mesh (100m x 100m)
1096        points, vertices, boundary = rectangular(2, 2, 100, 100)
1097
1098        #Create shallow water domain
1099        domain = Domain(points, vertices, boundary)
1100        domain.default_order = 2
1101
1102        domain.filename = 'datatest'
1103
1104        prjfile = domain.filename + '_elevation.prj'
1105        ascfile = domain.filename + '_elevation.asc'
1106        swwfile = domain.filename + '.sww'
1107
1108        domain.set_datadir('.')
1109        domain.format = 'sww'
1110        domain.smooth = True
1111        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1112
1113        #
1114        domain.set_quantity('elevation', lambda x,y: -x-y)
1115        domain.set_quantity('stage', 0)       
1116
1117        B = Transmissive_boundary(domain)
1118        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1119
1120
1121        #
1122        sww = get_dataobject(domain)
1123        sww.store_connectivity()
1124        sww.store_timestep('stage')
1125
1126        domain.evolve_to_end(finaltime = 0.01)
1127        sww.store_timestep('stage')
1128
1129        cellsize = 10  #10m grid
1130
1131       
1132        #Check contents
1133        #Get NetCDF
1134
1135        fid = NetCDFFile(sww.filename, 'r')
1136
1137        # Get the variables
1138        x = fid.variables['x'][:]
1139        y = fid.variables['y'][:]
1140        z = fid.variables['elevation'][:]
1141        time = fid.variables['time'][:]
1142        stage = fid.variables['stage'][:]
1143
1144
1145        #Export to ascii/prj files
1146        sww2dem(domain.filename,
1147                quantity = 'elevation',
1148                cellsize = cellsize,
1149                easting_min = 308530,
1150                easting_max = 308570,
1151                northing_min = 6189050, 
1152                northing_max = 6189100,                               
1153                verbose = False,
1154                format = 'asc')
1155       
1156        fid.close()
1157
1158       
1159        #Check prj (meta data)
1160        prjid = open(prjfile)
1161        lines = prjid.readlines()
1162        prjid.close()
1163
1164        L = lines[0].strip().split()
1165        assert L[0].strip().lower() == 'projection'
1166        assert L[1].strip().lower() == 'utm'
1167
1168        L = lines[1].strip().split()
1169        assert L[0].strip().lower() == 'zone'
1170        assert L[1].strip().lower() == '56'
1171
1172        L = lines[2].strip().split()
1173        assert L[0].strip().lower() == 'datum'
1174        assert L[1].strip().lower() == 'wgs84'
1175
1176        L = lines[3].strip().split()
1177        assert L[0].strip().lower() == 'zunits'
1178        assert L[1].strip().lower() == 'no'
1179
1180        L = lines[4].strip().split()
1181        assert L[0].strip().lower() == 'units'
1182        assert L[1].strip().lower() == 'meters'
1183
1184        L = lines[5].strip().split()
1185        assert L[0].strip().lower() == 'spheroid'
1186        assert L[1].strip().lower() == 'wgs84'
1187
1188        L = lines[6].strip().split()
1189        assert L[0].strip().lower() == 'xshift'
1190        assert L[1].strip().lower() == '500000'
1191
1192        L = lines[7].strip().split()
1193        assert L[0].strip().lower() == 'yshift'
1194        assert L[1].strip().lower() == '10000000'
1195
1196        L = lines[8].strip().split()
1197        assert L[0].strip().lower() == 'parameters'
1198
1199
1200        #Check asc file
1201        ascid = open(ascfile)
1202        lines = ascid.readlines()
1203        ascid.close()
1204
1205        L = lines[0].strip().split()
1206        assert L[0].strip().lower() == 'ncols'
1207        assert L[1].strip().lower() == '5'
1208
1209        L = lines[1].strip().split()
1210        assert L[0].strip().lower() == 'nrows'
1211        assert L[1].strip().lower() == '6'
1212
1213        L = lines[2].strip().split()
1214        assert L[0].strip().lower() == 'xllcorner'
1215        assert allclose(float(L[1].strip().lower()), 308530)
1216
1217        L = lines[3].strip().split()
1218        assert L[0].strip().lower() == 'yllcorner'
1219        assert allclose(float(L[1].strip().lower()), 6189050)
1220
1221        L = lines[4].strip().split()
1222        assert L[0].strip().lower() == 'cellsize'
1223        assert allclose(float(L[1].strip().lower()), cellsize)
1224
1225        L = lines[5].strip().split()
1226        assert L[0].strip() == 'NODATA_value'
1227        assert L[1].strip().lower() == '-9999'
1228
1229        #Check grid values
1230        for i, line in enumerate(lines[6:]):
1231            for j, value in enumerate( line.split() ):
1232                #assert float(value) == -(10-i+j)*cellsize               
1233                assert float(value) == -(10-i+j+3)*cellsize
1234       
1235
1236
1237        #Cleanup
1238        os.remove(prjfile)
1239        os.remove(ascfile)
1240        os.remove(swwfile)
1241
1242
1243
1244    def test_sww2dem_asc_stage_reduction(self):
1245        """Test that sww information can be converted correctly to asc/prj
1246        format readable by e.g. ArcView
1247
1248        This tests the reduction of quantity stage using min
1249        """
1250
1251        import time, os
1252        from Numeric import array, zeros, allclose, Float, concatenate
1253        from Scientific.IO.NetCDF import NetCDFFile
1254
1255        #Setup
1256        self.domain.filename = 'datatest'
1257
1258        prjfile = self.domain.filename + '_stage.prj'
1259        ascfile = self.domain.filename + '_stage.asc'
1260        swwfile = self.domain.filename + '.sww'
1261
1262        self.domain.set_datadir('.')
1263        self.domain.format = 'sww'
1264        self.domain.smooth = True
1265        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1266
1267        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1268
1269
1270        sww = get_dataobject(self.domain)
1271        sww.store_connectivity()
1272        sww.store_timestep('stage')
1273
1274        self.domain.evolve_to_end(finaltime = 0.01)
1275        sww.store_timestep('stage')
1276
1277        cellsize = 0.25
1278        #Check contents
1279        #Get NetCDF
1280
1281        fid = NetCDFFile(sww.filename, 'r')
1282
1283        # Get the variables
1284        x = fid.variables['x'][:]
1285        y = fid.variables['y'][:]
1286        z = fid.variables['elevation'][:]
1287        time = fid.variables['time'][:]
1288        stage = fid.variables['stage'][:]
1289
1290
1291        #Export to ascii/prj files
1292        sww2dem(self.domain.filename,
1293                quantity = 'stage',
1294                cellsize = cellsize,
1295                reduction = min,
1296                format = 'asc')
1297
1298
1299        #Check asc file
1300        ascid = open(ascfile)
1301        lines = ascid.readlines()
1302        ascid.close()
1303
1304        L = lines[0].strip().split()
1305        assert L[0].strip().lower() == 'ncols'
1306        assert L[1].strip().lower() == '5'
1307
1308        L = lines[1].strip().split()
1309        assert L[0].strip().lower() == 'nrows'
1310        assert L[1].strip().lower() == '5'
1311
1312        L = lines[2].strip().split()
1313        assert L[0].strip().lower() == 'xllcorner'
1314        assert allclose(float(L[1].strip().lower()), 308500)
1315
1316        L = lines[3].strip().split()
1317        assert L[0].strip().lower() == 'yllcorner'
1318        assert allclose(float(L[1].strip().lower()), 6189000)
1319
1320        L = lines[4].strip().split()
1321        assert L[0].strip().lower() == 'cellsize'
1322        assert allclose(float(L[1].strip().lower()), cellsize)
1323
1324        L = lines[5].strip().split()
1325        assert L[0].strip() == 'NODATA_value'
1326        assert L[1].strip().lower() == '-9999'
1327
1328
1329        #Check grid values (where applicable)
1330        for j in range(5):
1331            if j%2 == 0:
1332                L = lines[6+j].strip().split()
1333                jj = 4-j
1334                for i in range(5):
1335                    if i%2 == 0:
1336                        index = jj/2 + i/2*3
1337                        val0 = stage[0,index]
1338                        val1 = stage[1,index]
1339
1340                        #print i, j, index, ':', L[i], val0, val1
1341                        assert allclose(float(L[i]), min(val0, val1))
1342
1343
1344        fid.close()
1345
1346        #Cleanup
1347        os.remove(prjfile)
1348        os.remove(ascfile)
1349        #os.remove(swwfile)
1350
1351
1352
1353    def test_sww2dem_asc_derived_quantity(self):
1354        """Test that sww information can be converted correctly to asc/prj
1355        format readable by e.g. ArcView
1356
1357        This tests the use of derived quantities
1358        """
1359
1360        import time, os
1361        from Numeric import array, zeros, allclose, Float, concatenate
1362        from Scientific.IO.NetCDF import NetCDFFile
1363
1364        #Setup
1365        self.domain.filename = 'datatest'
1366
1367        prjfile = self.domain.filename + '_depth.prj'
1368        ascfile = self.domain.filename + '_depth.asc'
1369        swwfile = self.domain.filename + '.sww'
1370
1371        self.domain.set_datadir('.')
1372        self.domain.format = 'sww'
1373        self.domain.smooth = True
1374        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1375        self.domain.set_quantity('stage', 0.0)
1376
1377        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1378
1379
1380        sww = get_dataobject(self.domain)
1381        sww.store_connectivity()
1382        sww.store_timestep('stage')
1383
1384        self.domain.evolve_to_end(finaltime = 0.01)
1385        sww.store_timestep('stage')
1386
1387        cellsize = 0.25
1388        #Check contents
1389        #Get NetCDF
1390
1391        fid = NetCDFFile(sww.filename, 'r')
1392
1393        # Get the variables
1394        x = fid.variables['x'][:]
1395        y = fid.variables['y'][:]
1396        z = fid.variables['elevation'][:]
1397        time = fid.variables['time'][:]
1398        stage = fid.variables['stage'][:]
1399
1400
1401        #Export to ascii/prj files
1402        sww2dem(self.domain.filename,
1403                basename_out = 'datatest_depth',
1404                quantity = 'stage - elevation',
1405                cellsize = cellsize,
1406                reduction = min,
1407                format = 'asc',
1408                verbose = False)
1409
1410
1411        #Check asc file
1412        ascid = open(ascfile)
1413        lines = ascid.readlines()
1414        ascid.close()
1415
1416        L = lines[0].strip().split()
1417        assert L[0].strip().lower() == 'ncols'
1418        assert L[1].strip().lower() == '5'
1419
1420        L = lines[1].strip().split()
1421        assert L[0].strip().lower() == 'nrows'
1422        assert L[1].strip().lower() == '5'
1423
1424        L = lines[2].strip().split()
1425        assert L[0].strip().lower() == 'xllcorner'
1426        assert allclose(float(L[1].strip().lower()), 308500)
1427
1428        L = lines[3].strip().split()
1429        assert L[0].strip().lower() == 'yllcorner'
1430        assert allclose(float(L[1].strip().lower()), 6189000)
1431
1432        L = lines[4].strip().split()
1433        assert L[0].strip().lower() == 'cellsize'
1434        assert allclose(float(L[1].strip().lower()), cellsize)
1435
1436        L = lines[5].strip().split()
1437        assert L[0].strip() == 'NODATA_value'
1438        assert L[1].strip().lower() == '-9999'
1439
1440
1441        #Check grid values (where applicable)
1442        for j in range(5):
1443            if j%2 == 0:
1444                L = lines[6+j].strip().split()
1445                jj = 4-j
1446                for i in range(5):
1447                    if i%2 == 0:
1448                        index = jj/2 + i/2*3
1449                        val0 = stage[0,index] - z[index]
1450                        val1 = stage[1,index] - z[index]
1451
1452                        #print i, j, index, ':', L[i], val0, val1
1453                        assert allclose(float(L[i]), min(val0, val1))
1454
1455
1456        fid.close()
1457
1458        #Cleanup
1459        os.remove(prjfile)
1460        os.remove(ascfile)
1461        #os.remove(swwfile)
1462
1463
1464
1465
1466
1467    def test_sww2dem_asc_missing_points(self):
1468        """Test that sww information can be converted correctly to asc/prj
1469        format readable by e.g. ArcView
1470
1471        This test includes the writing of missing values
1472        """
1473
1474        import time, os
1475        from Numeric import array, zeros, allclose, Float, concatenate
1476        from Scientific.IO.NetCDF import NetCDFFile
1477
1478        #Setup mesh not coinciding with rectangle.
1479        #This will cause missing values to occur in gridded data
1480
1481
1482        points = [                        [1.0, 1.0],
1483                              [0.5, 0.5], [1.0, 0.5],
1484                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1485
1486        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1487
1488        #Create shallow water domain
1489        domain = Domain(points, vertices)
1490        domain.default_order=2
1491
1492
1493        #Set some field values
1494        domain.set_quantity('elevation', lambda x,y: -x-y)
1495        domain.set_quantity('friction', 0.03)
1496
1497
1498        ######################
1499        # Boundary conditions
1500        B = Transmissive_boundary(domain)
1501        domain.set_boundary( {'exterior': B} )
1502
1503
1504        ######################
1505        #Initial condition - with jumps
1506
1507        bed = domain.quantities['elevation'].vertex_values
1508        stage = zeros(bed.shape, Float)
1509
1510        h = 0.3
1511        for i in range(stage.shape[0]):
1512            if i % 2 == 0:
1513                stage[i,:] = bed[i,:] + h
1514            else:
1515                stage[i,:] = bed[i,:]
1516
1517        domain.set_quantity('stage', stage)
1518        domain.distribute_to_vertices_and_edges()
1519
1520        domain.filename = 'datatest'
1521
1522        prjfile = domain.filename + '_elevation.prj'
1523        ascfile = domain.filename + '_elevation.asc'
1524        swwfile = domain.filename + '.sww'
1525
1526        domain.set_datadir('.')
1527        domain.format = 'sww'
1528        domain.smooth = True
1529
1530        domain.geo_reference = Geo_reference(56,308500,6189000)
1531
1532        sww = get_dataobject(domain)
1533        sww.store_connectivity()
1534        sww.store_timestep('stage')
1535
1536        cellsize = 0.25
1537        #Check contents
1538        #Get NetCDF
1539
1540        fid = NetCDFFile(swwfile, 'r')
1541
1542        # Get the variables
1543        x = fid.variables['x'][:]
1544        y = fid.variables['y'][:]
1545        z = fid.variables['elevation'][:]
1546        time = fid.variables['time'][:]
1547
1548        try:
1549            geo_reference = Geo_reference(NetCDFObject=fid)
1550        except AttributeError, e:
1551            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1552
1553        #Export to ascii/prj files
1554        sww2dem(domain.filename,
1555                quantity = 'elevation',
1556                cellsize = cellsize,
1557                verbose = False,
1558                format = 'asc')
1559
1560
1561        #Check asc file
1562        ascid = open(ascfile)
1563        lines = ascid.readlines()
1564        ascid.close()
1565
1566        L = lines[0].strip().split()
1567        assert L[0].strip().lower() == 'ncols'
1568        assert L[1].strip().lower() == '5'
1569
1570        L = lines[1].strip().split()
1571        assert L[0].strip().lower() == 'nrows'
1572        assert L[1].strip().lower() == '5'
1573
1574        L = lines[2].strip().split()
1575        assert L[0].strip().lower() == 'xllcorner'
1576        assert allclose(float(L[1].strip().lower()), 308500)
1577
1578        L = lines[3].strip().split()
1579        assert L[0].strip().lower() == 'yllcorner'
1580        assert allclose(float(L[1].strip().lower()), 6189000)
1581
1582        L = lines[4].strip().split()
1583        assert L[0].strip().lower() == 'cellsize'
1584        assert allclose(float(L[1].strip().lower()), cellsize)
1585
1586        L = lines[5].strip().split()
1587        assert L[0].strip() == 'NODATA_value'
1588        assert L[1].strip().lower() == '-9999'
1589
1590
1591        #Check grid values
1592        for j in range(5):
1593            L = lines[6+j].strip().split()
1594            y = (4-j) * cellsize
1595            for i in range(5):
1596                if i+j >= 4:
1597                    assert allclose(float(L[i]), -i*cellsize - y)
1598                else:
1599                    #Missing values
1600                    assert allclose(float(L[i]), -9999)
1601
1602
1603
1604        fid.close()
1605
1606        #Cleanup
1607        os.remove(prjfile)
1608        os.remove(ascfile)
1609        os.remove(swwfile)
1610
1611    def test_sww2ers_simple(self):
1612        """Test that sww information can be converted correctly to asc/prj
1613        format readable by e.g. ArcView
1614        """
1615
1616        import time, os
1617        from Numeric import array, zeros, allclose, Float, concatenate
1618        from Scientific.IO.NetCDF import NetCDFFile
1619
1620        #Setup
1621        self.domain.filename = 'datatest'
1622
1623        headerfile = self.domain.filename + '.ers'
1624        swwfile = self.domain.filename + '.sww'
1625
1626        self.domain.set_datadir('.')
1627        self.domain.format = 'sww'
1628        self.domain.smooth = True
1629        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1630
1631        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1632
1633        sww = get_dataobject(self.domain)
1634        sww.store_connectivity()
1635        sww.store_timestep('stage')
1636
1637        self.domain.evolve_to_end(finaltime = 0.01)
1638        sww.store_timestep('stage')
1639
1640        cellsize = 0.25
1641        #Check contents
1642        #Get NetCDF
1643
1644        fid = NetCDFFile(sww.filename, 'r')
1645
1646        # Get the variables
1647        x = fid.variables['x'][:]
1648        y = fid.variables['y'][:]
1649        z = fid.variables['elevation'][:]
1650        time = fid.variables['time'][:]
1651        stage = fid.variables['stage'][:]
1652
1653
1654        #Export to ers files
1655        #sww2ers(self.domain.filename,
1656        #        quantity = 'elevation',
1657        #        cellsize = cellsize,
1658        #        verbose = False)
1659               
1660        sww2dem(self.domain.filename,
1661                quantity = 'elevation',
1662                cellsize = cellsize,
1663                verbose = False,
1664                format = 'ers')
1665
1666        #Check header data
1667        from ermapper_grids import read_ermapper_header, read_ermapper_data
1668       
1669        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
1670        #print header
1671        assert header['projection'].lower() == '"utm-56"'
1672        assert header['datum'].lower() == '"wgs84"'
1673        assert header['units'].lower() == '"meters"'   
1674        assert header['value'].lower() == '"elevation"'         
1675        assert header['xdimension'] == '0.25'
1676        assert header['ydimension'] == '0.25'   
1677        assert float(header['eastings']) == 308500.0   #xllcorner
1678        assert float(header['northings']) == 6189000.0 #yllcorner       
1679        assert int(header['nroflines']) == 5
1680        assert int(header['nrofcellsperline']) == 5     
1681        assert int(header['nullcellvalue']) == 0   #?           
1682        #FIXME - there is more in the header                   
1683
1684               
1685        #Check grid data               
1686        grid = read_ermapper_data(self.domain.filename + '_elevation') 
1687       
1688        #FIXME (Ole): Why is this the desired reference grid for -x-y?
1689        ref_grid = [0,      0,     0,     0,     0,
1690                    -1,    -1.25, -1.5,  -1.75, -2.0,
1691                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
1692                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1693                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
1694                                         
1695                                         
1696        assert allclose(grid, ref_grid)
1697
1698        fid.close()
1699       
1700        #Cleanup
1701        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1702        #Done (Ole) - it was because sww2ers didn't close it's sww file
1703        os.remove(sww.filename)
1704
1705
1706
1707    def test_ferret2sww(self):
1708        """Test that georeferencing etc works when converting from
1709        ferret format (lat/lon) to sww format (UTM)
1710        """
1711        from Scientific.IO.NetCDF import NetCDFFile
1712
1713        #The test file has
1714        # LON = 150.66667, 150.83334, 151, 151.16667
1715        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1716        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1717        #
1718        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1719        # Fourth value (index==3) is -6.50198 cm
1720
1721
1722        from coordinate_transforms.redfearn import redfearn
1723
1724        fid = NetCDFFile('small_ha.nc')
1725        first_value = fid.variables['HA'][:][0,0,0]
1726        fourth_value = fid.variables['HA'][:][0,0,3]
1727
1728
1729        #Call conversion (with zero origin)
1730        ferret2sww('small', verbose=False,
1731                   origin = (56, 0, 0))
1732
1733
1734        #Work out the UTM coordinates for first point
1735        zone, e, n = redfearn(-34.5, 150.66667)
1736        #print zone, e, n
1737
1738        #Read output file 'small.sww'
1739        fid = NetCDFFile('small.sww')
1740
1741        x = fid.variables['x'][:]
1742        y = fid.variables['y'][:]
1743
1744        #Check that first coordinate is correctly represented
1745        assert allclose(x[0], e)
1746        assert allclose(y[0], n)
1747
1748        #Check first value
1749        stage = fid.variables['stage'][:]
1750        xmomentum = fid.variables['xmomentum'][:]
1751        ymomentum = fid.variables['ymomentum'][:]
1752
1753        #print ymomentum
1754
1755        assert allclose(stage[0,0], first_value/100)  #Meters
1756
1757        #Check fourth value
1758        assert allclose(stage[0,3], fourth_value/100)  #Meters
1759
1760        fid.close()
1761
1762        #Cleanup
1763        import os
1764        os.remove('small.sww')
1765
1766
1767
1768    def test_ferret2sww_2(self):
1769        """Test that georeferencing etc works when converting from
1770        ferret format (lat/lon) to sww format (UTM)
1771        """
1772        from Scientific.IO.NetCDF import NetCDFFile
1773
1774        #The test file has
1775        # LON = 150.66667, 150.83334, 151, 151.16667
1776        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1777        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1778        #
1779        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1780        # Fourth value (index==3) is -6.50198 cm
1781
1782
1783        from coordinate_transforms.redfearn import redfearn
1784
1785        fid = NetCDFFile('small_ha.nc')
1786
1787        #Pick a coordinate and a value
1788
1789        time_index = 1
1790        lat_index = 0
1791        lon_index = 2
1792
1793        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
1794        test_time = fid.variables['TIME'][:][time_index]
1795        test_lat = fid.variables['LAT'][:][lat_index]
1796        test_lon = fid.variables['LON'][:][lon_index]
1797
1798        linear_point_index = lat_index*4 + lon_index
1799        fid.close()
1800
1801        #Call conversion (with zero origin)
1802        ferret2sww('small', verbose=False,
1803                   origin = (56, 0, 0))
1804
1805
1806        #Work out the UTM coordinates for test point
1807        zone, e, n = redfearn(test_lat, test_lon)
1808
1809        #Read output file 'small.sww'
1810        fid = NetCDFFile('small.sww')
1811
1812        x = fid.variables['x'][:]
1813        y = fid.variables['y'][:]
1814
1815        #Check that test coordinate is correctly represented
1816        assert allclose(x[linear_point_index], e)
1817        assert allclose(y[linear_point_index], n)
1818
1819        #Check test value
1820        stage = fid.variables['stage'][:]
1821
1822        assert allclose(stage[time_index, linear_point_index], test_value/100)
1823
1824        fid.close()
1825
1826        #Cleanup
1827        import os
1828        os.remove('small.sww')
1829
1830
1831
1832    def test_ferret2sww3(self):
1833        """Elevation included
1834        """
1835        from Scientific.IO.NetCDF import NetCDFFile
1836
1837        #The test file has
1838        # LON = 150.66667, 150.83334, 151, 151.16667
1839        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1840        # ELEVATION = [-1 -2 -3 -4
1841        #             -5 -6 -7 -8
1842        #              ...
1843        #              ...      -16]
1844        # where the top left corner is -1m,
1845        # and the ll corner is -13.0m
1846        #
1847        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1848        # Fourth value (index==3) is -6.50198 cm
1849
1850        from coordinate_transforms.redfearn import redfearn
1851        import os
1852        fid1 = NetCDFFile('test_ha.nc','w')
1853        fid2 = NetCDFFile('test_ua.nc','w')
1854        fid3 = NetCDFFile('test_va.nc','w')
1855        fid4 = NetCDFFile('test_e.nc','w')
1856
1857        h1_list = [150.66667,150.83334,151.]
1858        h2_list = [-34.5,-34.33333]
1859
1860        long_name = 'LON'
1861        lat_name = 'LAT'
1862
1863        nx = 3
1864        ny = 2
1865
1866        for fid in [fid1,fid2,fid3]:
1867            fid.createDimension(long_name,nx)
1868            fid.createVariable(long_name,'d',(long_name,))
1869            fid.variables[long_name].point_spacing='uneven'
1870            fid.variables[long_name].units='degrees_east'
1871            fid.variables[long_name].assignValue(h1_list)
1872
1873            fid.createDimension(lat_name,ny)
1874            fid.createVariable(lat_name,'d',(lat_name,))
1875            fid.variables[lat_name].point_spacing='uneven'
1876            fid.variables[lat_name].units='degrees_north'
1877            fid.variables[lat_name].assignValue(h2_list)
1878
1879            fid.createDimension('TIME',2)
1880            fid.createVariable('TIME','d',('TIME',))
1881            fid.variables['TIME'].point_spacing='uneven'
1882            fid.variables['TIME'].units='seconds'
1883            fid.variables['TIME'].assignValue([0.,1.])
1884            if fid == fid3: break
1885
1886
1887        for fid in [fid4]:
1888            fid.createDimension(long_name,nx)
1889            fid.createVariable(long_name,'d',(long_name,))
1890            fid.variables[long_name].point_spacing='uneven'
1891            fid.variables[long_name].units='degrees_east'
1892            fid.variables[long_name].assignValue(h1_list)
1893
1894            fid.createDimension(lat_name,ny)
1895            fid.createVariable(lat_name,'d',(lat_name,))
1896            fid.variables[lat_name].point_spacing='uneven'
1897            fid.variables[lat_name].units='degrees_north'
1898            fid.variables[lat_name].assignValue(h2_list)
1899
1900        name = {}
1901        name[fid1]='HA'
1902        name[fid2]='UA'
1903        name[fid3]='VA'
1904        name[fid4]='ELEVATION'
1905
1906        units = {}
1907        units[fid1]='cm'
1908        units[fid2]='cm/s'
1909        units[fid3]='cm/s'
1910        units[fid4]='m'
1911
1912        values = {}
1913        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
1914        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
1915        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
1916        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
1917
1918        for fid in [fid1,fid2,fid3]:
1919          fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
1920          fid.variables[name[fid]].point_spacing='uneven'
1921          fid.variables[name[fid]].units=units[fid]
1922          fid.variables[name[fid]].assignValue(values[fid])
1923          fid.variables[name[fid]].missing_value = -99999999.
1924          if fid == fid3: break
1925
1926        for fid in [fid4]:
1927            fid.createVariable(name[fid],'d',(lat_name,long_name))
1928            fid.variables[name[fid]].point_spacing='uneven'
1929            fid.variables[name[fid]].units=units[fid]
1930            fid.variables[name[fid]].assignValue(values[fid])
1931            fid.variables[name[fid]].missing_value = -99999999.
1932
1933
1934        fid1.sync(); fid1.close()
1935        fid2.sync(); fid2.close()
1936        fid3.sync(); fid3.close()
1937        fid4.sync(); fid4.close()
1938
1939        fid1 = NetCDFFile('test_ha.nc','r')
1940        fid2 = NetCDFFile('test_e.nc','r')
1941        fid3 = NetCDFFile('test_va.nc','r')
1942
1943
1944        first_amp = fid1.variables['HA'][:][0,0,0]
1945        third_amp = fid1.variables['HA'][:][0,0,2]
1946        first_elevation = fid2.variables['ELEVATION'][0,0]
1947        third_elevation= fid2.variables['ELEVATION'][:][0,2]
1948        first_speed = fid3.variables['VA'][0,0,0]
1949        third_speed = fid3.variables['VA'][:][0,0,2]
1950
1951        fid1.close()
1952        fid2.close()
1953        fid3.close()
1954
1955        #Call conversion (with zero origin)
1956        ferret2sww('test', verbose=False,
1957                   origin = (56, 0, 0))
1958
1959        os.remove('test_va.nc')
1960        os.remove('test_ua.nc')
1961        os.remove('test_ha.nc')
1962        os.remove('test_e.nc')
1963
1964        #Read output file 'test.sww'
1965        fid = NetCDFFile('test.sww')
1966
1967
1968        #Check first value
1969        elevation = fid.variables['elevation'][:]
1970        stage = fid.variables['stage'][:]
1971        xmomentum = fid.variables['xmomentum'][:]
1972        ymomentum = fid.variables['ymomentum'][:]
1973
1974        #print ymomentum
1975        first_height = first_amp/100 - first_elevation
1976        third_height = third_amp/100 - third_elevation
1977        first_momentum=first_speed*first_height/100
1978        third_momentum=third_speed*third_height/100
1979
1980        assert allclose(ymomentum[0][0],first_momentum)  #Meters
1981        assert allclose(ymomentum[0][2],third_momentum)  #Meters
1982
1983        fid.close()
1984
1985        #Cleanup
1986        os.remove('test.sww')
1987
1988
1989
1990
1991    def test_ferret2sww_nz_origin(self):
1992        from Scientific.IO.NetCDF import NetCDFFile
1993        from coordinate_transforms.redfearn import redfearn
1994
1995        #Call conversion (with nonzero origin)
1996        ferret2sww('small', verbose=False,
1997                   origin = (56, 100000, 200000))
1998
1999
2000        #Work out the UTM coordinates for first point
2001        zone, e, n = redfearn(-34.5, 150.66667)
2002
2003        #Read output file 'small.sww'
2004        fid = NetCDFFile('small.sww', 'r')
2005
2006        x = fid.variables['x'][:]
2007        y = fid.variables['y'][:]
2008
2009        #Check that first coordinate is correctly represented
2010        assert allclose(x[0], e-100000)
2011        assert allclose(y[0], n-200000)
2012
2013        fid.close()
2014
2015        #Cleanup
2016        import os
2017        os.remove('small.sww')
2018
2019
2020
2021    def test_sww_extent(self):
2022        """Not a test, rather a look at the sww format
2023        """
2024
2025        import time, os
2026        from Numeric import array, zeros, allclose, Float, concatenate
2027        from Scientific.IO.NetCDF import NetCDFFile
2028
2029        self.domain.filename = 'datatest' + str(id(self))
2030        self.domain.format = 'sww'
2031        self.domain.smooth = True
2032        self.domain.reduction = mean
2033        self.domain.set_datadir('.')
2034
2035
2036        sww = get_dataobject(self.domain)
2037        sww.store_connectivity()
2038        sww.store_timestep('stage')
2039        self.domain.time = 2.
2040
2041        #Modify stage at second timestep
2042        stage = self.domain.quantities['stage'].vertex_values
2043        self.domain.set_quantity('stage', stage/2)
2044
2045        sww.store_timestep('stage')
2046
2047        file_and_extension_name = self.domain.filename + ".sww"
2048        #print "file_and_extension_name",file_and_extension_name
2049        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2050               extent_sww(file_and_extension_name )
2051
2052        assert allclose(xmin, 0.0)
2053        assert allclose(xmax, 1.0)
2054        assert allclose(ymin, 0.0)
2055        assert allclose(ymax, 1.0)
2056        assert allclose(stagemin, -0.85)
2057        assert allclose(stagemax, 0.15)
2058
2059
2060        #Cleanup
2061        os.remove(sww.filename)
2062
2063
2064
2065    def test_sww2domain(self):
2066        ################################################
2067        #Create a test domain, and evolve and save it.
2068        ################################################
2069        from mesh_factory import rectangular
2070        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2071             Constant_height, Time_boundary, Transmissive_boundary
2072        from Numeric import array
2073
2074        #Create basic mesh
2075
2076        yiel=0.01
2077        points, vertices, boundary = rectangular(10,10)
2078
2079        #Create shallow water domain
2080        domain = Domain(points, vertices, boundary)
2081        domain.geo_reference = Geo_reference(56,11,11)
2082        domain.smooth = False
2083        domain.visualise = False
2084        domain.store = True
2085        domain.filename = 'bedslope'
2086        domain.default_order=2
2087        #Bed-slope and friction
2088        domain.set_quantity('elevation', lambda x,y: -x/3)
2089        domain.set_quantity('friction', 0.1)
2090        # Boundary conditions
2091        from math import sin, pi
2092        Br = Reflective_boundary(domain)
2093        Bt = Transmissive_boundary(domain)
2094        Bd = Dirichlet_boundary([0.2,0.,0.])
2095        Bw = Time_boundary(domain=domain,
2096                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2097
2098        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2099        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2100
2101        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2102        #Initial condition
2103        h = 0.05
2104        elevation = domain.quantities['elevation'].vertex_values
2105        domain.set_quantity('stage', elevation + h)
2106        #elevation = domain.get_quantity('elevation')
2107        #domain.set_quantity('stage', elevation + h)
2108
2109        domain.check_integrity()
2110        #Evolution
2111        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2112        #    domain.write_time()
2113            pass
2114
2115
2116        ##########################################
2117        #Import the example's file as a new domain
2118        ##########################################
2119        from data_manager import sww2domain
2120        from Numeric import allclose
2121        import os
2122
2123        filename = domain.datadir+os.sep+domain.filename+'.sww'
2124        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2125        #points, vertices, boundary = rectangular(15,15)
2126        #domain2.boundary = boundary
2127        ###################
2128        ##NOW TEST IT!!!
2129        ###################
2130
2131        bits = ['vertex_coordinates']
2132        for quantity in ['elevation']+domain.quantities_to_be_stored:
2133            bits.append('quantities["%s"].get_integral()'%quantity)
2134            bits.append('get_quantity("%s")'%quantity)
2135
2136        for bit in bits:
2137            #print 'testing that domain.'+bit+' has been restored'
2138            #print bit
2139        #print 'done'
2140            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2141
2142        ######################################
2143        #Now evolve them both, just to be sure
2144        ######################################x
2145        visualise = False
2146        #visualise = True
2147        domain.visualise = visualise
2148        domain.time = 0.
2149        from time import sleep
2150
2151        final = .1
2152        domain.set_quantity('friction', 0.1)
2153        domain.store = False
2154        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2155
2156
2157        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2158            if visualise: sleep(1.)
2159            #domain.write_time()
2160            pass
2161
2162        final = final - (domain2.starttime-domain.starttime)
2163        #BUT since domain1 gets time hacked back to 0:
2164        final = final + (domain2.starttime-domain.starttime)
2165
2166        domain2.smooth = False
2167        domain2.visualise = visualise
2168        domain2.store = False
2169        domain2.default_order=2
2170        domain2.set_quantity('friction', 0.1)
2171        #Bed-slope and friction
2172        # Boundary conditions
2173        Bd2=Dirichlet_boundary([0.2,0.,0.])
2174        domain2.boundary = domain.boundary
2175        #print 'domain2.boundary'
2176        #print domain2.boundary
2177        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2178        #domain2.set_boundary({'exterior': Bd})
2179
2180        domain2.check_integrity()
2181
2182        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2183            if visualise: sleep(1.)
2184            #domain2.write_time()
2185            pass
2186
2187        ###################
2188        ##NOW TEST IT!!!
2189        ##################
2190
2191        bits = [ 'vertex_coordinates']
2192
2193        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2194            bits.append('quantities["%s"].get_integral()'%quantity)
2195            bits.append('get_quantity("%s")'%quantity)
2196
2197        for bit in bits:
2198            #print bit
2199            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2200
2201
2202    def test_sww2domain2(self):
2203        ##################################################################
2204        #Same as previous test, but this checks how NaNs are handled.
2205        ##################################################################
2206
2207
2208        from mesh_factory import rectangular
2209        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2210             Constant_height, Time_boundary, Transmissive_boundary
2211        from Numeric import array
2212
2213        #Create basic mesh
2214        points, vertices, boundary = rectangular(2,2)
2215
2216        #Create shallow water domain
2217        domain = Domain(points, vertices, boundary)
2218        domain.smooth = False
2219        domain.visualise = False
2220        domain.store = True
2221        domain.filename = 'bedslope'
2222        domain.default_order=2
2223        domain.quantities_to_be_stored=['stage']
2224
2225        domain.set_quantity('elevation', lambda x,y: -x/3)
2226        domain.set_quantity('friction', 0.1)
2227
2228        from math import sin, pi
2229        Br = Reflective_boundary(domain)
2230        Bt = Transmissive_boundary(domain)
2231        Bd = Dirichlet_boundary([0.2,0.,0.])
2232        Bw = Time_boundary(domain=domain,
2233                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2234
2235        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2236
2237        h = 0.05
2238        elevation = domain.quantities['elevation'].vertex_values
2239        domain.set_quantity('stage', elevation + h)
2240
2241        domain.check_integrity()
2242
2243        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2244            pass
2245            #domain.write_time()
2246
2247
2248
2249        ##################################
2250        #Import the 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
2258        #Fail because NaNs are present
2259        try:
2260            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2261            assert True == False
2262        except:
2263            #Now import it, filling NaNs to be 0
2264            filler = 0
2265            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2266        bits = [ 'geo_reference.get_xllcorner()',
2267                'geo_reference.get_yllcorner()',
2268                'vertex_coordinates']
2269
2270        for quantity in ['elevation']+domain.quantities_to_be_stored:
2271            bits.append('quantities["%s"].get_integral()'%quantity)
2272            bits.append('get_quantity("%s")'%quantity)
2273
2274        for bit in bits:
2275        #    print 'testing that domain.'+bit+' has been restored'
2276            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2277
2278        assert max(max(domain2.get_quantity('xmomentum')))==filler
2279        assert min(min(domain2.get_quantity('xmomentum')))==filler
2280        assert max(max(domain2.get_quantity('ymomentum')))==filler
2281        assert min(min(domain2.get_quantity('ymomentum')))==filler
2282
2283        #print 'passed'
2284
2285        #cleanup
2286        #import os
2287        #os.remove(domain.datadir+'/'+domain.filename+'.sww')
2288
2289
2290    #def test_weed(self):
2291        from data_manager import weed
2292
2293        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2294        volumes1 = [[0,1,2],[3,4,5]]
2295        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2296        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2297
2298        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2299
2300        assert len(points2)==len(coordinates2)
2301        for i in range(len(coordinates2)):
2302            coordinate = tuple(coordinates2[i])
2303            assert points2.has_key(coordinate)
2304            points2[coordinate]=i
2305
2306        for triangle in volumes1:
2307            for coordinate in triangle:
2308                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2309                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2310
2311
2312    #FIXME This fails - smooth makes the comparism too hard for allclose
2313    def ztest_sww2domain3(self):
2314        ################################################
2315        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2316        ################################################
2317        from mesh_factory import rectangular
2318        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2319             Constant_height, Time_boundary, Transmissive_boundary
2320        from Numeric import array
2321        #Create basic mesh
2322
2323        yiel=0.01
2324        points, vertices, boundary = rectangular(10,10)
2325
2326        #Create shallow water domain
2327        domain = Domain(points, vertices, boundary)
2328        domain.geo_reference = Geo_reference(56,11,11)
2329        domain.smooth = True
2330        domain.visualise = False
2331        domain.store = True
2332        domain.filename = 'bedslope'
2333        domain.default_order=2
2334        #Bed-slope and friction
2335        domain.set_quantity('elevation', lambda x,y: -x/3)
2336        domain.set_quantity('friction', 0.1)
2337        # Boundary conditions
2338        from math import sin, pi
2339        Br = Reflective_boundary(domain)
2340        Bt = Transmissive_boundary(domain)
2341        Bd = Dirichlet_boundary([0.2,0.,0.])
2342        Bw = Time_boundary(domain=domain,
2343                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2344
2345        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2346
2347        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2348        #Initial condition
2349        h = 0.05
2350        elevation = domain.quantities['elevation'].vertex_values
2351        domain.set_quantity('stage', elevation + h)
2352        #elevation = domain.get_quantity('elevation')
2353        #domain.set_quantity('stage', elevation + h)
2354
2355        domain.check_integrity()
2356        #Evolution
2357        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2358        #    domain.write_time()
2359            pass
2360
2361
2362        ##########################################
2363        #Import the example's file as a new domain
2364        ##########################################
2365        from data_manager import sww2domain
2366        from Numeric import allclose
2367        import os
2368
2369        filename = domain.datadir+os.sep+domain.filename+'.sww'
2370        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2371        #points, vertices, boundary = rectangular(15,15)
2372        #domain2.boundary = boundary
2373        ###################
2374        ##NOW TEST IT!!!
2375        ###################
2376
2377        #FIXME smooth domain so that they can be compared
2378
2379
2380        bits = []#'vertex_coordinates']
2381        for quantity in ['elevation']+domain.quantities_to_be_stored:
2382            bits.append('quantities["%s"].get_integral()'%quantity)
2383            #bits.append('get_quantity("%s")'%quantity)
2384
2385        for bit in bits:
2386            #print 'testing that domain.'+bit+' has been restored'
2387            #print bit
2388            #print 'done'
2389            #print ('domain.'+bit), eval('domain.'+bit)
2390            #print ('domain2.'+bit), eval('domain2.'+bit)
2391            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2392            pass
2393
2394        ######################################
2395        #Now evolve them both, just to be sure
2396        ######################################x
2397        visualise = False
2398        visualise = True
2399        domain.visualise = visualise
2400        domain.time = 0.
2401        from time import sleep
2402
2403        final = .5
2404        domain.set_quantity('friction', 0.1)
2405        domain.store = False
2406        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2407
2408        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2409            if visualise: sleep(.03)
2410            #domain.write_time()
2411            pass
2412
2413        domain2.smooth = True
2414        domain2.visualise = visualise
2415        domain2.store = False
2416        domain2.default_order=2
2417        domain2.set_quantity('friction', 0.1)
2418        #Bed-slope and friction
2419        # Boundary conditions
2420        Bd2=Dirichlet_boundary([0.2,0.,0.])
2421        Br2 = Reflective_boundary(domain2)
2422        domain2.boundary = domain.boundary
2423        #print 'domain2.boundary'
2424        #print domain2.boundary
2425        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2426        #domain2.boundary = domain.boundary
2427        #domain2.set_boundary({'exterior': Bd})
2428
2429        domain2.check_integrity()
2430
2431        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2432            if visualise: sleep(.03)
2433            #domain2.write_time()
2434            pass
2435
2436        ###################
2437        ##NOW TEST IT!!!
2438        ##################
2439
2440        bits = [ 'vertex_coordinates']
2441
2442        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2443            #bits.append('quantities["%s"].get_integral()'%quantity)
2444            bits.append('get_quantity("%s")'%quantity)
2445
2446        for bit in bits:
2447            print bit
2448            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2449
2450
2451    def test_decimate_dem(self):
2452        """Test decimation of dem file
2453        """
2454
2455        import os
2456        from Numeric import ones, allclose, Float, arange
2457        from Scientific.IO.NetCDF import NetCDFFile
2458
2459        #Write test dem file
2460        root = 'decdemtest'
2461
2462        filename = root + '.dem'
2463        fid = NetCDFFile(filename, 'w')
2464
2465        fid.institution = 'Geoscience Australia'
2466        fid.description = 'NetCDF DEM format for compact and portable ' +\
2467                          'storage of spatial point data'
2468
2469        nrows = 15
2470        ncols = 18
2471
2472        fid.ncols = ncols
2473        fid.nrows = nrows
2474        fid.xllcorner = 2000.5
2475        fid.yllcorner = 3000.5
2476        fid.cellsize = 25
2477        fid.NODATA_value = -9999
2478
2479        fid.zone = 56
2480        fid.false_easting = 0.0
2481        fid.false_northing = 0.0
2482        fid.projection = 'UTM'
2483        fid.datum = 'WGS84'
2484        fid.units = 'METERS'
2485
2486        fid.createDimension('number_of_points', nrows*ncols)
2487
2488        fid.createVariable('elevation', Float, ('number_of_points',))
2489
2490        elevation = fid.variables['elevation']
2491
2492        elevation[:] = (arange(nrows*ncols))
2493
2494        fid.close()
2495
2496        #generate the elevation values expected in the decimated file
2497        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2498                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2499                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2500                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2501                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2502                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2503                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2504                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2505                         (144+145+146+162+163+164+180+181+182) / 9.0,
2506                         (148+149+150+166+167+168+184+185+186) / 9.0,
2507                         (152+153+154+170+171+172+188+189+190) / 9.0,
2508                         (156+157+158+174+175+176+192+193+194) / 9.0,
2509                         (216+217+218+234+235+236+252+253+254) / 9.0,
2510                         (220+221+222+238+239+240+256+257+258) / 9.0,
2511                         (224+225+226+242+243+244+260+261+262) / 9.0,
2512                         (228+229+230+246+247+248+264+265+266) / 9.0]
2513
2514        #generate a stencil for computing the decimated values
2515        stencil = ones((3,3), Float) / 9.0
2516
2517        decimate_dem(root, stencil=stencil, cellsize_new=100)
2518
2519        #Open decimated NetCDF file
2520        fid = NetCDFFile(root + '_100.dem', 'r')
2521
2522        # Get decimated elevation
2523        elevation = fid.variables['elevation']
2524
2525        #Check values
2526        assert allclose(elevation, ref_elevation)
2527
2528        #Cleanup
2529        fid.close()
2530
2531        os.remove(root + '.dem')
2532        os.remove(root + '_100.dem')
2533
2534    def test_decimate_dem_NODATA(self):
2535        """Test decimation of dem file that includes NODATA values
2536        """
2537
2538        import os
2539        from Numeric import ones, allclose, Float, arange, reshape
2540        from Scientific.IO.NetCDF import NetCDFFile
2541
2542        #Write test dem file
2543        root = 'decdemtest'
2544
2545        filename = root + '.dem'
2546        fid = NetCDFFile(filename, 'w')
2547
2548        fid.institution = 'Geoscience Australia'
2549        fid.description = 'NetCDF DEM format for compact and portable ' +\
2550                          'storage of spatial point data'
2551
2552        nrows = 15
2553        ncols = 18
2554        NODATA_value = -9999
2555
2556        fid.ncols = ncols
2557        fid.nrows = nrows
2558        fid.xllcorner = 2000.5
2559        fid.yllcorner = 3000.5
2560        fid.cellsize = 25
2561        fid.NODATA_value = NODATA_value
2562
2563        fid.zone = 56
2564        fid.false_easting = 0.0
2565        fid.false_northing = 0.0
2566        fid.projection = 'UTM'
2567        fid.datum = 'WGS84'
2568        fid.units = 'METERS'
2569
2570        fid.createDimension('number_of_points', nrows*ncols)
2571
2572        fid.createVariable('elevation', Float, ('number_of_points',))
2573
2574        elevation = fid.variables['elevation']
2575
2576        #generate initial elevation values
2577        elevation_tmp = (arange(nrows*ncols))
2578        #add some NODATA values
2579        elevation_tmp[0]   = NODATA_value
2580        elevation_tmp[95]  = NODATA_value
2581        elevation_tmp[188] = NODATA_value
2582        elevation_tmp[189] = NODATA_value
2583        elevation_tmp[190] = NODATA_value
2584        elevation_tmp[209] = NODATA_value
2585        elevation_tmp[252] = NODATA_value
2586
2587        elevation[:] = elevation_tmp
2588
2589        fid.close()
2590
2591        #generate the elevation values expected in the decimated file
2592        ref_elevation = [NODATA_value,
2593                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2594                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2595                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2596                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2597                         NODATA_value,
2598                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2599                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2600                         (144+145+146+162+163+164+180+181+182) / 9.0,
2601                         (148+149+150+166+167+168+184+185+186) / 9.0,
2602                         NODATA_value,
2603                         (156+157+158+174+175+176+192+193+194) / 9.0,
2604                         NODATA_value,
2605                         (220+221+222+238+239+240+256+257+258) / 9.0,
2606                         (224+225+226+242+243+244+260+261+262) / 9.0,
2607                         (228+229+230+246+247+248+264+265+266) / 9.0]
2608
2609        #generate a stencil for computing the decimated values
2610        stencil = ones((3,3), Float) / 9.0
2611
2612        decimate_dem(root, stencil=stencil, cellsize_new=100)
2613
2614        #Open decimated NetCDF file
2615        fid = NetCDFFile(root + '_100.dem', 'r')
2616
2617        # Get decimated elevation
2618        elevation = fid.variables['elevation']
2619
2620        #Check values
2621        assert allclose(elevation, ref_elevation)
2622
2623        #Cleanup
2624        fid.close()
2625
2626        os.remove(root + '.dem')
2627        os.remove(root + '_100.dem')
2628
2629    def xxxtestz_sww2ers_real(self):
2630        """Test that sww information can be converted correctly to asc/prj
2631        format readable by e.g. ArcView
2632        """
2633
2634        import time, os
2635        from Numeric import array, zeros, allclose, Float, concatenate
2636        from Scientific.IO.NetCDF import NetCDFFile
2637
2638        # the memory optimised least squares
2639        #  cellsize = 20,   # this one seems to hang
2640        #  cellsize = 200000, # Ran 1 test in 269.703s
2641                                #Ran 1 test in 267.344s
2642        #  cellsize = 20000,  # Ran 1 test in 460.922s
2643        #  cellsize = 2000   #Ran 1 test in 5340.250s
2644        #  cellsize = 200   #this one seems to hang, building matirx A
2645
2646        # not optimised
2647        # seems to hang
2648        #  cellsize = 2000   # Ran 1 test in 5334.563s
2649        #Export to ascii/prj files
2650        sww2dem('karratha_100m',
2651                quantity = 'depth',
2652                cellsize = 200000,
2653                verbose = True)
2654
2655
2656
2657
2658#-------------------------------------------------------------
2659if __name__ == "__main__":
2660    suite = unittest.makeSuite(Test_Data_Manager,'test')
2661    #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest')
2662    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')   
2663    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
2664    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
2665    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
2666    runner = unittest.TextTestRunner()
2667    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.