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

Last change on this file since 8634 was 8634, checked in by steve, 11 years ago

Found bug causing segmentation fault in calc_grid_values_ext.c
File size: 56.0 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 example has rows with a large number of zeros
457
458        ncols         2001
459        nrows         2
460        xllcorner     308500
461        yllcorner     6189000
462        cellsize      1.000000
463        NODATA_value  -9999
464        0.0 ....
465        """
466
467        import time, os
468        from Scientific.IO.NetCDF import NetCDFFile
469
470        #Setup
471
472        from mesh_factory import rectangular_cross
473
474        #Create basic mesh (100m x 100m)
475        points, vertices, boundary = rectangular_cross(20, 1, 20.0, 1.0)
476
477        #Create shallow water domain
478        domain = Domain(points, vertices, boundary)
479        domain.default_order = 1
480
481        domain.set_name('datatest')
482
483        prjfile = domain.get_name() + '_elevation.prj'
484        ascfile = domain.get_name() + '_elevation.asc'
485        swwfile = domain.get_name() + '.sww'
486
487        domain.set_datadir('.')
488        domain.format = 'sww'
489        domain.smooth = True
490        domain.geo_reference = Geo_reference(56, 308500, 6189000)
491
492        #
493        domain.set_quantity('elevation', 0)
494        domain.set_quantity('stage', 0)
495
496        B = Transmissive_boundary(domain)
497        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
498
499
500        #
501        sww = SWW_file(domain)
502        sww.store_connectivity()
503        sww.store_timestep()
504       
505        domain.tight_slope_limiters = 1
506        domain.evolve_to_end(finaltime = 0.01)
507        sww.store_timestep()
508
509        cellsize = 1.0  #0.1 grid
510
511
512        #Check contents
513        #Get NetCDF
514
515        fid = NetCDFFile(sww.filename, netcdf_mode_r)
516
517        # Get the variables
518        x = fid.variables['x'][:]
519        y = fid.variables['y'][:]
520        z = fid.variables['elevation'][:]
521        time = fid.variables['time'][:]
522        stage = fid.variables['stage'][:]
523
524
525        #Export to ascii/prj files
526        sww2dem(domain.get_name()+'.sww',
527                domain.get_name() + '_elevation.asc',
528                quantity = 'elevation',
529                cellsize = cellsize,
530                number_of_decimal_places = 9,
531                verbose = self.verbose,
532                block_size=2)
533
534
535        #Check prj (meta data)
536        prjid = open(prjfile)
537        lines = prjid.readlines()
538        prjid.close()
539
540
541        L = lines[0].strip().split()
542        assert L[0].strip().lower() == 'projection'
543        assert L[1].strip().lower() == 'utm'
544
545        L = lines[1].strip().split()
546        assert L[0].strip().lower() == 'zone'
547        assert L[1].strip().lower() == '56'
548
549        L = lines[2].strip().split()
550        assert L[0].strip().lower() == 'datum'
551        assert L[1].strip().lower() == 'wgs84'
552
553        L = lines[3].strip().split()
554        assert L[0].strip().lower() == 'zunits'
555        assert L[1].strip().lower() == 'no'
556
557        L = lines[4].strip().split()
558        assert L[0].strip().lower() == 'units'
559        assert L[1].strip().lower() == 'meters'
560
561        L = lines[5].strip().split()
562        assert L[0].strip().lower() == 'spheroid'
563        assert L[1].strip().lower() == 'wgs84'
564
565        L = lines[6].strip().split()
566        assert L[0].strip().lower() == 'xshift'
567        assert L[1].strip().lower() == '500000'
568
569        L = lines[7].strip().split()
570        assert L[0].strip().lower() == 'yshift'
571        assert L[1].strip().lower() == '10000000'
572
573        L = lines[8].strip().split()
574        assert L[0].strip().lower() == 'parameters'
575
576
577        #Check asc file
578        ascid = open(ascfile)
579        lines = ascid.readlines()
580        ascid.close()
581
582
583
584        L = lines[0].strip().split()
585        assert L[0].strip().lower() == 'ncols'
586        assert L[1].strip().lower() == '21'
587
588        L = lines[1].strip().split()
589        assert L[0].strip().lower() == 'nrows'
590        assert L[1].strip().lower() == '2'
591
592        L = lines[2].strip().split()
593        assert L[0].strip().lower() == 'xllcorner'
594        assert num.allclose(float(L[1].strip().lower()), 308500)
595
596        L = lines[3].strip().split()
597        assert L[0].strip().lower() == 'yllcorner'
598        assert num.allclose(float(L[1].strip().lower()), 6189000)
599
600        L = lines[4].strip().split()
601        assert L[0].strip().lower() == 'cellsize'
602        assert num.allclose(float(L[1].strip().lower()), cellsize)
603
604        L = lines[5].strip().split()
605        assert L[0].strip() == 'NODATA_value'
606        assert L[1].strip().lower() == '-9999'
607
608        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
609        for i, line in enumerate(lines[6:]):
610            for j, value in enumerate( line.split() ):
611                #print value
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 mesh can be restricted by bounding box
633
634        Original extent is 100m x 100m:
635
636        Eastings:   308500 -  308600
637        Northings: 6189000 - 6189100
638
639        Bounding box changes this to the 50m x 50m square defined by
640
641        Eastings:   308530 -  308570
642        Northings: 6189050 - 6189100
643
644        The cropped values should be
645
646         -130 -140 -150 -160 -170
647         -120 -130 -140 -150 -160
648         -110 -120 -130 -140 -150
649         -100 -110 -120 -130 -140
650          -90 -100 -110 -120 -130
651          -80  -90 -100 -110 -120
652
653        and the new lower reference point should be
654        Eastings:   308530
655        Northings: 6189050
656
657        Original dataset is the same as in test_sww2dem_larger()
658
659        """
660
661        import time, os
662        from Scientific.IO.NetCDF import NetCDFFile
663
664        #Setup
665
666        from mesh_factory import rectangular
667
668        #Create basic mesh (100m x 100m)
669        points, vertices, boundary = rectangular(2, 2, 100, 100)
670
671        #Create shallow water domain
672        domain = Domain(points, vertices, boundary)
673        domain.default_order = 2
674
675        domain.set_name('datatest')
676
677        prjfile = domain.get_name() + '_elevation.prj'
678        ascfile = domain.get_name() + '_elevation.asc'
679        swwfile = domain.get_name() + '.sww'
680
681        domain.set_datadir('.')
682        domain.format = 'sww'
683        domain.smooth = True
684        domain.geo_reference = Geo_reference(56, 308500, 6189000)
685
686        #
687        domain.set_quantity('elevation', lambda x,y: -x-y)
688        domain.set_quantity('stage', 0)
689
690        B = Transmissive_boundary(domain)
691        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
692
693
694        #
695        sww = SWW_file(domain)
696        sww.store_connectivity()
697        sww.store_timestep()
698
699        #domain.tight_slope_limiters = 1
700        domain.evolve_to_end(finaltime = 0.01)
701        sww.store_timestep()
702
703        cellsize = 10  #10m grid
704
705
706        #Check contents
707        #Get NetCDF
708
709        fid = NetCDFFile(sww.filename, netcdf_mode_r)
710
711        # Get the variables
712        x = fid.variables['x'][:]
713        y = fid.variables['y'][:]
714        z = fid.variables['elevation'][:]
715        time = fid.variables['time'][:]
716        stage = fid.variables['stage'][:]
717
718
719        # Export to ascii/prj files
720        sww2dem(domain.get_name() + '.sww',
721                domain.get_name() + '_elevation.asc',
722                quantity = 'elevation',
723                cellsize = cellsize,
724                number_of_decimal_places = 9,
725                easting_min = 308530,
726                easting_max = 308570,
727                northing_min = 6189050,
728                northing_max = 6189100,
729                verbose = self.verbose)
730
731        fid.close()
732
733
734        # Check prj (meta data)
735        prjid = open(prjfile)
736        lines = prjid.readlines()
737        prjid.close()
738
739        L = lines[0].strip().split()
740        assert L[0].strip().lower() == 'projection'
741        assert L[1].strip().lower() == 'utm'
742
743        L = lines[1].strip().split()
744        assert L[0].strip().lower() == 'zone'
745        assert L[1].strip().lower() == '56'
746
747        L = lines[2].strip().split()
748        assert L[0].strip().lower() == 'datum'
749        assert L[1].strip().lower() == 'wgs84'
750
751        L = lines[3].strip().split()
752        assert L[0].strip().lower() == 'zunits'
753        assert L[1].strip().lower() == 'no'
754
755        L = lines[4].strip().split()
756        assert L[0].strip().lower() == 'units'
757        assert L[1].strip().lower() == 'meters'
758
759        L = lines[5].strip().split()
760        assert L[0].strip().lower() == 'spheroid'
761        assert L[1].strip().lower() == 'wgs84'
762
763        L = lines[6].strip().split()
764        assert L[0].strip().lower() == 'xshift'
765        assert L[1].strip().lower() == '500000'
766
767        L = lines[7].strip().split()
768        assert L[0].strip().lower() == 'yshift'
769        assert L[1].strip().lower() == '10000000'
770
771        L = lines[8].strip().split()
772        assert L[0].strip().lower() == 'parameters'
773
774
775        #Check asc file
776        ascid = open(ascfile)
777        lines = ascid.readlines()
778        ascid.close()
779
780        L = lines[0].strip().split()
781        assert L[0].strip().lower() == 'ncols'
782        assert L[1].strip().lower() == '5'
783
784        L = lines[1].strip().split()
785        assert L[0].strip().lower() == 'nrows'
786        assert L[1].strip().lower() == '6'
787
788        L = lines[2].strip().split()
789        assert L[0].strip().lower() == 'xllcorner'
790        assert num.allclose(float(L[1].strip().lower()), 308530)
791
792        L = lines[3].strip().split()
793        assert L[0].strip().lower() == 'yllcorner'
794        assert num.allclose(float(L[1].strip().lower()), 6189050)
795
796        L = lines[4].strip().split()
797        assert L[0].strip().lower() == 'cellsize'
798        assert num.allclose(float(L[1].strip().lower()), cellsize)
799
800        L = lines[5].strip().split()
801        assert L[0].strip() == 'NODATA_value'
802        assert L[1].strip().lower() == '-9999'
803
804        #Check grid values
805        for i, line in enumerate(lines[6:]):
806            for j, value in enumerate( line.split() ):
807                #assert float(value) == -(10-i+j)*cellsize
808                assert float(value) == -(10-i+j+3)*cellsize
809
810
811
812        #Cleanup
813        os.remove(prjfile)
814        os.remove(ascfile)
815        os.remove(swwfile)
816
817
818
819    def test_sww2dem_asc_stage_reduction(self):
820        """Test that sww information can be converted correctly to asc/prj
821        format readable by e.g. ArcView
822
823        This tests the reduction of quantity stage using min
824        """
825
826        import time, os
827        from Scientific.IO.NetCDF import NetCDFFile
828
829        #Setup
830        self.domain.set_name('datatest')
831
832        prjfile = self.domain.get_name() + '_stage.prj'
833        ascfile = self.domain.get_name() + '_stage.asc'
834        swwfile = self.domain.get_name() + '.sww'
835
836        self.domain.set_datadir('.')
837        self.domain.format = 'sww'
838        self.domain.smooth = True
839        self.domain.set_quantity('elevation', lambda x,y: -x-y)
840
841        self.domain.geo_reference = Geo_reference(56,308500,6189000)
842
843
844        sww = SWW_file(self.domain)
845        sww.store_connectivity()
846        sww.store_timestep()
847
848        #self.domain.tight_slope_limiters = 1
849        self.domain.evolve_to_end(finaltime = 0.01)
850        sww.store_timestep()
851
852        cellsize = 0.25
853        #Check contents
854        #Get NetCDF
855
856        fid = NetCDFFile(sww.filename, netcdf_mode_r)
857
858        # Get the variables
859        x = fid.variables['x'][:]
860        y = fid.variables['y'][:]
861        z = fid.variables['elevation'][:]
862        time = fid.variables['time'][:]
863        stage = fid.variables['stage'][:]
864
865
866        #Export to ascii/prj files
867        sww2dem(self.domain.get_name() + '.sww',
868                self.domain.get_name() + '_stage.asc',
869                quantity = 'stage',
870                cellsize = cellsize,
871                number_of_decimal_places = 9,
872                reduction = min)
873
874
875        #Check asc file
876        ascid = open(ascfile)
877        lines = ascid.readlines()
878        ascid.close()
879
880        L = lines[0].strip().split()
881        assert L[0].strip().lower() == 'ncols'
882        assert L[1].strip().lower() == '5'
883
884        L = lines[1].strip().split()
885        assert L[0].strip().lower() == 'nrows'
886        assert L[1].strip().lower() == '5'
887
888        L = lines[2].strip().split()
889        assert L[0].strip().lower() == 'xllcorner'
890        assert num.allclose(float(L[1].strip().lower()), 308500)
891
892        L = lines[3].strip().split()
893        assert L[0].strip().lower() == 'yllcorner'
894        assert num.allclose(float(L[1].strip().lower()), 6189000)
895
896        L = lines[4].strip().split()
897        assert L[0].strip().lower() == 'cellsize'
898        assert num.allclose(float(L[1].strip().lower()), cellsize)
899
900        L = lines[5].strip().split()
901        assert L[0].strip() == 'NODATA_value'
902        assert L[1].strip().lower() == '-9999'
903
904
905        #Check grid values (where applicable)
906        for j in range(5):
907            if j%2 == 0:
908                L = lines[6+j].strip().split()
909                jj = 4-j
910                for i in range(5):
911                    if i%2 == 0:
912                        index = jj/2 + i/2*3
913                        val0 = stage[0,index]
914                        val1 = stage[1,index]
915
916                        #print i, j, index, ':', L[i], val0, val1
917                        assert num.allclose(float(L[i]), min(val0, val1))
918
919
920        fid.close()
921
922        #Cleanup
923        os.remove(prjfile)
924        os.remove(ascfile)
925        os.remove(swwfile)
926
927    def test_sww2dem_asc_stage_time(self):
928        """Test that sww information can be converted correctly to asc/prj
929        format readable by e.g. ArcView
930
931        This tests the reduction of quantity stage using min
932        """
933
934        import time, os
935        from Scientific.IO.NetCDF import NetCDFFile
936
937        #Setup
938        self.domain.set_name('datatest')
939
940        prjfile = self.domain.get_name() + '_stage.prj'
941        ascfile = self.domain.get_name() + '_stage.asc'
942        swwfile = self.domain.get_name() + '.sww'
943
944        self.domain.set_datadir('.')
945        self.domain.format = 'sww'
946        self.domain.smooth = True
947        self.domain.set_quantity('elevation', lambda x,y: -x-y)
948
949        self.domain.geo_reference = Geo_reference(56,308500,6189000)
950
951        sww = SWW_file(self.domain)
952        sww.store_connectivity()
953        sww.store_timestep()
954
955        #self.domain.tight_slope_limiters = 1
956        self.domain.evolve_to_end(finaltime = 0.01)
957        sww.store_timestep()
958
959        cellsize = 0.25
960        #Check contents
961        #Get NetCDF
962
963        fid = NetCDFFile(sww.filename, netcdf_mode_r)
964
965        # Get the variables
966        x = fid.variables['x'][:]
967        y = fid.variables['y'][:]
968        z = fid.variables['elevation'][:]
969        time = fid.variables['time'][:]
970        stage = fid.variables['stage'][:]
971
972        #Export to ascii/prj files
973        sww2dem(self.domain.get_name() + '.sww',
974                self.domain.get_name() + '_stage.asc',
975                quantity = 'stage',
976                cellsize = cellsize,
977                number_of_decimal_places = 9,
978                reduction = 1,
979                verbose=self.verbose)
980
981
982        #Check asc file
983        ascid = open(ascfile)
984        lines = ascid.readlines()
985        ascid.close()
986
987        L = lines[0].strip().split()
988        assert L[0].strip().lower() == 'ncols'
989        assert L[1].strip().lower() == '5'
990
991        L = lines[1].strip().split()
992        assert L[0].strip().lower() == 'nrows'
993        assert L[1].strip().lower() == '5'
994
995        L = lines[2].strip().split()
996        assert L[0].strip().lower() == 'xllcorner'
997        assert num.allclose(float(L[1].strip().lower()), 308500)
998
999        L = lines[3].strip().split()
1000        assert L[0].strip().lower() == 'yllcorner'
1001        assert num.allclose(float(L[1].strip().lower()), 6189000)
1002
1003        L = lines[4].strip().split()
1004        assert L[0].strip().lower() == 'cellsize'
1005        assert num.allclose(float(L[1].strip().lower()), cellsize)
1006
1007        L = lines[5].strip().split()
1008        assert L[0].strip() == 'NODATA_value'
1009        assert L[1].strip().lower() == '-9999'
1010
1011        #Check grid values (where applicable)
1012        for j in range(5):
1013            if j%2 == 0:
1014                L = lines[6+j].strip().split()
1015                jj = 4-j
1016                for i in range(5):
1017                    if i%2 == 0:
1018                        index = jj/2 + i/2*3
1019                       
1020                        val = stage[1,index]
1021                   
1022                        assert num.allclose(float(L[i]), val)
1023
1024        fid.close()
1025
1026        #Cleanup
1027        os.remove(prjfile)
1028        os.remove(ascfile)
1029        os.remove(swwfile)
1030
1031
1032    def test_sww2dem_asc_derived_quantity(self):
1033        """Test that sww information can be converted correctly to asc/prj
1034        format readable by e.g. ArcView
1035
1036        This tests the use of derived quantities
1037        """
1038
1039        import time, os
1040        from Scientific.IO.NetCDF import NetCDFFile
1041
1042        #Setup
1043        self.domain.set_name('datatest')
1044
1045        prjfile = self.domain.get_name() + '_depth.prj'
1046        ascfile = self.domain.get_name() + '_depth.asc'
1047        swwfile = self.domain.get_name() + '.sww'
1048
1049        self.domain.set_datadir('.')
1050        self.domain.format = 'sww'
1051        self.domain.smooth = True
1052        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1053        self.domain.set_quantity('stage', 0.0)
1054
1055        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1056
1057
1058        sww = SWW_file(self.domain)
1059        sww.store_connectivity()
1060        sww.store_timestep()
1061
1062        #self.domain.tight_slope_limiters = 1
1063        self.domain.evolve_to_end(finaltime = 0.01)
1064        sww.store_timestep()
1065
1066        cellsize = 0.25
1067        #Check contents
1068        #Get NetCDF
1069
1070        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1071
1072        # Get the variables
1073        x = fid.variables['x'][:]
1074        y = fid.variables['y'][:]
1075        z = fid.variables['elevation'][:]
1076        time = fid.variables['time'][:]
1077        stage = fid.variables['stage'][:]
1078
1079
1080        #Export to ascii/prj files
1081        sww2dem(self.domain.get_name()+'.sww',
1082                name_out = 'datatest_depth.asc',
1083                quantity = 'stage - elevation',
1084                cellsize = cellsize,
1085                number_of_decimal_places = 9,
1086                reduction = min,
1087                verbose = self.verbose)
1088
1089
1090        #Check asc file
1091        ascid = open(ascfile)
1092        lines = ascid.readlines()
1093        ascid.close()
1094
1095        L = lines[0].strip().split()
1096        assert L[0].strip().lower() == 'ncols'
1097        assert L[1].strip().lower() == '5'
1098
1099        L = lines[1].strip().split()
1100        assert L[0].strip().lower() == 'nrows'
1101        assert L[1].strip().lower() == '5'
1102
1103        L = lines[2].strip().split()
1104        assert L[0].strip().lower() == 'xllcorner'
1105        assert num.allclose(float(L[1].strip().lower()), 308500)
1106
1107        L = lines[3].strip().split()
1108        assert L[0].strip().lower() == 'yllcorner'
1109        assert num.allclose(float(L[1].strip().lower()), 6189000)
1110
1111        L = lines[4].strip().split()
1112        assert L[0].strip().lower() == 'cellsize'
1113        assert num.allclose(float(L[1].strip().lower()), cellsize)
1114
1115        L = lines[5].strip().split()
1116        assert L[0].strip() == 'NODATA_value'
1117        assert L[1].strip().lower() == '-9999'
1118
1119
1120        #Check grid values (where applicable)
1121        for j in range(5):
1122            if j%2 == 0:
1123                L = lines[6+j].strip().split()
1124                jj = 4-j
1125                for i in range(5):
1126                    if i%2 == 0:
1127                        index = jj/2 + i/2*3
1128                        val0 = stage[0,index] - z[index]
1129                        val1 = stage[1,index] - z[index]
1130
1131                        #print i, j, index, ':', L[i], val0, val1
1132                        assert num.allclose(float(L[i]), min(val0, val1))
1133
1134
1135        fid.close()
1136
1137        #Cleanup
1138        os.remove(prjfile)
1139        os.remove(ascfile)
1140        os.remove(swwfile)
1141
1142
1143
1144
1145
1146    def test_sww2dem_asc_missing_points(self):
1147        """Test that sww information can be converted correctly to asc/prj
1148        format readable by e.g. ArcView
1149
1150        This test includes the writing of missing values
1151        """
1152
1153        import time, os
1154        from Scientific.IO.NetCDF import NetCDFFile
1155
1156        #Setup mesh not coinciding with rectangle.
1157        #This will cause missing values to occur in gridded data
1158
1159
1160        points = [                        [1.0, 1.0],
1161                              [0.5, 0.5], [1.0, 0.5],
1162                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1163
1164        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1165
1166        #Create shallow water domain
1167        domain = Domain(points, vertices)
1168        domain.default_order=2
1169
1170
1171        #Set some field values
1172        domain.set_quantity('elevation', lambda x,y: -x-y)
1173        domain.set_quantity('friction', 0.03)
1174
1175
1176        ######################
1177        # Boundary conditions
1178        B = Transmissive_boundary(domain)
1179        domain.set_boundary( {'exterior': B} )
1180
1181
1182        ######################
1183        #Initial condition - with jumps
1184
1185        bed = domain.quantities['elevation'].vertex_values
1186        stage = num.zeros(bed.shape, num.float)
1187
1188        h = 0.3
1189        for i in range(stage.shape[0]):
1190            if i % 2 == 0:
1191                stage[i,:] = bed[i,:] + h
1192            else:
1193                stage[i,:] = bed[i,:]
1194
1195        domain.set_quantity('stage', stage)
1196        domain.distribute_to_vertices_and_edges()
1197
1198        domain.set_name('datatest')
1199
1200        prjfile = domain.get_name() + '_elevation.prj'
1201        ascfile = domain.get_name() + '_elevation.asc'
1202        swwfile = domain.get_name() + '.sww'
1203
1204        domain.set_datadir('.')
1205        domain.format = 'sww'
1206        domain.smooth = True
1207
1208        domain.geo_reference = Geo_reference(56,308500,6189000)
1209
1210        sww = SWW_file(domain)
1211        sww.store_connectivity()
1212        sww.store_timestep()
1213
1214        cellsize = 0.25
1215        #Check contents
1216        #Get NetCDF
1217
1218        fid = NetCDFFile(swwfile, netcdf_mode_r)
1219
1220        # Get the variables
1221        x = fid.variables['x'][:]
1222        y = fid.variables['y'][:]
1223        z = fid.variables['elevation'][:]
1224        time = fid.variables['time'][:]
1225
1226        try:
1227            geo_reference = Geo_reference(NetCDFObject=fid)
1228        except AttributeError, e:
1229            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1230
1231        #Export to ascii/prj files
1232        sww2dem(domain.get_name()+'.sww',
1233                domain.get_name()+'_elevation.asc',
1234                quantity = 'elevation',
1235                cellsize = cellsize,
1236                number_of_decimal_places = 9,
1237                verbose = self.verbose)
1238
1239
1240        #Check asc file
1241        ascid = open(ascfile)
1242        lines = ascid.readlines()
1243        ascid.close()
1244
1245        L = lines[0].strip().split()
1246        assert L[0].strip().lower() == 'ncols'
1247        assert L[1].strip().lower() == '5'
1248
1249        L = lines[1].strip().split()
1250        assert L[0].strip().lower() == 'nrows'
1251        assert L[1].strip().lower() == '5'
1252
1253        L = lines[2].strip().split()
1254        assert L[0].strip().lower() == 'xllcorner'
1255        assert num.allclose(float(L[1].strip().lower()), 308500)
1256
1257        L = lines[3].strip().split()
1258        assert L[0].strip().lower() == 'yllcorner'
1259        assert num.allclose(float(L[1].strip().lower()), 6189000)
1260
1261        L = lines[4].strip().split()
1262        assert L[0].strip().lower() == 'cellsize'
1263        assert num.allclose(float(L[1].strip().lower()), cellsize)
1264
1265        L = lines[5].strip().split()
1266        assert L[0].strip() == 'NODATA_value'
1267        assert L[1].strip().lower() == '-9999'
1268
1269        #Check grid values
1270        for j in range(5):
1271            L = lines[6+j].strip().split()
1272            assert len(L) == 5
1273            y = (4-j) * cellsize
1274
1275            for i in range(5):
1276                #print i
1277                if i+j >= 4:
1278                    assert num.allclose(float(L[i]), -i*cellsize - y)
1279                else:
1280                    #Missing values
1281                    assert num.allclose(float(L[i]), -9999)
1282
1283
1284
1285        fid.close()
1286
1287        #Cleanup
1288        os.remove(prjfile)
1289        os.remove(ascfile)
1290        os.remove(swwfile)
1291
1292
1293
1294    def test_sww2ers_simple(self):
1295        """Test that sww information can be converted correctly to ers
1296        format
1297        """
1298
1299        import time, os
1300        from Scientific.IO.NetCDF import NetCDFFile
1301
1302
1303        NODATA_value = 1758323
1304
1305        #Setup
1306        self.domain.set_name('datatest')
1307
1308        headerfile = self.domain.get_name() + '.ers'
1309        swwfile = self.domain.get_name() + '.sww'
1310
1311        self.domain.set_datadir('.')
1312        self.domain.format = 'sww'
1313        self.domain.smooth = True
1314        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1315
1316        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1317
1318        sww = SWW_file(self.domain)
1319        sww.store_connectivity()
1320        sww.store_timestep()
1321
1322        #self.domain.tight_slope_limiters = 1
1323        self.domain.evolve_to_end(finaltime = 0.01)
1324        sww.store_timestep()
1325
1326        cellsize = 0.25
1327        #Check contents
1328        #Get NetCDF
1329
1330        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1331
1332        # Get the variables
1333        x = fid.variables['x'][:]
1334        y = fid.variables['y'][:]
1335        z = fid.variables['elevation'][:]
1336        time = fid.variables['time'][:]
1337        stage = fid.variables['stage'][:]
1338
1339
1340        #Export to ers files
1341        outname = self.domain.get_name() + '_elevation.ers'
1342        sww2dem(self.domain.get_name() + '.sww',
1343                outname,
1344                quantity = 'elevation',
1345                cellsize = cellsize,
1346                number_of_decimal_places = 9,
1347                NODATA_value = NODATA_value,
1348                verbose = self.verbose)
1349
1350        #Check header data
1351        from ermapper_grids import read_ermapper_header, read_ermapper_data
1352
1353        header = read_ermapper_header(outname)
1354
1355        assert header['projection'].lower() == '"utm-56"'
1356        assert header['datum'].lower() == '"wgs84"'
1357        assert header['units'].lower() == '"meters"'
1358        assert header['value'].lower() == '"elevation"'
1359        assert header['xdimension'] == '0.25'
1360        assert header['ydimension'] == '0.25'
1361        assert float(header['eastings']) == 308500.0   #xllcorner
1362        assert float(header['northings']) == 6189000.0 #yllcorner
1363        assert int(header['nroflines']) == 5
1364        assert int(header['nrofcellsperline']) == 5
1365        assert int(header['nullcellvalue']) == NODATA_value
1366        #FIXME - there is more in the header
1367
1368
1369        #Check grid data
1370        grid = read_ermapper_data(self.domain.get_name() + '_elevation')
1371
1372
1373        ref_grid = [-1,    -1.25, -1.5,  -1.75, -2.0,
1374                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
1375                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
1376                    -0.25, -0.5,  -0.75, -1.0,  -1.25,
1377                    -0.0,  -0.25, -0.5,  -0.75, -1.0]
1378
1379
1380        #print grid.reshape((5,5))
1381        assert num.allclose(grid, ref_grid)
1382
1383        fid.close()
1384
1385        #Cleanup
1386        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
1387        #Done (Ole) - it was because sww2ers didn't close it's sww file
1388        os.remove(sww.filename)
1389        os.remove(self.domain.get_name() + '_elevation')
1390        os.remove(self.domain.get_name() + '_elevation.ers')
1391       
1392    def test_export_grid_parallel(self):
1393        """Test that sww information can be converted correctly to asc/prj
1394        format readable by e.g. ArcView
1395        """
1396
1397        import time, os
1398        from Scientific.IO.NetCDF import NetCDFFile
1399
1400        base_name = 'tegp'
1401        #Setup
1402        self.domain.set_name(base_name+'_P0_8')
1403        swwfile = self.domain.get_name() + '.sww'
1404
1405        self.domain.set_datadir('.')
1406        self.domain.format = 'sww'
1407        self.domain.smooth = True
1408        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1409        self.domain.set_quantity('stage', 1.0)
1410
1411        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1412
1413        sww = SWW_file(self.domain)
1414        sww.store_connectivity()
1415        sww.store_timestep()
1416        self.domain.evolve_to_end(finaltime = 0.0001)
1417        #Setup
1418        self.domain.set_name(base_name+'_P1_8')
1419        swwfile2 = self.domain.get_name() + '.sww'
1420        sww = SWW_file(self.domain)
1421        sww.store_connectivity()
1422        sww.store_timestep()
1423        self.domain.evolve_to_end(finaltime = 0.0002)
1424        sww.store_timestep()
1425
1426        cellsize = 0.25
1427        #Check contents
1428        #Get NetCDF
1429
1430        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1431
1432        # Get the variables
1433        x = fid.variables['x'][:]
1434        y = fid.variables['y'][:]
1435        z = fid.variables['elevation'][:]
1436        time = fid.variables['time'][:]
1437        stage = fid.variables['stage'][:]
1438
1439        fid.close()
1440
1441        #Export to ascii/prj files
1442        extra_name_out = 'yeah'
1443        sww2dem_batch(base_name,
1444                    quantities = ['elevation', 'depth'],
1445                    extra_name_out = extra_name_out,
1446                    cellsize = cellsize,
1447                    verbose = self.verbose,
1448                    format = 'asc')
1449
1450        prjfile = base_name + '_P0_8_elevation_yeah.prj'
1451        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
1452        #Check asc file
1453        ascid = open(ascfile)
1454        lines = ascid.readlines()
1455        ascid.close()
1456        #Check grid values
1457        for j in range(5):
1458            L = lines[6+j].strip().split()
1459            y = (4-j) * cellsize
1460            for i in range(5):
1461                #print " -i*cellsize - y",  -i*cellsize - y
1462                #print "float(L[i])", float(L[i])
1463                assert num.allclose(float(L[i]), -i*cellsize - y)               
1464        #Cleanup
1465        os.remove(prjfile)
1466        os.remove(ascfile)
1467
1468        prjfile = base_name + '_P1_8_elevation_yeah.prj'
1469        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
1470        #Check asc file
1471        ascid = open(ascfile)
1472        lines = ascid.readlines()
1473        ascid.close()
1474        #Check grid values
1475        for j in range(5):
1476            L = lines[6+j].strip().split()
1477            y = (4-j) * cellsize
1478            for i in range(5):
1479                #print " -i*cellsize - y",  -i*cellsize - y
1480                #print "float(L[i])", float(L[i])
1481                assert num.allclose(float(L[i]), -i*cellsize - y)               
1482        #Cleanup
1483        os.remove(prjfile)
1484        os.remove(ascfile)
1485        os.remove(swwfile)
1486
1487        #Check asc file
1488        ascfile = base_name + '_P0_8_depth_yeah.asc'
1489        prjfile = base_name + '_P0_8_depth_yeah.prj'
1490        ascid = open(ascfile)
1491        lines = ascid.readlines()
1492        ascid.close()
1493        #Check grid values
1494        for j in range(5):
1495            L = lines[6+j].strip().split()
1496            y = (4-j) * cellsize
1497            for i in range(5):
1498                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1499        #Cleanup
1500        os.remove(prjfile)
1501        os.remove(ascfile)
1502
1503        #Check asc file
1504        ascfile = base_name + '_P1_8_depth_yeah.asc'
1505        prjfile = base_name + '_P1_8_depth_yeah.prj'
1506        ascid = open(ascfile)
1507        lines = ascid.readlines()
1508        ascid.close()
1509        #Check grid values
1510        for j in range(5):
1511            L = lines[6+j].strip().split()
1512            y = (4-j) * cellsize
1513            for i in range(5):
1514                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1515        #Cleanup
1516        os.remove(prjfile)
1517        os.remove(ascfile)
1518        os.remove(swwfile2)
1519
1520
1521    def test_export_grid(self):
1522        """
1523        test_export_grid(self):
1524        Test that sww information can be converted correctly to asc/prj
1525        format readable by e.g. ArcView
1526        """
1527
1528        import time, os
1529        from Scientific.IO.NetCDF import NetCDFFile
1530
1531        try:
1532            os.remove('teg*.sww')
1533        except:
1534            pass
1535
1536        #Setup
1537        self.domain.set_name('teg')
1538
1539        prjfile = self.domain.get_name() + '_elevation.prj'
1540        ascfile = self.domain.get_name() + '_elevation.asc'
1541        swwfile = self.domain.get_name() + '.sww'
1542
1543        self.domain.set_datadir('.')
1544        self.domain.smooth = True
1545        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1546        self.domain.set_quantity('stage', 1.0)
1547
1548        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1549
1550        sww = SWW_file(self.domain)
1551        sww.store_connectivity()
1552        sww.store_timestep()
1553        self.domain.evolve_to_end(finaltime = 0.01)
1554        sww.store_timestep()
1555
1556        cellsize = 0.25
1557        #Check contents
1558        #Get NetCDF
1559
1560        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1561
1562        # Get the variables
1563        x = fid.variables['x'][:]
1564        y = fid.variables['y'][:]
1565        z = fid.variables['elevation'][:]
1566        time = fid.variables['time'][:]
1567        stage = fid.variables['stage'][:]
1568
1569        fid.close()
1570
1571        #Export to ascii/prj files
1572        sww2dem_batch(self.domain.get_name(),
1573                quantities = 'elevation',
1574                cellsize = cellsize,
1575                verbose = self.verbose,
1576                format = 'asc')
1577
1578        #Check asc file
1579        ascid = open(ascfile)
1580        lines = ascid.readlines()
1581        ascid.close()
1582
1583        L = lines[2].strip().split()
1584        assert L[0].strip().lower() == 'xllcorner'
1585        assert num.allclose(float(L[1].strip().lower()), 308500)
1586
1587        L = lines[3].strip().split()
1588        assert L[0].strip().lower() == 'yllcorner'
1589        assert num.allclose(float(L[1].strip().lower()), 6189000)
1590
1591        #Check grid values
1592        #print '==='
1593        for j in range(5):
1594            L = lines[6+j].strip().split()
1595            y = (4-j) * cellsize
1596            for i in range(5):
1597                #print float(L[i])
1598                assert num.allclose(float(L[i]), -i*cellsize - y)
1599               
1600        #Cleanup
1601        os.remove(prjfile)
1602        os.remove(ascfile)
1603        os.remove(swwfile)
1604
1605    def test_export_gridII(self):
1606        """
1607        test_export_gridII(self):
1608        Test that sww information can be converted correctly to asc/prj
1609        format readable by e.g. ArcView
1610        """
1611
1612        import time, os
1613        from Scientific.IO.NetCDF import NetCDFFile
1614
1615        try:
1616            os.remove('teg*.sww')
1617        except:
1618            pass
1619
1620        #Setup
1621        self.domain.set_name('tegII')
1622
1623        swwfile = self.domain.get_name() + '.sww'
1624
1625        self.domain.set_datadir('.')
1626        self.domain.smooth = True
1627        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1628        self.domain.set_quantity('stage', 1.0)
1629
1630        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1631
1632        sww = SWW_file(self.domain)
1633        sww.store_connectivity()
1634        sww.store_timestep()
1635        self.domain.evolve_to_end(finaltime = 0.01)
1636        sww.store_timestep()
1637
1638        cellsize = 0.25
1639        #Check contents
1640        #Get NetCDF
1641
1642        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1643
1644        # Get the variables
1645        x = fid.variables['x'][:]
1646        y = fid.variables['y'][:]
1647        z = fid.variables['elevation'][:]
1648        time = fid.variables['time'][:]
1649        stage = fid.variables['stage'][:]
1650        xmomentum = fid.variables['xmomentum'][:]
1651        ymomentum = fid.variables['ymomentum'][:]       
1652
1653        #print 'stage', stage
1654        #print 'xmom', xmomentum
1655        #print 'ymom', ymomentum       
1656
1657        fid.close()
1658
1659        #Export to ascii/prj files
1660        if True:
1661            sww2dem_batch(self.domain.get_name(),
1662                        quantities = ['elevation', 'depth'],
1663                        cellsize = cellsize,
1664                        verbose = self.verbose,
1665                        format = 'asc')
1666
1667        else:
1668            sww2dem_batch(self.domain.get_name(),
1669                quantities = ['depth'],
1670                cellsize = cellsize,
1671                verbose = self.verbose,
1672                format = 'asc')
1673
1674
1675            export_grid(self.domain.get_name(),
1676                quantities = ['elevation'],
1677                cellsize = cellsize,
1678                verbose = self.verbose,
1679                format = 'asc')
1680
1681        prjfile = self.domain.get_name() + '_elevation.prj'
1682        ascfile = self.domain.get_name() + '_elevation.asc'
1683       
1684        #Check asc file
1685        ascid = open(ascfile)
1686        lines = ascid.readlines()
1687        ascid.close()
1688
1689        L = lines[2].strip().split()
1690        assert L[0].strip().lower() == 'xllcorner'
1691        assert num.allclose(float(L[1].strip().lower()), 308500)
1692
1693        L = lines[3].strip().split()
1694        assert L[0].strip().lower() == 'yllcorner'
1695        assert num.allclose(float(L[1].strip().lower()), 6189000)
1696
1697        #print "ascfile", ascfile
1698        #Check grid values
1699        for j in range(5):
1700            L = lines[6+j].strip().split()
1701            y = (4-j) * cellsize
1702            for i in range(5):
1703                #print " -i*cellsize - y",  -i*cellsize - y
1704                #print "float(L[i])", float(L[i])
1705                assert num.allclose(float(L[i]), -i*cellsize - y)
1706
1707        #Cleanup
1708        os.remove(prjfile)
1709        os.remove(ascfile)
1710       
1711        #Check asc file
1712        ascfile = self.domain.get_name() + '_depth.asc'
1713        prjfile = self.domain.get_name() + '_depth.prj'
1714        ascid = open(ascfile)
1715        lines = ascid.readlines()
1716        ascid.close()
1717
1718        L = lines[2].strip().split()
1719        assert L[0].strip().lower() == 'xllcorner'
1720        assert num.allclose(float(L[1].strip().lower()), 308500)
1721
1722        L = lines[3].strip().split()
1723        assert L[0].strip().lower() == 'yllcorner'
1724        assert num.allclose(float(L[1].strip().lower()), 6189000)
1725
1726        #Check grid values
1727        for j in range(5):
1728            L = lines[6+j].strip().split()
1729            y = (4-j) * cellsize
1730            for i in range(5):
1731                #print " -i*cellsize - y",  -i*cellsize - y
1732                #print "float(L[i])", float(L[i])               
1733                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1734
1735        #Cleanup
1736        os.remove(prjfile)
1737        os.remove(ascfile)
1738        os.remove(swwfile)
1739
1740
1741    def test_export_gridIII(self):
1742        """
1743        test_export_gridIII
1744        Test that sww information can be converted correctly to asc/prj
1745        format readable by e.g. ArcView
1746        """
1747
1748        import time, os
1749        from Scientific.IO.NetCDF import NetCDFFile
1750
1751        try:
1752            os.remove('teg*.sww')
1753        except:
1754            pass
1755
1756        #Setup
1757       
1758        self.domain.set_name('tegIII')
1759
1760        swwfile = self.domain.get_name() + '.sww'
1761
1762        self.domain.set_datadir('.')
1763        self.domain.format = 'sww'
1764        self.domain.smooth = True
1765        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1766        self.domain.set_quantity('stage', 1.0)
1767
1768        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1769       
1770        sww = SWW_file(self.domain)
1771        sww.store_connectivity()
1772        sww.store_timestep() #'stage')
1773        self.domain.evolve_to_end(finaltime = 0.01)
1774        sww.store_timestep() #'stage')
1775
1776        cellsize = 0.25
1777        #Check contents
1778        #Get NetCDF
1779
1780        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1781
1782        # Get the variables
1783        x = fid.variables['x'][:]
1784        y = fid.variables['y'][:]
1785        z = fid.variables['elevation'][:]
1786        time = fid.variables['time'][:]
1787        stage = fid.variables['stage'][:]
1788
1789        fid.close()
1790
1791        #Export to ascii/prj files
1792        extra_name_out = 'yeah'
1793        if True:
1794            sww2dem_batch(self.domain.get_name(),
1795                        quantities = ['elevation', 'depth'],
1796                        extra_name_out = extra_name_out,
1797                        cellsize = cellsize,
1798                        verbose = self.verbose,
1799                        format = 'asc')
1800
1801        else:
1802            sww2dem_batch(self.domain.get_name(),
1803                quantities = ['depth'],
1804                cellsize = cellsize,
1805                verbose = self.verbose,
1806                format = 'asc')
1807
1808
1809            sww2dem_batch(self.domain.get_name(),
1810                quantities = ['elevation'],
1811                cellsize = cellsize,
1812                verbose = self.verbose,
1813                format = 'asc')
1814
1815        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
1816        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
1817       
1818        #Check asc file
1819        ascid = open(ascfile)
1820        lines = ascid.readlines()
1821        ascid.close()
1822
1823        L = lines[2].strip().split()
1824        assert L[0].strip().lower() == 'xllcorner'
1825        assert num.allclose(float(L[1].strip().lower()), 308500)
1826
1827        L = lines[3].strip().split()
1828        assert L[0].strip().lower() == 'yllcorner'
1829        assert num.allclose(float(L[1].strip().lower()), 6189000)
1830
1831        #print "ascfile", ascfile
1832        #Check grid values
1833        for j in range(5):
1834            L = lines[6+j].strip().split()
1835            y = (4-j) * cellsize
1836            for i in range(5):
1837                #print " -i*cellsize - y",  -i*cellsize - y
1838                #print "float(L[i])", float(L[i])
1839                assert num.allclose(float(L[i]), -i*cellsize - y)
1840               
1841        #Cleanup
1842        os.remove(prjfile)
1843        os.remove(ascfile)
1844       
1845        #Check asc file
1846        ascfile = self.domain.get_name() + '_depth_yeah.asc'
1847        prjfile = self.domain.get_name() + '_depth_yeah.prj'
1848        ascid = open(ascfile)
1849        lines = ascid.readlines()
1850        ascid.close()
1851
1852        L = lines[2].strip().split()
1853        assert L[0].strip().lower() == 'xllcorner'
1854        assert num.allclose(float(L[1].strip().lower()), 308500)
1855
1856        L = lines[3].strip().split()
1857        assert L[0].strip().lower() == 'yllcorner'
1858        assert num.allclose(float(L[1].strip().lower()), 6189000)
1859
1860        #Check grid values
1861        for j in range(5):
1862            L = lines[6+j].strip().split()
1863            y = (4-j) * cellsize
1864            for i in range(5):
1865                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1866
1867        #Cleanup
1868        os.remove(prjfile)
1869        os.remove(ascfile)
1870        os.remove(swwfile)
1871
1872    def test_export_grid_bad(self):
1873        """Test that Bad input throws exception error
1874        """
1875
1876        try:
1877            sww2dem_batch('a_small_round-egg',
1878                        quantities = ['elevation', 'depth'],
1879                        cellsize = 99,
1880                        verbose = self.verbose,
1881                        format = 'asc')
1882        except IOError:
1883            pass
1884        else:
1885            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
1886       
1887
1888#################################################################################
1889
1890if __name__ == "__main__":
1891    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
1892    suite = unittest.makeSuite(Test_Sww2Dem, 'test_')
1893    runner = unittest.TextTestRunner(verbosity=2)
1894    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.