source: trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py @ 7776

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

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

File size: 87.2 KB
Line 
1#!/usr/bin/env python
2#
3
4"""
5Set of tests for the now-defunct data manager module.
6
7These could be split up into their correct modules.
8"""
9
10import unittest
11import copy
12import numpy as num
13import sys
14               
15import tempfile
16import os
17import shutil
18from struct import pack, unpack
19
20from Scientific.IO.NetCDF import NetCDFFile
21
22
23from anuga.anuga_exceptions import ANUGAError
24from anuga.file.sww import SWW_file
25from anuga.coordinate_transforms.geo_reference import Geo_reference                       
26from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
27from anuga.abstract_2d_finite_volumes.util import file_function
28from anuga.utilities.system_tools import get_pathname_from_package, \
29                                            get_revision_number
30from anuga.utilities.file_utils import del_dir
31from anuga.utilities.numerical_tools import ensure_numeric, mean
32from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
33from anuga.config import netcdf_float, epsilon, g
34
35from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array
36from anuga.file.sts import create_sts_boundary
37
38
39# import all the boundaries - some are generic, some are shallow water
40from boundaries import Reflective_boundary, \
41            Field_boundary, Transmissive_momentum_set_stage_boundary, \
42            Transmissive_stage_zero_momentum_boundary
43from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
44     import Transmissive_boundary, Dirichlet_boundary, \
45            Time_boundary, File_boundary, AWI_boundary
46
47# This is needed to run the tests of local functions
48from anuga.file_conversion.urs2sts import urs2sts
49from anuga.coordinate_transforms.redfearn import redfearn
50from anuga.coordinate_transforms.geo_reference import Geo_reference, \
51     DEFAULT_ZONE
52from anuga.geospatial_data.geospatial_data import Geospatial_data
53
54from shallow_water_domain import Domain
55
56# use helper methods from other unit test
57from anuga.file.test_mux import Test_Mux
58
59
60class Test_Data_Manager(Test_Mux):
61    # Class variable
62    verbose = False
63
64    def set_verbose(self):
65        Test_Data_Manager.verbose = True
66       
67    def setUp(self):
68        import time
69        from mesh_factory import rectangular
70       
71        self.verbose = Test_Data_Manager.verbose
72        # Create basic mesh
73        points, vertices, boundary = rectangular(2, 2)
74
75        # Create shallow water domain
76        domain = Domain(points, vertices, boundary)
77        domain.default_order = 2
78
79        # Set some field values
80        domain.set_quantity('elevation', lambda x,y: -x)
81        domain.set_quantity('friction', 0.03)
82
83
84        ######################
85        # Boundary conditions
86        B = Transmissive_boundary(domain)
87        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
88
89
90        ######################
91        #Initial condition - with jumps
92        bed = domain.quantities['elevation'].vertex_values
93        stage = num.zeros(bed.shape, num.float)
94
95        h = 0.3
96        for i in range(stage.shape[0]):
97            if i % 2 == 0:
98                stage[i,:] = bed[i,:] + h
99            else:
100                stage[i,:] = bed[i,:]
101
102        domain.set_quantity('stage', stage)
103
104
105        domain.distribute_to_vertices_and_edges()               
106        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
107
108
109        self.domain = domain
110
111        C = domain.get_vertex_coordinates()
112        self.X = C[:,0:6:2].copy()
113        self.Y = C[:,1:6:2].copy()
114
115        self.F = bed
116
117        #Write A testfile (not realistic. Values aren't realistic)
118        self.test_MOST_file = 'most_small'
119
120        longitudes = [150.66667, 150.83334, 151., 151.16667]
121        latitudes = [-34.5, -34.33333, -34.16667, -34]
122
123        long_name = 'LON'
124        lat_name = 'LAT'
125
126        nx = 4
127        ny = 4
128        six = 6
129
130
131        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
132            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
133
134            fid.createDimension(long_name,nx)
135            fid.createVariable(long_name,netcdf_float,(long_name,))
136            fid.variables[long_name].point_spacing='uneven'
137            fid.variables[long_name].units='degrees_east'
138            fid.variables[long_name].assignValue(longitudes)
139
140            fid.createDimension(lat_name,ny)
141            fid.createVariable(lat_name,netcdf_float,(lat_name,))
142            fid.variables[lat_name].point_spacing='uneven'
143            fid.variables[lat_name].units='degrees_north'
144            fid.variables[lat_name].assignValue(latitudes)
145
146            fid.createDimension('TIME',six)
147            fid.createVariable('TIME',netcdf_float,('TIME',))
148            fid.variables['TIME'].point_spacing='uneven'
149            fid.variables['TIME'].units='seconds'
150            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
151
152
153            name = ext[1:3].upper()
154            if name == 'E.': name = 'ELEVATION'
155            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
156            fid.variables[name].units='CENTIMETERS'
157            fid.variables[name].missing_value=-1.e+034
158
159            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
160                                              [-0.1214216, 0, 0, 0],
161                                              [0, 0, 0, 0],
162                                              [0, 0, 0, 0]],
163                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
164                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
165                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
166                                              [0, 2.291054e-005, 2.291054e-005, 0]],
167                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
168                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
169                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
170                                              [0, 0.0001374632, 0.0001374632, 0]],
171                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
172                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
173                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
174                                              [0, 0.0002520159, 0.0002520159, 0]],
175                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
176                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
177                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
178                                              [0, 0.0003665686, 0.0003665686, 0]],
179                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
180                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
181                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
182                                              [0, 0.0004811212, 0.0004811212, 0]]])
183
184
185            fid.close()
186
187
188
189
190    def tearDown(self):
191        import os
192        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
193            #print 'Trying to remove', self.test_MOST_file + ext
194            os.remove(self.test_MOST_file + ext)
195
196    def test_sww_constant(self):
197        """Test that constant sww information can be written correctly
198        (non smooth)
199        """
200        self.domain.set_name('datatest' + str(id(self)))
201        self.domain.format = 'sww' #Remove??
202        self.domain.smooth = False
203       
204        sww = SWW_file(self.domain)
205        sww.store_connectivity()
206
207        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
208
209        # Get the variables
210        x = fid.variables['x']
211        y = fid.variables['y']
212        z = fid.variables['elevation']
213        V = fid.variables['volumes']
214
215        assert num.allclose (x[:], self.X.flatten())
216        assert num.allclose (y[:], self.Y.flatten())
217        assert num.allclose (z[:], self.F.flatten())
218
219        P = len(self.domain)
220        for k in range(P):
221            assert V[k, 0] == 3*k
222            assert V[k, 1] == 3*k+1
223            assert V[k, 2] == 3*k+2
224
225        fid.close()
226        os.remove(sww.filename)
227
228    def test_sww_header(self):
229        """Test that constant sww information can be written correctly
230        (non smooth)
231        """
232        self.domain.set_name('datatest' + str(id(self)))
233        self.domain.format = 'sww' #Remove??
234        self.domain.smooth = False
235
236        sww = SWW_file(self.domain)
237        sww.store_connectivity()
238
239        # Check contents
240        # Get NetCDF
241        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
242
243        # Get the variables
244        sww_revision = fid.revision_number
245        try:
246            revision_number = get_revision_number()
247        except:
248            revision_number = None
249           
250        assert str(revision_number) == sww_revision
251        fid.close()
252
253        #print "sww.filename", sww.filename
254        os.remove(sww.filename)
255
256    def test_sww_range(self):
257        """Test that constant sww information can be written correctly
258        Use non-smooth to be able to compare to quantity values.
259        """
260
261        self.domain.set_name('datatest' + str(id(self)))
262        self.domain.format = 'sww'
263        self.domain.set_store_vertices_uniquely(True)
264       
265        sww = SWW_file(self.domain)       
266
267        dqs = self.domain.get_quantity('stage')
268        dqx = self.domain.get_quantity('xmomentum')
269        dqy = self.domain.get_quantity('ymomentum')       
270        xmom_min = ymom_min = stage_min = sys.maxint
271        xmom_max = ymom_max = stage_max = -stage_min       
272        for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
273            wmax = max(dqs.get_values().flatten())
274            if wmax > stage_max: stage_max = wmax
275            wmin = min(dqs.get_values().flatten())
276            if wmin < stage_min: stage_min = wmin           
277           
278            uhmax = max(dqx.get_values().flatten())
279            if uhmax > xmom_max: xmom_max = uhmax
280            uhmin = min(dqx.get_values().flatten())
281            if uhmin < xmom_min: xmom_min = uhmin                       
282           
283            vhmax = max(dqy.get_values().flatten())
284            if vhmax > ymom_max: ymom_max = vhmax
285            vhmin = min(dqy.get_values().flatten())
286            if vhmin < ymom_min: ymom_min = vhmin                                   
287           
288           
289           
290        # Get NetCDF
291        fid = NetCDFFile(sww.filename, netcdf_mode_r) # Open existing file for append
292
293        # Get the variables
294        range = fid.variables['stage_range'][:]
295        assert num.allclose(range,[stage_min, stage_max])
296
297        range = fid.variables['xmomentum_range'][:]
298        #print range
299        assert num.allclose(range, [xmom_min, xmom_max])
300       
301        range = fid.variables['ymomentum_range'][:]
302        #print range
303        assert num.allclose(range, [ymom_min, ymom_max])       
304
305
306       
307        fid.close()
308        os.remove(sww.filename)
309
310    def test_sww_extrema(self):
311        """Test that extrema of quantities can be retrieved at every vertex
312        Extrema are updated at every *internal* timestep
313        """
314
315        domain = self.domain
316       
317        domain.set_name('extrema_test' + str(id(self)))
318        domain.format = 'sww'
319        domain.smooth = True
320
321        assert domain.quantities_to_be_monitored is None
322        assert domain.monitor_polygon is None
323        assert domain.monitor_time_interval is None       
324       
325        domain.set_quantities_to_be_monitored(['xmomentum',
326                                               'ymomentum',
327                                               'stage-elevation'])
328
329        assert domain.monitor_polygon is None
330        assert domain.monitor_time_interval is None
331
332
333        domain.set_quantities_to_be_monitored(['xmomentum',
334                                               'ymomentum',
335                                               'stage-elevation'],
336                                              polygon=domain.get_boundary_polygon(),
337                                              time_interval=[0,1])
338       
339       
340        assert len(domain.quantities_to_be_monitored) == 3
341        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
342        assert domain.quantities_to_be_monitored.has_key('xmomentum')               
343        assert domain.quantities_to_be_monitored.has_key('ymomentum')       
344
345       
346        #domain.protect_against_isolated_degenerate_timesteps = True
347        #domain.tight_slope_limiters = 1
348        domain.tight_slope_limiters = 0 # Backwards compatibility
349        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
350       
351       
352        sww = SWW_file(domain)
353
354        for t in domain.evolve(yieldstep = 1, finaltime = 1):
355            pass
356            #print domain.timestepping_statistics()
357            domain.quantity_statistics(precision = '%.8f') # Silent
358
359           
360        # Get NetCDF
361        fid = NetCDFFile(sww.filename, netcdf_mode_r) # Open existing file for append
362
363        # Get the variables
364        extrema = fid.variables['stage-elevation.extrema'][:]
365        assert num.allclose(extrema, [0.00, 0.30])
366
367        loc = fid.variables['stage-elevation.min_location'][:]
368        assert num.allclose(loc, [0.16666667, 0.33333333])
369
370        loc = fid.variables['stage-elevation.max_location'][:]       
371        assert num.allclose(loc, [0.8333333, 0.16666667])       
372
373        time = fid.variables['stage-elevation.max_time'][:]
374        assert num.allclose(time, 0.0)               
375
376        extrema = fid.variables['xmomentum.extrema'][:]
377        assert num.allclose(extrema,[-0.06062178, 0.47873023]) or \
378            num.allclose(extrema, [-0.06062178, 0.47847986]) or \
379            num.allclose(extrema, [-0.06062178, 0.47848481]) or \
380            num.allclose(extrema, [-0.06062178, 0.47763887]) # 18/09/09
381       
382        extrema = fid.variables['ymomentum.extrema'][:]
383        assert num.allclose(extrema,[0.00, 0.0625786]) or num.allclose(extrema,[0.00, 0.06062178])
384
385        time_interval = fid.variables['extrema.time_interval'][:]
386        assert num.allclose(time_interval, [0,1])
387       
388        polygon = fid.variables['extrema.polygon'][:]       
389        assert num.allclose(polygon, domain.get_boundary_polygon())
390       
391        fid.close()
392        #print "sww.filename", sww.filename
393        os.remove(sww.filename)
394
395       
396       
397    def test_sww_constant_smooth(self):
398        """Test that constant sww information can be written correctly
399        (non smooth)
400        """
401        self.domain.set_name('datatest' + str(id(self)))
402        self.domain.format = 'sww'
403        self.domain.smooth = True
404
405        sww = SWW_file(self.domain)
406        sww.store_connectivity()
407
408        # Check contents
409        # Get NetCDF
410        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
411
412        # Get the variables
413        X = fid.variables['x'][:]
414        Y = fid.variables['y'][:]
415        Z = fid.variables['elevation'][:]
416        V = fid.variables['volumes']
417
418        assert num.allclose([X[0], Y[0]], num.array([0.0, 0.0]))
419        assert num.allclose([X[1], Y[1]], num.array([0.0, 0.5]))
420        assert num.allclose([X[2], Y[2]], num.array([0.0, 1.0]))
421        assert num.allclose([X[4], Y[4]], num.array([0.5, 0.5]))
422        assert num.allclose([X[7], Y[7]], num.array([1.0, 0.5]))
423
424        assert Z[4] == -0.5
425
426        assert V[2,0] == 4
427        assert V[2,1] == 5
428        assert V[2,2] == 1
429        assert V[4,0] == 6
430        assert V[4,1] == 7
431        assert V[4,2] == 3
432
433        fid.close()
434        os.remove(sww.filename)
435       
436
437    def test_sww_variable(self):
438        """Test that sww information can be written correctly
439        """
440        self.domain.set_name('datatest' + str(id(self)))
441        self.domain.format = 'sww'
442        self.domain.smooth = True
443        self.domain.reduction = mean
444
445        sww = SWW_file(self.domain)
446        sww.store_connectivity()
447        sww.store_timestep()
448
449        # Check contents
450        # Get NetCDF
451        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
452
453
454        # Get the variables
455        time = fid.variables['time']
456        stage = fid.variables['stage']
457
458        Q = self.domain.quantities['stage']
459        Q0 = Q.vertex_values[:,0]
460        Q1 = Q.vertex_values[:,1]
461        Q2 = Q.vertex_values[:,2]
462
463        A = stage[0,:]
464        #print A[0], (Q2[0,0] + Q1[1,0])/2
465        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
466        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
467        assert num.allclose(A[2], Q0[3])
468        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
469
470        #Center point
471        assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
472                                   Q0[5] + Q2[6] + Q1[7])/6)
473       
474        fid.close()
475        os.remove(sww.filename)
476
477
478    def test_sww_variable2(self):
479        """Test that sww information can be written correctly
480        multiple timesteps. Use average as reduction operator
481        """
482
483        import time, os
484        from Scientific.IO.NetCDF import NetCDFFile
485
486        self.domain.set_name('datatest' + str(id(self)))
487        self.domain.format = 'sww'
488        self.domain.smooth = True
489
490        self.domain.reduction = mean
491
492        sww = SWW_file(self.domain)
493        sww.store_connectivity()
494        sww.store_timestep()
495        #self.domain.tight_slope_limiters = 1
496        self.domain.evolve_to_end(finaltime = 0.01)
497        sww.store_timestep()
498
499
500        # Check contents
501        # Get NetCDF
502        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
503
504        # Get the variables
505        x = fid.variables['x']
506        y = fid.variables['y']
507        z = fid.variables['elevation']
508        time = fid.variables['time']
509        stage = fid.variables['stage']
510
511        #Check values
512        Q = self.domain.quantities['stage']
513        Q0 = Q.vertex_values[:,0]
514        Q1 = Q.vertex_values[:,1]
515        Q2 = Q.vertex_values[:,2]
516
517        A = stage[1,:]
518        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
519        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
520        assert num.allclose(A[2], Q0[3])
521        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
522
523        #Center point
524        assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
525                                   Q0[5] + Q2[6] + Q1[7])/6)
526
527
528        fid.close()
529
530        #Cleanup
531        os.remove(sww.filename)
532
533    def no_test_sww_variable3(self):
534        """Test that sww information can be written correctly
535        multiple timesteps using a different reduction operator (min)
536        """
537
538        # Different reduction in sww files has been made obsolete.
539       
540        import time, os
541        from Scientific.IO.NetCDF import NetCDFFile
542
543        self.domain.set_name('datatest' + str(id(self)))
544        self.domain.format = 'sww'
545        self.domain.smooth = True
546        self.domain.reduction = min
547
548        sww = SWW_file(self.domain)
549        sww.store_connectivity()
550        sww.store_timestep()
551        #self.domain.tight_slope_limiters = 1
552        self.domain.evolve_to_end(finaltime = 0.01)
553        sww.store_timestep()
554
555
556        #Check contents
557        #Get NetCDF
558        fid = NetCDFFile(sww.filename, netcdf_mode_r)
559
560        # Get the variables
561        x = fid.variables['x']
562        y = fid.variables['y']
563        z = fid.variables['elevation']
564        time = fid.variables['time']
565        stage = fid.variables['stage']
566
567        #Check values
568        Q = self.domain.quantities['stage']
569        Q0 = Q.vertex_values[:,0]
570        Q1 = Q.vertex_values[:,1]
571        Q2 = Q.vertex_values[:,2]
572
573        A = stage[1,:]
574        assert num.allclose(A[0], min(Q2[0], Q1[1]))
575        assert num.allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
576        assert num.allclose(A[2], Q0[3])
577        assert num.allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
578
579        #Center point
580        assert num.allclose(A[4], min(Q1[0], Q2[1], Q0[2],
581                                      Q0[5], Q2[6], Q1[7]))
582
583
584        fid.close()
585
586        #Cleanup
587        os.remove(sww.filename)
588
589
590    def test_sync(self):
591        """test_sync - Test info stored at each timestep is as expected (incl initial condition)
592        """
593
594        import time, os, config
595        from Scientific.IO.NetCDF import NetCDFFile
596
597        self.domain.set_name('synctest')
598        self.domain.format = 'sww'
599        self.domain.smooth = False
600        self.domain.store = True
601
602        self.domain.tight_slope_limiters = True
603        self.domain.use_centroid_velocities = True       
604       
605        # In this case tight_slope_limiters as default
606        # in conjunction with protection
607        # against isolated degenerate timesteps works.
608        #self.domain.tight_slope_limiters = 1
609        #self.domain.protect_against_isolated_degenerate_timesteps = True
610
611        #print 'tight_sl', self.domain.tight_slope_limiters
612       
613
614        #Evolution
615        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
616           
617            #########self.domain.write_time(track_speeds=True)
618            stage = self.domain.quantities['stage'].vertex_values
619
620            #Get NetCDF
621            fid = NetCDFFile(self.domain.writer.filename, netcdf_mode_r)
622            stage_file = fid.variables['stage']
623           
624            if t == 0.0:
625                assert num.allclose(stage, self.initial_stage)
626                assert num.allclose(stage_file[:], stage.flatten())
627            else:
628                assert not num.allclose(stage, self.initial_stage)
629                assert not num.allclose(stage_file[:], stage.flatten())
630
631            fid.close()
632
633        os.remove(self.domain.writer.filename)
634
635
636    def test_sww_minimum_storable_height(self):
637        """Test that sww information can be written correctly
638        multiple timesteps using a different reduction operator (min)
639        """
640
641        import time, os
642        from Scientific.IO.NetCDF import NetCDFFile
643
644        self.domain.set_name('datatest' + str(id(self)))
645        self.domain.format = 'sww'
646        self.domain.smooth = True
647        self.domain.reduction = min
648        self.domain.minimum_storable_height = 100
649
650        sww = SWW_file(self.domain)
651        sww.store_connectivity()
652        sww.store_timestep()
653
654        #self.domain.tight_slope_limiters = 1
655        self.domain.evolve_to_end(finaltime = 0.01)
656        sww.store_timestep()
657
658
659        #Check contents
660        #Get NetCDF
661        fid = NetCDFFile(sww.filename, netcdf_mode_r)
662
663
664        # Get the variables
665        x = fid.variables['x']
666        y = fid.variables['y']
667        z = fid.variables['elevation']
668        time = fid.variables['time']
669        stage = fid.variables['stage']
670        xmomentum = fid.variables['xmomentum']
671        ymomentum = fid.variables['ymomentum']       
672
673        #Check values
674        Q = self.domain.quantities['stage']
675        Q0 = Q.vertex_values[:,0]
676        Q1 = Q.vertex_values[:,1]
677        Q2 = Q.vertex_values[:,2]
678
679        A = stage[1,:]
680        assert num.allclose(stage[1,:], z[:])
681
682
683        assert num.allclose(xmomentum, 0.0)
684        assert num.allclose(ymomentum, 0.0)       
685       
686        fid.close()
687
688        #Cleanup
689        os.remove(sww.filename)
690
691
692    def Not_a_test_sww_DSG(self):
693        """Not a test, rather a look at the sww format
694        """
695
696        import time, os
697        from Scientific.IO.NetCDF import NetCDFFile
698
699        self.domain.set_name('datatest' + str(id(self)))
700        self.domain.format = 'sww'
701        self.domain.smooth = True
702        self.domain.reduction = mean
703
704        sww = SWW_file(self.domain)
705        sww.store_connectivity()
706        sww.store_timestep()
707
708        #Check contents
709        #Get NetCDF
710        fid = NetCDFFile(sww.filename, netcdf_mode_r)
711
712        # Get the variables
713        x = fid.variables['x']
714        y = fid.variables['y']
715        z = fid.variables['elevation']
716
717        volumes = fid.variables['volumes']
718        time = fid.variables['time']
719
720        # 2D
721        stage = fid.variables['stage']
722
723        X = x[:]
724        Y = y[:]
725        Z = z[:]
726        V = volumes[:]
727        T = time[:]
728        S = stage[:,:]
729
730#         print "****************************"
731#         print "X ",X
732#         print "****************************"
733#         print "Y ",Y
734#         print "****************************"
735#         print "Z ",Z
736#         print "****************************"
737#         print "V ",V
738#         print "****************************"
739#         print "Time ",T
740#         print "****************************"
741#         print "Stage ",S
742#         print "****************************"
743
744
745        fid.close()
746
747        #Cleanup
748        os.remove(sww.filename)
749
750
751
752    def test_export_grid(self):
753        """
754        test_export_grid(self):
755        Test that sww information can be converted correctly to asc/prj
756        format readable by e.g. ArcView
757        """
758
759        import time, os
760        from Scientific.IO.NetCDF import NetCDFFile
761
762        try:
763            os.remove('teg*.sww')
764        except:
765            pass
766
767        #Setup
768        self.domain.set_name('teg')
769
770        prjfile = self.domain.get_name() + '_elevation.prj'
771        ascfile = self.domain.get_name() + '_elevation.asc'
772        swwfile = self.domain.get_name() + '.sww'
773
774        self.domain.set_datadir('.')
775        self.domain.smooth = True
776        self.domain.set_quantity('elevation', lambda x,y: -x-y)
777        self.domain.set_quantity('stage', 1.0)
778
779        self.domain.geo_reference = Geo_reference(56,308500,6189000)
780
781        sww = SWW_file(self.domain)
782        sww.store_connectivity()
783        sww.store_timestep()
784        self.domain.evolve_to_end(finaltime = 0.01)
785        sww.store_timestep()
786
787        cellsize = 0.25
788        #Check contents
789        #Get NetCDF
790
791        fid = NetCDFFile(sww.filename, netcdf_mode_r)
792
793        # Get the variables
794        x = fid.variables['x'][:]
795        y = fid.variables['y'][:]
796        z = fid.variables['elevation'][:]
797        time = fid.variables['time'][:]
798        stage = fid.variables['stage'][:]
799
800        fid.close()
801
802        #Export to ascii/prj files
803        export_grid(self.domain.get_name(),
804                quantities = 'elevation',
805                cellsize = cellsize,
806                verbose = self.verbose,
807                format = 'asc')
808
809        #Check asc file
810        ascid = open(ascfile)
811        lines = ascid.readlines()
812        ascid.close()
813
814        L = lines[2].strip().split()
815        assert L[0].strip().lower() == 'xllcorner'
816        assert num.allclose(float(L[1].strip().lower()), 308500)
817
818        L = lines[3].strip().split()
819        assert L[0].strip().lower() == 'yllcorner'
820        assert num.allclose(float(L[1].strip().lower()), 6189000)
821
822        #Check grid values
823        for j in range(5):
824            L = lines[6+j].strip().split()
825            y = (4-j) * cellsize
826            for i in range(5):
827                assert num.allclose(float(L[i]), -i*cellsize - y)
828               
829        #Cleanup
830        os.remove(prjfile)
831        os.remove(ascfile)
832        os.remove(swwfile)
833
834    def test_export_gridII(self):
835        """
836        test_export_gridII(self):
837        Test that sww information can be converted correctly to asc/prj
838        format readable by e.g. ArcView
839        """
840
841        import time, os
842        from Scientific.IO.NetCDF import NetCDFFile
843
844        try:
845            os.remove('teg*.sww')
846        except:
847            pass
848
849        #Setup
850        self.domain.set_name('tegII')
851
852        swwfile = self.domain.get_name() + '.sww'
853
854        self.domain.set_datadir('.')
855        self.domain.smooth = True
856        self.domain.set_quantity('elevation', lambda x,y: -x-y)
857        self.domain.set_quantity('stage', 1.0)
858
859        self.domain.geo_reference = Geo_reference(56,308500,6189000)
860
861        sww = SWW_file(self.domain)
862        sww.store_connectivity()
863        sww.store_timestep()
864        self.domain.evolve_to_end(finaltime = 0.01)
865        sww.store_timestep()
866
867        cellsize = 0.25
868        #Check contents
869        #Get NetCDF
870
871        fid = NetCDFFile(sww.filename, netcdf_mode_r)
872
873        # Get the variables
874        x = fid.variables['x'][:]
875        y = fid.variables['y'][:]
876        z = fid.variables['elevation'][:]
877        time = fid.variables['time'][:]
878        stage = fid.variables['stage'][:]
879        xmomentum = fid.variables['xmomentum'][:]
880        ymomentum = fid.variables['ymomentum'][:]       
881
882        #print 'stage', stage
883        #print 'xmom', xmomentum
884        #print 'ymom', ymomentum       
885
886        fid.close()
887
888        #Export to ascii/prj files
889        if True:
890            export_grid(self.domain.get_name(),
891                        quantities = ['elevation', 'depth'],
892                        cellsize = cellsize,
893                        verbose = self.verbose,
894                        format = 'asc')
895
896        else:
897            export_grid(self.domain.get_name(),
898                quantities = ['depth'],
899                cellsize = cellsize,
900                verbose = self.verbose,
901                format = 'asc')
902
903
904            export_grid(self.domain.get_name(),
905                quantities = ['elevation'],
906                cellsize = cellsize,
907                verbose = self.verbose,
908                format = 'asc')
909
910        prjfile = self.domain.get_name() + '_elevation.prj'
911        ascfile = self.domain.get_name() + '_elevation.asc'
912       
913        #Check asc file
914        ascid = open(ascfile)
915        lines = ascid.readlines()
916        ascid.close()
917
918        L = lines[2].strip().split()
919        assert L[0].strip().lower() == 'xllcorner'
920        assert num.allclose(float(L[1].strip().lower()), 308500)
921
922        L = lines[3].strip().split()
923        assert L[0].strip().lower() == 'yllcorner'
924        assert num.allclose(float(L[1].strip().lower()), 6189000)
925
926        #print "ascfile", ascfile
927        #Check grid values
928        for j in range(5):
929            L = lines[6+j].strip().split()
930            y = (4-j) * cellsize
931            for i in range(5):
932                #print " -i*cellsize - y",  -i*cellsize - y
933                #print "float(L[i])", float(L[i])
934                assert num.allclose(float(L[i]), -i*cellsize - y)
935
936        #Cleanup
937        os.remove(prjfile)
938        os.remove(ascfile)
939       
940        #Check asc file
941        ascfile = self.domain.get_name() + '_depth.asc'
942        prjfile = self.domain.get_name() + '_depth.prj'
943        ascid = open(ascfile)
944        lines = ascid.readlines()
945        ascid.close()
946
947        L = lines[2].strip().split()
948        assert L[0].strip().lower() == 'xllcorner'
949        assert num.allclose(float(L[1].strip().lower()), 308500)
950
951        L = lines[3].strip().split()
952        assert L[0].strip().lower() == 'yllcorner'
953        assert num.allclose(float(L[1].strip().lower()), 6189000)
954
955        #Check grid values
956        for j in range(5):
957            L = lines[6+j].strip().split()
958            y = (4-j) * cellsize
959            for i in range(5):
960                #print " -i*cellsize - y",  -i*cellsize - y
961                #print "float(L[i])", float(L[i])               
962                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
963
964        #Cleanup
965        os.remove(prjfile)
966        os.remove(ascfile)
967        os.remove(swwfile)
968
969
970    def test_export_gridIII(self):
971        """
972        test_export_gridIII
973        Test that sww information can be converted correctly to asc/prj
974        format readable by e.g. ArcView
975        """
976
977        import time, os
978        from Scientific.IO.NetCDF import NetCDFFile
979
980        try:
981            os.remove('teg*.sww')
982        except:
983            pass
984
985        #Setup
986       
987        self.domain.set_name('tegIII')
988
989        swwfile = self.domain.get_name() + '.sww'
990
991        self.domain.set_datadir('.')
992        self.domain.format = 'sww'
993        self.domain.smooth = True
994        self.domain.set_quantity('elevation', lambda x,y: -x-y)
995        self.domain.set_quantity('stage', 1.0)
996
997        self.domain.geo_reference = Geo_reference(56,308500,6189000)
998       
999        sww = SWW_file(self.domain)
1000        sww.store_connectivity()
1001        sww.store_timestep() #'stage')
1002        self.domain.evolve_to_end(finaltime = 0.01)
1003        sww.store_timestep() #'stage')
1004
1005        cellsize = 0.25
1006        #Check contents
1007        #Get NetCDF
1008
1009        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1010
1011        # Get the variables
1012        x = fid.variables['x'][:]
1013        y = fid.variables['y'][:]
1014        z = fid.variables['elevation'][:]
1015        time = fid.variables['time'][:]
1016        stage = fid.variables['stage'][:]
1017
1018        fid.close()
1019
1020        #Export to ascii/prj files
1021        extra_name_out = 'yeah'
1022        if True:
1023            export_grid(self.domain.get_name(),
1024                        quantities = ['elevation', 'depth'],
1025                        extra_name_out = extra_name_out,
1026                        cellsize = cellsize,
1027                        verbose = self.verbose,
1028                        format = 'asc')
1029
1030        else:
1031            export_grid(self.domain.get_name(),
1032                quantities = ['depth'],
1033                cellsize = cellsize,
1034                verbose = self.verbose,
1035                format = 'asc')
1036
1037
1038            export_grid(self.domain.get_name(),
1039                quantities = ['elevation'],
1040                cellsize = cellsize,
1041                verbose = self.verbose,
1042                format = 'asc')
1043
1044        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
1045        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
1046       
1047        #Check asc file
1048        ascid = open(ascfile)
1049        lines = ascid.readlines()
1050        ascid.close()
1051
1052        L = lines[2].strip().split()
1053        assert L[0].strip().lower() == 'xllcorner'
1054        assert num.allclose(float(L[1].strip().lower()), 308500)
1055
1056        L = lines[3].strip().split()
1057        assert L[0].strip().lower() == 'yllcorner'
1058        assert num.allclose(float(L[1].strip().lower()), 6189000)
1059
1060        #print "ascfile", ascfile
1061        #Check grid values
1062        for j in range(5):
1063            L = lines[6+j].strip().split()
1064            y = (4-j) * cellsize
1065            for i in range(5):
1066                #print " -i*cellsize - y",  -i*cellsize - y
1067                #print "float(L[i])", float(L[i])
1068                assert num.allclose(float(L[i]), -i*cellsize - y)
1069               
1070        #Cleanup
1071        os.remove(prjfile)
1072        os.remove(ascfile)
1073       
1074        #Check asc file
1075        ascfile = self.domain.get_name() + '_depth_yeah.asc'
1076        prjfile = self.domain.get_name() + '_depth_yeah.prj'
1077        ascid = open(ascfile)
1078        lines = ascid.readlines()
1079        ascid.close()
1080
1081        L = lines[2].strip().split()
1082        assert L[0].strip().lower() == 'xllcorner'
1083        assert num.allclose(float(L[1].strip().lower()), 308500)
1084
1085        L = lines[3].strip().split()
1086        assert L[0].strip().lower() == 'yllcorner'
1087        assert num.allclose(float(L[1].strip().lower()), 6189000)
1088
1089        #Check grid values
1090        for j in range(5):
1091            L = lines[6+j].strip().split()
1092            y = (4-j) * cellsize
1093            for i in range(5):
1094                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1095
1096        #Cleanup
1097        os.remove(prjfile)
1098        os.remove(ascfile)
1099        os.remove(swwfile)
1100
1101    def test_export_grid_bad(self):
1102        """Test that sww information can be converted correctly to asc/prj
1103        format readable by e.g. ArcView
1104        """
1105
1106        try:
1107            export_grid('a_small_round-egg',
1108                        quantities = ['elevation', 'depth'],
1109                        cellsize = 99,
1110                        verbose = self.verbose,
1111                        format = 'asc')
1112        except IOError:
1113            pass
1114        else:
1115            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
1116
1117    def test_export_grid_parallel(self):
1118        """Test that sww information can be converted correctly to asc/prj
1119        format readable by e.g. ArcView
1120        """
1121
1122        import time, os
1123        from Scientific.IO.NetCDF import NetCDFFile
1124
1125        base_name = 'tegp'
1126        #Setup
1127        self.domain.set_name(base_name+'_P0_8')
1128        swwfile = self.domain.get_name() + '.sww'
1129
1130        self.domain.set_datadir('.')
1131        self.domain.format = 'sww'
1132        self.domain.smooth = True
1133        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1134        self.domain.set_quantity('stage', 1.0)
1135
1136        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1137
1138        sww = SWW_file(self.domain)
1139        sww.store_connectivity()
1140        sww.store_timestep()
1141        self.domain.evolve_to_end(finaltime = 0.0001)
1142        #Setup
1143        self.domain.set_name(base_name+'_P1_8')
1144        swwfile2 = self.domain.get_name() + '.sww'
1145        sww = SWW_file(self.domain)
1146        sww.store_connectivity()
1147        sww.store_timestep()
1148        self.domain.evolve_to_end(finaltime = 0.0002)
1149        sww.store_timestep()
1150
1151        cellsize = 0.25
1152        #Check contents
1153        #Get NetCDF
1154
1155        fid = NetCDFFile(sww.filename, netcdf_mode_r)
1156
1157        # Get the variables
1158        x = fid.variables['x'][:]
1159        y = fid.variables['y'][:]
1160        z = fid.variables['elevation'][:]
1161        time = fid.variables['time'][:]
1162        stage = fid.variables['stage'][:]
1163
1164        fid.close()
1165
1166        #Export to ascii/prj files
1167        extra_name_out = 'yeah'
1168        export_grid(base_name,
1169                    quantities = ['elevation', 'depth'],
1170                    extra_name_out = extra_name_out,
1171                    cellsize = cellsize,
1172                    verbose = self.verbose,
1173                    format = 'asc')
1174
1175        prjfile = base_name + '_P0_8_elevation_yeah.prj'
1176        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
1177        #Check asc file
1178        ascid = open(ascfile)
1179        lines = ascid.readlines()
1180        ascid.close()
1181        #Check grid values
1182        for j in range(5):
1183            L = lines[6+j].strip().split()
1184            y = (4-j) * cellsize
1185            for i in range(5):
1186                #print " -i*cellsize - y",  -i*cellsize - y
1187                #print "float(L[i])", float(L[i])
1188                assert num.allclose(float(L[i]), -i*cellsize - y)               
1189        #Cleanup
1190        os.remove(prjfile)
1191        os.remove(ascfile)
1192
1193        prjfile = base_name + '_P1_8_elevation_yeah.prj'
1194        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
1195        #Check asc file
1196        ascid = open(ascfile)
1197        lines = ascid.readlines()
1198        ascid.close()
1199        #Check grid values
1200        for j in range(5):
1201            L = lines[6+j].strip().split()
1202            y = (4-j) * cellsize
1203            for i in range(5):
1204                #print " -i*cellsize - y",  -i*cellsize - y
1205                #print "float(L[i])", float(L[i])
1206                assert num.allclose(float(L[i]), -i*cellsize - y)               
1207        #Cleanup
1208        os.remove(prjfile)
1209        os.remove(ascfile)
1210        os.remove(swwfile)
1211
1212        #Check asc file
1213        ascfile = base_name + '_P0_8_depth_yeah.asc'
1214        prjfile = base_name + '_P0_8_depth_yeah.prj'
1215        ascid = open(ascfile)
1216        lines = ascid.readlines()
1217        ascid.close()
1218        #Check grid values
1219        for j in range(5):
1220            L = lines[6+j].strip().split()
1221            y = (4-j) * cellsize
1222            for i in range(5):
1223                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1224        #Cleanup
1225        os.remove(prjfile)
1226        os.remove(ascfile)
1227
1228        #Check asc file
1229        ascfile = base_name + '_P1_8_depth_yeah.asc'
1230        prjfile = base_name + '_P1_8_depth_yeah.prj'
1231        ascid = open(ascfile)
1232        lines = ascid.readlines()
1233        ascid.close()
1234        #Check grid values
1235        for j in range(5):
1236            L = lines[6+j].strip().split()
1237            y = (4-j) * cellsize
1238            for i in range(5):
1239                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
1240        #Cleanup
1241        os.remove(prjfile)
1242        os.remove(ascfile)
1243        os.remove(swwfile2)
1244
1245
1246
1247    def DISABLEDtest_sww2domain2(self):
1248        ##################################################################
1249        #Same as previous test, but this checks how NaNs are handled.
1250        ##################################################################
1251
1252        #FIXME: See ticket 223
1253
1254        from mesh_factory import rectangular
1255
1256        #Create basic mesh
1257        points, vertices, boundary = rectangular(2,2)
1258
1259        #Create shallow water domain
1260        domain = Domain(points, vertices, boundary)
1261        domain.smooth = False
1262        domain.store = True
1263        domain.set_name('test_file')
1264        domain.set_datadir('.')
1265        domain.default_order=2
1266
1267        domain.set_quantity('elevation', lambda x,y: -x/3)
1268        domain.set_quantity('friction', 0.1)
1269
1270        from math import sin, pi
1271        Br = Reflective_boundary(domain)
1272        Bt = Transmissive_boundary(domain)
1273        Bd = Dirichlet_boundary([0.2,0.,0.])
1274        Bw = Time_boundary(domain=domain,
1275                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
1276
1277        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1278
1279        h = 0.05
1280        elevation = domain.quantities['elevation'].vertex_values
1281        domain.set_quantity('stage', elevation + h)
1282
1283        domain.check_integrity()
1284
1285        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
1286            pass
1287            #domain.write_time()
1288
1289
1290
1291        ##################################
1292        #Import the file as a new domain
1293        ##################################
1294        from data_manager import load_sww_as_domain
1295        import os
1296
1297        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
1298
1299        # Fail because NaNs are present
1300        #domain2 = sww2domain(filename,
1301        #                     boundary,
1302        #                     fail_if_NaN=True,
1303        #                     verbose=self.verbose)       
1304        try:
1305            domain2 = load_sww_as_domain(filename,
1306                                 boundary,
1307                                 fail_if_NaN=True,
1308                                 verbose=self.verbose)
1309        except DataDomainError:
1310            # Now import it, filling NaNs to be -9999
1311            filler = -9999
1312            domain2 = load_sww_as_domain(filename,
1313                                 None,
1314                                 fail_if_NaN=False,
1315                                 NaN_filler=filler,
1316                                 verbose=self.verbose)
1317        else:
1318            raise Exception, 'should have failed'
1319
1320           
1321        # Now import it, filling NaNs to be 0
1322        filler = -9999
1323        domain2 = load_sww_as_domain(filename,
1324                             None,
1325                             fail_if_NaN=False,
1326                             NaN_filler=filler,
1327                             verbose=self.verbose)           
1328                             
1329        import sys; sys.exit() 
1330           
1331        # Clean up
1332        os.remove(filename)
1333
1334
1335        bits = ['geo_reference.get_xllcorner()',
1336                'geo_reference.get_yllcorner()',
1337                'vertex_coordinates']
1338
1339        for quantity in domain.quantities_to_be_stored:
1340            bits.append('get_quantity("%s").get_integral()' %quantity)
1341            bits.append('get_quantity("%s").get_values()' %quantity)
1342
1343        for bit in bits:
1344        #    print 'testing that domain.'+bit+' has been restored'
1345            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
1346
1347        print 
1348        print
1349        print domain2.get_quantity('xmomentum').get_values()
1350        print
1351        print domain2.get_quantity('stage').get_values()
1352        print
1353             
1354        print 'filler', filler
1355        print max(domain2.get_quantity('xmomentum').get_values().flat)
1356       
1357        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
1358        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
1359        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
1360        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
1361
1362
1363
1364    #FIXME This fails - smooth makes the comparism too hard for allclose
1365    def ztest_sww2domain3(self):
1366        ################################################
1367        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
1368        ################################################
1369        from mesh_factory import rectangular
1370        #Create basic mesh
1371
1372        yiel=0.01
1373        points, vertices, boundary = rectangular(10,10)
1374
1375        #Create shallow water domain
1376        domain = Domain(points, vertices, boundary)
1377        domain.geo_reference = Geo_reference(56,11,11)
1378        domain.smooth = True
1379        domain.store = True
1380        domain.set_name('bedslope')
1381        domain.default_order=2
1382        #Bed-slope and friction
1383        domain.set_quantity('elevation', lambda x,y: -x/3)
1384        domain.set_quantity('friction', 0.1)
1385        # Boundary conditions
1386        from math import sin, pi
1387        Br = Reflective_boundary(domain)
1388        Bt = Transmissive_boundary(domain)
1389        Bd = Dirichlet_boundary([0.2,0.,0.])
1390        Bw = Time_boundary(domain=domain,
1391                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
1392
1393        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
1394
1395        domain.quantities_to_be_stored['xmomentum'] = 2
1396        domain.quantities_to_be_stored['ymomentum'] = 2       
1397        #Initial condition
1398        h = 0.05
1399        elevation = domain.quantities['elevation'].vertex_values
1400        domain.set_quantity('stage', elevation + h)
1401
1402
1403        domain.check_integrity()
1404        #Evolution
1405        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
1406        #    domain.write_time()
1407            pass
1408
1409
1410        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
1411        domain2 = load_sww_as_domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
1412        #points, vertices, boundary = rectangular(15,15)
1413        #domain2.boundary = boundary
1414        ###################
1415        ##NOW TEST IT!!!
1416        ###################
1417
1418        os.remove(domain.get_name() + '.sww')
1419
1420        #FIXME smooth domain so that they can be compared
1421
1422
1423        bits = []
1424        for quantity in domain.quantities_to_be_stored:
1425            bits.append('quantities["%s"].get_integral()'%quantity)
1426
1427
1428        for bit in bits:
1429            #print 'testing that domain.'+bit+' has been restored'
1430            #print bit
1431            #print 'done'
1432            #print ('domain.'+bit), eval('domain.'+bit)
1433            #print ('domain2.'+bit), eval('domain2.'+bit)
1434            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
1435            pass
1436
1437        ######################################
1438        #Now evolve them both, just to be sure
1439        ######################################x
1440        domain.time = 0.
1441        from time import sleep
1442
1443        final = .5
1444        domain.set_quantity('friction', 0.1)
1445        domain.store = False
1446        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
1447
1448        for t in domain.evolve(yieldstep = yiel, finaltime = final):
1449            #domain.write_time()
1450            pass
1451
1452        domain2.smooth = True
1453        domain2.store = False
1454        domain2.default_order=2
1455        domain2.set_quantity('friction', 0.1)
1456        #Bed-slope and friction
1457        # Boundary conditions
1458        Bd2=Dirichlet_boundary([0.2,0.,0.])
1459        Br2 = Reflective_boundary(domain2)
1460        domain2.boundary = domain.boundary
1461        #print 'domain2.boundary'
1462        #print domain2.boundary
1463        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
1464        #domain2.boundary = domain.boundary
1465        #domain2.set_boundary({'exterior': Bd})
1466
1467        domain2.check_integrity()
1468
1469        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
1470            #domain2.write_time()
1471            pass
1472
1473        ###################
1474        ##NOW TEST IT!!!
1475        ##################
1476
1477        print '><><><><>>'
1478        bits = [ 'vertex_coordinates']
1479
1480        for quantity in ['elevation','xmomentum','ymomentum']:
1481            #bits.append('quantities["%s"].get_integral()'%quantity)
1482            bits.append('get_quantity("%s").get_values()' %quantity)
1483
1484        for bit in bits:
1485            print bit
1486            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
1487
1488
1489    def test_decimate_dem(self):
1490        """Test decimation of dem file
1491        """
1492
1493        import os
1494        from Scientific.IO.NetCDF import NetCDFFile
1495
1496        #Write test dem file
1497        root = 'decdemtest'
1498
1499        filename = root + '.dem'
1500        fid = NetCDFFile(filename, netcdf_mode_w)
1501
1502        fid.institution = 'Geoscience Australia'
1503        fid.description = 'NetCDF DEM format for compact and portable ' +\
1504                          'storage of spatial point data'
1505
1506        nrows = 15
1507        ncols = 18
1508
1509        fid.ncols = ncols
1510        fid.nrows = nrows
1511        fid.xllcorner = 2000.5
1512        fid.yllcorner = 3000.5
1513        fid.cellsize = 25
1514        fid.NODATA_value = -9999
1515
1516        fid.zone = 56
1517        fid.false_easting = 0.0
1518        fid.false_northing = 0.0
1519        fid.projection = 'UTM'
1520        fid.datum = 'WGS84'
1521        fid.units = 'METERS'
1522
1523        fid.createDimension('number_of_points', nrows*ncols)
1524
1525        fid.createVariable('elevation', netcdf_float, ('number_of_points',))
1526
1527        elevation = fid.variables['elevation']
1528
1529        elevation[:] = (num.arange(nrows*ncols))
1530
1531        fid.close()
1532
1533        #generate the elevation values expected in the decimated file
1534        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
1535                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
1536                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
1537                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
1538                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
1539                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
1540                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
1541                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
1542                         (144+145+146+162+163+164+180+181+182) / 9.0,
1543                         (148+149+150+166+167+168+184+185+186) / 9.0,
1544                         (152+153+154+170+171+172+188+189+190) / 9.0,
1545                         (156+157+158+174+175+176+192+193+194) / 9.0,
1546                         (216+217+218+234+235+236+252+253+254) / 9.0,
1547                         (220+221+222+238+239+240+256+257+258) / 9.0,
1548                         (224+225+226+242+243+244+260+261+262) / 9.0,
1549                         (228+229+230+246+247+248+264+265+266) / 9.0]
1550
1551        # generate a stencil for computing the decimated values
1552        stencil = num.ones((3,3), num.float) / 9.0
1553
1554        decimate_dem(root, stencil=stencil, cellsize_new=100)
1555
1556        # Open decimated NetCDF file
1557        fid = NetCDFFile(root + '_100.dem', netcdf_mode_r)
1558
1559        # Get decimated elevation
1560        elevation = fid.variables['elevation']
1561
1562        # Check values
1563        assert num.allclose(elevation, ref_elevation)
1564
1565        # Cleanup
1566        fid.close()
1567
1568        os.remove(root + '.dem')
1569        os.remove(root + '_100.dem')
1570
1571    def test_decimate_dem_NODATA(self):
1572        """Test decimation of dem file that includes NODATA values
1573        """
1574
1575        import os
1576        from Scientific.IO.NetCDF import NetCDFFile
1577
1578        # Write test dem file
1579        root = 'decdemtest'
1580
1581        filename = root + '.dem'
1582        fid = NetCDFFile(filename, netcdf_mode_w)
1583
1584        fid.institution = 'Geoscience Australia'
1585        fid.description = 'NetCDF DEM format for compact and portable ' +\
1586                          'storage of spatial point data'
1587
1588        nrows = 15
1589        ncols = 18
1590        NODATA_value = -9999
1591
1592        fid.ncols = ncols
1593        fid.nrows = nrows
1594        fid.xllcorner = 2000.5
1595        fid.yllcorner = 3000.5
1596        fid.cellsize = 25
1597        fid.NODATA_value = NODATA_value
1598
1599        fid.zone = 56
1600        fid.false_easting = 0.0
1601        fid.false_northing = 0.0
1602        fid.projection = 'UTM'
1603        fid.datum = 'WGS84'
1604        fid.units = 'METERS'
1605
1606        fid.createDimension('number_of_points', nrows*ncols)
1607
1608        fid.createVariable('elevation', netcdf_float, ('number_of_points',))
1609
1610        elevation = fid.variables['elevation']
1611
1612        # Generate initial elevation values
1613        elevation_tmp = (num.arange(nrows*ncols))
1614
1615        # Add some NODATA values
1616        elevation_tmp[0]   = NODATA_value
1617        elevation_tmp[95]  = NODATA_value
1618        elevation_tmp[188] = NODATA_value
1619        elevation_tmp[189] = NODATA_value
1620        elevation_tmp[190] = NODATA_value
1621        elevation_tmp[209] = NODATA_value
1622        elevation_tmp[252] = NODATA_value
1623
1624        elevation[:] = elevation_tmp
1625
1626        fid.close()
1627
1628        # Generate the elevation values expected in the decimated file
1629        ref_elevation = [NODATA_value,
1630                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
1631                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
1632                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
1633                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
1634                         NODATA_value,
1635                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
1636                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
1637                         (144+145+146+162+163+164+180+181+182) / 9.0,
1638                         (148+149+150+166+167+168+184+185+186) / 9.0,
1639                         NODATA_value,
1640                         (156+157+158+174+175+176+192+193+194) / 9.0,
1641                         NODATA_value,
1642                         (220+221+222+238+239+240+256+257+258) / 9.0,
1643                         (224+225+226+242+243+244+260+261+262) / 9.0,
1644                         (228+229+230+246+247+248+264+265+266) / 9.0]
1645
1646        # Generate a stencil for computing the decimated values
1647        stencil = num.ones((3,3), num.float) / 9.0
1648
1649        decimate_dem(root, stencil=stencil, cellsize_new=100)
1650
1651        # Open decimated NetCDF file
1652        fid = NetCDFFile(root + '_100.dem', netcdf_mode_r)
1653
1654        # Get decimated elevation
1655        elevation = fid.variables['elevation']
1656
1657        # Check values
1658        assert num.allclose(elevation, ref_elevation)
1659
1660        # Cleanup
1661        fid.close()
1662
1663        os.remove(root + '.dem')
1664        os.remove(root + '_100.dem')     
1665
1666       
1667    def test_file_boundary_stsIV_sinewave_ordering(self):
1668        """test_file_boundary_stsIV_sinewave_ordering(self):
1669        Read correct points from ordering file and apply sts to boundary
1670        This one uses a sine wave and compares to time boundary
1671        """
1672       
1673        from anuga.shallow_water import Domain
1674        from anuga.shallow_water import Reflective_boundary
1675        from anuga.shallow_water import Dirichlet_boundary
1676        from anuga.shallow_water import File_boundary
1677        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1678
1679        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
1680        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
1681        tide = 0.35
1682        time_step_count = 50
1683        time_step = 0.1
1684        times_ref = num.arange(0, time_step_count*time_step, time_step)
1685       
1686        n=len(lat_long_points)
1687        first_tstep=num.ones(n,num.int)
1688        last_tstep=(time_step_count)*num.ones(n,num.int)
1689       
1690        gauge_depth=20*num.ones(n,num.float)
1691       
1692        ha1=num.ones((n,time_step_count),num.float)
1693        ua1=3.*num.ones((n,time_step_count),num.float)
1694        va1=2.*num.ones((n,time_step_count),num.float)
1695        for i in range(n):
1696            ha1[i]=num.sin(times_ref)
1697       
1698       
1699        base_name, files = self.write_mux2(lat_long_points,
1700                                           time_step_count, time_step,
1701                                           first_tstep, last_tstep,
1702                                           depth=gauge_depth,
1703                                           ha=ha1,
1704                                           ua=ua1,
1705                                           va=va1)
1706
1707        # Write order file
1708        file_handle, order_base_name = tempfile.mkstemp("")
1709        os.close(file_handle)
1710        os.remove(order_base_name)
1711        d=","
1712        order_file=order_base_name+'order.txt'
1713        fid=open(order_file,'w')
1714       
1715        # Write Header
1716        header='index, longitude, latitude\n'
1717        fid.write(header)
1718        indices=[3,0,1]
1719        for i in indices:
1720            line=str(i)+d+str(lat_long_points[i][1])+d+\
1721                str(lat_long_points[i][0])+"\n"
1722            fid.write(line)
1723        fid.close()
1724
1725        sts_file=base_name
1726        urs2sts(base_name, basename_out=sts_file,
1727                ordering_filename=order_file,
1728                mean_stage=tide,
1729                verbose=False)
1730        self.delete_mux(files)
1731       
1732       
1733       
1734        # Now read the sts file and check that values have been stored correctly.
1735        fid = NetCDFFile(sts_file + '.sts')
1736
1737        # Check the time vector
1738        times = fid.variables['time'][:]
1739       
1740        #print times
1741
1742        # Check sts quantities
1743        stage = fid.variables['stage'][:]
1744        xmomentum = fid.variables['xmomentum'][:]
1745        ymomentum = fid.variables['ymomentum'][:]
1746        elevation = fid.variables['elevation'][:]
1747
1748        #print stage
1749        #print xmomentum
1750        #print ymomentum
1751        #print elevation
1752       
1753       
1754
1755        # Create beginnings of boundary polygon based on sts_boundary
1756        boundary_polygon = create_sts_boundary(base_name)
1757       
1758        os.remove(order_file)
1759
1760        # Append the remaining part of the boundary polygon to be defined by
1761        # the user
1762        bounding_polygon_utm=[]
1763        for point in bounding_polygon:
1764            zone,easting,northing=redfearn(point[0],point[1])
1765            bounding_polygon_utm.append([easting,northing])
1766
1767        boundary_polygon.append(bounding_polygon_utm[3])
1768        boundary_polygon.append(bounding_polygon_utm[4])
1769
1770        #print 'boundary_polygon', boundary_polygon
1771       
1772        plot=False
1773        if plot:
1774            from pylab import plot,show,axis
1775            boundary_polygon=ensure_numeric(boundary_polygon)
1776            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
1777            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
1778            plot(boundary_polygon[:,0], boundary_polygon[:,1])
1779            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1])
1780            show()
1781
1782        assert num.allclose(bounding_polygon_utm,boundary_polygon)
1783
1784
1785        extent_res=1000000
1786        meshname = 'urs_test_mesh' + '.tsh'
1787        interior_regions=None
1788        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1789       
1790        # have to change boundary tags from last example because now bounding
1791        # polygon starts in different place.
1792        create_mesh_from_regions(boundary_polygon,
1793                                 boundary_tags=boundary_tags,
1794                                 maximum_triangle_area=extent_res,
1795                                 filename=meshname,
1796                                 interior_regions=interior_regions,
1797                                 verbose=False)
1798       
1799        domain_fbound = Domain(meshname)
1800        domain_fbound.set_quantity('stage', tide)
1801        Bf = File_boundary(sts_file+'.sts', 
1802                           domain_fbound, 
1803                           boundary_polygon=boundary_polygon)
1804        Br = Reflective_boundary(domain_fbound)
1805
1806        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1807        finaltime=time_step*(time_step_count-1)
1808        yieldstep=time_step
1809        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
1810   
1811        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1812                                                   finaltime=finaltime, 
1813                                                   skip_initial_step=False)):
1814            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
1815   
1816       
1817        domain_time = Domain(meshname)
1818        domain_time.set_quantity('stage', tide)
1819        Br = Reflective_boundary(domain_time)
1820        Bw = Time_boundary(domain=domain_time,
1821                         f=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
1822        domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
1823       
1824        temp_time=num.zeros(int(finaltime/yieldstep)+1,num.float)
1825        for i, t in enumerate(domain_time.evolve(yieldstep=yieldstep,
1826                                                   finaltime=finaltime, 
1827                                                   skip_initial_step=False)):
1828            temp_time[i]=domain_time.quantities['stage'].centroid_values[2]
1829
1830
1831
1832        #print temp_fbound
1833        #print temp_time
1834
1835        #print domain_fbound.quantities['stage'].vertex_values
1836        #print domain_time.quantities['stage'].vertex_values
1837       
1838        assert num.allclose(temp_fbound, temp_time)               
1839        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
1840                            domain_time.quantities['stage'].vertex_values)
1841                       
1842        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
1843                            domain_time.quantities['xmomentum'].vertex_values)                       
1844                       
1845        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
1846                            domain_time.quantities['ymomentum'].vertex_values)                                               
1847       
1848
1849        try:
1850            os.remove(sts_file+'.sts')
1851        except:
1852            # Windoze can't remove this file for some reason
1853            pass
1854       
1855        os.remove(meshname)
1856       
1857       
1858
1859       
1860       
1861    def test_file_boundary_sts_time_limit(self):
1862        """test_file_boundary_stsIV_sinewave_ordering(self):
1863        Read correct points from ordering file and apply sts to boundary
1864        This one uses a sine wave and compares to time boundary
1865       
1866        This one tests that times used can be limited by upper limit
1867        """
1868       
1869        from anuga.shallow_water import Domain
1870        from anuga.shallow_water import Reflective_boundary
1871        from anuga.shallow_water import Dirichlet_boundary
1872        from anuga.shallow_water import File_boundary
1873        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1874
1875        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
1876        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
1877        tide = 0.35
1878        time_step_count = 50
1879        time_step = 0.1
1880        times_ref = num.arange(0, time_step_count*time_step, time_step)
1881       
1882        n=len(lat_long_points)
1883        first_tstep=num.ones(n,num.int)
1884        last_tstep=(time_step_count)*num.ones(n,num.int)
1885       
1886        gauge_depth=20*num.ones(n,num.float)
1887       
1888        ha1=num.ones((n,time_step_count),num.float)
1889        ua1=3.*num.ones((n,time_step_count),num.float)
1890        va1=2.*num.ones((n,time_step_count),num.float)
1891        for i in range(n):
1892            ha1[i]=num.sin(times_ref)
1893       
1894       
1895        base_name, files = self.write_mux2(lat_long_points,
1896                                           time_step_count, time_step,
1897                                           first_tstep, last_tstep,
1898                                           depth=gauge_depth,
1899                                           ha=ha1,
1900                                           ua=ua1,
1901                                           va=va1)
1902
1903        # Write order file
1904        file_handle, order_base_name = tempfile.mkstemp("")
1905        os.close(file_handle)
1906        os.remove(order_base_name)
1907        d=","
1908        order_file=order_base_name+'order.txt'
1909        fid=open(order_file,'w')
1910       
1911        # Write Header
1912        header='index, longitude, latitude\n'
1913        fid.write(header)
1914        indices=[3,0,1]
1915        for i in indices:
1916            line=str(i)+d+str(lat_long_points[i][1])+d+\
1917                str(lat_long_points[i][0])+"\n"
1918            fid.write(line)
1919        fid.close()
1920
1921        sts_file=base_name
1922        urs2sts(base_name, basename_out=sts_file,
1923                ordering_filename=order_file,
1924                mean_stage=tide,
1925                verbose=False)
1926        self.delete_mux(files)
1927       
1928       
1929       
1930        # Now read the sts file and check that values have been stored correctly.
1931        fid = NetCDFFile(sts_file + '.sts')
1932
1933        # Check the time vector
1934        times = fid.variables['time'][:]
1935        starttime = fid.starttime[0]
1936        #print times
1937        #print starttime
1938
1939        # Check sts quantities
1940        stage = fid.variables['stage'][:]
1941        xmomentum = fid.variables['xmomentum'][:]
1942        ymomentum = fid.variables['ymomentum'][:]
1943        elevation = fid.variables['elevation'][:]
1944
1945       
1946
1947        # Create beginnings of boundary polygon based on sts_boundary
1948        boundary_polygon = create_sts_boundary(base_name)
1949       
1950        os.remove(order_file)
1951
1952        # Append the remaining part of the boundary polygon to be defined by
1953        # the user
1954        bounding_polygon_utm=[]
1955        for point in bounding_polygon:
1956            zone,easting,northing=redfearn(point[0],point[1])
1957            bounding_polygon_utm.append([easting,northing])
1958
1959        boundary_polygon.append(bounding_polygon_utm[3])
1960        boundary_polygon.append(bounding_polygon_utm[4])
1961
1962        #print 'boundary_polygon', boundary_polygon
1963       
1964
1965        assert num.allclose(bounding_polygon_utm,boundary_polygon)
1966
1967
1968        extent_res=1000000
1969        meshname = 'urs_test_mesh' + '.tsh'
1970        interior_regions=None
1971        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1972       
1973        # have to change boundary tags from last example because now bounding
1974        # polygon starts in different place.
1975        create_mesh_from_regions(boundary_polygon,
1976                                 boundary_tags=boundary_tags,
1977                                 maximum_triangle_area=extent_res,
1978                                 filename=meshname,
1979                                 interior_regions=interior_regions,
1980                                 verbose=False)
1981       
1982        domain_fbound = Domain(meshname)
1983        domain_fbound.set_quantity('stage', tide)
1984       
1985       
1986        Bf = File_boundary(sts_file+'.sts', 
1987                           domain_fbound, 
1988                           boundary_polygon=boundary_polygon)
1989        time_vec = Bf.F.get_time()
1990        assert num.allclose(Bf.F.starttime, starttime)
1991        assert num.allclose(time_vec, times_ref)                                   
1992       
1993        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
1994            Bf = File_boundary(sts_file+'.sts', 
1995                               domain_fbound, 
1996                               time_limit=time_limit+starttime,
1997                               boundary_polygon=boundary_polygon)
1998       
1999            time_vec = Bf.F.get_time()
2000            assert num.allclose(Bf.F.starttime, starttime)           
2001            assert num.alltrue(time_vec < time_limit)
2002           
2003           
2004        try:   
2005            Bf = File_boundary(sts_file+'.sts', 
2006                               domain_fbound, 
2007                               time_limit=-1+starttime,
2008                               boundary_polygon=boundary_polygon)           
2009            time_vec = Bf.F.get_time()   
2010            print time_vec   
2011        except AssertionError:
2012            pass
2013        else:
2014            raise Exception, 'Should have raised Exception here'
2015
2016    def test_lon_lat2grid(self):
2017        lonlatdep = [
2018            [ 113.06700134  ,  -26.06669998 ,   1.        ] ,
2019            [ 113.06700134  ,  -26.33329964 ,   3.        ] ,
2020            [ 113.19999695  ,  -26.06669998 ,   2.        ] ,
2021            [ 113.19999695  ,  -26.33329964 ,   4.        ] ]
2022           
2023        long, lat, quantity = lon_lat2grid(lonlatdep)
2024
2025        for i, result in enumerate(lat):
2026            assert lonlatdep [i][1] == result
2027        assert len(lat) == 2 
2028
2029        for i, result in enumerate(long):
2030            assert lonlatdep [i*2][0] == result
2031        assert len(long) == 2
2032
2033        for i,q in enumerate(quantity):
2034            assert q == i+1
2035           
2036    def test_lon_lat2grid_bad(self):
2037        lonlatdep  = [
2038            [ -26.06669998,  113.06700134,    1.        ],
2039            [ -26.06669998 , 113.19999695 ,   2.        ],
2040            [ -26.06669998 , 113.33300018,    3.        ],
2041            [ -26.06669998 , 113.43299866   , 4.        ],
2042            [ -26.20000076 , 113.06700134,    5.        ],
2043            [ -26.20000076 , 113.19999695 ,   6.        ],
2044            [ -26.20000076 , 113.33300018  ,  7.        ],
2045            [ -26.20000076 , 113.43299866   , 8.        ],
2046            [ -26.33329964 , 113.06700134,    9.        ],
2047            [ -26.33329964 , 113.19999695 ,   10.        ],
2048            [ -26.33329964 , 113.33300018  ,  11.        ],
2049            [ -26.33329964 , 113.43299866 ,   12.        ],
2050            [ -26.43330002 , 113.06700134 ,   13        ],
2051            [ -26.43330002 , 113.19999695 ,   14.        ],
2052            [ -26.43330002 , 113.33300018,    15.        ],
2053            [ -26.43330002 , 113.43299866,    16.        ]]
2054        try:
2055            long, lat, quantity = lon_lat2grid(lonlatdep)
2056        except AssertionError:
2057            pass
2058        else:
2059            msg = 'Should have raised exception'
2060            raise msg
2061       
2062    def test_lon_lat2gridII(self):
2063        lonlatdep = [
2064            [ 113.06700134  ,  -26.06669998 ,   1.        ] ,
2065            [ 113.06700134  ,  -26.33329964 ,   2.        ] ,
2066            [ 113.19999695  ,  -26.06669998 ,   3.        ] ,
2067            [ 113.19999695  ,  -26.344329964 ,   4.        ] ]
2068        try:
2069            long, lat, quantity = lon_lat2grid(lonlatdep)
2070        except AssertionError:
2071            pass
2072        else:
2073            msg = 'Should have raised exception'
2074            raise msg
2075       
2076    #### END TESTS FOR URS 2 SWW  ###
2077
2078
2079    def test_triangulation(self):
2080        #
2081       
2082       
2083        filename = tempfile.mktemp("_data_manager.sww")
2084        outfile = NetCDFFile(filename, netcdf_mode_w)
2085        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
2086        volumes = (0,1,2)
2087        elevation = [0,1,2]
2088        new_origin = None
2089        new_origin = Geo_reference(56, 0, 0)
2090        times = [0, 10]
2091        number_of_volumes = len(volumes)
2092        number_of_points = len(points_utm)
2093        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
2094        sww.store_header(outfile, times, number_of_volumes,
2095                         number_of_points, description='fully sick testing',
2096                         verbose=self.verbose,sww_precision=netcdf_float)
2097        sww.store_triangulation(outfile, points_utm, volumes,
2098                                elevation,  new_origin=new_origin,
2099                                verbose=self.verbose)       
2100        outfile.close()
2101        fid = NetCDFFile(filename)
2102
2103        x = fid.variables['x'][:]
2104        y = fid.variables['y'][:]
2105        fid.close()
2106
2107        assert num.allclose(num.array(map(None, x,y)), points_utm)
2108        os.remove(filename)
2109
2110       
2111    def test_triangulationII(self):
2112        #
2113       
2114       
2115        filename = tempfile.mktemp("_data_manager.sww")
2116        outfile = NetCDFFile(filename, netcdf_mode_w)
2117        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
2118        volumes = (0,1,2)
2119        elevation = [0,1,2]
2120        new_origin = None
2121        #new_origin = Geo_reference(56, 0, 0)
2122        times = [0, 10]
2123        number_of_volumes = len(volumes)
2124        number_of_points = len(points_utm)
2125        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
2126        sww.store_header(outfile, times, number_of_volumes,
2127                         number_of_points, description='fully sick testing',
2128                         verbose=self.verbose,sww_precision=netcdf_float)
2129        sww.store_triangulation(outfile, points_utm, volumes,
2130                                new_origin=new_origin,
2131                                verbose=self.verbose)
2132        sww.store_static_quantities(outfile, elevation=elevation)                               
2133                               
2134        outfile.close()
2135        fid = NetCDFFile(filename)
2136
2137        x = fid.variables['x'][:]
2138        y = fid.variables['y'][:]
2139        results_georef = Geo_reference()
2140        results_georef.read_NetCDF(fid)
2141        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
2142        fid.close()
2143
2144        assert num.allclose(num.array(map(None, x,y)), points_utm)
2145        os.remove(filename)
2146
2147       
2148    def test_triangulation_new_origin(self):
2149        #
2150       
2151       
2152        filename = tempfile.mktemp("_data_manager.sww")
2153        outfile = NetCDFFile(filename, netcdf_mode_w)
2154        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
2155        volumes = (0,1,2)
2156        elevation = [0,1,2]
2157        new_origin = None
2158        new_origin = Geo_reference(56, 1, 554354)
2159        points_utm = new_origin.change_points_geo_ref(points_utm)
2160        times = [0, 10]
2161        number_of_volumes = len(volumes)
2162        number_of_points = len(points_utm)
2163        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
2164        sww.store_header(outfile, times, number_of_volumes,
2165                         number_of_points, description='fully sick testing',
2166                         verbose=self.verbose,sww_precision=netcdf_float)
2167        sww.store_triangulation(outfile, points_utm, volumes,
2168                                elevation,  new_origin=new_origin,
2169                                verbose=self.verbose)
2170        outfile.close()
2171        fid = NetCDFFile(filename)
2172
2173        x = fid.variables['x'][:]
2174        y = fid.variables['y'][:]
2175        results_georef = Geo_reference()
2176        results_georef.read_NetCDF(fid)
2177        assert results_georef == new_origin
2178        fid.close()
2179
2180        absolute = Geo_reference(56, 0,0)
2181        assert num.allclose(num.array(
2182            absolute.change_points_geo_ref(map(None, x,y),
2183                                           new_origin)),points_utm)
2184       
2185        os.remove(filename)
2186       
2187    def test_triangulation_points_georeference(self):
2188        #
2189       
2190       
2191        filename = tempfile.mktemp("_data_manager.sww")
2192        outfile = NetCDFFile(filename, netcdf_mode_w)
2193        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
2194        volumes = (0,1,2)
2195        elevation = [0,1,2]
2196        new_origin = None
2197        points_georeference = Geo_reference(56, 1, 554354)
2198        points_utm = points_georeference.change_points_geo_ref(points_utm)
2199        times = [0, 10]
2200        number_of_volumes = len(volumes)
2201        number_of_points = len(points_utm)
2202        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
2203        sww.store_header(outfile, times, number_of_volumes,
2204                         number_of_points, description='fully sick testing',
2205                         verbose=self.verbose,sww_precision=netcdf_float)
2206        sww.store_triangulation(outfile, points_utm, volumes,
2207                                elevation,  new_origin=new_origin,
2208                                points_georeference=points_georeference,
2209                                verbose=self.verbose)       
2210        outfile.close()
2211        fid = NetCDFFile(filename)
2212
2213        x = fid.variables['x'][:]
2214        y = fid.variables['y'][:]
2215        results_georef = Geo_reference()
2216        results_georef.read_NetCDF(fid)
2217        assert results_georef == points_georeference
2218        fid.close()
2219
2220        assert num.allclose(num.array(map(None, x,y)), points_utm)
2221        os.remove(filename)
2222       
2223    def test_triangulation_2_geo_refs(self):
2224        #
2225       
2226       
2227        filename = tempfile.mktemp("_data_manager.sww")
2228        outfile = NetCDFFile(filename, netcdf_mode_w)
2229        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
2230        volumes = (0,1,2)
2231        elevation = [0,1,2]
2232        new_origin = Geo_reference(56, 1, 1)
2233        points_georeference = Geo_reference(56, 0, 0)
2234        points_utm = points_georeference.change_points_geo_ref(points_utm)
2235        times = [0, 10]
2236        number_of_volumes = len(volumes)
2237        number_of_points = len(points_utm)
2238        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
2239        sww.store_header(outfile, times, number_of_volumes,
2240                         number_of_points, description='fully sick testing',
2241                         verbose=self.verbose,sww_precision=netcdf_float)
2242        sww.store_triangulation(outfile, points_utm, volumes,
2243                                elevation,  new_origin=new_origin,
2244                                points_georeference=points_georeference,
2245                                verbose=self.verbose)       
2246        outfile.close()
2247        fid = NetCDFFile(filename)
2248
2249        x = fid.variables['x'][:]
2250        y = fid.variables['y'][:]
2251        results_georef = Geo_reference()
2252        results_georef.read_NetCDF(fid)
2253        assert results_georef == new_origin
2254        fid.close()
2255
2256
2257        absolute = Geo_reference(56, 0,0)
2258        assert num.allclose(num.array(
2259            absolute.change_points_geo_ref(map(None, x,y),
2260                                           new_origin)),points_utm)
2261        os.remove(filename)
2262       
2263
2264       
2265    def test_get_all_swwfiles(self):
2266        try:
2267            swwfiles = get_all_swwfiles('','test.txt')  #Invalid
2268        except IOError:
2269            pass
2270        else:
2271            raise 'Should have raised exception' 
2272       
2273    def test_get_all_swwfiles1(self):
2274       
2275        temp_dir = tempfile.mkdtemp('','sww_test')
2276        filename0 = tempfile.mktemp('.sww','test',temp_dir)
2277        filename1 = tempfile.mktemp('.sww','test',temp_dir)
2278        filename2 = tempfile.mktemp('.sww','test',temp_dir)
2279        filename3 = tempfile.mktemp('.sww','test',temp_dir)
2280       
2281        #print'filename', filename0,filename1,filename2,filename3
2282       
2283        fid0 = open(filename0, 'w')
2284        fid1 = open(filename1, 'w')
2285        fid2 = open(filename2, 'w')
2286        fid3 = open(filename3, 'w')
2287
2288        fid0.write('hello')
2289        fid1.write('hello')
2290        fid2.write('hello')
2291        fid3.write('hello')
2292       
2293        fid0.close()
2294        fid1.close()
2295        fid2.close()
2296        fid3.close()
2297       
2298       
2299        dir, name0 = os.path.split(filename0)
2300        #print 'dir',dir,name0
2301       
2302        iterate=get_all_swwfiles(dir,'test')
2303       
2304        del_dir(temp_dir)
2305#        removeall(temp_dir)
2306
2307        _, name0 = os.path.split(filename0) 
2308        #print'name0',name0[:-4],iterate[0]   
2309        _, name1 = os.path.split(filename1)       
2310        _, name2 = os.path.split(filename2)       
2311        _, name3 = os.path.split(filename3)       
2312
2313        assert name0[:-4] in iterate
2314        assert name1[:-4] in iterate
2315        assert name2[:-4] in iterate
2316        assert name3[:-4] in iterate
2317       
2318        assert len(iterate)==4
2319
2320 
2321    def test_points2polygon(self): 
2322        att_dict = {}
2323        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[0.0, 0.0]])
2324        att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
2325        att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
2326       
2327        fileName = tempfile.mktemp(".csv")
2328       
2329        G = Geospatial_data(pointlist, att_dict)
2330       
2331        G.export_points_file(fileName)
2332       
2333        polygon = points2polygon(fileName)
2334       
2335        # This test may fail if the order changes
2336        assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])
2337       
2338   
2339    def test_csv2polygons(self):
2340        """test_csv2polygons
2341        """
2342       
2343        path = get_pathname_from_package('anuga.shallow_water')               
2344        testfile = os.path.join(path, 'polygon_values_example.csv')               
2345
2346        polygons, values = csv2polygons(testfile, 
2347                                        value_name='floors')
2348
2349        assert len(polygons) == 7, 'Must have seven polygons'
2350        assert len(values) == 7, 'Must have seven values'
2351           
2352        # Known floor values
2353        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
2354       
2355        # Known polygon values
2356        known_polys = {'1': [[422681.61,871117.55],
2357                             [422691.02,871117.60],
2358                             [422690.87,871084.23],
2359                             [422649.36,871081.85],
2360                             [422649.36,871080.39],
2361                             [422631.86,871079.50],
2362                             [422631.72,871086.75],
2363                             [422636.75,871087.20],
2364                             [422636.75,871091.50],
2365                             [422649.66,871092.09],
2366                             [422649.83,871084.91],
2367                             [422652.94,871084.90],
2368                             [422652.84,871092.39],
2369                             [422681.83,871093.73],
2370                             [422681.61,871117.55]],
2371                       '2': [[422664.22,870785.46],
2372                             [422672.48,870780.14],
2373                             [422668.17,870772.62],
2374                             [422660.35,870777.17],
2375                             [422664.22,870785.46]],
2376                       '3': [[422661.30,871215.06],
2377                             [422667.50,871215.70],
2378                             [422668.30,871204.86],
2379                             [422662.21,871204.33],
2380                             [422661.30,871215.06]],
2381                       '4': [[422473.44,871191.22],
2382                             [422478.33,871192.26],
2383                             [422479.52,871186.03],
2384                             [422474.78,871185.14],
2385                             [422473.44,871191.22]],
2386                       '5': [[422369.69,871049.29],
2387                             [422378.63,871053.58],
2388                             [422383.91,871044.51],
2389                             [422374.97,871040.32],
2390                             [422369.69,871049.29]],
2391                       '8': [[422730.56,871203.13],
2392                             [422734.10,871204.90],
2393                             [422735.26,871202.18],
2394                             [422731.87,871200.58],
2395                             [422730.56,871203.13]],
2396                       '9': [[422659.85,871213.80],
2397                             [422660.91,871210.97],
2398                             [422655.42,871208.85],
2399                             [422654.36,871211.68],
2400                             [422659.85,871213.80]]
2401                       }       
2402       
2403
2404       
2405               
2406        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
2407            assert id in polygons.keys()
2408            assert id in values.keys()           
2409
2410            assert int(values[id]) == int(floors[id])
2411            assert len(polygons[id]) == len(known_polys[id])
2412            assert num.allclose(polygons[id], known_polys[id])
2413
2414
2415    def test_csv2polygons_with_clipping(self):
2416        """test_csv2polygons with optional clipping
2417        """
2418        #FIXME(Ole): Not Done!!
2419       
2420        path = get_pathname_from_package('anuga.shallow_water')               
2421        testfile = os.path.join(path, 'polygon_values_example.csv')               
2422
2423        polygons, values = csv2polygons(testfile, 
2424                                        value_name='floors',
2425                                        clipping_polygons=None)
2426
2427        assert len(polygons) == 7, 'Must have seven polygons'
2428        assert len(values) == 7, 'Must have seven values'
2429           
2430        # Known floor values
2431        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
2432       
2433        # Known polygon values
2434        known_polys = {'1': [[422681.61,871117.55],
2435                             [422691.02,871117.60],
2436                             [422690.87,871084.23],
2437                             [422649.36,871081.85],
2438                             [422649.36,871080.39],
2439                             [422631.86,871079.50],
2440                             [422631.72,871086.75],
2441                             [422636.75,871087.20],
2442                             [422636.75,871091.50],
2443                             [422649.66,871092.09],
2444                             [422649.83,871084.91],
2445                             [422652.94,871084.90],
2446                             [422652.84,871092.39],
2447                             [422681.83,871093.73],
2448                             [422681.61,871117.55]],
2449                       '2': [[422664.22,870785.46],
2450                             [422672.48,870780.14],
2451                             [422668.17,870772.62],
2452                             [422660.35,870777.17],
2453                             [422664.22,870785.46]],
2454                       '3': [[422661.30,871215.06],
2455                             [422667.50,871215.70],
2456                             [422668.30,871204.86],
2457                             [422662.21,871204.33],
2458                             [422661.30,871215.06]],
2459                       '4': [[422473.44,871191.22],
2460                             [422478.33,871192.26],
2461                             [422479.52,871186.03],
2462                             [422474.78,871185.14],
2463                             [422473.44,871191.22]],
2464                       '5': [[422369.69,871049.29],
2465                             [422378.63,871053.58],
2466                             [422383.91,871044.51],
2467                             [422374.97,871040.32],
2468                             [422369.69,871049.29]],
2469                       '8': [[422730.56,871203.13],
2470                             [422734.10,871204.90],
2471                             [422735.26,871202.18],
2472                             [422731.87,871200.58],
2473                             [422730.56,871203.13]],
2474                       '9': [[422659.85,871213.80],
2475                             [422660.91,871210.97],
2476                             [422655.42,871208.85],
2477                             [422654.36,871211.68],
2478                             [422659.85,871213.80]]
2479                       }       
2480       
2481
2482       
2483               
2484        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
2485            assert id in polygons.keys()
2486            assert id in values.keys()           
2487
2488            assert int(values[id]) == int(floors[id])
2489            assert len(polygons[id]) == len(known_polys[id])
2490            assert num.allclose(polygons[id], known_polys[id])
2491
2492
2493
2494
2495   
2496    def test_csv2building_polygons(self):
2497        """test_csv2building_polygons
2498        """
2499       
2500        path = get_pathname_from_package('anuga.shallow_water')               
2501        testfile = os.path.join(path, 'polygon_values_example.csv')               
2502
2503        polygons, values = csv2building_polygons(testfile, 
2504                                                 floor_height=3)
2505
2506        assert len(polygons) == 7, 'Must have seven polygons'
2507        assert len(values) == 7, 'Must have seven values'
2508           
2509        # Known floor values
2510        floors = {'1': 6, '2': 0, '3': 3, '4': 6, '5': 0, '8': 3, '9': 3}
2511       
2512               
2513        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
2514            assert id in polygons.keys()
2515            assert id in values.keys()           
2516
2517            assert float(values[id]) == float(floors[id])
2518
2519
2520#-------------------------------------------------------------
2521
2522if __name__ == "__main__":
2523    #suite = unittest.makeSuite(Test_Data_Manager, 'test_sww2domain2')
2524    suite = unittest.makeSuite(Test_Data_Manager, 'test_sww')
2525   
2526   
2527   
2528    # FIXME(Ole): When Ross has implemented logging, we can
2529    # probably get rid of all this:
2530    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
2531        Test_Data_Manager.verbose=True
2532        saveout = sys.stdout   
2533        filename = ".temp_verbose"
2534        fid = open(filename, 'w')
2535        sys.stdout = fid
2536    else:
2537        pass
2538    runner = unittest.TextTestRunner() #verbosity=2)
2539    runner.run(suite)
2540
2541    # Cleaning up
2542    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
2543        sys.stdout = saveout
2544        fid.close() 
2545        os.remove(filename)
2546
2547
Note: See TracBrowser for help on using the repository browser.