source: trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py @ 7776

Last change on this file since 7776 was 7776, checked in by hudson, 14 years ago

Removed redundant data_manager class. Unit tests are running, but may fail.

File size: 41.1 KB
Line 
1import unittest, os
2import numpy as num
3
4from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
5from anuga.shallow_water.shallow_water_domain import Domain
6from anuga.config import netcdf_mode_r
7
8from anuga.coordinate_transforms.geo_reference import Geo_reference, \
9     DEFAULT_ZONE
10
11from anuga.file.sww import SWW_file
12     
13# boundary functions
14from anuga.shallow_water.boundaries import Reflective_boundary, \
15            Field_boundary, Transmissive_momentum_set_stage_boundary, \
16            Transmissive_stage_zero_momentum_boundary
17from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
18     import Transmissive_boundary, Dirichlet_boundary, \
19            Time_boundary, File_boundary, AWI_boundary
20
21# local modules
22from sww2dem import sww2dem
23
24class Test_Sww2Dem(unittest.TestCase):
25    def setUp(self):
26        import time
27        from mesh_factory import rectangular
28       
29        # Create basic mesh
30        points, vertices, boundary = rectangular(2, 2)
31
32        # Create shallow water domain
33        domain = Domain(points, vertices, boundary)
34        domain.default_order = 2
35
36        # Set some field values
37        domain.set_quantity('elevation', lambda x,y: -x)
38        domain.set_quantity('friction', 0.03)
39
40
41        ######################
42        # Boundary conditions
43        B = Transmissive_boundary(domain)
44        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
45
46
47        ######################
48        #Initial condition - with jumps
49        bed = domain.quantities['elevation'].vertex_values
50        stage = num.zeros(bed.shape, num.float)
51
52        h = 0.3
53        for i in range(stage.shape[0]):
54            if i % 2 == 0:
55                stage[i,:] = bed[i,:] + h
56            else:
57                stage[i,:] = bed[i,:]
58
59        domain.set_quantity('stage', stage)
60
61
62        domain.distribute_to_vertices_and_edges()               
63        self.domain = domain
64
65        C = domain.get_vertex_coordinates()
66        self.X = C[:,0:6:2].copy()
67        self.Y = C[:,1:6:2].copy()
68
69        self.F = bed
70        self.verbose = False
71
72
73    def tearDown(self):
74        pass
75
76
77    def test_sww2dem_asc_elevation_depth(self):
78        """test_sww2dem_asc_elevation_depth
79       
80        Test that sww information can be converted correctly to asc/prj
81        format readable by e.g. ArcView
82       
83        Also check geo_reference is correct
84        """
85
86        import time, os
87        from Scientific.IO.NetCDF import NetCDFFile
88
89        # Setup
90        self.domain.set_name('datatest')
91
92        prjfile = self.domain.get_name() + '_elevation.prj'
93        ascfile = self.domain.get_name() + '_elevation.asc'
94        swwfile = self.domain.get_name() + '.sww'
95
96        self.domain.set_datadir('.')
97        self.domain.format = 'sww'
98        self.domain.smooth = True
99        self.domain.set_quantity('elevation', lambda x,y: -x-y)
100        self.domain.set_quantity('stage', 1.0)
101
102        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
103
104        sww = SWW_file(self.domain)
105        sww.store_connectivity()
106        sww.store_timestep()
107
108
109        self.domain.evolve_to_end(finaltime = 0.01)
110        sww.store_timestep()
111
112        cellsize = 0.25
113        #Check contents
114        #Get NetCDF
115
116        fid = NetCDFFile(sww.filename, netcdf_mode_r)
117
118        # Get the variables
119        x = fid.variables['x'][:]
120        y = fid.variables['y'][:]
121        z = fid.variables['elevation'][:]
122        time = fid.variables['time'][:]
123        stage = fid.variables['stage'][:]
124       
125        # Check georeferencig: zone, xllcorner and yllcorner
126        assert fid.zone[0] == 56
127        assert fid.xllcorner[0] == 308500
128        assert fid.yllcorner[0] == 6189000
129               
130
131        fid.close()
132
133        #Export to ascii/prj files
134        sww2dem(self.domain.get_name(),
135                quantity = 'elevation',
136                cellsize = cellsize,
137                number_of_decimal_places = 9,
138                verbose = self.verbose,
139                format = 'asc')
140
141        #Check prj (meta data)
142        prjid = open(prjfile)
143        lines = prjid.readlines()
144        prjid.close()
145
146        L = lines[0].strip().split()
147        assert L[0].strip().lower() == 'projection'
148        assert L[1].strip().lower() == 'utm'
149
150        L = lines[1].strip().split()
151        assert L[0].strip().lower() == 'zone'
152        assert L[1].strip().lower() == '56'
153
154        L = lines[2].strip().split()
155        assert L[0].strip().lower() == 'datum'
156        assert L[1].strip().lower() == 'wgs84'
157
158        L = lines[3].strip().split()
159        assert L[0].strip().lower() == 'zunits'
160        assert L[1].strip().lower() == 'no'
161
162        L = lines[4].strip().split()
163        assert L[0].strip().lower() == 'units'
164        assert L[1].strip().lower() == 'meters'
165
166        L = lines[5].strip().split()
167        assert L[0].strip().lower() == 'spheroid'
168        assert L[1].strip().lower() == 'wgs84'
169
170        L = lines[6].strip().split()
171        assert L[0].strip().lower() == 'xshift'
172        assert L[1].strip().lower() == '500000'
173
174        L = lines[7].strip().split()
175        assert L[0].strip().lower() == 'yshift'
176        assert L[1].strip().lower() == '10000000'
177
178        L = lines[8].strip().split()
179        assert L[0].strip().lower() == 'parameters'
180
181
182        #Check asc file
183        ascid = open(ascfile)
184        lines = ascid.readlines()
185        ascid.close()
186
187        L = lines[0].strip().split()
188        assert L[0].strip().lower() == 'ncols'
189        assert L[1].strip().lower() == '5'
190
191        L = lines[1].strip().split()
192        assert L[0].strip().lower() == 'nrows'
193        assert L[1].strip().lower() == '5'
194
195        L = lines[2].strip().split()
196        assert L[0].strip().lower() == 'xllcorner'
197        assert num.allclose(float(L[1].strip().lower()), 308500)
198
199        L = lines[3].strip().split()
200        assert L[0].strip().lower() == 'yllcorner'
201        assert num.allclose(float(L[1].strip().lower()), 6189000)
202
203        L = lines[4].strip().split()
204        assert L[0].strip().lower() == 'cellsize'
205        assert num.allclose(float(L[1].strip().lower()), cellsize)
206
207        L = lines[5].strip().split()
208        assert L[0].strip() == 'NODATA_value'
209        assert L[1].strip().lower() == '-9999'
210
211        #Check grid values
212        for j in range(5):
213            L = lines[6+j].strip().split()
214            y = (4-j) * cellsize
215            for i in range(5):
216                assert num.allclose(float(L[i]), -i*cellsize - y)
217               
218        #Cleanup
219        os.remove(prjfile)
220        os.remove(ascfile)
221
222        #Export to ascii/prj files
223        sww2dem(self.domain.get_name(),
224                quantity = 'depth',
225                cellsize = cellsize,
226                number_of_decimal_places = 9,
227                verbose = self.verbose,
228                format = 'asc')
229       
230        #Check asc file
231        ascfile = self.domain.get_name() + '_depth.asc'
232        prjfile = self.domain.get_name() + '_depth.prj'
233        ascid = open(ascfile)
234        lines = ascid.readlines()
235        ascid.close()
236
237        L = lines[0].strip().split()
238        assert L[0].strip().lower() == 'ncols'
239        assert L[1].strip().lower() == '5'
240
241        L = lines[1].strip().split()
242        assert L[0].strip().lower() == 'nrows'
243        assert L[1].strip().lower() == '5'
244
245        L = lines[2].strip().split()
246        assert L[0].strip().lower() == 'xllcorner'
247        assert num.allclose(float(L[1].strip().lower()), 308500)
248
249        L = lines[3].strip().split()
250        assert L[0].strip().lower() == 'yllcorner'
251        assert num.allclose(float(L[1].strip().lower()), 6189000)
252
253        L = lines[4].strip().split()
254        assert L[0].strip().lower() == 'cellsize'
255        assert num.allclose(float(L[1].strip().lower()), cellsize)
256
257        L = lines[5].strip().split()
258        assert L[0].strip() == 'NODATA_value'
259        assert L[1].strip().lower() == '-9999'
260
261        #Check grid values
262        for j in range(5):
263            L = lines[6+j].strip().split()
264            y = (4-j) * cellsize
265            for i in range(5):
266                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
267
268
269        #Cleanup
270        os.remove(prjfile)
271        os.remove(ascfile)
272        os.remove(swwfile)
273
274
275
276    def test_sww2dem_larger(self):
277        """Test that sww information can be converted correctly to asc/prj
278        format readable by e.g. ArcView. Here:
279
280        ncols         11
281        nrows         11
282        xllcorner     308500
283        yllcorner     6189000
284        cellsize      10.000000
285        NODATA_value  -9999
286        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
287         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
288         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
289         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
290         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
291         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
292         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
293         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
294         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
295         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
296           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
297
298        """
299
300        import time, os
301        from Scientific.IO.NetCDF import NetCDFFile
302
303        #Create basic mesh (100m x 100m)
304        points, vertices, boundary = rectangular(2, 2, 100, 100)
305
306        #Create shallow water domain
307        domain = Domain(points, vertices, boundary)
308        domain.default_order = 2
309
310        domain.set_name('datatest')
311
312        prjfile = domain.get_name() + '_elevation.prj'
313        ascfile = domain.get_name() + '_elevation.asc'
314        swwfile = domain.get_name() + '.sww'
315
316        domain.set_datadir('.')
317        domain.format = 'sww'
318        domain.smooth = True
319        domain.geo_reference = Geo_reference(56, 308500, 6189000)
320
321        #
322        domain.set_quantity('elevation', lambda x,y: -x-y)
323        domain.set_quantity('stage', 0)
324
325        B = Transmissive_boundary(domain)
326        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
327
328
329        #
330        sww = SWW_file(domain)
331        sww.store_connectivity()
332        sww.store_timestep()
333       
334        domain.tight_slope_limiters = 1
335        domain.evolve_to_end(finaltime = 0.01)
336        sww.store_timestep()
337
338        cellsize = 10  #10m grid
339
340
341        #Check contents
342        #Get NetCDF
343
344        fid = NetCDFFile(sww.filename, netcdf_mode_r)
345
346        # Get the variables
347        x = fid.variables['x'][:]
348        y = fid.variables['y'][:]
349        z = fid.variables['elevation'][:]
350        time = fid.variables['time'][:]
351        stage = fid.variables['stage'][:]
352
353
354        #Export to ascii/prj files
355        sww2dem(domain.get_name(),
356                quantity = 'elevation',
357                cellsize = cellsize,
358                number_of_decimal_places = 9,
359                verbose = self.verbose,
360                format = 'asc',
361                block_size=2)
362
363
364        #Check prj (meta data)
365        prjid = open(prjfile)
366        lines = prjid.readlines()
367        prjid.close()
368
369        L = lines[0].strip().split()
370        assert L[0].strip().lower() == 'projection'
371        assert L[1].strip().lower() == 'utm'
372
373        L = lines[1].strip().split()
374        assert L[0].strip().lower() == 'zone'
375        assert L[1].strip().lower() == '56'
376
377        L = lines[2].strip().split()
378        assert L[0].strip().lower() == 'datum'
379        assert L[1].strip().lower() == 'wgs84'
380
381        L = lines[3].strip().split()
382        assert L[0].strip().lower() == 'zunits'
383        assert L[1].strip().lower() == 'no'
384
385        L = lines[4].strip().split()
386        assert L[0].strip().lower() == 'units'
387        assert L[1].strip().lower() == 'meters'
388
389        L = lines[5].strip().split()
390        assert L[0].strip().lower() == 'spheroid'
391        assert L[1].strip().lower() == 'wgs84'
392
393        L = lines[6].strip().split()
394        assert L[0].strip().lower() == 'xshift'
395        assert L[1].strip().lower() == '500000'
396
397        L = lines[7].strip().split()
398        assert L[0].strip().lower() == 'yshift'
399        assert L[1].strip().lower() == '10000000'
400
401        L = lines[8].strip().split()
402        assert L[0].strip().lower() == 'parameters'
403
404
405        #Check asc file
406        ascid = open(ascfile)
407        lines = ascid.readlines()
408        ascid.close()
409
410        L = lines[0].strip().split()
411        assert L[0].strip().lower() == 'ncols'
412        assert L[1].strip().lower() == '11'
413
414        L = lines[1].strip().split()
415        assert L[0].strip().lower() == 'nrows'
416        assert L[1].strip().lower() == '11'
417
418        L = lines[2].strip().split()
419        assert L[0].strip().lower() == 'xllcorner'
420        assert num.allclose(float(L[1].strip().lower()), 308500)
421
422        L = lines[3].strip().split()
423        assert L[0].strip().lower() == 'yllcorner'
424        assert num.allclose(float(L[1].strip().lower()), 6189000)
425
426        L = lines[4].strip().split()
427        assert L[0].strip().lower() == 'cellsize'
428        assert num.allclose(float(L[1].strip().lower()), cellsize)
429
430        L = lines[5].strip().split()
431        assert L[0].strip() == 'NODATA_value'
432        assert L[1].strip().lower() == '-9999'
433
434        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
435        for i, line in enumerate(lines[6:]):
436            for j, value in enumerate( line.split() ):
437                assert num.allclose(float(value), -(10-i+j)*cellsize,
438                                    atol=1.0e-12, rtol=1.0e-12)
439
440                # Note: Equality can be obtained in this case,
441                # but it is better to use allclose.
442                #assert float(value) == -(10-i+j)*cellsize
443
444
445        fid.close()
446
447        #Cleanup
448        os.remove(prjfile)
449        os.remove(ascfile)
450        os.remove(swwfile)
451
452
453
454    def test_sww2dem_larger_zero(self):
455        """Test that sww information can be converted correctly to asc/prj
456        format readable by e.g. Arcview. This example has rows with a
457        large number of zeros
458
459        ncols         2001
460        nrows         2
461        xllcorner     308500
462        yllcorner     6189000
463        cellsize      1.000000
464        NODATA_value  -9999
465        0.0 ....
466        """
467
468        import time, os
469        from Scientific.IO.NetCDF import NetCDFFile
470
471        #Setup
472
473        from mesh_factory import rectangular_cross
474
475        #Create basic mesh (100m x 100m)
476        points, vertices, boundary = rectangular_cross(2000, 1, 2000.0, 1.0)
477
478        #Create shallow water domain
479        domain = Domain(points, vertices, boundary)
480        domain.default_order = 1
481
482        domain.set_name('datatest')
483
484        prjfile = domain.get_name() + '_elevation.prj'
485        ascfile = domain.get_name() + '_elevation.asc'
486        swwfile = domain.get_name() + '.sww'
487
488        domain.set_datadir('.')
489        domain.format = 'sww'
490        domain.smooth = True
491        domain.geo_reference = Geo_reference(56, 308500, 6189000)
492
493        #
494        domain.set_quantity('elevation', 0)
495        domain.set_quantity('stage', 0)
496
497        B = Transmissive_boundary(domain)
498        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
499
500
501        #
502        sww = SWW_file(domain)
503        sww.store_connectivity()
504        sww.store_timestep()
505       
506        domain.tight_slope_limiters = 1
507        domain.evolve_to_end(finaltime = 0.01)
508        sww.store_timestep()
509
510        cellsize = 1.0  #0.1 grid
511
512
513        #Check contents
514        #Get NetCDF
515
516        fid = NetCDFFile(sww.filename, netcdf_mode_r)
517
518        # Get the variables
519        x = fid.variables['x'][:]
520        y = fid.variables['y'][:]
521        z = fid.variables['elevation'][:]
522        time = fid.variables['time'][:]
523        stage = fid.variables['stage'][:]
524
525
526        #Export to ascii/prj files
527        sww2dem(domain.get_name(),
528                quantity = 'elevation',
529                cellsize = cellsize,
530                number_of_decimal_places = 9,
531                verbose = self.verbose,
532                format = 'asc',
533                block_size=2)
534
535
536        #Check prj (meta data)
537        prjid = open(prjfile)
538        lines = prjid.readlines()
539        prjid.close()
540
541
542        L = lines[0].strip().split()
543        assert L[0].strip().lower() == 'projection'
544        assert L[1].strip().lower() == 'utm'
545
546        L = lines[1].strip().split()
547        assert L[0].strip().lower() == 'zone'
548        assert L[1].strip().lower() == '56'
549
550        L = lines[2].strip().split()
551        assert L[0].strip().lower() == 'datum'
552        assert L[1].strip().lower() == 'wgs84'
553
554        L = lines[3].strip().split()
555        assert L[0].strip().lower() == 'zunits'
556        assert L[1].strip().lower() == 'no'
557
558        L = lines[4].strip().split()
559        assert L[0].strip().lower() == 'units'
560        assert L[1].strip().lower() == 'meters'
561
562        L = lines[5].strip().split()
563        assert L[0].strip().lower() == 'spheroid'
564        assert L[1].strip().lower() == 'wgs84'
565
566        L = lines[6].strip().split()
567        assert L[0].strip().lower() == 'xshift'
568        assert L[1].strip().lower() == '500000'
569
570        L = lines[7].strip().split()
571        assert L[0].strip().lower() == 'yshift'
572        assert L[1].strip().lower() == '10000000'
573
574        L = lines[8].strip().split()
575        assert L[0].strip().lower() == 'parameters'
576
577
578        #Check asc file
579        ascid = open(ascfile)
580        lines = ascid.readlines()
581        ascid.close()
582
583
584
585        L = lines[0].strip().split()
586        assert L[0].strip().lower() == 'ncols'
587        assert L[1].strip().lower() == '2001'
588
589        L = lines[1].strip().split()
590        assert L[0].strip().lower() == 'nrows'
591        assert L[1].strip().lower() == '2'
592
593        L = lines[2].strip().split()
594        assert L[0].strip().lower() == 'xllcorner'
595        assert num.allclose(float(L[1].strip().lower()), 308500)
596
597        L = lines[3].strip().split()
598        assert L[0].strip().lower() == 'yllcorner'
599        assert num.allclose(float(L[1].strip().lower()), 6189000)
600
601        L = lines[4].strip().split()
602        assert L[0].strip().lower() == 'cellsize'
603        assert num.allclose(float(L[1].strip().lower()), cellsize)
604
605        L = lines[5].strip().split()
606        assert L[0].strip() == 'NODATA_value'
607        assert L[1].strip().lower() == '-9999'
608
609        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
610        for i, line in enumerate(lines[6:]):
611            for j, value in enumerate( line.split() ):
612                assert num.allclose(float(value), 0.0,
613                                    atol=1.0e-12, rtol=1.0e-12)
614
615                # Note: Equality can be obtained in this case,
616                # but it is better to use allclose.
617                #assert float(value) == -(10-i+j)*cellsize
618
619
620        fid.close()
621
622
623        #Cleanup
624        os.remove(prjfile)
625        os.remove(ascfile)
626        os.remove(swwfile)
627
628
629
630
631    def test_sww2dem_boundingbox(self):
632        """Test that sww information can be converted correctly to asc/prj
633        format readable by e.g. ArcView.
634        This will test that mesh can be restricted by bounding box
635
636        Original extent is 100m x 100m:
637
638        Eastings:   308500 -  308600
639        Northings: 6189000 - 6189100
640
641        Bounding box changes this to the 50m x 50m square defined by
642
643        Eastings:   308530 -  308570
644        Northings: 6189050 - 6189100
645
646        The cropped values should be
647
648         -130 -140 -150 -160 -170
649         -120 -130 -140 -150 -160
650         -110 -120 -130 -140 -150
651         -100 -110 -120 -130 -140
652          -90 -100 -110 -120 -130
653          -80  -90 -100 -110 -120
654
655        and the new lower reference point should be
656        Eastings:   308530
657        Northings: 6189050
658
659        Original dataset is the same as in test_sww2dem_larger()
660
661        """
662
663        import time, os
664        from Scientific.IO.NetCDF import NetCDFFile
665
666        #Setup
667
668        from mesh_factory import rectangular
669
670        #Create basic mesh (100m x 100m)
671        points, vertices, boundary = rectangular(2, 2, 100, 100)
672
673        #Create shallow water domain
674        domain = Domain(points, vertices, boundary)
675        domain.default_order = 2
676
677        domain.set_name('datatest')
678
679        prjfile = domain.get_name() + '_elevation.prj'
680        ascfile = domain.get_name() + '_elevation.asc'
681        swwfile = domain.get_name() + '.sww'
682
683        domain.set_datadir('.')
684        domain.format = 'sww'
685        domain.smooth = True
686        domain.geo_reference = Geo_reference(56, 308500, 6189000)
687
688        #
689        domain.set_quantity('elevation', lambda x,y: -x-y)
690        domain.set_quantity('stage', 0)
691
692        B = Transmissive_boundary(domain)
693        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
694
695
696        #
697        sww = SWW_file(domain)
698        sww.store_connectivity()
699        sww.store_timestep()
700
701        #domain.tight_slope_limiters = 1
702        domain.evolve_to_end(finaltime = 0.01)
703        sww.store_timestep()
704
705        cellsize = 10  #10m grid
706
707
708        #Check contents
709        #Get NetCDF
710
711        fid = NetCDFFile(sww.filename, netcdf_mode_r)
712
713        # Get the variables
714        x = fid.variables['x'][:]
715        y = fid.variables['y'][:]
716        z = fid.variables['elevation'][:]
717        time = fid.variables['time'][:]
718        stage = fid.variables['stage'][:]
719
720
721        # Export to ascii/prj files
722        sww2dem(domain.get_name(),
723                quantity = 'elevation',
724                cellsize = cellsize,
725                number_of_decimal_places = 9,
726                easting_min = 308530,
727                easting_max = 308570,
728                northing_min = 6189050,
729                northing_max = 6189100,
730                verbose = self.verbose,
731                format = 'asc')
732
733        fid.close()
734
735
736        # Check prj (meta data)
737        prjid = open(prjfile)
738        lines = prjid.readlines()
739        prjid.close()
740
741        L = lines[0].strip().split()
742        assert L[0].strip().lower() == 'projection'
743        assert L[1].strip().lower() == 'utm'
744
745        L = lines[1].strip().split()
746        assert L[0].strip().lower() == 'zone'
747        assert L[1].strip().lower() == '56'
748
749        L = lines[2].strip().split()
750        assert L[0].strip().lower() == 'datum'
751        assert L[1].strip().lower() == 'wgs84'
752
753        L = lines[3].strip().split()
754        assert L[0].strip().lower() == 'zunits'
755        assert L[1].strip().lower() == 'no'
756
757        L = lines[4].strip().split()
758        assert L[0].strip().lower() == 'units'
759        assert L[1].strip().lower() == 'meters'
760
761        L = lines[5].strip().split()
762        assert L[0].strip().lower() == 'spheroid'
763        assert L[1].strip().lower() == 'wgs84'
764
765        L = lines[6].strip().split()
766        assert L[0].strip().lower() == 'xshift'
767        assert L[1].strip().lower() == '500000'
768
769        L = lines[7].strip().split()
770        assert L[0].strip().lower() == 'yshift'
771        assert L[1].strip().lower() == '10000000'
772
773        L = lines[8].strip().split()
774        assert L[0].strip().lower() == 'parameters'
775
776
777        #Check asc file
778        ascid = open(ascfile)
779        lines = ascid.readlines()
780        ascid.close()
781
782        L = lines[0].strip().split()
783        assert L[0].strip().lower() == 'ncols'
784        assert L[1].strip().lower() == '5'
785
786        L = lines[1].strip().split()
787        assert L[0].strip().lower() == 'nrows'
788        assert L[1].strip().lower() == '6'
789
790        L = lines[2].strip().split()
791        assert L[0].strip().lower() == 'xllcorner'
792        assert num.allclose(float(L[1].strip().lower()), 308530)
793
794        L = lines[3].strip().split()
795        assert L[0].strip().lower() == 'yllcorner'
796        assert num.allclose(float(L[1].strip().lower()), 6189050)
797
798        L = lines[4].strip().split()
799        assert L[0].strip().lower() == 'cellsize'
800        assert num.allclose(float(L[1].strip().lower()), cellsize)
801
802        L = lines[5].strip().split()
803        assert L[0].strip() == 'NODATA_value'
804        assert L[1].strip().lower() == '-9999'
805
806        #Check grid values
807        for i, line in enumerate(lines[6:]):
808            for j, value in enumerate( line.split() ):
809                #assert float(value) == -(10-i+j)*cellsize
810                assert float(value) == -(10-i+j+3)*cellsize
811
812
813
814        #Cleanup
815        os.remove(prjfile)
816        os.remove(ascfile)
817        os.remove(swwfile)
818
819
820
821    def test_sww2dem_asc_stage_reduction(self):
822        """Test that sww information can be converted correctly to asc/prj
823        format readable by e.g. ArcView
824
825        This tests the reduction of quantity stage using min
826        """
827
828        import time, os
829        from Scientific.IO.NetCDF import NetCDFFile
830
831        #Setup
832        self.domain.set_name('datatest')
833
834        prjfile = self.domain.get_name() + '_stage.prj'
835        ascfile = self.domain.get_name() + '_stage.asc'
836        swwfile = self.domain.get_name() + '.sww'
837
838        self.domain.set_datadir('.')
839        self.domain.format = 'sww'
840        self.domain.smooth = True
841        self.domain.set_quantity('elevation', lambda x,y: -x-y)
842
843        self.domain.geo_reference = Geo_reference(56,308500,6189000)
844
845
846        sww = SWW_file(self.domain)
847        sww.store_connectivity()
848        sww.store_timestep()
849
850        #self.domain.tight_slope_limiters = 1
851        self.domain.evolve_to_end(finaltime = 0.01)
852        sww.store_timestep()
853
854        cellsize = 0.25
855        #Check contents
856        #Get NetCDF
857
858        fid = NetCDFFile(sww.filename, netcdf_mode_r)
859
860        # Get the variables
861        x = fid.variables['x'][:]
862        y = fid.variables['y'][:]
863        z = fid.variables['elevation'][:]
864        time = fid.variables['time'][:]
865        stage = fid.variables['stage'][:]
866
867
868        #Export to ascii/prj files
869        sww2dem(self.domain.get_name(),
870                quantity = 'stage',
871                cellsize = cellsize,
872                number_of_decimal_places = 9,
873                reduction = min,
874                format = 'asc',
875                verbose=self.verbose)
876
877
878        #Check asc file
879        ascid = open(ascfile)
880        lines = ascid.readlines()
881        ascid.close()
882
883        L = lines[0].strip().split()
884        assert L[0].strip().lower() == 'ncols'
885        assert L[1].strip().lower() == '5'
886
887        L = lines[1].strip().split()
888        assert L[0].strip().lower() == 'nrows'
889        assert L[1].strip().lower() == '5'
890
891        L = lines[2].strip().split()
892        assert L[0].strip().lower() == 'xllcorner'
893        assert num.allclose(float(L[1].strip().lower()), 308500)
894
895        L = lines[3].strip().split()
896        assert L[0].strip().lower() == 'yllcorner'
897        assert num.allclose(float(L[1].strip().lower()), 6189000)
898
899        L = lines[4].strip().split()
900        assert L[0].strip().lower() == 'cellsize'
901        assert num.allclose(float(L[1].strip().lower()), cellsize)
902
903        L = lines[5].strip().split()
904        assert L[0].strip() == 'NODATA_value'
905        assert L[1].strip().lower() == '-9999'
906
907
908        #Check grid values (where applicable)
909        for j in range(5):
910            if j%2 == 0:
911                L = lines[6+j].strip().split()
912                jj = 4-j
913                for i in range(5):
914                    if i%2 == 0:
915                        index = jj/2 + i/2*3
916                        val0 = stage[0,index]
917                        val1 = stage[1,index]
918
919                        #print i, j, index, ':', L[i], val0, val1
920                        assert num.allclose(float(L[i]), min(val0, val1))
921
922
923        fid.close()
924
925        #Cleanup
926        os.remove(prjfile)
927        os.remove(ascfile)
928        os.remove(swwfile)
929
930    def test_sww2dem_asc_stage_time(self):
931        """Test that sww information can be converted correctly to asc/prj
932        format readable by e.g. ArcView
933
934        This tests the reduction of quantity stage using min
935        """
936
937        import time, os
938        from Scientific.IO.NetCDF import NetCDFFile
939
940        #Setup
941        self.domain.set_name('datatest')
942
943        prjfile = self.domain.get_name() + '_stage.prj'
944        ascfile = self.domain.get_name() + '_stage.asc'
945        swwfile = self.domain.get_name() + '.sww'
946
947        self.domain.set_datadir('.')
948        self.domain.format = 'sww'
949        self.domain.smooth = True
950        self.domain.set_quantity('elevation', lambda x,y: -x-y)
951
952        self.domain.geo_reference = Geo_reference(56,308500,6189000)
953
954        sww = SWW_file(self.domain)
955        sww.store_connectivity()
956        sww.store_timestep()
957
958        #self.domain.tight_slope_limiters = 1
959        self.domain.evolve_to_end(finaltime = 0.01)
960        sww.store_timestep()
961
962        cellsize = 0.25
963        #Check contents
964        #Get NetCDF
965
966        fid = NetCDFFile(sww.filename, netcdf_mode_r)
967
968        # Get the variables
969        x = fid.variables['x'][:]
970        y = fid.variables['y'][:]
971        z = fid.variables['elevation'][:]
972        time = fid.variables['time'][:]
973        stage = fid.variables['stage'][:]
974
975        #Export to ascii/prj files
976        sww2dem(self.domain.get_name(),
977                quantity = 'stage',
978                cellsize = cellsize,
979                number_of_decimal_places = 9,
980                reduction = 1,
981                format = 'asc',
982                verbose=self.verbose)
983
984
985        #Check asc file
986        ascid = open(ascfile)
987        lines = ascid.readlines()
988        ascid.close()
989
990        L = lines[0].strip().split()
991        assert L[0].strip().lower() == 'ncols'
992        assert L[1].strip().lower() == '5'
993
994        L = lines[1].strip().split()
995        assert L[0].strip().lower() == 'nrows'
996        assert L[1].strip().lower() == '5'
997
998        L = lines[2].strip().split()
999        assert L[0].strip().lower() == 'xllcorner'
1000        assert num.allclose(float(L[1].strip().lower()), 308500)
1001
1002        L = lines[3].strip().split()
1003        assert L[0].strip().lower() == 'yllcorner'
1004        assert num.allclose(float(L[1].strip().lower()), 6189000)
1005
1006        L = lines[4].strip().split()
1007        assert L[0].strip().lower() == 'cellsize'
1008        assert num.allclose(float(L[1].strip().lower()), cellsize)
1009
1010        L = lines[5].strip().split()
1011        assert L[0].strip() == 'NODATA_value'
1012        assert L[1].strip().lower() == '-9999'
1013
1014        #Check grid values (where applicable)
1015        for j in range(5):
1016            if j%2 == 0:
1017                L = lines[6+j].strip().split()
1018                jj = 4-j
1019                for i in range(5):
1020                    if i%2 == 0:
1021                        index = jj/2 + i/2*3
1022                       
1023                        val = stage[1,index]
1024                   
1025                        assert num.allclose(float(L[i]), val)
1026
1027        fid.close()
1028
1029        #Cleanup
1030        os.remove(prjfile)
1031        os.remove(ascfile)
1032        os.remove(swwfile)
1033
1034
1035    def test_sww2dem_asc_derived_quantity(self):
1036        """Test that sww information can be converted correctly to asc/prj
1037        format readable by e.g. ArcView
1038
1039        This tests the use of derived quantities
1040        """
1041
1042        import time, os
1043        from Scientific.IO.NetCDF import NetCDFFile
1044
1045        #Setup
1046        self.domain.set_name('datatest')
1047
1048        prjfile = self.domain.get_name() + '_depth.prj'
1049        ascfile = self.domain.get_name() + '_depth.asc'
1050        swwfile = self.domain.get_name() + '.sww'
1051
1052        self.domain.set_datadir('.')
1053        self.domain.format = 'sww'
1054        self.domain.smooth = True
1055        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1056        self.domain.set_quantity('stage', 0.0)
1057
1058        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1059
1060
1061        sww = SWW_file(self.domain)
1062        sww.store_connectivity()
1063        sww.store_timestep()
1064
1065        #self.domain.tight_slope_limiters = 1
1066        self.domain.evolve_to_end(finaltime = 0.01)
1067        sww.store_timestep()
1068
1069        cellsize = 0.25
1070        #Check contents
1071        #Get NetCDF
1072
1073        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1074
1075        # Get the variables
1076        x = fid.variables['x'][:]
1077        y = fid.variables['y'][:]
1078        z = fid.variables['elevation'][:]
1079        time = fid.variables['time'][:]
1080        stage = fid.variables['stage'][:]
1081
1082
1083        #Export to ascii/prj files
1084        sww2dem(self.domain.get_name(),
1085                basename_out = 'datatest_depth',
1086                quantity = 'stage - elevation',
1087                cellsize = cellsize,
1088                number_of_decimal_places = 9,
1089                reduction = min,
1090                format = 'asc',
1091                verbose = self.verbose)
1092
1093
1094        #Check asc file
1095        ascid = open(ascfile)
1096        lines = ascid.readlines()
1097        ascid.close()
1098
1099        L = lines[0].strip().split()
1100        assert L[0].strip().lower() == 'ncols'
1101        assert L[1].strip().lower() == '5'
1102
1103        L = lines[1].strip().split()
1104        assert L[0].strip().lower() == 'nrows'
1105        assert L[1].strip().lower() == '5'
1106
1107        L = lines[2].strip().split()
1108        assert L[0].strip().lower() == 'xllcorner'
1109        assert num.allclose(float(L[1].strip().lower()), 308500)
1110
1111        L = lines[3].strip().split()
1112        assert L[0].strip().lower() == 'yllcorner'
1113        assert num.allclose(float(L[1].strip().lower()), 6189000)
1114
1115        L = lines[4].strip().split()
1116        assert L[0].strip().lower() == 'cellsize'
1117        assert num.allclose(float(L[1].strip().lower()), cellsize)
1118
1119        L = lines[5].strip().split()
1120        assert L[0].strip() == 'NODATA_value'
1121        assert L[1].strip().lower() == '-9999'
1122
1123
1124        #Check grid values (where applicable)
1125        for j in range(5):
1126            if j%2 == 0:
1127                L = lines[6+j].strip().split()
1128                jj = 4-j
1129                for i in range(5):
1130                    if i%2 == 0:
1131                        index = jj/2 + i/2*3
1132                        val0 = stage[0,index] - z[index]
1133                        val1 = stage[1,index] - z[index]
1134
1135                        #print i, j, index, ':', L[i], val0, val1
1136                        assert num.allclose(float(L[i]), min(val0, val1))
1137
1138
1139        fid.close()
1140
1141        #Cleanup
1142        os.remove(prjfile)
1143        os.remove(ascfile)
1144        os.remove(swwfile)
1145
1146
1147
1148
1149
1150    def test_sww2dem_asc_missing_points(self):
1151        """Test that sww information can be converted correctly to asc/prj
1152        format readable by e.g. ArcView
1153
1154        This test includes the writing of missing values
1155        """
1156
1157        import time, os
1158        from Scientific.IO.NetCDF import NetCDFFile
1159
1160        #Setup mesh not coinciding with rectangle.
1161        #This will cause missing values to occur in gridded data
1162
1163
1164        points = [                        [1.0, 1.0],
1165                              [0.5, 0.5], [1.0, 0.5],
1166                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1167
1168        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1169
1170        #Create shallow water domain
1171        domain = Domain(points, vertices)
1172        domain.default_order=2
1173
1174
1175        #Set some field values
1176        domain.set_quantity('elevation', lambda x,y: -x-y)
1177        domain.set_quantity('friction', 0.03)
1178
1179
1180        ######################
1181        # Boundary conditions
1182        B = Transmissive_boundary(domain)
1183        domain.set_boundary( {'exterior': B} )
1184
1185
1186        ######################
1187        #Initial condition - with jumps
1188
1189        bed = domain.quantities['elevation'].vertex_values
1190        stage = num.zeros(bed.shape, num.float)
1191
1192        h = 0.3
1193        for i in range(stage.shape[0]):
1194            if i % 2 == 0:
1195                stage[i,:] = bed[i,:] + h
1196            else:
1197                stage[i,:] = bed[i,:]
1198
1199        domain.set_quantity('stage', stage)
1200        domain.distribute_to_vertices_and_edges()
1201
1202        domain.set_name('datatest')
1203
1204        prjfile = domain.get_name() + '_elevation.prj'
1205        ascfile = domain.get_name() + '_elevation.asc'
1206        swwfile = domain.get_name() + '.sww'
1207
1208        domain.set_datadir('.')
1209        domain.format = 'sww'
1210        domain.smooth = True
1211
1212        domain.geo_reference = Geo_reference(56,308500,6189000)
1213
1214        sww = SWW_file(domain)
1215        sww.store_connectivity()
1216        sww.store_timestep()
1217
1218        cellsize = 0.25
1219        #Check contents
1220        #Get NetCDF
1221
1222        fid = NetCDFFile(swwfile, netcdf_mode_r)
1223
1224        # Get the variables
1225        x = fid.variables['x'][:]
1226        y = fid.variables['y'][:]
1227        z = fid.variables['elevation'][:]
1228        time = fid.variables['time'][:]
1229
1230        try:
1231            geo_reference = Geo_reference(NetCDFObject=fid)
1232        except AttributeError, e:
1233            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1234
1235        #Export to ascii/prj files
1236        sww2dem(domain.get_name(),
1237                quantity = 'elevation',
1238                cellsize = cellsize,
1239                number_of_decimal_places = 9,
1240                verbose = self.verbose,
1241                format = 'asc')
1242
1243
1244        #Check asc file
1245        ascid = open(ascfile)
1246        lines = ascid.readlines()
1247        ascid.close()
1248
1249        L = lines[0].strip().split()
1250        assert L[0].strip().lower() == 'ncols'
1251        assert L[1].strip().lower() == '5'
1252
1253        L = lines[1].strip().split()
1254        assert L[0].strip().lower() == 'nrows'
1255        assert L[1].strip().lower() == '5'
1256
1257        L = lines[2].strip().split()
1258        assert L[0].strip().lower() == 'xllcorner'
1259        assert num.allclose(float(L[1].strip().lower()), 308500)
1260
1261        L = lines[3].strip().split()
1262        assert L[0].strip().lower() == 'yllcorner'
1263        assert num.allclose(float(L[1].strip().lower()), 6189000)
1264
1265        L = lines[4].strip().split()
1266        assert L[0].strip().lower() == 'cellsize'
1267        assert num.allclose(float(L[1].strip().lower()), cellsize)
1268
1269        L = lines[5].strip().split()
1270        assert L[0].strip() == 'NODATA_value'
1271        assert L[1].strip().lower() == '-9999'
1272
1273        #Check grid values
1274        for j in range(5):
1275            L = lines[6+j].strip().split()
1276            assert len(L) == 5
1277            y = (4-j) * cellsize
1278
1279            for i in range(5):
1280                #print i
1281                if i+j >= 4:
1282                    assert num.allclose(float(L[i]), -i*cellsize - y)
1283                else:
1284                    #Missing values
1285                    assert num.allclose(float(L[i]), -9999)
1286
1287
1288
1289        fid.close()
1290
1291        #Cleanup
1292        os.remove(prjfile)
1293        os.remove(ascfile)
1294        os.remove(swwfile)
1295
1296
1297
1298    def test_sww2ers_simple(self):
1299        """Test that sww information can be converted correctly to asc/prj
1300        format readable by e.g. ArcView
1301        """
1302
1303        import time, os
1304        from Scientific.IO.NetCDF import NetCDFFile
1305
1306
1307        NODATA_value = 1758323
1308
1309        #Setup
1310        self.domain.set_name('datatest')
1311
1312        headerfile = self.domain.get_name() + '.ers'
1313        swwfile = self.domain.get_name() + '.sww'
1314
1315        self.domain.set_datadir('.')
1316        self.domain.format = 'sww'
1317        self.domain.smooth = True
1318        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1319
1320        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1321
1322        sww = SWW_file(self.domain)
1323        sww.store_connectivity()
1324        sww.store_timestep()
1325
1326        #self.domain.tight_slope_limiters = 1
1327        self.domain.evolve_to_end(finaltime = 0.01)
1328        sww.store_timestep()
1329
1330        cellsize = 0.25
1331        #Check contents
1332        #Get NetCDF
1333
1334        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1335
1336        # Get the variables
1337        x = fid.variables['x'][:]
1338        y = fid.variables['y'][:]
1339        z = fid.variables['elevation'][:]
1340        time = fid.variables['time'][:]
1341        stage = fid.variables['stage'][:]
1342
1343
1344        #Export to ers files
1345        sww2dem(self.domain.get_name(),
1346                quantity = 'elevation',
1347                cellsize = cellsize,
1348                number_of_decimal_places = 9,
1349                NODATA_value = NODATA_value,
1350                verbose = self.verbose,
1351                format = 'ers')
1352
1353        #Check header data
1354        from ermapper_grids import read_ermapper_header, read_ermapper_data
1355
1356        header = read_ermapper_header(self.domain.get_name() + '_elevation.ers')
1357        #print header
1358        assert header['projection'].lower() == '"utm-56"'
1359        assert header['datum'].lower() == '"wgs84"'
1360        assert header['units'].lower() == '"meters"'
1361        assert header['value'].lower() == '"elevation"'
1362        assert header['xdimension'] == '0.25'
1363        assert header['ydimension'] == '0.25'
1364        assert float(header['eastings']) == 308500.0   #xllcorner
1365        assert float(header['northings']) == 6189000.0 #yllcorner
1366        assert int(header['nroflines']) == 5
1367        assert int(header['nrofcellsperline']) == 5
1368        assert int(header['nullcellvalue']) == NODATA_value
1369        #FIXME - there is more in the header
1370
1371
1372        #Check grid data
1373        grid = read_ermapper_data(self.domain.get_name() + '_elevation')
1374
1375        #FIXME (Ole): Why is this the desired reference grid for -x-y?
1376        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
1377                    -1,    -1.25, -1.5,  -1.75, -2.0,
1378                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
1379                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1380                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
1381
1382
1383        #print grid
1384        assert num.allclose(grid, ref_grid)
1385
1386        fid.close()
1387
1388        #Cleanup
1389        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1390        #Done (Ole) - it was because sww2ers didn't close it's sww file
1391        os.remove(sww.filename)
1392        os.remove(self.domain.get_name() + '_elevation')
1393        os.remove(self.domain.get_name() + '_elevation.ers')
1394
1395#################################################################################
1396
1397if __name__ == "__main__":
1398    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
1399    suite = unittest.makeSuite(Test_Sww2Dem, 'test')
1400    runner = unittest.TextTestRunner(verbosity=1)
1401    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.