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

Last change on this file since 7866 was 7866, checked in by hudson, 13 years ago

More swb tests passing. Cleaned up some pylint errors.
File size: 56.3 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, sww2dem_batch
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()+'.sww',
135                self.domain.get_name()+'_elevation.asc',
136                quantity = 'elevation',
137                cellsize = cellsize,
138                number_of_decimal_places = 9,
139                verbose = self.verbose)
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        ascfile = self.domain.get_name() + '_depth.asc'
223        prjfile = self.domain.get_name() + '_depth.prj'
224
225        #Export to ascii/prj files
226        sww2dem(self.domain.get_name()+'.sww',
227                ascfile,
228                quantity = 'depth',
229                cellsize = cellsize,
230                number_of_decimal_places = 9,
231                verbose = self.verbose)
232       
233        #Check asc file
234        ascid = open(ascfile)
235        lines = ascid.readlines()
236        ascid.close()
237
238        L = lines[0].strip().split()
239        assert L[0].strip().lower() == 'ncols'
240        assert L[1].strip().lower() == '5'
241
242        L = lines[1].strip().split()
243        assert L[0].strip().lower() == 'nrows'
244        assert L[1].strip().lower() == '5'
245
246        L = lines[2].strip().split()
247        assert L[0].strip().lower() == 'xllcorner'
248        assert num.allclose(float(L[1].strip().lower()), 308500)
249
250        L = lines[3].strip().split()
251        assert L[0].strip().lower() == 'yllcorner'
252        assert num.allclose(float(L[1].strip().lower()), 6189000)
253
254        L = lines[4].strip().split()
255        assert L[0].strip().lower() == 'cellsize'
256        assert num.allclose(float(L[1].strip().lower()), cellsize)
257
258        L = lines[5].strip().split()
259        assert L[0].strip() == 'NODATA_value'
260        assert L[1].strip().lower() == '-9999'
261
262        #Check grid values
263        for j in range(5):
264            L = lines[6+j].strip().split()
265            y = (4-j) * cellsize
266            for i in range(5):
267                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
268
269
270        #Cleanup
271        os.remove(prjfile)
272        os.remove(ascfile)
273        os.remove(swwfile)
274
275
276
277    def test_sww2dem_larger(self):
278        """Test that sww information can be converted correctly to asc/prj
279        format readable by e.g. ArcView. Here:
280
281        ncols         11
282        nrows         11
283        xllcorner     308500
284        yllcorner     6189000
285        cellsize      10.000000
286        NODATA_value  -9999
287        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
288         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
289         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
290         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
291         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
292         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
293         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
294         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
295         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
296         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
297           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
298
299        """
300
301        import time, os
302        from Scientific.IO.NetCDF import NetCDFFile
303
304        #Create basic mesh (100m x 100m)
305        points, vertices, boundary = rectangular(2, 2, 100, 100)
306
307        #Create shallow water domain
308        domain = Domain(points, vertices, boundary)
309        domain.default_order = 2
310
311        domain.set_name('datatest')
312
313        prjfile = domain.get_name() + '_elevation.prj'
314        ascfile = domain.get_name() + '_elevation.asc'
315        swwfile = domain.get_name() + '.sww'
316
317        domain.set_datadir('.')
318        domain.format = 'sww'
319        domain.smooth = True
320        domain.geo_reference = Geo_reference(56, 308500, 6189000)
321
322        #
323        domain.set_quantity('elevation', lambda x,y: -x-y)
324        domain.set_quantity('stage', 0)
325
326        B = Transmissive_boundary(domain)
327        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
328
329
330        #
331        sww = SWW_file(domain)
332        sww.store_connectivity()
333        sww.store_timestep()
334       
335        domain.tight_slope_limiters = 1
336        domain.evolve_to_end(finaltime = 0.01)
337        sww.store_timestep()
338
339        cellsize = 10  #10m grid
340
341
342        #Check contents
343        #Get NetCDF
344
345        fid = NetCDFFile(sww.filename, netcdf_mode_r)
346
347        # Get the variables
348        x = fid.variables['x'][:]
349        y = fid.variables['y'][:]
350        z = fid.variables['elevation'][:]
351        time = fid.variables['time'][:]
352        stage = fid.variables['stage'][:]
353
354
355        #Export to ascii/prj files
356        sww2dem(domain.get_name() + '.sww',
357                domain.get_name() + '_elevation.asc',
358                quantity = 'elevation',
359                cellsize = cellsize,
360                number_of_decimal_places = 9,
361                verbose = self.verbose,
362                block_size=2)
363
364
365        #Check prj (meta data)
366        prjid = open(prjfile)
367        lines = prjid.readlines()
368        prjid.close()
369
370        L = lines[0].strip().split()
371        assert L[0].strip().lower() == 'projection'
372        assert L[1].strip().lower() == 'utm'
373
374        L = lines[1].strip().split()
375        assert L[0].strip().lower() == 'zone'
376        assert L[1].strip().lower() == '56'
377
378        L = lines[2].strip().split()
379        assert L[0].strip().lower() == 'datum'
380        assert L[1].strip().lower() == 'wgs84'
381
382        L = lines[3].strip().split()
383        assert L[0].strip().lower() == 'zunits'
384        assert L[1].strip().lower() == 'no'
385
386        L = lines[4].strip().split()
387        assert L[0].strip().lower() == 'units'
388        assert L[1].strip().lower() == 'meters'
389
390        L = lines[5].strip().split()
391        assert L[0].strip().lower() == 'spheroid'
392        assert L[1].strip().lower() == 'wgs84'
393
394        L = lines[6].strip().split()
395        assert L[0].strip().lower() == 'xshift'
396        assert L[1].strip().lower() == '500000'
397
398        L = lines[7].strip().split()
399        assert L[0].strip().lower() == 'yshift'
400        assert L[1].strip().lower() == '10000000'
401
402        L = lines[8].strip().split()
403        assert L[0].strip().lower() == 'parameters'
404
405
406        #Check asc file
407        ascid = open(ascfile)
408        lines = ascid.readlines()
409        ascid.close()
410
411        L = lines[0].strip().split()
412        assert L[0].strip().lower() == 'ncols'
413        assert L[1].strip().lower() == '11'
414
415        L = lines[1].strip().split()
416        assert L[0].strip().lower() == 'nrows'
417        assert L[1].strip().lower() == '11'
418
419        L = lines[2].strip().split()
420        assert L[0].strip().lower() == 'xllcorner'
421        assert num.allclose(float(L[1].strip().lower()), 308500)
422
423        L = lines[3].strip().split()
424        assert L[0].strip().lower() == 'yllcorner'
425        assert num.allclose(float(L[1].strip().lower()), 6189000)
426
427        L = lines[4].strip().split()
428        assert L[0].strip().lower() == 'cellsize'
429        assert num.allclose(float(L[1].strip().lower()), cellsize)
430
431        L = lines[5].strip().split()
432        assert L[0].strip() == 'NODATA_value'
433        assert L[1].strip().lower() == '-9999'
434
435        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
436        for i, line in enumerate(lines[6:]):
437            for j, value in enumerate( line.split() ):
438                assert num.allclose(float(value), -(10-i+j)*cellsize,
439                                    atol=1.0e-12, rtol=1.0e-12)
440
441                # Note: Equality can be obtained in this case,
442                # but it is better to use allclose.
443                #assert float(value) == -(10-i+j)*cellsize
444
445
446        fid.close()
447
448        #Cleanup
449        os.remove(prjfile)
450        os.remove(ascfile)
451        os.remove(swwfile)
452
453
454
455    def test_sww2dem_larger_zero(self):
456        """Test that sww information can be converted correctly to asc/prj
457        format readable by e.g. Arcview. This example has rows with a
458        large number of zeros
459
460        ncols         2001
461        nrows         2
462        xllcorner     308500
463        yllcorner     6189000
464        cellsize      1.000000
465        NODATA_value  -9999
466        0.0 ....
467        """
468
469        import time, os
470        from Scientific.IO.NetCDF import NetCDFFile
471
472        #Setup
473
474        from mesh_factory import rectangular_cross
475
476        #Create basic mesh (100m x 100m)
477        points, vertices, boundary = rectangular_cross(2000, 1, 2000.0, 1.0)
478
479        #Create shallow water domain
480        domain = Domain(points, vertices, boundary)
481        domain.default_order = 1
482
483        domain.set_name('datatest')
484
485        prjfile = domain.get_name() + '_elevation.prj'
486        ascfile = domain.get_name() + '_elevation.asc'
487        swwfile = domain.get_name() + '.sww'
488
489        domain.set_datadir('.')
490        domain.format = 'sww'
491        domain.smooth = True
492        domain.geo_reference = Geo_reference(56, 308500, 6189000)
493
494        #
495        domain.set_quantity('elevation', 0)
496        domain.set_quantity('stage', 0)
497
498        B = Transmissive_boundary(domain)
499        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
500
501
502        #
503        sww = SWW_file(domain)
504        sww.store_connectivity()
505        sww.store_timestep()
506       
507        domain.tight_slope_limiters = 1
508        domain.evolve_to_end(finaltime = 0.01)
509        sww.store_timestep()
510
511        cellsize = 1.0  #0.1 grid
512
513
514        #Check contents
515        #Get NetCDF
516
517        fid = NetCDFFile(sww.filename, netcdf_mode_r)
518
519        # Get the variables
520        x = fid.variables['x'][:]
521        y = fid.variables['y'][:]
522        z = fid.variables['elevation'][:]
523        time = fid.variables['time'][:]
524        stage = fid.variables['stage'][:]
525
526
527        #Export to ascii/prj files
528        sww2dem(domain.get_name()+'.sww',
529                domain.get_name() + '_elevation.asc',
530                quantity = 'elevation',
531                cellsize = cellsize,
532                number_of_decimal_places = 9,
533                verbose = self.verbose,
534                block_size=2)
535
536
537        #Check prj (meta data)
538        prjid = open(prjfile)
539        lines = prjid.readlines()
540        prjid.close()
541
542
543        L = lines[0].strip().split()
544        assert L[0].strip().lower() == 'projection'
545        assert L[1].strip().lower() == 'utm'
546
547        L = lines[1].strip().split()
548        assert L[0].strip().lower() == 'zone'
549        assert L[1].strip().lower() == '56'
550
551        L = lines[2].strip().split()
552        assert L[0].strip().lower() == 'datum'
553        assert L[1].strip().lower() == 'wgs84'
554
555        L = lines[3].strip().split()
556        assert L[0].strip().lower() == 'zunits'
557        assert L[1].strip().lower() == 'no'
558
559        L = lines[4].strip().split()
560        assert L[0].strip().lower() == 'units'
561        assert L[1].strip().lower() == 'meters'
562
563        L = lines[5].strip().split()
564        assert L[0].strip().lower() == 'spheroid'
565        assert L[1].strip().lower() == 'wgs84'
566
567        L = lines[6].strip().split()
568        assert L[0].strip().lower() == 'xshift'
569        assert L[1].strip().lower() == '500000'
570
571        L = lines[7].strip().split()
572        assert L[0].strip().lower() == 'yshift'
573        assert L[1].strip().lower() == '10000000'
574
575        L = lines[8].strip().split()
576        assert L[0].strip().lower() == 'parameters'
577
578
579        #Check asc file
580        ascid = open(ascfile)
581        lines = ascid.readlines()
582        ascid.close()
583
584
585
586        L = lines[0].strip().split()
587        assert L[0].strip().lower() == 'ncols'
588        assert L[1].strip().lower() == '2001'
589
590        L = lines[1].strip().split()
591        assert L[0].strip().lower() == 'nrows'
592        assert L[1].strip().lower() == '2'
593
594        L = lines[2].strip().split()
595        assert L[0].strip().lower() == 'xllcorner'
596        assert num.allclose(float(L[1].strip().lower()), 308500)
597
598        L = lines[3].strip().split()
599        assert L[0].strip().lower() == 'yllcorner'
600        assert num.allclose(float(L[1].strip().lower()), 6189000)
601
602        L = lines[4].strip().split()
603        assert L[0].strip().lower() == 'cellsize'
604        assert num.allclose(float(L[1].strip().lower()), cellsize)
605
606        L = lines[5].strip().split()
607        assert L[0].strip() == 'NODATA_value'
608        assert L[1].strip().lower() == '-9999'
609
610        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
611        for i, line in enumerate(lines[6:]):
612            for j, value in enumerate( line.split() ):
613                assert num.allclose(float(value), 0.0,
614                                    atol=1.0e-12, rtol=1.0e-12)
615
616                # Note: Equality can be obtained in this case,
617                # but it is better to use allclose.
618                #assert float(value) == -(10-i+j)*cellsize
619
620
621        fid.close()
622
623
624        #Cleanup
625        os.remove(prjfile)
626        os.remove(ascfile)
627        os.remove(swwfile)
628
629
630
631
632    def test_sww2dem_boundingbox(self):
633        """Test that sww information can be converted correctly to asc/prj
634        format readable by e.g. ArcView.
635        This will test that mesh can be restricted by bounding box
636
637        Original extent is 100m x 100m:
638
639        Eastings:   308500 -  308600
640        Northings: 6189000 - 6189100
641
642        Bounding box changes this to the 50m x 50m square defined by
643
644        Eastings:   308530 -  308570
645        Northings: 6189050 - 6189100
646
647        The cropped values should be
648
649         -130 -140 -150 -160 -170
650         -120 -130 -140 -150 -160
651         -110 -120 -130 -140 -150
652         -100 -110 -120 -130 -140
653          -90 -100 -110 -120 -130
654          -80  -90 -100 -110 -120
655
656        and the new lower reference point should be
657        Eastings:   308530
658        Northings: 6189050
659
660        Original dataset is the same as in test_sww2dem_larger()
661
662        """
663
664        import time, os
665        from Scientific.IO.NetCDF import NetCDFFile
666
667        #Setup
668
669        from mesh_factory import rectangular
670
671        #Create basic mesh (100m x 100m)
672        points, vertices, boundary = rectangular(2, 2, 100, 100)
673
674        #Create shallow water domain
675        domain = Domain(points, vertices, boundary)
676        domain.default_order = 2
677
678        domain.set_name('datatest')
679
680        prjfile = domain.get_name() + '_elevation.prj'
681        ascfile = domain.get_name() + '_elevation.asc'
682        swwfile = domain.get_name() + '.sww'
683
684        domain.set_datadir('.')
685        domain.format = 'sww'
686        domain.smooth = True
687        domain.geo_reference = Geo_reference(56, 308500, 6189000)
688
689        #
690        domain.set_quantity('elevation', lambda x,y: -x-y)
691        domain.set_quantity('stage', 0)
692
693        B = Transmissive_boundary(domain)
694        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
695
696
697        #
698        sww = SWW_file(domain)
699        sww.store_connectivity()
700        sww.store_timestep()
701
702        #domain.tight_slope_limiters = 1
703        domain.evolve_to_end(finaltime = 0.01)
704        sww.store_timestep()
705
706        cellsize = 10  #10m grid
707
708
709        #Check contents
710        #Get NetCDF
711
712        fid = NetCDFFile(sww.filename, netcdf_mode_r)
713
714        # Get the variables
715        x = fid.variables['x'][:]
716        y = fid.variables['y'][:]
717        z = fid.variables['elevation'][:]
718        time = fid.variables['time'][:]
719        stage = fid.variables['stage'][:]
720
721
722        # Export to ascii/prj files
723        sww2dem(domain.get_name() + '.sww',
724                domain.get_name() + '_elevation.asc',
725                quantity = 'elevation',
726                cellsize = cellsize,
727                number_of_decimal_places = 9,
728                easting_min = 308530,
729                easting_max = 308570,
730                northing_min = 6189050,
731                northing_max = 6189100,
732                verbose = self.verbose)
733
734        fid.close()
735
736
737        # Check prj (meta data)
738        prjid = open(prjfile)
739        lines = prjid.readlines()
740        prjid.close()
741
742        L = lines[0].strip().split()
743        assert L[0].strip().lower() == 'projection'
744        assert L[1].strip().lower() == 'utm'
745
746        L = lines[1].strip().split()
747        assert L[0].strip().lower() == 'zone'
748        assert L[1].strip().lower() == '56'
749
750        L = lines[2].strip().split()
751        assert L[0].strip().lower() == 'datum'
752        assert L[1].strip().lower() == 'wgs84'
753
754        L = lines[3].strip().split()
755        assert L[0].strip().lower() == 'zunits'
756        assert L[1].strip().lower() == 'no'
757
758        L = lines[4].strip().split()
759        assert L[0].strip().lower() == 'units'
760        assert L[1].strip().lower() == 'meters'
761
762        L = lines[5].strip().split()
763        assert L[0].strip().lower() == 'spheroid'
764        assert L[1].strip().lower() == 'wgs84'
765
766        L = lines[6].strip().split()
767        assert L[0].strip().lower() == 'xshift'
768        assert L[1].strip().lower() == '500000'
769
770        L = lines[7].strip().split()
771        assert L[0].strip().lower() == 'yshift'
772        assert L[1].strip().lower() == '10000000'
773
774        L = lines[8].strip().split()
775        assert L[0].strip().lower() == 'parameters'
776
777
778        #Check asc file
779        ascid = open(ascfile)
780        lines = ascid.readlines()
781        ascid.close()
782
783        L = lines[0].strip().split()
784        assert L[0].strip().lower() == 'ncols'
785        assert L[1].strip().lower() == '5'
786
787        L = lines[1].strip().split()
788        assert L[0].strip().lower() == 'nrows'
789        assert L[1].strip().lower() == '6'
790
791        L = lines[2].strip().split()
792        assert L[0].strip().lower() == 'xllcorner'
793        assert num.allclose(float(L[1].strip().lower()), 308530)
794
795        L = lines[3].strip().split()
796        assert L[0].strip().lower() == 'yllcorner'
797        assert num.allclose(float(L[1].strip().lower()), 6189050)
798
799        L = lines[4].strip().split()
800        assert L[0].strip().lower() == 'cellsize'
801        assert num.allclose(float(L[1].strip().lower()), cellsize)
802
803        L = lines[5].strip().split()
804        assert L[0].strip() == 'NODATA_value'
805        assert L[1].strip().lower() == '-9999'
806
807        #Check grid values
808        for i, line in enumerate(lines[6:]):
809            for j, value in enumerate( line.split() ):
810                #assert float(value) == -(10-i+j)*cellsize
811                assert float(value) == -(10-i+j+3)*cellsize
812
813
814
815        #Cleanup
816        os.remove(prjfile)
817        os.remove(ascfile)
818        os.remove(swwfile)
819
820
821
822    def test_sww2dem_asc_stage_reduction(self):
823        """Test that sww information can be converted correctly to asc/prj
824        format readable by e.g. ArcView
825
826        This tests the reduction of quantity stage using min
827        """
828
829        import time, os
830        from Scientific.IO.NetCDF import NetCDFFile
831
832        #Setup
833        self.domain.set_name('datatest')
834
835        prjfile = self.domain.get_name() + '_stage.prj'
836        ascfile = self.domain.get_name() + '_stage.asc'
837        swwfile = self.domain.get_name() + '.sww'
838
839        self.domain.set_datadir('.')
840        self.domain.format = 'sww'
841        self.domain.smooth = True
842        self.domain.set_quantity('elevation', lambda x,y: -x-y)
843
844        self.domain.geo_reference = Geo_reference(56,308500,6189000)
845
846
847        sww = SWW_file(self.domain)
848        sww.store_connectivity()
849        sww.store_timestep()
850
851        #self.domain.tight_slope_limiters = 1
852        self.domain.evolve_to_end(finaltime = 0.01)
853        sww.store_timestep()
854
855        cellsize = 0.25
856        #Check contents
857        #Get NetCDF
858
859        fid = NetCDFFile(sww.filename, netcdf_mode_r)
860
861        # Get the variables
862        x = fid.variables['x'][:]
863        y = fid.variables['y'][:]
864        z = fid.variables['elevation'][:]
865        time = fid.variables['time'][:]
866        stage = fid.variables['stage'][:]
867
868
869        #Export to ascii/prj files
870        sww2dem(self.domain.get_name() + '.sww',
871                self.domain.get_name() + '_stage.asc',
872                quantity = 'stage',
873                cellsize = cellsize,
874                number_of_decimal_places = 9,
875                reduction = min)
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() + '.sww',
977                self.domain.get_name() + '_stage.asc',
978                quantity = 'stage',
979                cellsize = cellsize,
980                number_of_decimal_places = 9,
981                reduction = 1,
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()+'.sww',
1085                name_out = 'datatest_depth.asc',
1086                quantity = 'stage - elevation',
1087                cellsize = cellsize,
1088                number_of_decimal_places = 9,
1089                reduction = min,
1090                verbose = self.verbose)
1091
1092
1093        #Check asc file
1094        ascid = open(ascfile)
1095        lines = ascid.readlines()
1096        ascid.close()
1097
1098        L = lines[0].strip().split()
1099        assert L[0].strip().lower() == 'ncols'
1100        assert L[1].strip().lower() == '5'
1101
1102        L = lines[1].strip().split()
1103        assert L[0].strip().lower() == 'nrows'
1104        assert L[1].strip().lower() == '5'
1105
1106        L = lines[2].strip().split()
1107        assert L[0].strip().lower() == 'xllcorner'
1108        assert num.allclose(float(L[1].strip().lower()), 308500)
1109
1110        L = lines[3].strip().split()
1111        assert L[0].strip().lower() == 'yllcorner'
1112        assert num.allclose(float(L[1].strip().lower()), 6189000)
1113
1114        L = lines[4].strip().split()
1115        assert L[0].strip().lower() == 'cellsize'
1116        assert num.allclose(float(L[1].strip().lower()), cellsize)
1117
1118        L = lines[5].strip().split()
1119        assert L[0].strip() == 'NODATA_value'
1120        assert L[1].strip().lower() == '-9999'
1121
1122
1123        #Check grid values (where applicable)
1124        for j in range(5):
1125            if j%2 == 0:
1126                L = lines[6+j].strip().split()
1127                jj = 4-j
1128                for i in range(5):
1129                    if i%2 == 0:
1130                        index = jj/2 + i/2*3
1131                        val0 = stage[0,index] - z[index]
1132                        val1 = stage[1,index] - z[index]
1133
1134                        #print i, j, index, ':', L[i], val0, val1
1135                        assert num.allclose(float(L[i]), min(val0, val1))
1136
1137
1138        fid.close()
1139
1140        #Cleanup
1141        os.remove(prjfile)
1142        os.remove(ascfile)
1143        os.remove(swwfile)
1144
1145
1146
1147
1148
1149    def test_sww2dem_asc_missing_points(self):
1150        """Test that sww information can be converted correctly to asc/prj
1151        format readable by e.g. ArcView
1152
1153        This test includes the writing of missing values
1154        """
1155
1156        import time, os
1157        from Scientific.IO.NetCDF import NetCDFFile
1158
1159        #Setup mesh not coinciding with rectangle.
1160        #This will cause missing values to occur in gridded data
1161
1162
1163        points = [                        [1.0, 1.0],
1164                              [0.5, 0.5], [1.0, 0.5],
1165                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1166
1167        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1168
1169        #Create shallow water domain
1170        domain = Domain(points, vertices)
1171        domain.default_order=2
1172
1173
1174        #Set some field values
1175        domain.set_quantity('elevation', lambda x,y: -x-y)
1176        domain.set_quantity('friction', 0.03)
1177
1178
1179        ######################
1180        # Boundary conditions
1181        B = Transmissive_boundary(domain)
1182        domain.set_boundary( {'exterior': B} )
1183
1184
1185        ######################
1186        #Initial condition - with jumps
1187
1188        bed = domain.quantities['elevation'].vertex_values
1189        stage = num.zeros(bed.shape, num.float)
1190
1191        h = 0.3
1192        for i in range(stage.shape[0]):
1193            if i % 2 == 0:
1194                stage[i,:] = bed[i,:] + h
1195            else:
1196                stage[i,:] = bed[i,:]
1197
1198        domain.set_quantity('stage', stage)
1199        domain.distribute_to_vertices_and_edges()
1200
1201        domain.set_name('datatest')
1202
1203        prjfile = domain.get_name() + '_elevation.prj'
1204        ascfile = domain.get_name() + '_elevation.asc'
1205        swwfile = domain.get_name() + '.sww'
1206
1207        domain.set_datadir('.')
1208        domain.format = 'sww'
1209        domain.smooth = True
1210
1211        domain.geo_reference = Geo_reference(56,308500,6189000)
1212
1213        sww = SWW_file(domain)
1214        sww.store_connectivity()
1215        sww.store_timestep()
1216
1217        cellsize = 0.25
1218        #Check contents
1219        #Get NetCDF
1220
1221        fid = NetCDFFile(swwfile, netcdf_mode_r)
1222
1223        # Get the variables
1224        x = fid.variables['x'][:]
1225        y = fid.variables['y'][:]
1226        z = fid.variables['elevation'][:]
1227        time = fid.variables['time'][:]
1228
1229        try:
1230            geo_reference = Geo_reference(NetCDFObject=fid)
1231        except AttributeError, e:
1232            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1233
1234        #Export to ascii/prj files
1235        sww2dem(domain.get_name()+'.sww',
1236                domain.get_name()+'_elevation.asc',
1237                quantity = 'elevation',
1238                cellsize = cellsize,
1239                number_of_decimal_places = 9,
1240                verbose = self.verbose)
1241
1242
1243        #Check asc file
1244        ascid = open(ascfile)
1245        lines = ascid.readlines()
1246        ascid.close()
1247
1248        L = lines[0].strip().split()
1249        assert L[0].strip().lower() == 'ncols'
1250        assert L[1].strip().lower() == '5'
1251
1252        L = lines[1].strip().split()
1253        assert L[0].strip().lower() == 'nrows'
1254        assert L[1].strip().lower() == '5'
1255
1256        L = lines[2].strip().split()
1257        assert L[0].strip().lower() == 'xllcorner'
1258        assert num.allclose(float(L[1].strip().lower()), 308500)
1259
1260        L = lines[3].strip().split()
1261        assert L[0].strip().lower() == 'yllcorner'
1262        assert num.allclose(float(L[1].strip().lower()), 6189000)
1263
1264        L = lines[4].strip().split()
1265        assert L[0].strip().lower() == 'cellsize'
1266        assert num.allclose(float(L[1].strip().lower()), cellsize)
1267
1268        L = lines[5].strip().split()
1269        assert L[0].strip() == 'NODATA_value'
1270        assert L[1].strip().lower() == '-9999'
1271
1272        #Check grid values
1273        for j in range(5):
1274            L = lines[6+j].strip().split()
1275            assert len(L) == 5
1276            y = (4-j) * cellsize
1277
1278            for i in range(5):
1279                #print i
1280                if i+j >= 4:
1281                    assert num.allclose(float(L[i]), -i*cellsize - y)
1282                else:
1283                    #Missing values
1284                    assert num.allclose(float(L[i]), -9999)
1285
1286
1287
1288        fid.close()
1289
1290        #Cleanup
1291        os.remove(prjfile)
1292        os.remove(ascfile)
1293        os.remove(swwfile)
1294
1295
1296
1297    def test_sww2ers_simple(self):
1298        """Test that sww information can be converted correctly to asc/prj
1299        format readable by e.g. ArcView
1300        """
1301
1302        import time, os
1303        from Scientific.IO.NetCDF import NetCDFFile
1304
1305
1306        NODATA_value = 1758323
1307
1308        #Setup
1309        self.domain.set_name('datatest')
1310
1311        headerfile = self.domain.get_name() + '.ers'
1312        swwfile = self.domain.get_name() + '.sww'
1313
1314        self.domain.set_datadir('.')
1315        self.domain.format = 'sww'
1316        self.domain.smooth = True
1317        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1318
1319        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1320
1321        sww = SWW_file(self.domain)
1322        sww.store_connectivity()
1323        sww.store_timestep()
1324
1325        #self.domain.tight_slope_limiters = 1
1326        self.domain.evolve_to_end(finaltime = 0.01)
1327        sww.store_timestep()
1328
1329        cellsize = 0.25
1330        #Check contents
1331        #Get NetCDF
1332
1333        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1334
1335        # Get the variables
1336        x = fid.variables['x'][:]
1337        y = fid.variables['y'][:]
1338        z = fid.variables['elevation'][:]
1339        time = fid.variables['time'][:]
1340        stage = fid.variables['stage'][:]
1341
1342
1343        #Export to ers files
1344        outname = self.domain.get_name() + '_elevation.ers'
1345        sww2dem(self.domain.get_name() + '.sww',
1346                outname,
1347                quantity = 'elevation',
1348                cellsize = cellsize,
1349                number_of_decimal_places = 9,
1350                NODATA_value = NODATA_value,
1351                verbose = self.verbose)
1352
1353        #Check header data
1354        from ermapper_grids import read_ermapper_header, read_ermapper_data
1355
1356        header = read_ermapper_header(outname)
1357
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    def test_export_grid_parallel(self):
1396        """Test that sww information can be converted correctly to asc/prj
1397        format readable by e.g. ArcView
1398        """
1399
1400        import time, os
1401        from Scientific.IO.NetCDF import NetCDFFile
1402
1403        base_name = 'tegp'
1404        #Setup
1405        self.domain.set_name(base_name+'_P0_8')
1406        swwfile = self.domain.get_name() + '.sww'
1407
1408        self.domain.set_datadir('.')
1409        self.domain.format = 'sww'
1410        self.domain.smooth = True
1411        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1412        self.domain.set_quantity('stage', 1.0)
1413
1414        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1415
1416        sww = SWW_file(self.domain)
1417        sww.store_connectivity()
1418        sww.store_timestep()
1419        self.domain.evolve_to_end(finaltime = 0.0001)
1420        #Setup
1421        self.domain.set_name(base_name+'_P1_8')
1422        swwfile2 = self.domain.get_name() + '.sww'
1423        sww = SWW_file(self.domain)
1424        sww.store_connectivity()
1425        sww.store_timestep()
1426        self.domain.evolve_to_end(finaltime = 0.0002)
1427        sww.store_timestep()
1428
1429        cellsize = 0.25
1430        #Check contents
1431        #Get NetCDF
1432
1433        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1434
1435        # Get the variables
1436        x = fid.variables['x'][:]
1437        y = fid.variables['y'][:]
1438        z = fid.variables['elevation'][:]
1439        time = fid.variables['time'][:]
1440        stage = fid.variables['stage'][:]
1441
1442        fid.close()
1443
1444        #Export to ascii/prj files
1445        extra_name_out = 'yeah'
1446        sww2dem_batch(base_name,
1447                    quantities = ['elevation', 'depth'],
1448                    extra_name_out = extra_name_out,
1449                    cellsize = cellsize,
1450                    verbose = self.verbose,
1451                    format = 'asc')
1452
1453        prjfile = base_name + '_P0_8_elevation_yeah.prj'
1454        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
1455        #Check asc file
1456        ascid = open(ascfile)
1457        lines = ascid.readlines()
1458        ascid.close()
1459        #Check grid values
1460        for j in range(5):
1461            L = lines[6+j].strip().split()
1462            y = (4-j) * cellsize
1463            for i in range(5):
1464                #print " -i*cellsize - y",  -i*cellsize - y
1465                #print "float(L[i])", float(L[i])
1466                assert num.allclose(float(L[i]), -i*cellsize - y)               
1467        #Cleanup
1468        os.remove(prjfile)
1469        os.remove(ascfile)
1470
1471        prjfile = base_name + '_P1_8_elevation_yeah.prj'
1472        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
1473        #Check asc file
1474        ascid = open(ascfile)
1475        lines = ascid.readlines()
1476        ascid.close()
1477        #Check grid values
1478        for j in range(5):
1479            L = lines[6+j].strip().split()
1480            y = (4-j) * cellsize
1481            for i in range(5):
1482                #print " -i*cellsize - y",  -i*cellsize - y
1483                #print "float(L[i])", float(L[i])
1484                assert num.allclose(float(L[i]), -i*cellsize - y)               
1485        #Cleanup
1486        os.remove(prjfile)
1487        os.remove(ascfile)
1488        os.remove(swwfile)
1489
1490        #Check asc file
1491        ascfile = base_name + '_P0_8_depth_yeah.asc'
1492        prjfile = base_name + '_P0_8_depth_yeah.prj'
1493        ascid = open(ascfile)
1494        lines = ascid.readlines()
1495        ascid.close()
1496        #Check grid values
1497        for j in range(5):
1498            L = lines[6+j].strip().split()
1499            y = (4-j) * cellsize
1500            for i in range(5):
1501                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1502        #Cleanup
1503        os.remove(prjfile)
1504        os.remove(ascfile)
1505
1506        #Check asc file
1507        ascfile = base_name + '_P1_8_depth_yeah.asc'
1508        prjfile = base_name + '_P1_8_depth_yeah.prj'
1509        ascid = open(ascfile)
1510        lines = ascid.readlines()
1511        ascid.close()
1512        #Check grid values
1513        for j in range(5):
1514            L = lines[6+j].strip().split()
1515            y = (4-j) * cellsize
1516            for i in range(5):
1517                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1518        #Cleanup
1519        os.remove(prjfile)
1520        os.remove(ascfile)
1521        os.remove(swwfile2)
1522
1523
1524    def test_export_grid(self):
1525        """
1526        test_export_grid(self):
1527        Test that sww information can be converted correctly to asc/prj
1528        format readable by e.g. ArcView
1529        """
1530
1531        import time, os
1532        from Scientific.IO.NetCDF import NetCDFFile
1533
1534        try:
1535            os.remove('teg*.sww')
1536        except:
1537            pass
1538
1539        #Setup
1540        self.domain.set_name('teg')
1541
1542        prjfile = self.domain.get_name() + '_elevation.prj'
1543        ascfile = self.domain.get_name() + '_elevation.asc'
1544        swwfile = self.domain.get_name() + '.sww'
1545
1546        self.domain.set_datadir('.')
1547        self.domain.smooth = True
1548        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1549        self.domain.set_quantity('stage', 1.0)
1550
1551        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1552
1553        sww = SWW_file(self.domain)
1554        sww.store_connectivity()
1555        sww.store_timestep()
1556        self.domain.evolve_to_end(finaltime = 0.01)
1557        sww.store_timestep()
1558
1559        cellsize = 0.25
1560        #Check contents
1561        #Get NetCDF
1562
1563        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1564
1565        # Get the variables
1566        x = fid.variables['x'][:]
1567        y = fid.variables['y'][:]
1568        z = fid.variables['elevation'][:]
1569        time = fid.variables['time'][:]
1570        stage = fid.variables['stage'][:]
1571
1572        fid.close()
1573
1574        #Export to ascii/prj files
1575        sww2dem_batch(self.domain.get_name(),
1576                quantities = 'elevation',
1577                cellsize = cellsize,
1578                verbose = self.verbose,
1579                format = 'asc')
1580
1581        #Check asc file
1582        ascid = open(ascfile)
1583        lines = ascid.readlines()
1584        ascid.close()
1585
1586        L = lines[2].strip().split()
1587        assert L[0].strip().lower() == 'xllcorner'
1588        assert num.allclose(float(L[1].strip().lower()), 308500)
1589
1590        L = lines[3].strip().split()
1591        assert L[0].strip().lower() == 'yllcorner'
1592        assert num.allclose(float(L[1].strip().lower()), 6189000)
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                assert num.allclose(float(L[i]), -i*cellsize - y)
1600               
1601        #Cleanup
1602        os.remove(prjfile)
1603        os.remove(ascfile)
1604        os.remove(swwfile)
1605
1606    def test_export_gridII(self):
1607        """
1608        test_export_gridII(self):
1609        Test that sww information can be converted correctly to asc/prj
1610        format readable by e.g. ArcView
1611        """
1612
1613        import time, os
1614        from Scientific.IO.NetCDF import NetCDFFile
1615
1616        try:
1617            os.remove('teg*.sww')
1618        except:
1619            pass
1620
1621        #Setup
1622        self.domain.set_name('tegII')
1623
1624        swwfile = self.domain.get_name() + '.sww'
1625
1626        self.domain.set_datadir('.')
1627        self.domain.smooth = True
1628        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1629        self.domain.set_quantity('stage', 1.0)
1630
1631        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1632
1633        sww = SWW_file(self.domain)
1634        sww.store_connectivity()
1635        sww.store_timestep()
1636        self.domain.evolve_to_end(finaltime = 0.01)
1637        sww.store_timestep()
1638
1639        cellsize = 0.25
1640        #Check contents
1641        #Get NetCDF
1642
1643        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1644
1645        # Get the variables
1646        x = fid.variables['x'][:]
1647        y = fid.variables['y'][:]
1648        z = fid.variables['elevation'][:]
1649        time = fid.variables['time'][:]
1650        stage = fid.variables['stage'][:]
1651        xmomentum = fid.variables['xmomentum'][:]
1652        ymomentum = fid.variables['ymomentum'][:]       
1653
1654        #print 'stage', stage
1655        #print 'xmom', xmomentum
1656        #print 'ymom', ymomentum       
1657
1658        fid.close()
1659
1660        #Export to ascii/prj files
1661        if True:
1662            sww2dem_batch(self.domain.get_name(),
1663                        quantities = ['elevation', 'depth'],
1664                        cellsize = cellsize,
1665                        verbose = self.verbose,
1666                        format = 'asc')
1667
1668        else:
1669            sww2dem_batch(self.domain.get_name(),
1670                quantities = ['depth'],
1671                cellsize = cellsize,
1672                verbose = self.verbose,
1673                format = 'asc')
1674
1675
1676            export_grid(self.domain.get_name(),
1677                quantities = ['elevation'],
1678                cellsize = cellsize,
1679                verbose = self.verbose,
1680                format = 'asc')
1681
1682        prjfile = self.domain.get_name() + '_elevation.prj'
1683        ascfile = self.domain.get_name() + '_elevation.asc'
1684       
1685        #Check asc file
1686        ascid = open(ascfile)
1687        lines = ascid.readlines()
1688        ascid.close()
1689
1690        L = lines[2].strip().split()
1691        assert L[0].strip().lower() == 'xllcorner'
1692        assert num.allclose(float(L[1].strip().lower()), 308500)
1693
1694        L = lines[3].strip().split()
1695        assert L[0].strip().lower() == 'yllcorner'
1696        assert num.allclose(float(L[1].strip().lower()), 6189000)
1697
1698        #print "ascfile", ascfile
1699        #Check grid values
1700        for j in range(5):
1701            L = lines[6+j].strip().split()
1702            y = (4-j) * cellsize
1703            for i in range(5):
1704                #print " -i*cellsize - y",  -i*cellsize - y
1705                #print "float(L[i])", float(L[i])
1706                assert num.allclose(float(L[i]), -i*cellsize - y)
1707
1708        #Cleanup
1709        os.remove(prjfile)
1710        os.remove(ascfile)
1711       
1712        #Check asc file
1713        ascfile = self.domain.get_name() + '_depth.asc'
1714        prjfile = self.domain.get_name() + '_depth.prj'
1715        ascid = open(ascfile)
1716        lines = ascid.readlines()
1717        ascid.close()
1718
1719        L = lines[2].strip().split()
1720        assert L[0].strip().lower() == 'xllcorner'
1721        assert num.allclose(float(L[1].strip().lower()), 308500)
1722
1723        L = lines[3].strip().split()
1724        assert L[0].strip().lower() == 'yllcorner'
1725        assert num.allclose(float(L[1].strip().lower()), 6189000)
1726
1727        #Check grid values
1728        for j in range(5):
1729            L = lines[6+j].strip().split()
1730            y = (4-j) * cellsize
1731            for i in range(5):
1732                #print " -i*cellsize - y",  -i*cellsize - y
1733                #print "float(L[i])", float(L[i])               
1734                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1735
1736        #Cleanup
1737        os.remove(prjfile)
1738        os.remove(ascfile)
1739        os.remove(swwfile)
1740
1741
1742    def test_export_gridIII(self):
1743        """
1744        test_export_gridIII
1745        Test that sww information can be converted correctly to asc/prj
1746        format readable by e.g. ArcView
1747        """
1748
1749        import time, os
1750        from Scientific.IO.NetCDF import NetCDFFile
1751
1752        try:
1753            os.remove('teg*.sww')
1754        except:
1755            pass
1756
1757        #Setup
1758       
1759        self.domain.set_name('tegIII')
1760
1761        swwfile = self.domain.get_name() + '.sww'
1762
1763        self.domain.set_datadir('.')
1764        self.domain.format = 'sww'
1765        self.domain.smooth = True
1766        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1767        self.domain.set_quantity('stage', 1.0)
1768
1769        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1770       
1771        sww = SWW_file(self.domain)
1772        sww.store_connectivity()
1773        sww.store_timestep() #'stage')
1774        self.domain.evolve_to_end(finaltime = 0.01)
1775        sww.store_timestep() #'stage')
1776
1777        cellsize = 0.25
1778        #Check contents
1779        #Get NetCDF
1780
1781        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1782
1783        # Get the variables
1784        x = fid.variables['x'][:]
1785        y = fid.variables['y'][:]
1786        z = fid.variables['elevation'][:]
1787        time = fid.variables['time'][:]
1788        stage = fid.variables['stage'][:]
1789
1790        fid.close()
1791
1792        #Export to ascii/prj files
1793        extra_name_out = 'yeah'
1794        if True:
1795            sww2dem_batch(self.domain.get_name(),
1796                        quantities = ['elevation', 'depth'],
1797                        extra_name_out = extra_name_out,
1798                        cellsize = cellsize,
1799                        verbose = self.verbose,
1800                        format = 'asc')
1801
1802        else:
1803            sww2dem_batch(self.domain.get_name(),
1804                quantities = ['depth'],
1805                cellsize = cellsize,
1806                verbose = self.verbose,
1807                format = 'asc')
1808
1809
1810            sww2dem_batch(self.domain.get_name(),
1811                quantities = ['elevation'],
1812                cellsize = cellsize,
1813                verbose = self.verbose,
1814                format = 'asc')
1815
1816        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
1817        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
1818       
1819        #Check asc file
1820        ascid = open(ascfile)
1821        lines = ascid.readlines()
1822        ascid.close()
1823
1824        L = lines[2].strip().split()
1825        assert L[0].strip().lower() == 'xllcorner'
1826        assert num.allclose(float(L[1].strip().lower()), 308500)
1827
1828        L = lines[3].strip().split()
1829        assert L[0].strip().lower() == 'yllcorner'
1830        assert num.allclose(float(L[1].strip().lower()), 6189000)
1831
1832        #print "ascfile", ascfile
1833        #Check grid values
1834        for j in range(5):
1835            L = lines[6+j].strip().split()
1836            y = (4-j) * cellsize
1837            for i in range(5):
1838                #print " -i*cellsize - y",  -i*cellsize - y
1839                #print "float(L[i])", float(L[i])
1840                assert num.allclose(float(L[i]), -i*cellsize - y)
1841               
1842        #Cleanup
1843        os.remove(prjfile)
1844        os.remove(ascfile)
1845       
1846        #Check asc file
1847        ascfile = self.domain.get_name() + '_depth_yeah.asc'
1848        prjfile = self.domain.get_name() + '_depth_yeah.prj'
1849        ascid = open(ascfile)
1850        lines = ascid.readlines()
1851        ascid.close()
1852
1853        L = lines[2].strip().split()
1854        assert L[0].strip().lower() == 'xllcorner'
1855        assert num.allclose(float(L[1].strip().lower()), 308500)
1856
1857        L = lines[3].strip().split()
1858        assert L[0].strip().lower() == 'yllcorner'
1859        assert num.allclose(float(L[1].strip().lower()), 6189000)
1860
1861        #Check grid values
1862        for j in range(5):
1863            L = lines[6+j].strip().split()
1864            y = (4-j) * cellsize
1865            for i in range(5):
1866                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1867
1868        #Cleanup
1869        os.remove(prjfile)
1870        os.remove(ascfile)
1871        os.remove(swwfile)
1872
1873    def test_export_grid_bad(self):
1874        """Test that sww information can be converted correctly to asc/prj
1875        format readable by e.g. ArcView
1876        """
1877
1878        try:
1879            sww2dem_batch('a_small_round-egg',
1880                        quantities = ['elevation', 'depth'],
1881                        cellsize = 99,
1882                        verbose = self.verbose,
1883                        format = 'asc')
1884        except IOError:
1885            pass
1886        else:
1887            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
1888       
1889
1890#################################################################################
1891
1892if __name__ == "__main__":
1893    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
1894    suite = unittest.makeSuite(Test_Sww2Dem, 'test')
1895    runner = unittest.TextTestRunner(verbosity=1)
1896    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.