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

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

Some changes to allow netcdf4 use

File size: 30.4 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 anuga.file.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 get_all_swwfiles
31from anuga.utilities.file_utils import del_dir
32from anuga.utilities.numerical_tools import ensure_numeric, mean
33from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
34from anuga.config import netcdf_float, epsilon, g
35from anuga.pmesh.mesh_interface import create_mesh_from_regions
36from anuga.file_conversion.sww2dem import sww2dem_batch
37from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array, \
38                                load_csv_as_building_polygons, \
39                                load_csv_as_polygons
40from anuga.file.sts import create_sts_boundary
41from anuga.file.pts import load_pts_as_polygon
42from anuga.file.sww import Write_sww
43
44
45# import all the boundaries - some are generic, some are shallow water
46from boundaries import Reflective_boundary, \
47            Field_boundary, Transmissive_momentum_set_stage_boundary, \
48            Transmissive_stage_zero_momentum_boundary
49from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
50     import Transmissive_boundary, Dirichlet_boundary, \
51            Time_boundary, File_boundary, AWI_boundary
52
53# This is needed to run the tests of local functions
54from anuga.file_conversion.urs2sts import urs2sts
55from anuga.coordinate_transforms.redfearn import redfearn
56from anuga.coordinate_transforms.geo_reference import Geo_reference, \
57     DEFAULT_ZONE
58from anuga.geospatial_data.geospatial_data import Geospatial_data
59
60from shallow_water_domain import Domain
61
62# use helper methods from other unit test
63from anuga.file.test_mux import Test_Mux
64
65
66class Test_Data_Manager(Test_Mux):
67    # Class variable
68    verbose = False
69
70    def set_verbose(self):
71        Test_Data_Manager.verbose = True
72       
73    def setUp(self):
74        import time
75        from mesh_factory import rectangular
76       
77        self.verbose = Test_Data_Manager.verbose
78        # Create basic mesh
79        points, vertices, boundary = rectangular(2, 2)
80
81        # Create shallow water domain
82        domain = Domain(points, vertices, boundary)
83        domain.default_order = 2
84
85        # Set some field values
86        domain.set_quantity('elevation', lambda x,y: -x)
87        domain.set_quantity('friction', 0.03)
88
89
90        ######################
91        # Boundary conditions
92        B = Transmissive_boundary(domain)
93        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
94
95
96        ######################
97        #Initial condition - with jumps
98        bed = domain.quantities['elevation'].vertex_values
99        stage = num.zeros(bed.shape, num.float)
100
101        h = 0.3
102        for i in range(stage.shape[0]):
103            if i % 2 == 0:
104                stage[i,:] = bed[i,:] + h
105            else:
106                stage[i,:] = bed[i,:]
107
108        domain.set_quantity('stage', stage)
109
110
111        domain.distribute_to_vertices_and_edges()               
112        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
113
114
115        self.domain = domain
116
117        C = domain.get_vertex_coordinates()
118        self.X = C[:,0:6:2].copy()
119        self.Y = C[:,1:6:2].copy()
120
121        self.F = bed
122
123        #Write A testfile (not realistic. Values aren't realistic)
124        self.test_MOST_file = 'most_small'
125
126        longitudes = [150.66667, 150.83334, 151., 151.16667]
127        latitudes = [-34.5, -34.33333, -34.16667, -34]
128
129        long_name = 'LON'
130        lat_name = 'LAT'
131
132        nx = 4
133        ny = 4
134        six = 6
135
136
137        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
138            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
139
140            fid.createDimension(long_name,nx)
141            fid.createVariable(long_name,netcdf_float,(long_name,))
142            fid.variables[long_name].point_spacing='uneven'
143            fid.variables[long_name].units='degrees_east'
144            fid.variables[long_name][:] = longitudes
145
146            fid.createDimension(lat_name,ny)
147            fid.createVariable(lat_name,netcdf_float,(lat_name,))
148            fid.variables[lat_name].point_spacing='uneven'
149            fid.variables[lat_name].units='degrees_north'
150            fid.variables[lat_name][:] = latitudes
151
152            fid.createDimension('TIME',six)
153            fid.createVariable('TIME',netcdf_float,('TIME',))
154            fid.variables['TIME'].point_spacing='uneven'
155            fid.variables['TIME'].units='seconds'
156            fid.variables['TIME'][:] = [0.0, 0.1, 0.6, 1.1, 1.6, 2.1]
157
158
159            name = ext[1:3].upper()
160            if name == 'E.': name = 'ELEVATION'
161            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
162            fid.variables[name].units='CENTIMETERS'
163            fid.variables[name].missing_value=-1.e+034
164
165            fid.variables[name][:] = [[[0.3400644, 0, -46.63519, -6.50198],
166                                              [-0.1214216, 0, 0, 0],
167                                              [0, 0, 0, 0],
168                                              [0, 0, 0, 0]],
169                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
170                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
171                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
172                                              [0, 2.291054e-005, 2.291054e-005, 0]],
173                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
174                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
175                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
176                                              [0, 0.0001374632, 0.0001374632, 0]],
177                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
178                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
179                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
180                                              [0, 0.0002520159, 0.0002520159, 0]],
181                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
182                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
183                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
184                                              [0, 0.0003665686, 0.0003665686, 0]],
185                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
186                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
187                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
188                                              [0, 0.0004811212, 0.0004811212, 0]]]
189
190
191            fid.close()
192
193
194
195
196    def tearDown(self):
197        import os
198        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
199            #print 'Trying to remove', self.test_MOST_file + ext
200            os.remove(self.test_MOST_file + ext)
201
202    def test_sww_constant(self):
203        """Test that constant sww information can be written correctly
204        (non smooth)
205        """
206        self.domain.set_name('datatest' + str(id(self)))
207        self.domain.format = 'sww' #Remove??
208        self.domain.smooth = False
209       
210        sww = SWW_file(self.domain)
211        sww.store_connectivity()
212
213        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
214
215        # Get the variables
216        x = fid.variables['x']
217        y = fid.variables['y']
218        z = fid.variables['elevation']
219        V = fid.variables['volumes']
220
221        assert num.allclose (x[:], self.X.flatten())
222        assert num.allclose (y[:], self.Y.flatten())
223        assert num.allclose (z[:], self.F.flatten())
224
225        P = len(self.domain)
226        for k in range(P):
227            assert V[k, 0] == 3*k
228            assert V[k, 1] == 3*k+1
229            assert V[k, 2] == 3*k+2
230
231        fid.close()
232        os.remove(sww.filename)
233
234    def test_sww_header(self):
235        """Test that constant sww information can be written correctly
236        (non smooth)
237        """
238        self.domain.set_name('datatest' + str(id(self)))
239        self.domain.format = 'sww' #Remove??
240        self.domain.smooth = False
241
242        sww = SWW_file(self.domain)
243        sww.store_connectivity()
244
245        # Check contents
246        # Get NetCDF
247        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
248
249        # Get the variables
250        sww_revision = fid.revision_number
251        try:
252            revision_number = get_revision_number()
253        except:
254            revision_number = None
255           
256        assert str(revision_number) == sww_revision
257        fid.close()
258
259        #print "sww.filename", sww.filename
260        os.remove(sww.filename)
261
262    def test_sww_range(self):
263        """Test that constant sww information can be written correctly
264        Use non-smooth to be able to compare to quantity values.
265        """
266
267        self.domain.set_name('datatest' + str(id(self)))
268        self.domain.format = 'sww'
269        self.domain.set_store_vertices_uniquely(True)
270       
271        sww = SWW_file(self.domain)       
272
273        dqs = self.domain.get_quantity('stage')
274        dqx = self.domain.get_quantity('xmomentum')
275        dqy = self.domain.get_quantity('ymomentum')       
276        xmom_min = ymom_min = stage_min = sys.maxint
277        xmom_max = ymom_max = stage_max = -stage_min       
278        for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
279            wmax = max(dqs.get_values().flatten())
280            if wmax > stage_max: stage_max = wmax
281            wmin = min(dqs.get_values().flatten())
282            if wmin < stage_min: stage_min = wmin           
283           
284            uhmax = max(dqx.get_values().flatten())
285            if uhmax > xmom_max: xmom_max = uhmax
286            uhmin = min(dqx.get_values().flatten())
287            if uhmin < xmom_min: xmom_min = uhmin                       
288           
289            vhmax = max(dqy.get_values().flatten())
290            if vhmax > ymom_max: ymom_max = vhmax
291            vhmin = min(dqy.get_values().flatten())
292            if vhmin < ymom_min: ymom_min = vhmin                                   
293           
294           
295           
296        # Get NetCDF
297        fid = NetCDFFile(sww.filename, netcdf_mode_r) # Open existing file for append
298
299        # Get the variables
300        range = fid.variables['stage_range'][:]
301        assert num.allclose(range,[stage_min, stage_max])
302
303        range = fid.variables['xmomentum_range'][:]
304        #print range
305        assert num.allclose(range, [xmom_min, xmom_max])
306       
307        range = fid.variables['ymomentum_range'][:]
308        #print range
309        assert num.allclose(range, [ymom_min, ymom_max])       
310
311
312       
313        fid.close()
314        os.remove(sww.filename)
315
316    def test_sww_extrema(self):
317        """Test that extrema of quantities can be retrieved at every vertex
318        Extrema are updated at every *internal* timestep
319        """
320
321        domain = self.domain
322       
323        domain.set_name('extrema_test' + str(id(self)))
324        domain.format = 'sww'
325        domain.smooth = True
326
327        assert domain.quantities_to_be_monitored is None
328        assert domain.monitor_polygon is None
329        assert domain.monitor_time_interval is None       
330       
331        domain.set_quantities_to_be_monitored(['xmomentum',
332                                               'ymomentum',
333                                               'stage-elevation'])
334
335        assert domain.monitor_polygon is None
336        assert domain.monitor_time_interval is None
337
338
339        domain.set_quantities_to_be_monitored(['xmomentum',
340                                               'ymomentum',
341                                               'stage-elevation'],
342                                              polygon=domain.get_boundary_polygon(),
343                                              time_interval=[0,1])
344       
345       
346        assert len(domain.quantities_to_be_monitored) == 3
347        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
348        assert domain.quantities_to_be_monitored.has_key('xmomentum')               
349        assert domain.quantities_to_be_monitored.has_key('ymomentum')       
350
351       
352        #domain.protect_against_isolated_degenerate_timesteps = True
353        #domain.tight_slope_limiters = 1
354        domain.tight_slope_limiters = 0 # Backwards compatibility
355        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
356       
357       
358        sww = SWW_file(domain)
359
360        for t in domain.evolve(yieldstep = 1, finaltime = 1):
361            pass
362            #print domain.timestepping_statistics()
363            domain.quantity_statistics(precision = '%.8f') # Silent
364
365           
366        # Get NetCDF
367        fid = NetCDFFile(sww.filename, netcdf_mode_r) # Open existing file for append
368
369        # Get the variables
370        extrema = fid.variables['stage-elevation.extrema'][:]
371        assert num.allclose(extrema, [0.00, 0.30])
372
373        loc = fid.variables['stage-elevation.min_location'][:]
374        assert num.allclose(loc, [0.16666667, 0.33333333])
375
376        loc = fid.variables['stage-elevation.max_location'][:]       
377        assert num.allclose(loc, [0.8333333, 0.16666667])       
378
379        time = fid.variables['stage-elevation.max_time'][:]
380        assert num.allclose(time, 0.0)               
381
382        extrema = fid.variables['xmomentum.extrema'][:]
383        # Padarn Note: Had to add an extra possibility here (the final one [-0.06062178  0.47518688])
384        # to pass this assertion. Cannot see what would be causing this.
385        assert num.allclose(extrema,[-0.06062178, 0.47873023]) or \
386            num.allclose(extrema, [-0.06062178, 0.47847986]) or \
387            num.allclose(extrema, [-0.06062178, 0.47848481]) or \
388            num.allclose(extrema, [-0.06062178, 0.47763887]) or \
389            num.allclose(extrema, [-0.06062178, 0.46691909])or \
390            num.allclose(extrema, [-0.06062178, 0.47503704]) or \
391            num.allclose(extrema, [-0.06062178,  0.47518688])
392
393
394       
395        extrema = fid.variables['ymomentum.extrema'][:]
396        assert num.allclose(extrema,[0.00, 0.0625786]) or num.allclose(extrema,[0.00, 0.06062178])
397
398        time_interval = fid.variables['extrema.time_interval'][:]
399        assert num.allclose(time_interval, [0,1])
400       
401        polygon = fid.variables['extrema.polygon'][:]       
402        assert num.allclose(polygon, domain.get_boundary_polygon())
403       
404        fid.close()
405        #print "sww.filename", sww.filename
406        os.remove(sww.filename)
407
408       
409       
410    def test_sww_constant_smooth(self):
411        """Test that constant sww information can be written correctly
412        (non smooth)
413        """
414        self.domain.set_name('datatest' + str(id(self)))
415        self.domain.format = 'sww'
416        self.domain.smooth = True
417
418        sww = SWW_file(self.domain)
419        sww.store_connectivity()
420
421        # Check contents
422        # Get NetCDF
423        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
424
425        # Get the variables
426        X = fid.variables['x'][:]
427        Y = fid.variables['y'][:]
428        Z = fid.variables['elevation'][:]
429        V = fid.variables['volumes']
430
431        assert num.allclose([X[0], Y[0]], num.array([0.0, 0.0]))
432        assert num.allclose([X[1], Y[1]], num.array([0.0, 0.5]))
433        assert num.allclose([X[2], Y[2]], num.array([0.0, 1.0]))
434        assert num.allclose([X[4], Y[4]], num.array([0.5, 0.5]))
435        assert num.allclose([X[7], Y[7]], num.array([1.0, 0.5]))
436
437        assert Z[4] == -0.5
438
439        assert V[2,0] == 4
440        assert V[2,1] == 5
441        assert V[2,2] == 1
442        assert V[4,0] == 6
443        assert V[4,1] == 7
444        assert V[4,2] == 3
445
446        fid.close()
447        os.remove(sww.filename)
448       
449
450    def test_sww_variable(self):
451        """Test that sww information can be written correctly
452        """
453        self.domain.set_name('datatest' + str(id(self)))
454        self.domain.format = 'sww'
455        self.domain.smooth = True
456        self.domain.reduction = mean
457
458        sww = SWW_file(self.domain)
459        sww.store_connectivity()
460        sww.store_timestep()
461
462        # Check contents
463        # Get NetCDF
464        fid = NetCDFFile(sww.filename, netcdf_mode_r)  # Open existing file for append
465
466
467        # Get the variables
468        time = fid.variables['time']
469        stage = fid.variables['stage']
470
471        Q = self.domain.quantities['stage']
472        Q0 = Q.vertex_values[:,0]
473        Q1 = Q.vertex_values[:,1]
474        Q2 = Q.vertex_values[:,2]
475
476        A = stage[0,:]
477        #print A[0], (Q2[0,0] + Q1[1,0])/2
478        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
479        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
480        assert num.allclose(A[2], Q0[3])
481        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
482
483        #Center point
484        assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
485                                   Q0[5] + Q2[6] + Q1[7])/6)
486       
487        fid.close()
488        os.remove(sww.filename)
489
490
491    def test_sww_variable2(self):
492        """Test that sww information can be written correctly
493        multiple timesteps. Use average as reduction operator
494        """
495
496        import time, os
497
498        self.domain.set_name('datatest' + str(id(self)))
499        self.domain.format = 'sww'
500        self.domain.smooth = True
501
502        self.domain.reduction = mean
503
504        sww = SWW_file(self.domain)
505        sww.store_connectivity()
506        sww.store_timestep()
507        #self.domain.tight_slope_limiters = 1
508        self.domain.evolve_to_end(finaltime = 0.01)
509        sww.store_timestep()
510
511        # Check contents
512        # Get NetCDF
513        fid = NetCDFFile(sww.filename, netcdf_mode_r)
514
515        # Get the variables
516        x = fid.variables['x']
517        y = fid.variables['y']
518        z = fid.variables['elevation']
519        time = fid.variables['time']
520        stage = fid.variables['stage']
521
522        #Check values
523        Q = self.domain.quantities['stage']
524        Q0 = Q.vertex_values[:,0]
525        Q1 = Q.vertex_values[:,1]
526        Q2 = Q.vertex_values[:,2]
527
528        A = stage[1,:]
529        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
530        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
531        assert num.allclose(A[2], Q0[3])
532        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
533
534        #Center point
535        assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
536                                   Q0[5] + Q2[6] + Q1[7])/6)
537
538
539        fid.close()
540
541        #Cleanup
542        os.remove(sww.filename)
543
544
545    def test_sync(self):
546        """test_sync - Test info stored at each timestep is as expected (incl initial condition)
547        """
548
549        import time, os, config
550
551        self.domain.set_name('synctest')
552        self.domain.format = 'sww'
553        self.domain.smooth = False
554        self.domain.store = True
555
556        self.domain.tight_slope_limiters = True
557        self.domain.use_centroid_velocities = True       
558       
559        # In this case tight_slope_limiters as default
560        # in conjunction with protection
561        # against isolated degenerate timesteps works.
562        #self.domain.tight_slope_limiters = 1
563        #self.domain.protect_against_isolated_degenerate_timesteps = True
564
565        #print 'tight_sl', self.domain.tight_slope_limiters
566       
567
568        #Evolution
569        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
570           
571            #########self.domain.write_time(track_speeds=True)
572            stage = self.domain.quantities['stage'].vertex_values
573
574            #Get NetCDF
575            fid = NetCDFFile(self.domain.writer.filename, netcdf_mode_r)
576            stage_file = fid.variables['stage']
577           
578            if t == 0.0:
579                assert num.allclose(stage, self.initial_stage)
580                assert num.allclose(stage_file[:], stage.flatten())
581            else:
582                assert not num.allclose(stage, self.initial_stage)
583                assert not num.allclose(stage_file[:], stage.flatten())
584
585            fid.close()
586
587        os.remove(self.domain.writer.filename)
588
589
590    def test_sww_minimum_storable_height(self):
591        """Test that sww information can be written correctly
592        multiple timesteps using a different reduction operator (min)
593        """
594
595        import time, os
596
597        self.domain.set_name('datatest' + str(id(self)))
598        self.domain.format = 'sww'
599        self.domain.smooth = True
600        self.domain.reduction = min
601        self.domain.minimum_storable_height = 100
602
603        sww = SWW_file(self.domain)
604        sww.store_connectivity()
605        sww.store_timestep()
606
607        #self.domain.tight_slope_limiters = 1
608        self.domain.evolve_to_end(finaltime = 0.01)
609        sww.store_timestep()
610
611        #Check contents
612        #Get NetCDF
613        fid = NetCDFFile(sww.filename, netcdf_mode_r)
614
615        # Get the variables
616        x = fid.variables['x']
617        y = fid.variables['y']
618        z = fid.variables['elevation']
619        time = fid.variables['time']
620        stage = fid.variables['stage']
621        xmomentum = fid.variables['xmomentum']
622        ymomentum = fid.variables['ymomentum']       
623
624        #Check values
625        Q = self.domain.quantities['stage']
626        Q0 = Q.vertex_values[:,0]
627        Q1 = Q.vertex_values[:,1]
628        Q2 = Q.vertex_values[:,2]
629
630        A = stage[1,:]
631        assert num.allclose(stage[1,:], z[:])
632
633
634        assert num.allclose(xmomentum, 0.0)
635        assert num.allclose(ymomentum, 0.0)       
636       
637        fid.close()
638
639        #Cleanup
640        os.remove(sww.filename)
641
642
643    def test_file_boundary_stsIV_sinewave_ordering(self):
644        """test_file_boundary_stsIV_sinewave_ordering(self):
645        Read correct points from ordering file and apply sts to boundary
646        This one uses a sine wave and compares to time boundary
647        """
648
649        lat_long_points=[[6.01, 97.0], [6.02, 97.0], [6.05, 96.9], [6.0, 97.0]]
650        bounding_polygon=[[6.0, 97.0], [6.01, 97.0], [6.02,97.0], \
651                            [6.02,97.02], [6.00,97.02]]
652        tide = 0.35
653        time_step_count = 50
654        time_step = 0.1
655        times_ref = num.arange(0, time_step_count*time_step, time_step)
656       
657        n=len(lat_long_points)
658        first_tstep=num.ones(n,num.int)
659        last_tstep=(time_step_count)*num.ones(n,num.int)
660       
661        gauge_depth=20*num.ones(n,num.float)
662       
663        ha1=num.ones((n,time_step_count),num.float)
664        ua1=3.*num.ones((n,time_step_count),num.float)
665        va1=2.*num.ones((n,time_step_count),num.float)
666        for i in range(n):
667            ha1[i]=num.sin(times_ref)
668       
669       
670        base_name, files = self.write_mux2(lat_long_points,
671                                           time_step_count, time_step,
672                                           first_tstep, last_tstep,
673                                           depth=gauge_depth,
674                                           ha=ha1,
675                                           ua=ua1,
676                                           va=va1)
677
678        # Write order file
679        file_handle, order_base_name = tempfile.mkstemp("")
680        os.close(file_handle)
681        os.remove(order_base_name)
682        d=","
683        order_file=order_base_name+'order.txt'
684        fid=open(order_file,'w')
685       
686        # Write Header
687        header='index, longitude, latitude\n'
688        fid.write(header)
689        indices=[3,0,1]
690        for i in indices:
691            line=str(i)+d+str(lat_long_points[i][1])+d+\
692                str(lat_long_points[i][0])+"\n"
693            fid.write(line)
694        fid.close()
695
696        sts_file=base_name
697        urs2sts(base_name, basename_out=sts_file,
698                ordering_filename=order_file,
699                mean_stage=tide,
700                verbose=False)
701        self.delete_mux(files)
702       
703       
704       
705        # Now read the sts file and check that values have been stored correctly.
706        fid = NetCDFFile(sts_file + '.sts')
707
708        # Check the time vector
709        times = fid.variables['time'][:]
710       
711        #print times
712
713        # Check sts quantities
714        stage = fid.variables['stage'][:]
715        xmomentum = fid.variables['xmomentum'][:]
716        ymomentum = fid.variables['ymomentum'][:]
717        elevation = fid.variables['elevation'][:]       
718
719        # Create beginnings of boundary polygon based on sts_boundary
720        boundary_polygon = create_sts_boundary(base_name)
721       
722        os.remove(order_file)
723
724        # Append the remaining part of the boundary polygon to be defined by
725        # the user
726        bounding_polygon_utm=[]
727        for point in bounding_polygon:
728            zone,easting,northing=redfearn(point[0],point[1])
729            bounding_polygon_utm.append([easting,northing])
730
731        boundary_polygon.append(bounding_polygon_utm[3])
732        boundary_polygon.append(bounding_polygon_utm[4])
733
734        #print 'boundary_polygon', boundary_polygon
735       
736        plot=False
737        if plot:
738            from pylab import plot,show,axis
739            boundary_polygon=ensure_numeric(boundary_polygon)
740            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
741            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
742            plot(boundary_polygon[:,0], boundary_polygon[:,1])
743            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1])
744            show()
745
746        assert num.allclose(bounding_polygon_utm,boundary_polygon)
747
748
749        extent_res=1000000
750        meshname = 'urs_test_mesh' + '.tsh'
751        interior_regions=None
752        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
753       
754        # have to change boundary tags from last example because now bounding
755        # polygon starts in different place.
756        create_mesh_from_regions(boundary_polygon,
757                                 boundary_tags=boundary_tags,
758                                 maximum_triangle_area=extent_res,
759                                 filename=meshname,
760                                 interior_regions=interior_regions,
761                                 verbose=False)
762       
763        domain_fbound = Domain(meshname)
764        domain_fbound.set_quantity('stage', tide)
765        Bf = File_boundary(sts_file+'.sts', 
766                           domain_fbound, 
767                           boundary_polygon=boundary_polygon)
768        Br = Reflective_boundary(domain_fbound)
769
770        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
771        finaltime=time_step*(time_step_count-1)
772        yieldstep=time_step
773        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
774   
775        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
776                                                   finaltime=finaltime, 
777                                                   skip_initial_step=False)):
778            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
779   
780       
781        domain_time = Domain(meshname)
782        domain_time.set_quantity('stage', tide)
783        Br = Reflective_boundary(domain_time)
784        Bw = Time_boundary(domain=domain_time,
785                         function=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
786        domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
787       
788        temp_time=num.zeros(int(finaltime/yieldstep)+1,num.float)
789       
790        domain_time.set_starttime(domain_fbound.get_starttime())
791       
792        for i, t in enumerate(domain_time.evolve(yieldstep=yieldstep,
793                                                   finaltime=finaltime, 
794                                                   skip_initial_step=False)):
795            temp_time[i]=domain_time.quantities['stage'].centroid_values[2]
796       
797        assert num.allclose(temp_fbound, temp_time)               
798        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
799                            domain_time.quantities['stage'].vertex_values)
800                       
801        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
802                            domain_time.quantities['xmomentum'].vertex_values)                       
803                       
804        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
805                            domain_time.quantities['ymomentum'].vertex_values)                                               
806       
807
808        try:
809            os.remove(sts_file+'.sts')
810        except:
811            # Windoze can't remove this file for some reason
812            pass
813       
814        os.remove(meshname)
815       
816 
817    def test_points2polygon(self): 
818        att_dict = {}
819        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[0.0, 0.0]])
820        att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
821        att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
822       
823        fileName = tempfile.mktemp(".csv")
824       
825        G = Geospatial_data(pointlist, att_dict)
826       
827        G.export_points_file(fileName)
828       
829        polygon = load_pts_as_polygon(fileName)
830       
831        # This test may fail if the order changes
832        assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])
833       
834
835
836#-------------------------------------------------------------
837
838if __name__ == "__main__":
839    #suite = unittest.makeSuite(Test_Data_Manager, 'test_sww2domain2')
840    suite = unittest.makeSuite(Test_Data_Manager, 'test')
841   
842   
843   
844    # FIXME(Ole): When Ross has implemented logging, we can
845    # probably get rid of all this:
846    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
847        Test_Data_Manager.verbose=True
848        saveout = sys.stdout   
849        filename = ".temp_verbose"
850        fid = open(filename, 'w')
851        sys.stdout = fid
852    else:
853        pass
854    runner = unittest.TextTestRunner() #verbosity=2)
855    runner.run(suite)
856
857    # Cleaning up
858    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
859        sys.stdout = saveout
860        fid.close() 
861        os.remove(filename)
862
863
Note: See TracBrowser for help on using the repository browser.