source: inundation/pyvolution/test_data_manager.py @ 2007

Last change on this file since 2007 was 2007, checked in by duncan, 18 years ago

working on gippsland to sww
File size: 84.7 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_sww2dem_asc_elevation(self):
745        """Test that sww information can be converted correctly to asc/prj
746        format readable by e.g. ArcView
747        """
748
749        import time, os
750        from Numeric import array, zeros, allclose, Float, concatenate
751        from Scientific.IO.NetCDF import NetCDFFile
752
753        #Setup
754        self.domain.filename = 'datatest'
755
756        prjfile = self.domain.filename + '_elevation.prj'
757        ascfile = self.domain.filename + '_elevation.asc'
758        swwfile = self.domain.filename + '.sww'
759
760        self.domain.set_datadir('.')
761        self.domain.format = 'sww'
762        self.domain.smooth = True
763        self.domain.set_quantity('elevation', lambda x,y: -x-y)
764
765        self.domain.geo_reference = Geo_reference(56,308500,6189000)
766
767        sww = get_dataobject(self.domain)
768        sww.store_connectivity()
769        sww.store_timestep('stage')
770
771        self.domain.evolve_to_end(finaltime = 0.01)
772        sww.store_timestep('stage')
773
774        cellsize = 0.25
775        #Check contents
776        #Get NetCDF
777
778        fid = NetCDFFile(sww.filename, 'r')
779
780        # Get the variables
781        x = fid.variables['x'][:]
782        y = fid.variables['y'][:]
783        z = fid.variables['elevation'][:]
784        time = fid.variables['time'][:]
785        stage = fid.variables['stage'][:]
786
787
788        #Export to ascii/prj files
789        sww2dem(self.domain.filename,
790                quantity = 'elevation',
791                cellsize = cellsize,
792                verbose = False,
793                format = 'asc')
794               
795        #Check prj (meta data)
796        prjid = open(prjfile)
797        lines = prjid.readlines()
798        prjid.close()
799
800        L = lines[0].strip().split()
801        assert L[0].strip().lower() == 'projection'
802        assert L[1].strip().lower() == 'utm'
803
804        L = lines[1].strip().split()
805        assert L[0].strip().lower() == 'zone'
806        assert L[1].strip().lower() == '56'
807
808        L = lines[2].strip().split()
809        assert L[0].strip().lower() == 'datum'
810        assert L[1].strip().lower() == 'wgs84'
811
812        L = lines[3].strip().split()
813        assert L[0].strip().lower() == 'zunits'
814        assert L[1].strip().lower() == 'no'
815
816        L = lines[4].strip().split()
817        assert L[0].strip().lower() == 'units'
818        assert L[1].strip().lower() == 'meters'
819
820        L = lines[5].strip().split()
821        assert L[0].strip().lower() == 'spheroid'
822        assert L[1].strip().lower() == 'wgs84'
823
824        L = lines[6].strip().split()
825        assert L[0].strip().lower() == 'xshift'
826        assert L[1].strip().lower() == '500000'
827
828        L = lines[7].strip().split()
829        assert L[0].strip().lower() == 'yshift'
830        assert L[1].strip().lower() == '10000000'
831
832        L = lines[8].strip().split()
833        assert L[0].strip().lower() == 'parameters'
834
835
836        #Check asc file
837        ascid = open(ascfile)
838        lines = ascid.readlines()
839        ascid.close()
840
841        L = lines[0].strip().split()
842        assert L[0].strip().lower() == 'ncols'
843        assert L[1].strip().lower() == '5'
844
845        L = lines[1].strip().split()
846        assert L[0].strip().lower() == 'nrows'
847        assert L[1].strip().lower() == '5'
848
849        L = lines[2].strip().split()
850        assert L[0].strip().lower() == 'xllcorner'
851        assert allclose(float(L[1].strip().lower()), 308500)
852
853        L = lines[3].strip().split()
854        assert L[0].strip().lower() == 'yllcorner'
855        assert allclose(float(L[1].strip().lower()), 6189000)
856
857        L = lines[4].strip().split()
858        assert L[0].strip().lower() == 'cellsize'
859        assert allclose(float(L[1].strip().lower()), cellsize)
860
861        L = lines[5].strip().split()
862        assert L[0].strip() == 'NODATA_value'
863        assert L[1].strip().lower() == '-9999'
864
865        #Check grid values
866        for j in range(5):
867            L = lines[6+j].strip().split()
868            y = (4-j) * cellsize
869            for i in range(5):
870                assert allclose(float(L[i]), -i*cellsize - y)
871
872
873        fid.close()
874
875        #Cleanup
876        os.remove(prjfile)
877        os.remove(ascfile)
878        os.remove(swwfile)
879
880
881
882    def test_sww2dem_larger(self):
883        """Test that sww information can be converted correctly to asc/prj
884        format readable by e.g. ArcView. Here:
885
886        ncols         11
887        nrows         11
888        xllcorner     308500
889        yllcorner     6189000
890        cellsize      10.000000
891        NODATA_value  -9999
892        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
893         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
894         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
895         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
896         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
897         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
898         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
899         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
900         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
901         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
902           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
903       
904        """
905
906        import time, os
907        from Numeric import array, zeros, allclose, Float, concatenate
908        from Scientific.IO.NetCDF import NetCDFFile
909
910        #Setup
911
912        from mesh_factory import rectangular
913
914        #Create basic mesh (100m x 100m)
915        points, vertices, boundary = rectangular(2, 2, 100, 100)
916
917        #Create shallow water domain
918        domain = Domain(points, vertices, boundary)
919        domain.default_order = 2
920
921        domain.filename = 'datatest'
922
923        prjfile = domain.filename + '_elevation.prj'
924        ascfile = domain.filename + '_elevation.asc'
925        swwfile = domain.filename + '.sww'
926
927        domain.set_datadir('.')
928        domain.format = 'sww'
929        domain.smooth = True
930        domain.geo_reference = Geo_reference(56, 308500, 6189000)
931
932        #
933        domain.set_quantity('elevation', lambda x,y: -x-y)
934        domain.set_quantity('stage', 0)       
935
936        B = Transmissive_boundary(domain)
937        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
938
939
940        #
941        sww = get_dataobject(domain)
942        sww.store_connectivity()
943        sww.store_timestep('stage')
944
945        domain.evolve_to_end(finaltime = 0.01)
946        sww.store_timestep('stage')
947
948        cellsize = 10  #10m grid
949
950       
951        #Check contents
952        #Get NetCDF
953
954        fid = NetCDFFile(sww.filename, 'r')
955
956        # Get the variables
957        x = fid.variables['x'][:]
958        y = fid.variables['y'][:]
959        z = fid.variables['elevation'][:]
960        time = fid.variables['time'][:]
961        stage = fid.variables['stage'][:]
962
963
964        #Export to ascii/prj files
965        sww2dem(domain.filename,
966                quantity = 'elevation',
967                cellsize = cellsize,
968                verbose = False,
969                format = 'asc')
970       
971               
972        #Check prj (meta data)
973        prjid = open(prjfile)
974        lines = prjid.readlines()
975        prjid.close()
976
977        L = lines[0].strip().split()
978        assert L[0].strip().lower() == 'projection'
979        assert L[1].strip().lower() == 'utm'
980
981        L = lines[1].strip().split()
982        assert L[0].strip().lower() == 'zone'
983        assert L[1].strip().lower() == '56'
984
985        L = lines[2].strip().split()
986        assert L[0].strip().lower() == 'datum'
987        assert L[1].strip().lower() == 'wgs84'
988
989        L = lines[3].strip().split()
990        assert L[0].strip().lower() == 'zunits'
991        assert L[1].strip().lower() == 'no'
992
993        L = lines[4].strip().split()
994        assert L[0].strip().lower() == 'units'
995        assert L[1].strip().lower() == 'meters'
996
997        L = lines[5].strip().split()
998        assert L[0].strip().lower() == 'spheroid'
999        assert L[1].strip().lower() == 'wgs84'
1000
1001        L = lines[6].strip().split()
1002        assert L[0].strip().lower() == 'xshift'
1003        assert L[1].strip().lower() == '500000'
1004
1005        L = lines[7].strip().split()
1006        assert L[0].strip().lower() == 'yshift'
1007        assert L[1].strip().lower() == '10000000'
1008
1009        L = lines[8].strip().split()
1010        assert L[0].strip().lower() == 'parameters'
1011
1012
1013        #Check asc file
1014        ascid = open(ascfile)
1015        lines = ascid.readlines()
1016        ascid.close()
1017
1018        L = lines[0].strip().split()
1019        assert L[0].strip().lower() == 'ncols'
1020        assert L[1].strip().lower() == '11'
1021
1022        L = lines[1].strip().split()
1023        assert L[0].strip().lower() == 'nrows'
1024        assert L[1].strip().lower() == '11'
1025
1026        L = lines[2].strip().split()
1027        assert L[0].strip().lower() == 'xllcorner'
1028        assert allclose(float(L[1].strip().lower()), 308500)
1029
1030        L = lines[3].strip().split()
1031        assert L[0].strip().lower() == 'yllcorner'
1032        assert allclose(float(L[1].strip().lower()), 6189000)
1033
1034        L = lines[4].strip().split()
1035        assert L[0].strip().lower() == 'cellsize'
1036        assert allclose(float(L[1].strip().lower()), cellsize)
1037
1038        L = lines[5].strip().split()
1039        assert L[0].strip() == 'NODATA_value'
1040        assert L[1].strip().lower() == '-9999'
1041
1042        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1043        for i, line in enumerate(lines[6:]):
1044            for j, value in enumerate( line.split() ):
1045                #assert allclose(float(value), -(10-i+j)*cellsize)
1046                assert float(value) == -(10-i+j)*cellsize
1047       
1048
1049        fid.close()
1050
1051        #Cleanup
1052        os.remove(prjfile)
1053        os.remove(ascfile)
1054        os.remove(swwfile)
1055
1056
1057
1058    def test_sww2dem_boundingbox(self):
1059        """Test that sww information can be converted correctly to asc/prj
1060        format readable by e.g. ArcView.
1061        This will test that mesh can be restricted by bounding box
1062
1063        Original extent is 100m x 100m:
1064
1065        Eastings:   308500 -  308600
1066        Northings: 6189000 - 6189100
1067
1068        Bounding box changes this to the 50m x 50m square defined by
1069
1070        Eastings:   308530 -  308570
1071        Northings: 6189050 - 6189100
1072
1073        The cropped values should be
1074
1075         -130 -140 -150 -160 -170 
1076         -120 -130 -140 -150 -160 
1077         -110 -120 -130 -140 -150 
1078         -100 -110 -120 -130 -140 
1079          -90 -100 -110 -120 -130 
1080          -80  -90 -100 -110 -120
1081
1082        and the new lower reference point should be   
1083        Eastings:   308530
1084        Northings: 6189050
1085       
1086        Original dataset is the same as in test_sww2dem_larger()
1087       
1088        """
1089
1090        import time, os
1091        from Numeric import array, zeros, allclose, Float, concatenate
1092        from Scientific.IO.NetCDF import NetCDFFile
1093
1094        #Setup
1095
1096        from mesh_factory import rectangular
1097
1098        #Create basic mesh (100m x 100m)
1099        points, vertices, boundary = rectangular(2, 2, 100, 100)
1100
1101        #Create shallow water domain
1102        domain = Domain(points, vertices, boundary)
1103        domain.default_order = 2
1104
1105        domain.filename = 'datatest'
1106
1107        prjfile = domain.filename + '_elevation.prj'
1108        ascfile = domain.filename + '_elevation.asc'
1109        swwfile = domain.filename + '.sww'
1110
1111        domain.set_datadir('.')
1112        domain.format = 'sww'
1113        domain.smooth = True
1114        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1115
1116        #
1117        domain.set_quantity('elevation', lambda x,y: -x-y)
1118        domain.set_quantity('stage', 0)       
1119
1120        B = Transmissive_boundary(domain)
1121        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1122
1123
1124        #
1125        sww = get_dataobject(domain)
1126        sww.store_connectivity()
1127        sww.store_timestep('stage')
1128
1129        domain.evolve_to_end(finaltime = 0.01)
1130        sww.store_timestep('stage')
1131
1132        cellsize = 10  #10m grid
1133
1134       
1135        #Check contents
1136        #Get NetCDF
1137
1138        fid = NetCDFFile(sww.filename, 'r')
1139
1140        # Get the variables
1141        x = fid.variables['x'][:]
1142        y = fid.variables['y'][:]
1143        z = fid.variables['elevation'][:]
1144        time = fid.variables['time'][:]
1145        stage = fid.variables['stage'][:]
1146
1147
1148        #Export to ascii/prj files
1149        sww2dem(domain.filename,
1150                quantity = 'elevation',
1151                cellsize = cellsize,
1152                easting_min = 308530,
1153                easting_max = 308570,
1154                northing_min = 6189050, 
1155                northing_max = 6189100,                               
1156                verbose = False,
1157                format = 'asc')
1158       
1159        fid.close()
1160
1161       
1162        #Check prj (meta data)
1163        prjid = open(prjfile)
1164        lines = prjid.readlines()
1165        prjid.close()
1166
1167        L = lines[0].strip().split()
1168        assert L[0].strip().lower() == 'projection'
1169        assert L[1].strip().lower() == 'utm'
1170
1171        L = lines[1].strip().split()
1172        assert L[0].strip().lower() == 'zone'
1173        assert L[1].strip().lower() == '56'
1174
1175        L = lines[2].strip().split()
1176        assert L[0].strip().lower() == 'datum'
1177        assert L[1].strip().lower() == 'wgs84'
1178
1179        L = lines[3].strip().split()
1180        assert L[0].strip().lower() == 'zunits'
1181        assert L[1].strip().lower() == 'no'
1182
1183        L = lines[4].strip().split()
1184        assert L[0].strip().lower() == 'units'
1185        assert L[1].strip().lower() == 'meters'
1186
1187        L = lines[5].strip().split()
1188        assert L[0].strip().lower() == 'spheroid'
1189        assert L[1].strip().lower() == 'wgs84'
1190
1191        L = lines[6].strip().split()
1192        assert L[0].strip().lower() == 'xshift'
1193        assert L[1].strip().lower() == '500000'
1194
1195        L = lines[7].strip().split()
1196        assert L[0].strip().lower() == 'yshift'
1197        assert L[1].strip().lower() == '10000000'
1198
1199        L = lines[8].strip().split()
1200        assert L[0].strip().lower() == 'parameters'
1201
1202
1203        #Check asc file
1204        ascid = open(ascfile)
1205        lines = ascid.readlines()
1206        ascid.close()
1207
1208        L = lines[0].strip().split()
1209        assert L[0].strip().lower() == 'ncols'
1210        assert L[1].strip().lower() == '5'
1211
1212        L = lines[1].strip().split()
1213        assert L[0].strip().lower() == 'nrows'
1214        assert L[1].strip().lower() == '6'
1215
1216        L = lines[2].strip().split()
1217        assert L[0].strip().lower() == 'xllcorner'
1218        assert allclose(float(L[1].strip().lower()), 308530)
1219
1220        L = lines[3].strip().split()
1221        assert L[0].strip().lower() == 'yllcorner'
1222        assert allclose(float(L[1].strip().lower()), 6189050)
1223
1224        L = lines[4].strip().split()
1225        assert L[0].strip().lower() == 'cellsize'
1226        assert allclose(float(L[1].strip().lower()), cellsize)
1227
1228        L = lines[5].strip().split()
1229        assert L[0].strip() == 'NODATA_value'
1230        assert L[1].strip().lower() == '-9999'
1231
1232        #Check grid values
1233        for i, line in enumerate(lines[6:]):
1234            for j, value in enumerate( line.split() ):
1235                #assert float(value) == -(10-i+j)*cellsize               
1236                assert float(value) == -(10-i+j+3)*cellsize
1237       
1238
1239
1240        #Cleanup
1241        os.remove(prjfile)
1242        os.remove(ascfile)
1243        os.remove(swwfile)
1244
1245
1246
1247    def test_sww2dem_asc_stage_reduction(self):
1248        """Test that sww information can be converted correctly to asc/prj
1249        format readable by e.g. ArcView
1250
1251        This tests the reduction of quantity stage using min
1252        """
1253
1254        import time, os
1255        from Numeric import array, zeros, allclose, Float, concatenate
1256        from Scientific.IO.NetCDF import NetCDFFile
1257
1258        #Setup
1259        self.domain.filename = 'datatest'
1260
1261        prjfile = self.domain.filename + '_stage.prj'
1262        ascfile = self.domain.filename + '_stage.asc'
1263        swwfile = self.domain.filename + '.sww'
1264
1265        self.domain.set_datadir('.')
1266        self.domain.format = 'sww'
1267        self.domain.smooth = True
1268        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1269
1270        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1271
1272
1273        sww = get_dataobject(self.domain)
1274        sww.store_connectivity()
1275        sww.store_timestep('stage')
1276
1277        self.domain.evolve_to_end(finaltime = 0.01)
1278        sww.store_timestep('stage')
1279
1280        cellsize = 0.25
1281        #Check contents
1282        #Get NetCDF
1283
1284        fid = NetCDFFile(sww.filename, 'r')
1285
1286        # Get the variables
1287        x = fid.variables['x'][:]
1288        y = fid.variables['y'][:]
1289        z = fid.variables['elevation'][:]
1290        time = fid.variables['time'][:]
1291        stage = fid.variables['stage'][:]
1292
1293
1294        #Export to ascii/prj files
1295        sww2dem(self.domain.filename,
1296                quantity = 'stage',
1297                cellsize = cellsize,
1298                reduction = min,
1299                format = 'asc')
1300
1301
1302        #Check asc file
1303        ascid = open(ascfile)
1304        lines = ascid.readlines()
1305        ascid.close()
1306
1307        L = lines[0].strip().split()
1308        assert L[0].strip().lower() == 'ncols'
1309        assert L[1].strip().lower() == '5'
1310
1311        L = lines[1].strip().split()
1312        assert L[0].strip().lower() == 'nrows'
1313        assert L[1].strip().lower() == '5'
1314
1315        L = lines[2].strip().split()
1316        assert L[0].strip().lower() == 'xllcorner'
1317        assert allclose(float(L[1].strip().lower()), 308500)
1318
1319        L = lines[3].strip().split()
1320        assert L[0].strip().lower() == 'yllcorner'
1321        assert allclose(float(L[1].strip().lower()), 6189000)
1322
1323        L = lines[4].strip().split()
1324        assert L[0].strip().lower() == 'cellsize'
1325        assert allclose(float(L[1].strip().lower()), cellsize)
1326
1327        L = lines[5].strip().split()
1328        assert L[0].strip() == 'NODATA_value'
1329        assert L[1].strip().lower() == '-9999'
1330
1331
1332        #Check grid values (where applicable)
1333        for j in range(5):
1334            if j%2 == 0:
1335                L = lines[6+j].strip().split()
1336                jj = 4-j
1337                for i in range(5):
1338                    if i%2 == 0:
1339                        index = jj/2 + i/2*3
1340                        val0 = stage[0,index]
1341                        val1 = stage[1,index]
1342
1343                        #print i, j, index, ':', L[i], val0, val1
1344                        assert allclose(float(L[i]), min(val0, val1))
1345
1346
1347        fid.close()
1348
1349        #Cleanup
1350        os.remove(prjfile)
1351        os.remove(ascfile)
1352        #os.remove(swwfile)
1353
1354
1355
1356    def test_sww2dem_asc_derived_quantity(self):
1357        """Test that sww information can be converted correctly to asc/prj
1358        format readable by e.g. ArcView
1359
1360        This tests the use of derived quantities
1361        """
1362
1363        import time, os
1364        from Numeric import array, zeros, allclose, Float, concatenate
1365        from Scientific.IO.NetCDF import NetCDFFile
1366
1367        #Setup
1368        self.domain.filename = 'datatest'
1369
1370        prjfile = self.domain.filename + '_depth.prj'
1371        ascfile = self.domain.filename + '_depth.asc'
1372        swwfile = self.domain.filename + '.sww'
1373
1374        self.domain.set_datadir('.')
1375        self.domain.format = 'sww'
1376        self.domain.smooth = True
1377        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1378        self.domain.set_quantity('stage', 0.0)
1379
1380        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1381
1382
1383        sww = get_dataobject(self.domain)
1384        sww.store_connectivity()
1385        sww.store_timestep('stage')
1386
1387        self.domain.evolve_to_end(finaltime = 0.01)
1388        sww.store_timestep('stage')
1389
1390        cellsize = 0.25
1391        #Check contents
1392        #Get NetCDF
1393
1394        fid = NetCDFFile(sww.filename, 'r')
1395
1396        # Get the variables
1397        x = fid.variables['x'][:]
1398        y = fid.variables['y'][:]
1399        z = fid.variables['elevation'][:]
1400        time = fid.variables['time'][:]
1401        stage = fid.variables['stage'][:]
1402
1403
1404        #Export to ascii/prj files
1405        sww2dem(self.domain.filename,
1406                basename_out = 'datatest_depth',
1407                quantity = 'stage - elevation',
1408                cellsize = cellsize,
1409                reduction = min,
1410                format = 'asc',
1411                verbose = False)
1412
1413
1414        #Check asc file
1415        ascid = open(ascfile)
1416        lines = ascid.readlines()
1417        ascid.close()
1418
1419        L = lines[0].strip().split()
1420        assert L[0].strip().lower() == 'ncols'
1421        assert L[1].strip().lower() == '5'
1422
1423        L = lines[1].strip().split()
1424        assert L[0].strip().lower() == 'nrows'
1425        assert L[1].strip().lower() == '5'
1426
1427        L = lines[2].strip().split()
1428        assert L[0].strip().lower() == 'xllcorner'
1429        assert allclose(float(L[1].strip().lower()), 308500)
1430
1431        L = lines[3].strip().split()
1432        assert L[0].strip().lower() == 'yllcorner'
1433        assert allclose(float(L[1].strip().lower()), 6189000)
1434
1435        L = lines[4].strip().split()
1436        assert L[0].strip().lower() == 'cellsize'
1437        assert allclose(float(L[1].strip().lower()), cellsize)
1438
1439        L = lines[5].strip().split()
1440        assert L[0].strip() == 'NODATA_value'
1441        assert L[1].strip().lower() == '-9999'
1442
1443
1444        #Check grid values (where applicable)
1445        for j in range(5):
1446            if j%2 == 0:
1447                L = lines[6+j].strip().split()
1448                jj = 4-j
1449                for i in range(5):
1450                    if i%2 == 0:
1451                        index = jj/2 + i/2*3
1452                        val0 = stage[0,index] - z[index]
1453                        val1 = stage[1,index] - z[index]
1454
1455                        #print i, j, index, ':', L[i], val0, val1
1456                        assert allclose(float(L[i]), min(val0, val1))
1457
1458
1459        fid.close()
1460
1461        #Cleanup
1462        os.remove(prjfile)
1463        os.remove(ascfile)
1464        #os.remove(swwfile)
1465
1466
1467
1468
1469
1470    def test_sww2dem_asc_missing_points(self):
1471        """Test that sww information can be converted correctly to asc/prj
1472        format readable by e.g. ArcView
1473
1474        This test includes the writing of missing values
1475        """
1476
1477        import time, os
1478        from Numeric import array, zeros, allclose, Float, concatenate
1479        from Scientific.IO.NetCDF import NetCDFFile
1480
1481        #Setup mesh not coinciding with rectangle.
1482        #This will cause missing values to occur in gridded data
1483
1484
1485        points = [                        [1.0, 1.0],
1486                              [0.5, 0.5], [1.0, 0.5],
1487                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1488
1489        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1490
1491        #Create shallow water domain
1492        domain = Domain(points, vertices)
1493        domain.default_order=2
1494
1495
1496        #Set some field values
1497        domain.set_quantity('elevation', lambda x,y: -x-y)
1498        domain.set_quantity('friction', 0.03)
1499
1500
1501        ######################
1502        # Boundary conditions
1503        B = Transmissive_boundary(domain)
1504        domain.set_boundary( {'exterior': B} )
1505
1506
1507        ######################
1508        #Initial condition - with jumps
1509
1510        bed = domain.quantities['elevation'].vertex_values
1511        stage = zeros(bed.shape, Float)
1512
1513        h = 0.3
1514        for i in range(stage.shape[0]):
1515            if i % 2 == 0:
1516                stage[i,:] = bed[i,:] + h
1517            else:
1518                stage[i,:] = bed[i,:]
1519
1520        domain.set_quantity('stage', stage)
1521        domain.distribute_to_vertices_and_edges()
1522
1523        domain.filename = 'datatest'
1524
1525        prjfile = domain.filename + '_elevation.prj'
1526        ascfile = domain.filename + '_elevation.asc'
1527        swwfile = domain.filename + '.sww'
1528
1529        domain.set_datadir('.')
1530        domain.format = 'sww'
1531        domain.smooth = True
1532
1533        domain.geo_reference = Geo_reference(56,308500,6189000)
1534
1535        sww = get_dataobject(domain)
1536        sww.store_connectivity()
1537        sww.store_timestep('stage')
1538
1539        cellsize = 0.25
1540        #Check contents
1541        #Get NetCDF
1542
1543        fid = NetCDFFile(swwfile, 'r')
1544
1545        # Get the variables
1546        x = fid.variables['x'][:]
1547        y = fid.variables['y'][:]
1548        z = fid.variables['elevation'][:]
1549        time = fid.variables['time'][:]
1550
1551        try:
1552            geo_reference = Geo_reference(NetCDFObject=fid)
1553        except AttributeError, e:
1554            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1555
1556        #Export to ascii/prj files
1557        sww2dem(domain.filename,
1558                quantity = 'elevation',
1559                cellsize = cellsize,
1560                verbose = False,
1561                format = 'asc')
1562
1563
1564        #Check asc file
1565        ascid = open(ascfile)
1566        lines = ascid.readlines()
1567        ascid.close()
1568
1569        L = lines[0].strip().split()
1570        assert L[0].strip().lower() == 'ncols'
1571        assert L[1].strip().lower() == '5'
1572
1573        L = lines[1].strip().split()
1574        assert L[0].strip().lower() == 'nrows'
1575        assert L[1].strip().lower() == '5'
1576
1577        L = lines[2].strip().split()
1578        assert L[0].strip().lower() == 'xllcorner'
1579        assert allclose(float(L[1].strip().lower()), 308500)
1580
1581        L = lines[3].strip().split()
1582        assert L[0].strip().lower() == 'yllcorner'
1583        assert allclose(float(L[1].strip().lower()), 6189000)
1584
1585        L = lines[4].strip().split()
1586        assert L[0].strip().lower() == 'cellsize'
1587        assert allclose(float(L[1].strip().lower()), cellsize)
1588
1589        L = lines[5].strip().split()
1590        assert L[0].strip() == 'NODATA_value'
1591        assert L[1].strip().lower() == '-9999'
1592
1593
1594        #Check grid values
1595        for j in range(5):
1596            L = lines[6+j].strip().split()
1597            y = (4-j) * cellsize
1598            for i in range(5):
1599                if i+j >= 4:
1600                    assert allclose(float(L[i]), -i*cellsize - y)
1601                else:
1602                    #Missing values
1603                    assert allclose(float(L[i]), -9999)
1604
1605
1606
1607        fid.close()
1608
1609        #Cleanup
1610        os.remove(prjfile)
1611        os.remove(ascfile)
1612        os.remove(swwfile)
1613
1614    def test_sww2ers_simple(self):
1615        """Test that sww information can be converted correctly to asc/prj
1616        format readable by e.g. ArcView
1617        """
1618
1619        import time, os
1620        from Numeric import array, zeros, allclose, Float, concatenate
1621        from Scientific.IO.NetCDF import NetCDFFile
1622
1623        #Setup
1624        self.domain.filename = 'datatest'
1625
1626        headerfile = self.domain.filename + '.ers'
1627        swwfile = self.domain.filename + '.sww'
1628
1629        self.domain.set_datadir('.')
1630        self.domain.format = 'sww'
1631        self.domain.smooth = True
1632        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1633
1634        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1635
1636        sww = get_dataobject(self.domain)
1637        sww.store_connectivity()
1638        sww.store_timestep('stage')
1639
1640        self.domain.evolve_to_end(finaltime = 0.01)
1641        sww.store_timestep('stage')
1642
1643        cellsize = 0.25
1644        #Check contents
1645        #Get NetCDF
1646
1647        fid = NetCDFFile(sww.filename, 'r')
1648
1649        # Get the variables
1650        x = fid.variables['x'][:]
1651        y = fid.variables['y'][:]
1652        z = fid.variables['elevation'][:]
1653        time = fid.variables['time'][:]
1654        stage = fid.variables['stage'][:]
1655
1656
1657        #Export to ers files
1658        #sww2ers(self.domain.filename,
1659        #        quantity = 'elevation',
1660        #        cellsize = cellsize,
1661        #        verbose = False)
1662               
1663        sww2dem(self.domain.filename,
1664                quantity = 'elevation',
1665                cellsize = cellsize,
1666                verbose = False,
1667                format = 'ers')
1668
1669        #Check header data
1670        from ermapper_grids import read_ermapper_header, read_ermapper_data
1671       
1672        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
1673        #print header
1674        assert header['projection'].lower() == '"utm-56"'
1675        assert header['datum'].lower() == '"wgs84"'
1676        assert header['units'].lower() == '"meters"'   
1677        assert header['value'].lower() == '"elevation"'         
1678        assert header['xdimension'] == '0.25'
1679        assert header['ydimension'] == '0.25'   
1680        assert float(header['eastings']) == 308500.0   #xllcorner
1681        assert float(header['northings']) == 6189000.0 #yllcorner       
1682        assert int(header['nroflines']) == 5
1683        assert int(header['nrofcellsperline']) == 5     
1684        assert int(header['nullcellvalue']) == 0   #?           
1685        #FIXME - there is more in the header                   
1686
1687               
1688        #Check grid data               
1689        grid = read_ermapper_data(self.domain.filename + '_elevation') 
1690       
1691        #FIXME (Ole): Why is this the desired reference grid for -x-y?
1692        ref_grid = [0,      0,     0,     0,     0,
1693                    -1,    -1.25, -1.5,  -1.75, -2.0,
1694                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
1695                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1696                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
1697                                         
1698                                         
1699        assert allclose(grid, ref_grid)
1700
1701        fid.close()
1702       
1703        #Cleanup
1704        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1705        #Done (Ole) - it was because sww2ers didn't close it's sww file
1706        os.remove(sww.filename)
1707
1708
1709
1710    def test_ferret2sww(self):
1711        """Test that georeferencing etc works when converting from
1712        ferret format (lat/lon) to sww format (UTM)
1713        """
1714        from Scientific.IO.NetCDF import NetCDFFile
1715
1716        #The test file has
1717        # LON = 150.66667, 150.83334, 151, 151.16667
1718        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1719        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1720        #
1721        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1722        # Fourth value (index==3) is -6.50198 cm
1723
1724
1725        from coordinate_transforms.redfearn import redfearn
1726
1727        fid = NetCDFFile('small_ha.nc')
1728        first_value = fid.variables['HA'][:][0,0,0]
1729        fourth_value = fid.variables['HA'][:][0,0,3]
1730
1731
1732        #Call conversion (with zero origin)
1733        ferret2sww('small', verbose=False,
1734                   origin = (56, 0, 0))
1735
1736
1737        #Work out the UTM coordinates for first point
1738        zone, e, n = redfearn(-34.5, 150.66667)
1739        #print zone, e, n
1740
1741        #Read output file 'small.sww'
1742        fid = NetCDFFile('small.sww')
1743
1744        x = fid.variables['x'][:]
1745        y = fid.variables['y'][:]
1746
1747        #Check that first coordinate is correctly represented
1748        assert allclose(x[0], e)
1749        assert allclose(y[0], n)
1750
1751        #Check first value
1752        stage = fid.variables['stage'][:]
1753        xmomentum = fid.variables['xmomentum'][:]
1754        ymomentum = fid.variables['ymomentum'][:]
1755
1756        #print ymomentum
1757
1758        assert allclose(stage[0,0], first_value/100)  #Meters
1759
1760        #Check fourth value
1761        assert allclose(stage[0,3], fourth_value/100)  #Meters
1762
1763        fid.close()
1764
1765        #Cleanup
1766        import os
1767        os.remove('small.sww')
1768
1769
1770
1771    def test_ferret2sww_2(self):
1772        """Test that georeferencing etc works when converting from
1773        ferret format (lat/lon) to sww format (UTM)
1774        """
1775        from Scientific.IO.NetCDF import NetCDFFile
1776
1777        #The test file has
1778        # LON = 150.66667, 150.83334, 151, 151.16667
1779        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1780        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1781        #
1782        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1783        # Fourth value (index==3) is -6.50198 cm
1784
1785
1786        from coordinate_transforms.redfearn import redfearn
1787
1788        fid = NetCDFFile('small_ha.nc')
1789
1790        #Pick a coordinate and a value
1791
1792        time_index = 1
1793        lat_index = 0
1794        lon_index = 2
1795
1796        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
1797        test_time = fid.variables['TIME'][:][time_index]
1798        test_lat = fid.variables['LAT'][:][lat_index]
1799        test_lon = fid.variables['LON'][:][lon_index]
1800
1801        linear_point_index = lat_index*4 + lon_index
1802        fid.close()
1803
1804        #Call conversion (with zero origin)
1805        ferret2sww('small', verbose=False,
1806                   origin = (56, 0, 0))
1807
1808
1809        #Work out the UTM coordinates for test point
1810        zone, e, n = redfearn(test_lat, test_lon)
1811
1812        #Read output file 'small.sww'
1813        fid = NetCDFFile('small.sww')
1814
1815        x = fid.variables['x'][:]
1816        y = fid.variables['y'][:]
1817
1818        #Check that test coordinate is correctly represented
1819        assert allclose(x[linear_point_index], e)
1820        assert allclose(y[linear_point_index], n)
1821
1822        #Check test value
1823        stage = fid.variables['stage'][:]
1824
1825        assert allclose(stage[time_index, linear_point_index], test_value/100)
1826
1827        fid.close()
1828
1829        #Cleanup
1830        import os
1831        os.remove('small.sww')
1832
1833
1834
1835    def test_ferret2sww3(self):
1836        """Elevation included
1837        """
1838        from Scientific.IO.NetCDF import NetCDFFile
1839
1840        #The test file has
1841        # LON = 150.66667, 150.83334, 151, 151.16667
1842        # LAT = -34.5, -34.33333, -34.16667, -34 ;
1843        # ELEVATION = [-1 -2 -3 -4
1844        #             -5 -6 -7 -8
1845        #              ...
1846        #              ...      -16]
1847        # where the top left corner is -1m,
1848        # and the ll corner is -13.0m
1849        #
1850        # First value (index=0) in small_ha.nc is 0.3400644 cm,
1851        # Fourth value (index==3) is -6.50198 cm
1852
1853        from coordinate_transforms.redfearn import redfearn
1854        import os
1855        fid1 = NetCDFFile('test_ha.nc','w')
1856        fid2 = NetCDFFile('test_ua.nc','w')
1857        fid3 = NetCDFFile('test_va.nc','w')
1858        fid4 = NetCDFFile('test_e.nc','w')
1859
1860        h1_list = [150.66667,150.83334,151.]
1861        h2_list = [-34.5,-34.33333]
1862
1863        long_name = 'LON'
1864        lat_name = 'LAT'
1865
1866        nx = 3
1867        ny = 2
1868
1869        for fid in [fid1,fid2,fid3]:
1870            fid.createDimension(long_name,nx)
1871            fid.createVariable(long_name,'d',(long_name,))
1872            fid.variables[long_name].point_spacing='uneven'
1873            fid.variables[long_name].units='degrees_east'
1874            fid.variables[long_name].assignValue(h1_list)
1875
1876            fid.createDimension(lat_name,ny)
1877            fid.createVariable(lat_name,'d',(lat_name,))
1878            fid.variables[lat_name].point_spacing='uneven'
1879            fid.variables[lat_name].units='degrees_north'
1880            fid.variables[lat_name].assignValue(h2_list)
1881
1882            fid.createDimension('TIME',2)
1883            fid.createVariable('TIME','d',('TIME',))
1884            fid.variables['TIME'].point_spacing='uneven'
1885            fid.variables['TIME'].units='seconds'
1886            fid.variables['TIME'].assignValue([0.,1.])
1887            if fid == fid3: break
1888
1889
1890        for fid in [fid4]:
1891            fid.createDimension(long_name,nx)
1892            fid.createVariable(long_name,'d',(long_name,))
1893            fid.variables[long_name].point_spacing='uneven'
1894            fid.variables[long_name].units='degrees_east'
1895            fid.variables[long_name].assignValue(h1_list)
1896
1897            fid.createDimension(lat_name,ny)
1898            fid.createVariable(lat_name,'d',(lat_name,))
1899            fid.variables[lat_name].point_spacing='uneven'
1900            fid.variables[lat_name].units='degrees_north'
1901            fid.variables[lat_name].assignValue(h2_list)
1902
1903        name = {}
1904        name[fid1]='HA'
1905        name[fid2]='UA'
1906        name[fid3]='VA'
1907        name[fid4]='ELEVATION'
1908
1909        units = {}
1910        units[fid1]='cm'
1911        units[fid2]='cm/s'
1912        units[fid3]='cm/s'
1913        units[fid4]='m'
1914
1915        values = {}
1916        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
1917        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
1918        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
1919        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
1920
1921        for fid in [fid1,fid2,fid3]:
1922          fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
1923          fid.variables[name[fid]].point_spacing='uneven'
1924          fid.variables[name[fid]].units=units[fid]
1925          fid.variables[name[fid]].assignValue(values[fid])
1926          fid.variables[name[fid]].missing_value = -99999999.
1927          if fid == fid3: break
1928
1929        for fid in [fid4]:
1930            fid.createVariable(name[fid],'d',(lat_name,long_name))
1931            fid.variables[name[fid]].point_spacing='uneven'
1932            fid.variables[name[fid]].units=units[fid]
1933            fid.variables[name[fid]].assignValue(values[fid])
1934            fid.variables[name[fid]].missing_value = -99999999.
1935
1936
1937        fid1.sync(); fid1.close()
1938        fid2.sync(); fid2.close()
1939        fid3.sync(); fid3.close()
1940        fid4.sync(); fid4.close()
1941
1942        fid1 = NetCDFFile('test_ha.nc','r')
1943        fid2 = NetCDFFile('test_e.nc','r')
1944        fid3 = NetCDFFile('test_va.nc','r')
1945
1946
1947        first_amp = fid1.variables['HA'][:][0,0,0]
1948        third_amp = fid1.variables['HA'][:][0,0,2]
1949        first_elevation = fid2.variables['ELEVATION'][0,0]
1950        third_elevation= fid2.variables['ELEVATION'][:][0,2]
1951        first_speed = fid3.variables['VA'][0,0,0]
1952        third_speed = fid3.variables['VA'][:][0,0,2]
1953
1954        fid1.close()
1955        fid2.close()
1956        fid3.close()
1957
1958        #Call conversion (with zero origin)
1959        ferret2sww('test', verbose=False,
1960                   origin = (56, 0, 0))
1961
1962        os.remove('test_va.nc')
1963        os.remove('test_ua.nc')
1964        os.remove('test_ha.nc')
1965        os.remove('test_e.nc')
1966
1967        #Read output file 'test.sww'
1968        fid = NetCDFFile('test.sww')
1969
1970
1971        #Check first value
1972        elevation = fid.variables['elevation'][:]
1973        stage = fid.variables['stage'][:]
1974        xmomentum = fid.variables['xmomentum'][:]
1975        ymomentum = fid.variables['ymomentum'][:]
1976
1977        #print ymomentum
1978        first_height = first_amp/100 - first_elevation
1979        third_height = third_amp/100 - third_elevation
1980        first_momentum=first_speed*first_height/100
1981        third_momentum=third_speed*third_height/100
1982
1983        assert allclose(ymomentum[0][0],first_momentum)  #Meters
1984        assert allclose(ymomentum[0][2],third_momentum)  #Meters
1985
1986        fid.close()
1987
1988        #Cleanup
1989        os.remove('test.sww')
1990
1991
1992
1993
1994    def test_ferret2sww_nz_origin(self):
1995        from Scientific.IO.NetCDF import NetCDFFile
1996        from coordinate_transforms.redfearn import redfearn
1997
1998        #Call conversion (with nonzero origin)
1999        ferret2sww('small', verbose=False,
2000                   origin = (56, 100000, 200000))
2001
2002
2003        #Work out the UTM coordinates for first point
2004        zone, e, n = redfearn(-34.5, 150.66667)
2005
2006        #Read output file 'small.sww'
2007        fid = NetCDFFile('small.sww', 'r')
2008
2009        x = fid.variables['x'][:]
2010        y = fid.variables['y'][:]
2011
2012        #Check that first coordinate is correctly represented
2013        assert allclose(x[0], e-100000)
2014        assert allclose(y[0], n-200000)
2015
2016        fid.close()
2017
2018        #Cleanup
2019        os.remove('small.sww')
2020
2021
2022
2023    def test_sww_extent(self):
2024        """Not a test, rather a look at the sww format
2025        """
2026
2027        import time, os
2028        from Numeric import array, zeros, allclose, Float, concatenate
2029        from Scientific.IO.NetCDF import NetCDFFile
2030
2031        self.domain.filename = 'datatest' + str(id(self))
2032        self.domain.format = 'sww'
2033        self.domain.smooth = True
2034        self.domain.reduction = mean
2035        self.domain.set_datadir('.')
2036
2037
2038        sww = get_dataobject(self.domain)
2039        sww.store_connectivity()
2040        sww.store_timestep('stage')
2041        self.domain.time = 2.
2042
2043        #Modify stage at second timestep
2044        stage = self.domain.quantities['stage'].vertex_values
2045        self.domain.set_quantity('stage', stage/2)
2046
2047        sww.store_timestep('stage')
2048
2049        file_and_extension_name = self.domain.filename + ".sww"
2050        #print "file_and_extension_name",file_and_extension_name
2051        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2052               extent_sww(file_and_extension_name )
2053
2054        assert allclose(xmin, 0.0)
2055        assert allclose(xmax, 1.0)
2056        assert allclose(ymin, 0.0)
2057        assert allclose(ymax, 1.0)
2058        assert allclose(stagemin, -0.85)
2059        assert allclose(stagemax, 0.15)
2060
2061
2062        #Cleanup
2063        os.remove(sww.filename)
2064
2065
2066
2067    def test_sww2domain(self):
2068        ################################################
2069        #Create a test domain, and evolve and save it.
2070        ################################################
2071        from mesh_factory import rectangular
2072        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2073             Constant_height, Time_boundary, Transmissive_boundary
2074        from Numeric import array
2075
2076        #Create basic mesh
2077
2078        yiel=0.01
2079        points, vertices, boundary = rectangular(10,10)
2080
2081        #Create shallow water domain
2082        domain = Domain(points, vertices, boundary)
2083        domain.geo_reference = Geo_reference(56,11,11)
2084        domain.smooth = False
2085        domain.visualise = False
2086        domain.store = True
2087        domain.filename = 'bedslope'
2088        domain.default_order=2
2089        #Bed-slope and friction
2090        domain.set_quantity('elevation', lambda x,y: -x/3)
2091        domain.set_quantity('friction', 0.1)
2092        # Boundary conditions
2093        from math import sin, pi
2094        Br = Reflective_boundary(domain)
2095        Bt = Transmissive_boundary(domain)
2096        Bd = Dirichlet_boundary([0.2,0.,0.])
2097        Bw = Time_boundary(domain=domain,
2098                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2099
2100        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2101        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2102
2103        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2104        #Initial condition
2105        h = 0.05
2106        elevation = domain.quantities['elevation'].vertex_values
2107        domain.set_quantity('stage', elevation + h)
2108        #elevation = domain.get_quantity('elevation')
2109        #domain.set_quantity('stage', elevation + h)
2110
2111        domain.check_integrity()
2112        #Evolution
2113        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2114        #    domain.write_time()
2115            pass
2116
2117
2118        ##########################################
2119        #Import the example's file as a new domain
2120        ##########################################
2121        from data_manager import sww2domain
2122        from Numeric import allclose
2123        import os
2124
2125        filename = domain.datadir+os.sep+domain.filename+'.sww'
2126        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2127        #points, vertices, boundary = rectangular(15,15)
2128        #domain2.boundary = boundary
2129        ###################
2130        ##NOW TEST IT!!!
2131        ###################
2132
2133        bits = ['vertex_coordinates']
2134        for quantity in ['elevation']+domain.quantities_to_be_stored:
2135            bits.append('quantities["%s"].get_integral()'%quantity)
2136            bits.append('get_quantity("%s")'%quantity)
2137
2138        for bit in bits:
2139            #print 'testing that domain.'+bit+' has been restored'
2140            #print bit
2141        #print 'done'
2142            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2143
2144        ######################################
2145        #Now evolve them both, just to be sure
2146        ######################################x
2147        visualise = False
2148        #visualise = True
2149        domain.visualise = visualise
2150        domain.time = 0.
2151        from time import sleep
2152
2153        final = .1
2154        domain.set_quantity('friction', 0.1)
2155        domain.store = False
2156        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2157
2158
2159        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2160            if visualise: sleep(1.)
2161            #domain.write_time()
2162            pass
2163
2164        final = final - (domain2.starttime-domain.starttime)
2165        #BUT since domain1 gets time hacked back to 0:
2166        final = final + (domain2.starttime-domain.starttime)
2167
2168        domain2.smooth = False
2169        domain2.visualise = visualise
2170        domain2.store = False
2171        domain2.default_order=2
2172        domain2.set_quantity('friction', 0.1)
2173        #Bed-slope and friction
2174        # Boundary conditions
2175        Bd2=Dirichlet_boundary([0.2,0.,0.])
2176        domain2.boundary = domain.boundary
2177        #print 'domain2.boundary'
2178        #print domain2.boundary
2179        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2180        #domain2.set_boundary({'exterior': Bd})
2181
2182        domain2.check_integrity()
2183
2184        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2185            if visualise: sleep(1.)
2186            #domain2.write_time()
2187            pass
2188
2189        ###################
2190        ##NOW TEST IT!!!
2191        ##################
2192
2193        bits = [ 'vertex_coordinates']
2194
2195        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2196            bits.append('quantities["%s"].get_integral()'%quantity)
2197            bits.append('get_quantity("%s")'%quantity)
2198
2199        for bit in bits:
2200            #print bit
2201            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2202
2203
2204    def test_sww2domain2(self):
2205        ##################################################################
2206        #Same as previous test, but this checks how NaNs are handled.
2207        ##################################################################
2208
2209
2210        from mesh_factory import rectangular
2211        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2212             Constant_height, Time_boundary, Transmissive_boundary
2213        from Numeric import array
2214
2215        #Create basic mesh
2216        points, vertices, boundary = rectangular(2,2)
2217
2218        #Create shallow water domain
2219        domain = Domain(points, vertices, boundary)
2220        domain.smooth = False
2221        domain.visualise = False
2222        domain.store = True
2223        domain.filename = 'bedslope'
2224        domain.default_order=2
2225        domain.quantities_to_be_stored=['stage']
2226
2227        domain.set_quantity('elevation', lambda x,y: -x/3)
2228        domain.set_quantity('friction', 0.1)
2229
2230        from math import sin, pi
2231        Br = Reflective_boundary(domain)
2232        Bt = Transmissive_boundary(domain)
2233        Bd = Dirichlet_boundary([0.2,0.,0.])
2234        Bw = Time_boundary(domain=domain,
2235                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2236
2237        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2238
2239        h = 0.05
2240        elevation = domain.quantities['elevation'].vertex_values
2241        domain.set_quantity('stage', elevation + h)
2242
2243        domain.check_integrity()
2244
2245        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2246            pass
2247            #domain.write_time()
2248
2249
2250
2251        ##################################
2252        #Import the file as a new domain
2253        ##################################
2254        from data_manager import sww2domain
2255        from Numeric import allclose
2256        import os
2257
2258        filename = domain.datadir+os.sep+domain.filename+'.sww'
2259
2260        #Fail because NaNs are present
2261        try:
2262            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2263            assert True == False
2264        except:
2265            #Now import it, filling NaNs to be 0
2266            filler = 0
2267            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2268        bits = [ 'geo_reference.get_xllcorner()',
2269                'geo_reference.get_yllcorner()',
2270                'vertex_coordinates']
2271
2272        for quantity in ['elevation']+domain.quantities_to_be_stored:
2273            bits.append('quantities["%s"].get_integral()'%quantity)
2274            bits.append('get_quantity("%s")'%quantity)
2275
2276        for bit in bits:
2277        #    print 'testing that domain.'+bit+' has been restored'
2278            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2279
2280        assert max(max(domain2.get_quantity('xmomentum')))==filler
2281        assert min(min(domain2.get_quantity('xmomentum')))==filler
2282        assert max(max(domain2.get_quantity('ymomentum')))==filler
2283        assert min(min(domain2.get_quantity('ymomentum')))==filler
2284
2285        #print 'passed'
2286
2287        #cleanup
2288        #import os
2289        #os.remove(domain.datadir+'/'+domain.filename+'.sww')
2290
2291
2292    #def test_weed(self):
2293        from data_manager import weed
2294
2295        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2296        volumes1 = [[0,1,2],[3,4,5]]
2297        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2298        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2299
2300        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2301
2302        assert len(points2)==len(coordinates2)
2303        for i in range(len(coordinates2)):
2304            coordinate = tuple(coordinates2[i])
2305            assert points2.has_key(coordinate)
2306            points2[coordinate]=i
2307
2308        for triangle in volumes1:
2309            for coordinate in triangle:
2310                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2311                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2312
2313
2314    #FIXME This fails - smooth makes the comparism too hard for allclose
2315    def ztest_sww2domain3(self):
2316        ################################################
2317        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2318        ################################################
2319        from mesh_factory import rectangular
2320        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2321             Constant_height, Time_boundary, Transmissive_boundary
2322        from Numeric import array
2323        #Create basic mesh
2324
2325        yiel=0.01
2326        points, vertices, boundary = rectangular(10,10)
2327
2328        #Create shallow water domain
2329        domain = Domain(points, vertices, boundary)
2330        domain.geo_reference = Geo_reference(56,11,11)
2331        domain.smooth = True
2332        domain.visualise = False
2333        domain.store = True
2334        domain.filename = 'bedslope'
2335        domain.default_order=2
2336        #Bed-slope and friction
2337        domain.set_quantity('elevation', lambda x,y: -x/3)
2338        domain.set_quantity('friction', 0.1)
2339        # Boundary conditions
2340        from math import sin, pi
2341        Br = Reflective_boundary(domain)
2342        Bt = Transmissive_boundary(domain)
2343        Bd = Dirichlet_boundary([0.2,0.,0.])
2344        Bw = Time_boundary(domain=domain,
2345                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2346
2347        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2348
2349        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2350        #Initial condition
2351        h = 0.05
2352        elevation = domain.quantities['elevation'].vertex_values
2353        domain.set_quantity('stage', elevation + h)
2354        #elevation = domain.get_quantity('elevation')
2355        #domain.set_quantity('stage', elevation + h)
2356
2357        domain.check_integrity()
2358        #Evolution
2359        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2360        #    domain.write_time()
2361            pass
2362
2363
2364        ##########################################
2365        #Import the example's file as a new domain
2366        ##########################################
2367        from data_manager import sww2domain
2368        from Numeric import allclose
2369        import os
2370
2371        filename = domain.datadir+os.sep+domain.filename+'.sww'
2372        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2373        #points, vertices, boundary = rectangular(15,15)
2374        #domain2.boundary = boundary
2375        ###################
2376        ##NOW TEST IT!!!
2377        ###################
2378
2379        #FIXME smooth domain so that they can be compared
2380
2381
2382        bits = []#'vertex_coordinates']
2383        for quantity in ['elevation']+domain.quantities_to_be_stored:
2384            bits.append('quantities["%s"].get_integral()'%quantity)
2385            #bits.append('get_quantity("%s")'%quantity)
2386
2387        for bit in bits:
2388            #print 'testing that domain.'+bit+' has been restored'
2389            #print bit
2390            #print 'done'
2391            #print ('domain.'+bit), eval('domain.'+bit)
2392            #print ('domain2.'+bit), eval('domain2.'+bit)
2393            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2394            pass
2395
2396        ######################################
2397        #Now evolve them both, just to be sure
2398        ######################################x
2399        visualise = False
2400        visualise = True
2401        domain.visualise = visualise
2402        domain.time = 0.
2403        from time import sleep
2404
2405        final = .5
2406        domain.set_quantity('friction', 0.1)
2407        domain.store = False
2408        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2409
2410        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2411            if visualise: sleep(.03)
2412            #domain.write_time()
2413            pass
2414
2415        domain2.smooth = True
2416        domain2.visualise = visualise
2417        domain2.store = False
2418        domain2.default_order=2
2419        domain2.set_quantity('friction', 0.1)
2420        #Bed-slope and friction
2421        # Boundary conditions
2422        Bd2=Dirichlet_boundary([0.2,0.,0.])
2423        Br2 = Reflective_boundary(domain2)
2424        domain2.boundary = domain.boundary
2425        #print 'domain2.boundary'
2426        #print domain2.boundary
2427        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2428        #domain2.boundary = domain.boundary
2429        #domain2.set_boundary({'exterior': Bd})
2430
2431        domain2.check_integrity()
2432
2433        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2434            if visualise: sleep(.03)
2435            #domain2.write_time()
2436            pass
2437
2438        ###################
2439        ##NOW TEST IT!!!
2440        ##################
2441
2442        bits = [ 'vertex_coordinates']
2443
2444        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2445            #bits.append('quantities["%s"].get_integral()'%quantity)
2446            bits.append('get_quantity("%s")'%quantity)
2447
2448        for bit in bits:
2449            print bit
2450            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2451
2452
2453    def test_decimate_dem(self):
2454        """Test decimation of dem file
2455        """
2456
2457        import os
2458        from Numeric import ones, allclose, Float, arange
2459        from Scientific.IO.NetCDF import NetCDFFile
2460
2461        #Write test dem file
2462        root = 'decdemtest'
2463
2464        filename = root + '.dem'
2465        fid = NetCDFFile(filename, 'w')
2466
2467        fid.institution = 'Geoscience Australia'
2468        fid.description = 'NetCDF DEM format for compact and portable ' +\
2469                          'storage of spatial point data'
2470
2471        nrows = 15
2472        ncols = 18
2473
2474        fid.ncols = ncols
2475        fid.nrows = nrows
2476        fid.xllcorner = 2000.5
2477        fid.yllcorner = 3000.5
2478        fid.cellsize = 25
2479        fid.NODATA_value = -9999
2480
2481        fid.zone = 56
2482        fid.false_easting = 0.0
2483        fid.false_northing = 0.0
2484        fid.projection = 'UTM'
2485        fid.datum = 'WGS84'
2486        fid.units = 'METERS'
2487
2488        fid.createDimension('number_of_points', nrows*ncols)
2489
2490        fid.createVariable('elevation', Float, ('number_of_points',))
2491
2492        elevation = fid.variables['elevation']
2493
2494        elevation[:] = (arange(nrows*ncols))
2495
2496        fid.close()
2497
2498        #generate the elevation values expected in the decimated file
2499        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2500                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2501                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2502                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2503                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2504                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2505                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2506                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2507                         (144+145+146+162+163+164+180+181+182) / 9.0,
2508                         (148+149+150+166+167+168+184+185+186) / 9.0,
2509                         (152+153+154+170+171+172+188+189+190) / 9.0,
2510                         (156+157+158+174+175+176+192+193+194) / 9.0,
2511                         (216+217+218+234+235+236+252+253+254) / 9.0,
2512                         (220+221+222+238+239+240+256+257+258) / 9.0,
2513                         (224+225+226+242+243+244+260+261+262) / 9.0,
2514                         (228+229+230+246+247+248+264+265+266) / 9.0]
2515
2516        #generate a stencil for computing the decimated values
2517        stencil = ones((3,3), Float) / 9.0
2518
2519        decimate_dem(root, stencil=stencil, cellsize_new=100)
2520
2521        #Open decimated NetCDF file
2522        fid = NetCDFFile(root + '_100.dem', 'r')
2523
2524        # Get decimated elevation
2525        elevation = fid.variables['elevation']
2526
2527        #Check values
2528        assert allclose(elevation, ref_elevation)
2529
2530        #Cleanup
2531        fid.close()
2532
2533        os.remove(root + '.dem')
2534        os.remove(root + '_100.dem')
2535
2536    def test_decimate_dem_NODATA(self):
2537        """Test decimation of dem file that includes NODATA values
2538        """
2539
2540        import os
2541        from Numeric import ones, allclose, Float, arange, reshape
2542        from Scientific.IO.NetCDF import NetCDFFile
2543
2544        #Write test dem file
2545        root = 'decdemtest'
2546
2547        filename = root + '.dem'
2548        fid = NetCDFFile(filename, 'w')
2549
2550        fid.institution = 'Geoscience Australia'
2551        fid.description = 'NetCDF DEM format for compact and portable ' +\
2552                          'storage of spatial point data'
2553
2554        nrows = 15
2555        ncols = 18
2556        NODATA_value = -9999
2557
2558        fid.ncols = ncols
2559        fid.nrows = nrows
2560        fid.xllcorner = 2000.5
2561        fid.yllcorner = 3000.5
2562        fid.cellsize = 25
2563        fid.NODATA_value = NODATA_value
2564
2565        fid.zone = 56
2566        fid.false_easting = 0.0
2567        fid.false_northing = 0.0
2568        fid.projection = 'UTM'
2569        fid.datum = 'WGS84'
2570        fid.units = 'METERS'
2571
2572        fid.createDimension('number_of_points', nrows*ncols)
2573
2574        fid.createVariable('elevation', Float, ('number_of_points',))
2575
2576        elevation = fid.variables['elevation']
2577
2578        #generate initial elevation values
2579        elevation_tmp = (arange(nrows*ncols))
2580        #add some NODATA values
2581        elevation_tmp[0]   = NODATA_value
2582        elevation_tmp[95]  = NODATA_value
2583        elevation_tmp[188] = NODATA_value
2584        elevation_tmp[189] = NODATA_value
2585        elevation_tmp[190] = NODATA_value
2586        elevation_tmp[209] = NODATA_value
2587        elevation_tmp[252] = NODATA_value
2588
2589        elevation[:] = elevation_tmp
2590
2591        fid.close()
2592
2593        #generate the elevation values expected in the decimated file
2594        ref_elevation = [NODATA_value,
2595                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2596                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2597                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2598                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2599                         NODATA_value,
2600                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2601                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2602                         (144+145+146+162+163+164+180+181+182) / 9.0,
2603                         (148+149+150+166+167+168+184+185+186) / 9.0,
2604                         NODATA_value,
2605                         (156+157+158+174+175+176+192+193+194) / 9.0,
2606                         NODATA_value,
2607                         (220+221+222+238+239+240+256+257+258) / 9.0,
2608                         (224+225+226+242+243+244+260+261+262) / 9.0,
2609                         (228+229+230+246+247+248+264+265+266) / 9.0]
2610
2611        #generate a stencil for computing the decimated values
2612        stencil = ones((3,3), Float) / 9.0
2613
2614        decimate_dem(root, stencil=stencil, cellsize_new=100)
2615
2616        #Open decimated NetCDF file
2617        fid = NetCDFFile(root + '_100.dem', 'r')
2618
2619        # Get decimated elevation
2620        elevation = fid.variables['elevation']
2621
2622        #Check values
2623        assert allclose(elevation, ref_elevation)
2624
2625        #Cleanup
2626        fid.close()
2627
2628        os.remove(root + '.dem')
2629        os.remove(root + '_100.dem')
2630
2631    def xxxtestz_sww2ers_real(self):
2632        """Test that sww information can be converted correctly to asc/prj
2633        format readable by e.g. ArcView
2634        """
2635
2636        import time, os
2637        from Numeric import array, zeros, allclose, Float, concatenate
2638        from Scientific.IO.NetCDF import NetCDFFile
2639
2640        # the memory optimised least squares
2641        #  cellsize = 20,   # this one seems to hang
2642        #  cellsize = 200000, # Ran 1 test in 269.703s
2643                                #Ran 1 test in 267.344s
2644        #  cellsize = 20000,  # Ran 1 test in 460.922s
2645        #  cellsize = 2000   #Ran 1 test in 5340.250s
2646        #  cellsize = 200   #this one seems to hang, building matirx A
2647
2648        # not optimised
2649        # seems to hang
2650        #  cellsize = 2000   # Ran 1 test in 5334.563s
2651        #Export to ascii/prj files
2652        sww2dem('karratha_100m',
2653                quantity = 'depth',
2654                cellsize = 200000,
2655                verbose = True)
2656
2657    def test_read_asc(self):
2658        """Test conversion from dem in ascii format to native NetCDF xya format
2659        """
2660
2661        import time, os
2662        from Numeric import array, zeros, allclose, Float, concatenate
2663        from Scientific.IO.NetCDF import NetCDFFile
2664
2665        import data_manager
2666        #Write test asc file
2667        filename = tempfile.mktemp(".000")
2668        fid = open(filename, 'w')
2669        fid.write("""ncols         7
2670nrows         4
2671xllcorner     2000.5
2672yllcorner     3000.5
2673cellsize      25
2674NODATA_value  -9999
2675    97.921    99.285   125.588   180.830   258.645   342.872   415.836
2676   473.157   514.391   553.893   607.120   678.125   777.283   883.038
2677   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
2678   502.645   516.230   504.739   450.604   388.500   338.097   514.980
2679""")
2680        fid.close()
2681        bath_metadata, grid = \
2682                   data_manager._read_asc(filename, verbose=False)
2683        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
2684        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
2685        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
2686        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
2687        self.failUnless(grid[0][0]  == 97.921,  'Failed')
2688        self.failUnless(grid[3][6]  == 514.980,  'Failed')
2689
2690        os.remove(filename)
2691       
2692    def test_asc_csiro2sww(self):
2693        import tempfile
2694
2695        bath_dir = tempfile.mkdtemp()
2696        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
2697        #bath_dir = 'bath_data_manager_test'
2698        #print "os.getcwd( )",os.getcwd( )
2699        elevation_dir =  tempfile.mkdtemp()
2700        #elevation_dir = 'elev_expanded'
2701        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
2702        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
2703       
2704        fid = open(bath_dir_filename, 'w')
2705        fid.write(""" ncols             3
2706 nrows             2
2707 xllcorner    148.00000
2708 yllcorner    -38.00000
2709 cellsize       0.25
2710 nodata_value   -9999.0
2711    9000.000    60.000    3000.0
2712   -10.000   -20.000   -30.000
2713""")
2714        fid.close()
2715       
2716        fid = open(elevation_dir_filename1, 'w')
2717        fid.write(""" ncols             3
2718 nrows             2
2719 xllcorner    148.00000
2720 yllcorner    -38.00000
2721 cellsize       0.25
2722 nodata_value   -9999.0
2723    90.000    60.000    30.0
2724     0.000     0.000     0.000
2725""")
2726        fid.close()
2727
2728        fid = open(elevation_dir_filename2, 'w')
2729        fid.write(""" ncols             3
2730 nrows             2
2731 xllcorner    148.00000
2732 yllcorner    -38.00000
2733 cellsize       0.25
2734 nodata_value   -9999.0
2735    90.000    60.000    30.0
2736    10.000    10.000    10.000
2737""")
2738        fid.close()
2739       
2740        sww_file = 'test.sww'
2741        asc_csiro2sww(bath_dir,elevation_dir, sww_file)
2742
2743        # check the sww file
2744       
2745        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
2746        x = fid.variables['x'][:]
2747        y = fid.variables['y'][:]
2748        geo_ref = Geo_reference(NetCDFObject=fid)
2749        #print "geo_ref",geo_ref
2750        x_ref = geo_ref.get_xllcorner()
2751        y_ref = geo_ref.get_yllcorner()
2752        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
2753        assert allclose(x_ref, 587798.418) # (-38, 148)
2754        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
2755       
2756        #Zone:   55   
2757        #Easting:  588095.674  Northing: 5821451.722
2758        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
2759        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
2760
2761        #Zone:   55   
2762        #Easting:  632145.632  Northing: 5820863.269
2763        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
2764        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
2765
2766        #Zone:   55   
2767        #Easting:  609748.788  Northing: 5793447.860
2768        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
2769        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.860 - y_ref))
2770
2771        fid.close()
2772   
2773        #tidy up
2774        os.remove(bath_dir_filename)
2775        os.rmdir(bath_dir)
2776
2777        os.remove(elevation_dir_filename1)
2778        os.remove(elevation_dir_filename2)
2779        os.rmdir(elevation_dir)
2780
2781
2782        # remove sww file
2783        #co-ords
2784        #Zone:   55   
2785        #Easting:  587798.418  Northing: 5793713.242
2786        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 ''
2787        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494
2788
2789        #Zone:   55   
2790        #Easting:  587798.418  Northing: 5793713.242
2791        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 ''
2792        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494
2793
2794
2795        #Zone:   55   
2796        #Easting:  631699.669  Northing: 5793123.477
2797        #Latitude:   -38 0 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
2798
2799       
2800#-------------------------------------------------------------
2801if __name__ == "__main__":
2802    suite = unittest.makeSuite(Test_Data_Manager,'test')
2803    #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest')
2804    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')   
2805    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
2806    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
2807    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
2808    runner = unittest.TextTestRunner()
2809    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.