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

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

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.

File size: 41.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
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#################################################################################
1396
1397if __name__ == "__main__":
1398    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
1399    suite = unittest.makeSuite(Test_Sww2Dem, 'test')
1400    runner = unittest.TextTestRunner(verbosity=1)
1401    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.