source: trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py @ 7778

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

Cleaned up unit tests to use API.

File size: 79.1 KB
Line 
1import numpy as num
2import unittest
3import tempfile
4import os
5from Scientific.IO.NetCDF import NetCDFFile
6
7from anuga.utilities.system_tools import get_pathname_from_package
8from anuga.coordinate_transforms.geo_reference import Geo_reference     
9from anuga.coordinate_transforms.redfearn import redfearn
10from anuga.utilities.numerical_tools import ensure_numeric
11from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
12from anuga.file.sts import create_sts_boundary
13from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array
14
15from anuga.shallow_water.shallow_water_domain import Domain
16
17# boundary functions
18from anuga.shallow_water.boundaries import Reflective_boundary, \
19            Field_boundary, Transmissive_momentum_set_stage_boundary, \
20            Transmissive_stage_zero_momentum_boundary
21from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
22     import Transmissive_boundary, Dirichlet_boundary, \
23            Time_boundary, File_boundary, AWI_boundary
24
25from anuga.pmesh.mesh_interface import create_mesh_from_regions
26
27from urs2sts import urs2sts
28
29# Allow us to use helper methods from this test.
30from anuga.file.test_mux import Test_Mux
31
32class Test_Urs2Sts(Test_Mux):
33    """ A suite of tests to test urs2sts file conversion functions.
34        These tests are quite coarse-grained: converting a file
35        and checking that its headers and some of its contents
36        are correct.
37    """ 
38
39    def test_urs2sts0(self):
40        """
41        Test single source
42        """
43        tide=0
44        time_step_count = 3
45        time_step = 2
46        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
47        n=len(lat_long_points)
48        first_tstep=num.ones(n,num.int)
49        first_tstep[0]+=1
50        first_tstep[2]+=1
51        last_tstep=(time_step_count)*num.ones(n,num.int)
52        last_tstep[0]-=1
53
54        gauge_depth=20*num.ones(n,num.float)
55        ha=2*num.ones((n,time_step_count),num.float)
56        ha[0]=num.arange(0,time_step_count)
57        ha[1]=num.arange(time_step_count,2*time_step_count)
58        ha[2]=num.arange(2*time_step_count,3*time_step_count)
59        ha[3]=num.arange(3*time_step_count,4*time_step_count)
60        ua=5*num.ones((n,time_step_count),num.float)
61        va=-10*num.ones((n,time_step_count),num.float)
62
63        base_name, files = self.write_mux2(lat_long_points,
64                                      time_step_count, time_step,
65                                      first_tstep, last_tstep,
66                                      depth=gauge_depth,
67                                      ha=ha,
68                                      ua=ua,
69                                      va=va)
70
71        urs2sts(base_name,
72                basename_out=base_name, 
73                mean_stage=tide,verbose=False)
74
75        # now I want to check the sts file ...
76        sts_file = base_name + '.sts'
77
78        #Let's interigate the sww file
79        # Note, the sww info is not gridded.  It is point data.
80        fid = NetCDFFile(sts_file)
81
82        # Make x and y absolute
83        x = fid.variables['x'][:]
84        y = fid.variables['y'][:]
85
86        geo_reference = Geo_reference(NetCDFObject=fid)
87        points = geo_reference.get_absolute(map(None, x, y))
88        points = ensure_numeric(points)
89
90        x = points[:,0]
91        y = points[:,1]
92
93        #Check that first coordinate is correctly represented       
94        #Work out the UTM coordinates for first point
95        for i in range(4):
96            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1]) 
97            assert num.allclose([x[i],y[i]], [e,n])
98
99        #Check the time vector
100        times = fid.variables['time'][:]
101
102        times_actual = []
103        for i in range(time_step_count):
104            times_actual.append(time_step * i)
105
106        assert num.allclose(ensure_numeric(times),
107                            ensure_numeric(times_actual))
108
109        #Check first value
110        stage = fid.variables['stage'][:]
111        xmomentum = fid.variables['xmomentum'][:]
112        ymomentum = fid.variables['ymomentum'][:]
113        elevation = fid.variables['elevation'][:]
114
115        # Set original data used to write mux file to be zero when gauges are
116        #not recdoring
117        ha[0][0]=0.0
118        ha[0][time_step_count-1]=0.0;
119        ha[2][0]=0.0;
120        ua[0][0]=0.0
121        ua[0][time_step_count-1]=0.0;
122        ua[2][0]=0.0;
123        va[0][0]=0.0
124        va[0][time_step_count-1]=0.0;
125        va[2][0]=0.0;
126
127        assert num.allclose(num.transpose(ha),stage)  #Meters
128
129        #Check the momentums - ua
130        #momentum = velocity*(stage-elevation)
131        # elevation = - depth
132        #momentum = velocity_ua *(stage+depth)
133
134        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
135        for i in range(len(lat_long_points)):
136            depth[i]=gauge_depth[i]+tide+ha[i]
137        assert num.allclose(num.transpose(ua*depth),xmomentum) 
138
139        #Check the momentums - va
140        #momentum = velocity*(stage-elevation)
141        # elevation = - depth
142        #momentum = velocity_va *(stage+depth)
143
144        assert num.allclose(num.transpose(va*depth),ymomentum)
145
146        # check the elevation values.
147        # -ve since urs measures depth, sww meshers height,
148        assert num.allclose(-elevation, gauge_depth)  #Meters
149
150        fid.close()
151        self.delete_mux(files)
152        os.remove(sts_file)
153
154    def test_urs2sts_nonstandard_meridian(self):
155        """
156        Test single source using the meridian from zone 50 as a nonstandard meridian
157        """
158        tide=0
159        time_step_count = 3
160        time_step = 2
161        lat_long_points =[(-21.,114.5),(-21.,113.5),(-21.,114.), (-21.,115.)]
162        n=len(lat_long_points)
163        first_tstep=num.ones(n,num.int)
164        first_tstep[0]+=1
165        first_tstep[2]+=1
166        last_tstep=(time_step_count)*num.ones(n,num.int)
167        last_tstep[0]-=1
168
169        gauge_depth=20*num.ones(n,num.float)
170        ha=2*num.ones((n,time_step_count),num.float)
171        ha[0]=num.arange(0,time_step_count)
172        ha[1]=num.arange(time_step_count,2*time_step_count)
173        ha[2]=num.arange(2*time_step_count,3*time_step_count)
174        ha[3]=num.arange(3*time_step_count,4*time_step_count)
175        ua=5*num.ones((n,time_step_count),num.float)
176        va=-10*num.ones((n,time_step_count),num.float)
177
178        base_name, files = self.write_mux2(lat_long_points,
179                                           time_step_count, time_step,
180                                           first_tstep, last_tstep,
181                                           depth=gauge_depth,
182                                           ha=ha,
183                                           ua=ua,
184                                           va=va)
185
186        urs2sts(base_name,
187                basename_out=base_name, 
188                central_meridian=123,
189                mean_stage=tide,
190                verbose=False)
191
192        # now I want to check the sts file ...
193        sts_file = base_name + '.sts'
194
195        #Let's interigate the sww file
196        # Note, the sww info is not gridded.  It is point data.
197        fid = NetCDFFile(sts_file)
198
199        # Make x and y absolute
200        x = fid.variables['x'][:]
201        y = fid.variables['y'][:]
202
203        geo_reference = Geo_reference(NetCDFObject=fid)
204        points = geo_reference.get_absolute(map(None, x, y))
205        points = ensure_numeric(points)
206
207        x = points[:,0]
208        y = points[:,1]
209
210        # Check that all coordinate are correctly represented       
211        # Using the non standard projection (50)
212        for i in range(4):
213            zone, e, n = redfearn(lat_long_points[i][0],
214                                  lat_long_points[i][1],
215                                  central_meridian=123)
216            assert num.allclose([x[i],y[i]], [e,n])
217            assert zone==-1
218        self.delete_mux(files)
219       
220
221
222    def test_urs2sts_individual_sources(self):   
223        """Test that individual sources compare to actual urs output
224           Test that the first recording time is the smallest
225           over waveheight, easting and northing velocity
226        """
227       
228        # Get path where this test is run
229        path = get_pathname_from_package('anuga.shallow_water')       
230       
231        testdir = os.path.join(path, 'urs_test_data')
232        ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt')
233       
234        sources = ['1-z.grd','2-z.grd','3-z.grd']
235       
236        # Start times by source and station taken manually from urs header files
237        time_start_z = num.array([[10.0,11.5,13,14.5,17.7],
238                                  [9.8,11.2,12.7,14.2,17.4],
239                                  [9.5,10.9,12.4,13.9,17.1]])
240
241        time_start_e = time_start_n = time_start_z
242
243        # time step in urs output
244        delta_t = 0.1
245       
246        # Make sts file for each source
247        for k, source_filename in enumerate(sources):
248            source_number = k + 1 # Source numbering starts at 1
249           
250            urs_filenames = os.path.join(testdir, source_filename)
251            weights = [1.]
252            sts_name_out = 'test'
253           
254            urs2sts(urs_filenames,
255                    basename_out=sts_name_out,
256                    ordering_filename=ordering_filename,
257                    weights=weights,
258                    mean_stage=0.0,
259                    verbose=False)
260
261            # Read in sts file for this source file
262            fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r) # Open existing file for read
263            x = fid.variables['x'][:]+fid.xllcorner    # x-coordinates of vertices
264            y = fid.variables['y'][:]+fid.yllcorner    # y-coordinates of vertices
265            elevation = fid.variables['elevation'][:]
266            time=fid.variables['time'][:]+fid.starttime
267
268            # Get quantity data from sts file
269            quantity_names=['stage','xmomentum','ymomentum']
270            quantities = {}
271            for i, name in enumerate(quantity_names):
272                quantities[name] = fid.variables[name][:]
273           
274            # Make sure start time from sts file is the minimum starttime
275            # across all stations (for this source)
276            #print k, time_start_z[k,:]
277            starttime = min(time_start_z[k, :])
278            sts_starttime = fid.starttime[0]
279            msg = 'sts starttime for source %d was %f. Should have been %f'\
280                %(source_number, sts_starttime, starttime)
281            assert num.allclose(sts_starttime, starttime), msg             
282
283            # For each station, compare urs2sts output to known urs output
284            for j in range(len(x)):
285                index_start_urs = 0
286                index_end_urs = 0
287                index_start = 0
288                index_end = 0
289                count = 0
290
291                # read in urs test data for stage, e and n velocity
292                urs_file_name_z = 'z_'+str(source_number)+'_'+str(j)+'.csv'
293                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
294                urs_stage = dict['urs_stage']
295                urs_file_name_e = 'e_'+str(source_number)+'_'+str(j)+'.csv'
296                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
297                urs_e = dict['urs_e']
298                urs_file_name_n = 'n_'+str(source_number)+'_'+str(j)+'.csv'
299                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
300                urs_n = dict['urs_n']
301
302                # find start and end time for stage             
303                for i in range(len(urs_stage)):
304                    if urs_stage[i] == 0.0:
305                        index_start_urs_z = i+1
306                    if int(urs_stage[i]) == 99 and count <> 1:
307                        count +=1
308                        index_end_urs_z = i
309
310                if count == 0: index_end_urs_z = len(urs_stage)
311
312                start_times_z = time_start_z[source_number-1]
313
314                # start times for easting velocities should be the same as stage
315                start_times_e = time_start_e[source_number-1]
316                index_start_urs_e = index_start_urs_z
317                index_end_urs_e = index_end_urs_z
318
319                # start times for northing velocities should be the same as stage
320                start_times_n = time_start_n[source_number-1]
321                index_start_urs_n = index_start_urs_z
322                index_end_urs_n = index_end_urs_z
323               
324                # Check that actual start time matches header information for stage
325                msg = 'stage start time from urs file is not the same as the '
326                msg += 'header file for source %i and station %i' %(source_number,j)
327                assert num.allclose(index_start_urs_z,start_times_z[j]/delta_t), msg
328
329                msg = 'e velocity start time from urs file is not the same as the '
330                msg += 'header file for source %i and station %i' %(source_number,j)
331                assert num.allclose(index_start_urs_e,start_times_e[j]/delta_t), msg
332
333                msg = 'n velocity start time from urs file is not the same as the '
334                msg += 'header file for source %i and station %i' %(source_number,j)
335                assert num.allclose(index_start_urs_n,start_times_n[j]/delta_t), msg
336               
337                # get index for start and end time for sts quantities
338                index_start_stage = 0
339                index_end_stage = 0
340                count = 0
341                sts_stage = quantities['stage'][:,j]
342                for i in range(len(sts_stage)):
343                    if sts_stage[i] <> 0.0 and count <> 1:
344                        count += 1
345                        index_start_stage = i
346                    if int(sts_stage[i]) == 99 and count <> 1:
347                        count += 1
348                        index_end_stage = i
349
350                index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
351
352                sts_xmom = quantities['xmomentum'][:,j]
353                index_start_x = index_start_stage
354                index_end_x = index_start_x + len(urs_e[index_start_urs_e:index_end_urs_e])
355
356                sts_ymom = quantities['ymomentum'][:,j]
357                index_start_y = index_start_stage
358                index_end_y = index_start_y + len(urs_n[index_start_urs_n:index_end_urs_n])
359
360                # check that urs stage and sts stage are the same
361                msg = 'urs stage is not equal to sts stage for for source %i and station %i' %(source_number,j)
362                assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
363                                    sts_stage[index_start_stage:index_end_stage], 
364                                    rtol=1.0e-6, atol=1.0e-5 ), msg                               
365                               
366                # check that urs e velocity and sts xmomentum are the same
367                msg = 'urs e velocity is not equivalent to sts x momentum for for source %i and station %i' %(source_number,j)
368                assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]),
369                                sts_xmom[index_start_x:index_end_x], 
370                                rtol=1.0e-5, atol=1.0e-4 ), msg
371               
372                # check that urs n velocity and sts ymomentum are the same
373                #print 'urs n velocity', urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j])
374                #print 'sts momentum', sts_ymom[index_start_y:index_end_y]                                                             
375                msg = 'urs n velocity is not equivalent to sts y momentum for source %i and station %i' %(source_number,j)
376                assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
377                                -sts_ymom[index_start_y:index_end_y], 
378                                rtol=1.0e-5, atol=1.0e-4 ), msg
379                                               
380                               
381            fid.close()
382           
383        os.remove(sts_name_out+'.sts')
384
385    def test_urs2sts_combined_sources(self):   
386        """Test that combined sources compare to actual urs output
387           Test that the first recording time is the smallest
388           over waveheight, easting and northing velocity
389        """
390
391        # combined
392        time_start_z = num.array([9.5,10.9,12.4,13.9,17.1])
393        time_start_e = time_start_n = time_start_z
394         
395        # make sts file for combined sources
396        weights = [1., 2., 3.]
397       
398        path = get_pathname_from_package('anuga.shallow_water')       
399               
400        testdir = os.path.join(path, 'urs_test_data')       
401        ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt')
402
403        urs_filenames = [os.path.join(testdir,'1-z.grd'),
404                         os.path.join(testdir,'2-z.grd'),
405                         os.path.join(testdir,'3-z.grd')]
406        sts_name_out = 'test'
407       
408        urs2sts(urs_filenames,
409                basename_out=sts_name_out,
410                ordering_filename=ordering_filename,
411                weights=weights,
412                mean_stage=0.0,
413                verbose=False)
414       
415        # read in sts file for combined source
416        fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r)    # Open existing file for read
417        x = fid.variables['x'][:]+fid.xllcorner   # x-coordinates of vertices
418        y = fid.variables['y'][:]+fid.yllcorner   # y-coordinates of vertices
419        elevation = fid.variables['elevation'][:]
420        time=fid.variables['time'][:]+fid.starttime
421       
422       
423        # Check that stored permutation is as per default
424        permutation = range(len(x))
425        stored_permutation = fid.variables['permutation'][:]
426        msg = 'Permutation was not stored correctly. I got '
427        msg += str(stored_permutation)
428        assert num.allclose(stored_permutation, permutation), msg       
429
430        # Get quantity data from sts file
431        quantity_names=['stage','xmomentum','ymomentum']
432        quantities = {}
433        for i, name in enumerate(quantity_names):
434            quantities[name] = fid.variables[name][:]
435
436        # For each station, compare urs2sts output to known urs output   
437        delta_t = 0.1
438       
439        # Make sure start time from sts file is the minimum starttime
440        # across all stations (for this source)
441        starttime = min(time_start_z[:])
442        sts_starttime = fid.starttime[0]
443        msg = 'sts starttime was %f. Should have been %f'\
444            %(sts_starttime, starttime)
445        assert num.allclose(sts_starttime, starttime), msg
446   
447        #stations = [1,2,3]
448        #for j in stations:
449        for j in range(len(x)):
450            index_start_urs_z = 0
451            index_end_urs_z = 0
452            index_start_urs_e = 0
453            index_end_urs_e = 0
454            index_start_urs_n = 0
455            index_end_urs_n = 0
456            count = 0
457
458            # read in urs test data for stage, e and n velocity
459            urs_file_name_z = 'z_combined_'+str(j)+'.csv'
460            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
461            urs_stage = dict['urs_stage']
462            urs_file_name_e = 'e_combined_'+str(j)+'.csv'
463            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
464            urs_e = dict['urs_e']
465            urs_file_name_n = 'n_combined_'+str(j)+'.csv'
466            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
467            urs_n = dict['urs_n']
468
469            # find start and end time for stage         
470            for i in range(len(urs_stage)):
471                if urs_stage[i] == 0.0:
472                    index_start_urs_z = i+1
473                if int(urs_stage[i]) == 99 and count <> 1:
474                    count +=1
475                    index_end_urs_z = i
476
477            if count == 0: index_end_urs_z = len(urs_stage)
478
479            start_times_z = time_start_z[j]
480
481            start_times_e = time_start_e[j]
482            index_start_urs_e = index_start_urs_z
483
484            start_times_n = time_start_n[j]
485            index_start_urs_n = index_start_urs_z
486               
487            # Check that actual start time matches header information for stage
488            msg = 'stage start time from urs file is not the same as the '
489            msg += 'header file at station %i' %(j)
490            assert num.allclose(index_start_urs_z,start_times_z/delta_t), msg
491
492            msg = 'e velocity start time from urs file is not the same as the '
493            msg += 'header file at station %i' %(j)
494            assert num.allclose(index_start_urs_e,start_times_e/delta_t), msg
495
496            msg = 'n velocity start time from urs file is not the same as the '
497            msg += 'header file at station %i' %(j)
498            assert num.allclose(index_start_urs_n,start_times_n/delta_t), msg
499               
500            # get index for start and end time for sts quantities
501            index_start_stage = 0
502            index_end_stage = 0
503            index_start_x = 0
504            index_end_x = 0
505            index_start_y = 0
506            index_end_y = 0
507            count = 0
508            count1 = 0
509            sts_stage = quantities['stage'][:,j]
510            for i in range(len(sts_stage)):
511                if sts_stage[i] <> 0.0 and count <> 1:
512                    count += 1
513                    index_start_stage = i
514                if int(urs_stage[i]) == 99 and count <> 1:
515                    count +=1
516                    index_end_stage = i
517               
518            index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
519
520            index_start_x = index_start_stage
521            index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e])
522            sts_xmom = quantities['ymomentum'][:,j]
523
524            index_start_y = index_start_stage
525            index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n])
526            sts_ymom = quantities['ymomentum'][:,j]
527
528            # check that urs stage and sts stage are the same
529            msg = 'urs stage is not equal to sts stage for station %i' %j
530            #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z]
531            #print 'sts stage', sts_stage[index_start_stage:index_end_stage]
532            #print 'diff', max(urs_stage[index_start_urs_z:index_end_urs_z]-sts_stage[index_start_stage:index_end_stage])
533            #print 'index', index_start_stage, index_end_stage, len(sts_stage)
534            assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
535                            sts_stage[index_start_stage:index_end_stage], 
536                                rtol=1.0e-5, atol=1.0e-4 ), msg                               
537                               
538            # check that urs e velocity and sts xmomentum are the same         
539            msg = 'urs e velocity is not equivalent to sts xmomentum for station %i' %j
540            assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]),
541                            sts_xmom[index_start_x:index_end_x], 
542                            rtol=1.0e-5, atol=1.0e-4 ), msg
543               
544            # check that urs n velocity and sts ymomentum are the same                           
545            msg = 'urs n velocity is not equivalent to sts ymomentum for station %i' %j
546            assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
547                            sts_ymom[index_start_y:index_end_y], 
548                            rtol=1.0e-5, atol=1.0e-4 ), msg
549
550        fid.close()
551       
552        os.remove(sts_name_out+'.sts')
553       
554       
555       
556    def test_urs2sts_ordering(self):
557        """Test multiple sources with ordering file
558        """
559       
560        tide = 0.35
561        time_step_count = 6 # I made this a little longer (Ole)
562        time_step = 2
563        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
564        n=len(lat_long_points)
565        first_tstep=num.ones(n,num.int)
566        first_tstep[0]+=1
567        first_tstep[2]+=1
568        last_tstep=(time_step_count)*num.ones(n,num.int)
569        last_tstep[0]-=1
570
571        gauge_depth=20*num.ones(n,num.float)
572        ha=2*num.ones((n,time_step_count),num.float)
573        ha[0]=num.arange(0,time_step_count)
574        ha[1]=num.arange(time_step_count,2*time_step_count)
575        ha[2]=num.arange(2*time_step_count,3*time_step_count)
576        ha[3]=num.arange(3*time_step_count,4*time_step_count)
577        ua=5*num.ones((n,time_step_count),num.float)
578        va=-10*num.ones((n,time_step_count),num.float)
579
580        # Create two identical mux files to be combined by urs2sts
581        base_nameI, filesI = self.write_mux2(lat_long_points,
582                                             time_step_count, time_step,
583                                             first_tstep, last_tstep,
584                                             depth=gauge_depth,
585                                             ha=ha,
586                                             ua=ua,
587                                             va=va)
588
589        base_nameII, filesII = self.write_mux2(lat_long_points,
590                                               time_step_count, time_step,
591                                               first_tstep, last_tstep,
592                                               depth=gauge_depth,
593                                               ha=ha,
594                                               ua=ua,
595                                               va=va)
596
597                                               
598        # Create ordering file
599        permutation = [3,0,2]
600
601        _, ordering_filename = tempfile.mkstemp('')
602        order_fid = open(ordering_filename, 'w') 
603        order_fid.write('index, longitude, latitude\n')
604        for index in permutation:
605            order_fid.write('%d, %f, %f\n' %(index, 
606                                             lat_long_points[index][1], 
607                                             lat_long_points[index][0]))
608        order_fid.close()
609       
610           
611
612                                               
613        # Call urs2sts with multiple mux files
614        urs2sts([base_nameI, base_nameII], 
615                basename_out=base_nameI, 
616                ordering_filename=ordering_filename,
617                weights=[1.0, 1.0],
618                mean_stage=tide,
619                verbose=False)
620
621        # now I want to check the sts file ...
622        sts_file = base_nameI + '.sts'
623
624        #Let's interrogate the sts file
625        # Note, the sts info is not gridded.  It is point data.
626        fid = NetCDFFile(sts_file)
627       
628        # Check that original indices have been stored
629        stored_permutation = fid.variables['permutation'][:]
630        msg = 'Permutation was not stored correctly. I got '
631        msg += str(stored_permutation)
632        assert num.allclose(stored_permutation, permutation), msg
633       
634
635        # Make x and y absolute
636        x = fid.variables['x'][:]
637        y = fid.variables['y'][:]
638
639        geo_reference = Geo_reference(NetCDFObject=fid)
640        points = geo_reference.get_absolute(map(None, x, y))
641        points = ensure_numeric(points)
642
643        x = points[:,0]
644        y = points[:,1]
645
646        #print
647        #print x
648        #print y
649        for i, index in enumerate(permutation):
650            # Check that STS points are stored in the correct order
651           
652            # Work out the UTM coordinates sts point i
653            zone, e, n = redfearn(lat_long_points[index][0], 
654                                  lat_long_points[index][1])             
655
656            #print i, [x[i],y[i]], [e,n]
657            assert num.allclose([x[i],y[i]], [e,n])
658           
659                       
660        # Check the time vector
661        times = fid.variables['time'][:]
662
663        times_actual = []
664        for i in range(time_step_count):
665            times_actual.append(time_step * i)
666
667        assert num.allclose(ensure_numeric(times),
668                            ensure_numeric(times_actual))
669                       
670
671        # Check sts values
672        stage = fid.variables['stage'][:]
673        xmomentum = fid.variables['xmomentum'][:]
674        ymomentum = fid.variables['ymomentum'][:]
675        elevation = fid.variables['elevation'][:]
676
677        # Set original data used to write mux file to be zero when gauges are
678        # not recdoring
679       
680        ha[0][0]=0.0
681        ha[0][time_step_count-1]=0.0
682        ha[2][0]=0.0
683        ua[0][0]=0.0
684        ua[0][time_step_count-1]=0.0
685        ua[2][0]=0.0
686        va[0][0]=0.0
687        va[0][time_step_count-1]=0.0
688        va[2][0]=0.0;
689
690       
691        # The stage stored in the .sts file should be the sum of the stage
692        # in the two mux2 files because both have weights = 1. In this case
693        # the mux2 files are the same so stage == 2.0 * ha
694        #print 2.0*num.transpose(ha) - stage
695       
696        ha_permutation = num.take(ha, permutation, axis=0) 
697        ua_permutation = num.take(ua, permutation, axis=0)
698        va_permutation = num.take(va, permutation, axis=0)
699        gauge_depth_permutation = num.take(gauge_depth, permutation, axis=0)
700
701       
702        assert num.allclose(2.0*num.transpose(ha_permutation)+tide, stage)  # Meters
703
704        #Check the momentums - ua
705        #momentum = velocity*(stage-elevation)
706        # elevation = - depth
707        #momentum = velocity_ua *(stage+depth)
708
709        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
710        for i in range(len(lat_long_points)):
711            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
712            #2.0*ha necessary because using two files with weights=1 are used
713           
714        depth_permutation = num.take(depth, permutation, axis=0)                     
715       
716
717        # The xmomentum stored in the .sts file should be the sum of the ua
718        # in the two mux2 files multiplied by the depth.
719        assert num.allclose(2.0*num.transpose(ua_permutation*depth_permutation), xmomentum) 
720
721        #Check the momentums - va
722        #momentum = velocity*(stage-elevation)
723        # elevation = - depth
724        #momentum = velocity_va *(stage+depth)
725
726        # The ymomentum stored in the .sts file should be the sum of the va
727        # in the two mux2 files multiplied by the depth.
728        assert num.allclose(2.0*num.transpose(va_permutation*depth_permutation), ymomentum)
729
730        # check the elevation values.
731        # -ve since urs measures depth, sww meshers height,
732        assert num.allclose(-gauge_depth_permutation, elevation)  #Meters
733
734        fid.close()
735        self.delete_mux(filesI)
736        self.delete_mux(filesII)
737        os.remove(sts_file)
738       
739
740       
741       
742       
743    def Xtest_urs2sts_ordering_exception(self):
744        """Test that inconsistent lats and lons in ordering file are caught.
745        """
746       
747        tide=0
748        time_step_count = 3
749        time_step = 2
750        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
751        n=len(lat_long_points)
752        first_tstep=num.ones(n,num.int)
753        first_tstep[0]+=1
754        first_tstep[2]+=1
755        last_tstep=(time_step_count)*num.ones(n,num.int)
756        last_tstep[0]-=1
757
758        gauge_depth=20*num.ones(n,num.float)
759        ha=2*num.ones((n,time_step_count),num.float)
760        ha[0]=num.arange(0,time_step_count)
761        ha[1]=num.arange(time_step_count,2*time_step_count)
762        ha[2]=num.arange(2*time_step_count,3*time_step_count)
763        ha[3]=num.arange(3*time_step_count,4*time_step_count)
764        ua=5*num.ones((n,time_step_count),num.float)
765        va=-10*num.ones((n,time_step_count),num.float)
766
767        # Create two identical mux files to be combined by urs2sts
768        base_nameI, filesI = self.write_mux2(lat_long_points,
769                                             time_step_count, time_step,
770                                             first_tstep, last_tstep,
771                                             depth=gauge_depth,
772                                             ha=ha,
773                                             ua=ua,
774                                             va=va)
775
776        base_nameII, filesII = self.write_mux2(lat_long_points,
777                                               time_step_count, time_step,
778                                               first_tstep, last_tstep,
779                                               depth=gauge_depth,
780                                               ha=ha,
781                                               ua=ua,
782                                               va=va)
783
784                                               
785        # Create ordering file
786        permutation = [3,0,2]
787
788        # Do it wrongly and check that exception is being raised
789        _, ordering_filename = tempfile.mkstemp('')
790        order_fid = open(ordering_filename, 'w') 
791        order_fid.write('index, longitude, latitude\n')
792        for index in permutation:
793            order_fid.write('%d, %f, %f\n' %(index, 
794                                             lat_long_points[index][0], 
795                                             lat_long_points[index][1]))
796        order_fid.close()
797       
798        try:
799            urs2sts([base_nameI, base_nameII], 
800                    basename_out=base_nameI, 
801                    ordering_filename=ordering_filename,
802                    weights=[1.0, 1.0],
803                    mean_stage=tide,
804                    verbose=False) 
805            os.remove(ordering_filename)           
806        except:
807            pass
808        else:
809            msg = 'Should have caught wrong lat longs'
810            raise Exception, msg
811
812       
813        self.delete_mux(filesI)
814        self.delete_mux(filesII)
815
816       
817
818       
819    def test_urs2sts_ordering_different_sources(self):
820        """Test multiple sources with ordering file, different source files and weights.
821           This test also has more variable values than the previous ones
822        """
823       
824        tide = 1.5
825        time_step_count = 10
826        time_step = 0.2
827       
828        times_ref = num.arange(0, time_step_count*time_step, time_step)
829        #print 'time vector', times_ref
830       
831        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
832        n = len(lat_long_points)
833       
834        # Create non-trivial weights
835        #weights = [0.8, 1.5] # OK
836        #weights = [0.8, 10.5] # Fail (up to allclose tolerance)
837        #weights = [10.5, 10.5] # OK
838        #weights = [0.0, 10.5] # OK
839        #weights = [0.8, 0.] # OK               
840        #weights = [8, 0.1] # OK                       
841        #weights = [0.8, 10.0] # OK                               
842        #weights = [0.8, 10.6] # OK           
843        weights = [3.8, 7.6] # OK                   
844        #weights = [0.5, 0.5] # OK                           
845        #weights = [2., 2.] # OK                           
846        #weights = [0.0, 0.5] # OK                                         
847        #weights = [1.0, 1.0] # OK                                                 
848       
849       
850        # Create different timeseries starting and ending at different times
851        first_tstep=num.ones(n,num.int)
852        first_tstep[0]+=2   # Point 0 starts at 2
853        first_tstep[1]+=4   # Point 1 starts at 4       
854        first_tstep[2]+=3   # Point 2 starts at 3
855       
856        last_tstep=(time_step_count)*num.ones(n,num.int)
857        last_tstep[0]-=1    # Point 0 ends 1 step early
858        last_tstep[1]-=2    # Point 1 ends 2 steps early               
859        last_tstep[4]-=3    # Point 4 ends 3 steps early       
860       
861        #print
862        #print 'time_step_count', time_step_count
863        #print 'time_step', time_step
864        #print 'first_tstep', first_tstep
865        #print 'last_tstep', last_tstep               
866       
867       
868        # Create varying elevation data (positive values for seafloor)
869        gauge_depth=20*num.ones(n,num.float)
870        for i in range(n):
871            gauge_depth[i] += i**2
872           
873        #print 'gauge_depth', gauge_depth
874       
875        # Create data to be written to first mux file       
876        ha0=2*num.ones((n,time_step_count),num.float)
877        ha0[0]=num.arange(0,time_step_count)
878        ha0[1]=num.arange(time_step_count,2*time_step_count)
879        ha0[2]=num.arange(2*time_step_count,3*time_step_count)
880        ha0[3]=num.arange(3*time_step_count,4*time_step_count)
881        ua0=5*num.ones((n,time_step_count),num.float)
882        va0=-10*num.ones((n,time_step_count),num.float)
883
884        # Ensure data used to write mux file to be zero when gauges are
885        # not recording
886        for i in range(n):
887             # For each point
888             
889             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
890                 # For timesteps before and after recording range
891                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                 
892
893
894                 
895        #print
896        #print 'using varying start and end time'
897        #print 'ha0', ha0[4]
898        #print 'ua0', ua0
899        #print 'va0', va0       
900       
901        # Write first mux file to be combined by urs2sts
902        base_nameI, filesI = self.write_mux2(lat_long_points,
903                                             time_step_count, time_step,
904                                             first_tstep, last_tstep,
905                                             depth=gauge_depth,
906                                             ha=ha0,
907                                             ua=ua0,
908                                             va=va0)
909
910                                             
911                                             
912        # Create data to be written to second mux file       
913        ha1=num.ones((n,time_step_count),num.float)
914        ha1[0]=num.sin(times_ref)
915        ha1[1]=2*num.sin(times_ref - 3)
916        ha1[2]=5*num.sin(4*times_ref)
917        ha1[3]=num.sin(times_ref)
918        ha1[4]=num.sin(2*times_ref-0.7)
919               
920        ua1=num.zeros((n,time_step_count),num.float)
921        ua1[0]=3*num.cos(times_ref)       
922        ua1[1]=2*num.sin(times_ref-0.7)   
923        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
924        ua1[4]=2*num.ones(time_step_count)
925       
926        va1=num.zeros((n,time_step_count),num.float)
927        va1[0]=2*num.cos(times_ref-0.87)       
928        va1[1]=3*num.ones(time_step_count, num.int)       #array default#
929        va1[3]=2*num.sin(times_ref-0.71)       
930       
931       
932        # Ensure data used to write mux file to be zero when gauges are
933        # not recording
934        for i in range(n):
935             # For each point
936             
937             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
938                 # For timesteps before and after recording range
939                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                 
940
941
942        #print
943        #print 'using varying start and end time'
944        #print 'ha1', ha1[4]
945        #print 'ua1', ua1
946        #print 'va1', va1       
947                                             
948                                             
949        # Write second mux file to be combined by urs2sts                                             
950        base_nameII, filesII = self.write_mux2(lat_long_points,
951                                               time_step_count, time_step,
952                                               first_tstep, last_tstep,
953                                               depth=gauge_depth,
954                                               ha=ha1,
955                                               ua=ua1,
956                                               va=va1)
957
958                                               
959        # Create ordering file
960        permutation = [4,0,2]
961
962        _, ordering_filename = tempfile.mkstemp('')
963        order_fid = open(ordering_filename, 'w') 
964        order_fid.write('index, longitude, latitude\n')
965        for index in permutation:
966            order_fid.write('%d, %f, %f\n' %(index, 
967                                             lat_long_points[index][1], 
968                                             lat_long_points[index][0]))
969        order_fid.close()
970       
971           
972
973        #------------------------------------------------------------
974        # Now read the mux files one by one without weights and test
975       
976        # Call urs2sts with mux file #0
977        urs2sts([base_nameI], 
978                basename_out=base_nameI, 
979                ordering_filename=ordering_filename,
980                mean_stage=tide,
981                verbose=False)
982
983        # Now read the sts file and check that values have been stored correctly.
984        sts_file = base_nameI + '.sts'
985        fid = NetCDFFile(sts_file)
986       
987
988        # Check that original indices have been stored
989        stored_permutation = fid.variables['permutation'][:]
990        msg = 'Permutation was not stored correctly. I got '
991        msg += str(stored_permutation)
992        assert num.allclose(stored_permutation, permutation), msg
993       
994
995       
996       
997        # Make x and y absolute
998        x = fid.variables['x'][:]
999        y = fid.variables['y'][:]
1000
1001        geo_reference = Geo_reference(NetCDFObject=fid)
1002        points = geo_reference.get_absolute(map(None, x, y))
1003        points = ensure_numeric(points)
1004
1005        x = points[:,0]
1006        y = points[:,1]
1007
1008        for i, index in enumerate(permutation):
1009            # Check that STS points are stored in the correct order
1010           
1011            # Work out the UTM coordinates sts point i
1012            zone, e, n = redfearn(lat_long_points[index][0], 
1013                                  lat_long_points[index][1])             
1014
1015            #print i, [x[i],y[i]], [e,n]
1016            assert num.allclose([x[i],y[i]], [e,n])
1017           
1018                       
1019        # Check the time vector
1020        times = fid.variables['time'][:]
1021        assert num.allclose(ensure_numeric(times),
1022                            ensure_numeric(times_ref))
1023                       
1024
1025        # Check sts values for mux #0
1026        stage0 = fid.variables['stage'][:].copy()
1027        xmomentum0 = fid.variables['xmomentum'][:].copy()
1028        ymomentum0 = fid.variables['ymomentum'][:].copy()
1029        elevation0 = fid.variables['elevation'][:].copy()
1030
1031       
1032        #print 'stage', stage0
1033        #print 'xmomentum', xmomentum0
1034        #print 'ymomentum', ymomentum0       
1035        #print 'elevation', elevation0
1036       
1037        # The quantities stored in the .sts file should be the weighted sum of the
1038        # quantities written to the mux2 files subject to the permutation vector.
1039       
1040        ha_ref = num.take(ha0, permutation, axis=0)
1041        ua_ref = num.take(ua0, permutation, axis=0)       
1042        va_ref = num.take(va0, permutation, axis=0)               
1043
1044        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                     
1045       
1046        assert num.allclose(num.transpose(ha_ref)+tide, stage0)  # Meters
1047       
1048       
1049       
1050        #Check the momentums - ua
1051        #momentum = velocity*(stage-elevation)
1052        # elevation = - depth
1053        #momentum = velocity_ua *(stage+depth)
1054
1055        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
1056        for i in range(len(permutation)):
1057            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
1058
1059
1060        # The xmomentum stored in the .sts file should be the sum of the ua
1061        # in the two mux2 files multiplied by the depth.
1062        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum0) 
1063
1064        #Check the momentums - va
1065        #momentum = velocity*(stage-elevation)
1066        # elevation = - depth
1067        #momentum = velocity_va *(stage+depth)
1068
1069        # The ymomentum stored in the .sts file should be the sum of the va
1070        # in the two mux2 files multiplied by the depth.
1071       
1072       
1073        #print transpose(va_ref*depth_ref)
1074        #print ymomentum
1075        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum0)       
1076
1077        # check the elevation values.
1078        # -ve since urs measures depth, sww meshers height,
1079        assert num.allclose(-gauge_depth_ref, elevation0) 
1080
1081        fid.close()
1082        os.remove(sts_file)
1083       
1084       
1085
1086       
1087        # Call urs2sts with mux file #1
1088        urs2sts([base_nameII], 
1089                basename_out=base_nameI, 
1090                ordering_filename=ordering_filename,
1091                mean_stage=tide,
1092                verbose=False)
1093
1094        # Now read the sts file and check that values have been stored correctly.
1095        sts_file = base_nameI + '.sts'
1096        fid = NetCDFFile(sts_file)
1097       
1098       
1099        # Check that original indices have been stored
1100        stored_permutation = fid.variables['permutation'][:]
1101        msg = 'Permutation was not stored correctly. I got '
1102        msg += str(stored_permutation)
1103        assert num.allclose(stored_permutation, permutation), msg
1104       
1105        # Make x and y absolute
1106        x = fid.variables['x'][:]
1107        y = fid.variables['y'][:]
1108
1109        geo_reference = Geo_reference(NetCDFObject=fid)
1110        points = geo_reference.get_absolute(map(None, x, y))
1111        points = ensure_numeric(points)
1112
1113        x = points[:,0]
1114        y = points[:,1]
1115
1116        for i, index in enumerate(permutation):
1117            # Check that STS points are stored in the correct order
1118           
1119            # Work out the UTM coordinates sts point i
1120            zone, e, n = redfearn(lat_long_points[index][0], 
1121                                  lat_long_points[index][1])             
1122
1123            #print i, [x[i],y[i]], [e,n]
1124            assert num.allclose([x[i],y[i]], [e,n])
1125           
1126                       
1127        # Check the time vector
1128        times = fid.variables['time'][:]
1129        assert num.allclose(ensure_numeric(times),
1130                            ensure_numeric(times_ref))
1131                       
1132
1133        # Check sts values for mux #1
1134        stage1 = fid.variables['stage'][:].copy()
1135        xmomentum1 = fid.variables['xmomentum'][:].copy()
1136        ymomentum1 = fid.variables['ymomentum'][:].copy()
1137        elevation1 = fid.variables['elevation'][:].copy()
1138
1139       
1140        #print 'stage', stage1
1141        #print 'xmomentum', xmomentum1
1142        #print 'ymomentum', ymomentum1       
1143        #print 'elevation', elevation1
1144       
1145        # The quantities stored in the .sts file should be the weighted sum of the
1146        # quantities written to the mux2 files subject to the permutation vector.
1147       
1148        ha_ref = num.take(ha1, permutation, axis=0)
1149        ua_ref = num.take(ua1, permutation, axis=0)       
1150        va_ref = num.take(va1, permutation, axis=0)               
1151
1152        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
1153
1154
1155        #print
1156        #print stage1
1157        #print transpose(ha_ref)+tide - stage1
1158       
1159
1160        assert num.allclose(num.transpose(ha_ref)+tide, stage1)  # Meters
1161        #import sys; sys.exit()
1162
1163        #Check the momentums - ua
1164        #momentum = velocity*(stage-elevation)
1165        # elevation = - depth
1166        #momentum = velocity_ua *(stage+depth)
1167
1168        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
1169        for i in range(len(permutation)):
1170            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
1171
1172
1173        # The xmomentum stored in the .sts file should be the sum of the ua
1174        # in the two mux2 files multiplied by the depth.
1175        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum1) 
1176
1177        #Check the momentums - va
1178        #momentum = velocity*(stage-elevation)
1179        # elevation = - depth
1180        #momentum = velocity_va *(stage+depth)
1181
1182        # The ymomentum stored in the .sts file should be the sum of the va
1183        # in the two mux2 files multiplied by the depth.
1184       
1185       
1186        #print transpose(va_ref*depth_ref)
1187        #print ymomentum
1188        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum1)       
1189
1190        # check the elevation values.
1191        # -ve since urs measures depth, sww meshers height,
1192        assert num.allclose(-gauge_depth_ref, elevation1) 
1193
1194        fid.close()
1195        os.remove(sts_file)
1196       
1197        #----------------------
1198        # Then read the mux files together and test
1199       
1200                                               
1201        # Call urs2sts with multiple mux files
1202        urs2sts([base_nameI, base_nameII], 
1203                basename_out=base_nameI, 
1204                ordering_filename=ordering_filename,
1205                weights=weights,
1206                mean_stage=tide,
1207                verbose=False)
1208
1209        # Now read the sts file and check that values have been stored correctly.
1210        sts_file = base_nameI + '.sts'
1211        fid = NetCDFFile(sts_file)
1212
1213        # Make x and y absolute
1214        x = fid.variables['x'][:]
1215        y = fid.variables['y'][:]
1216
1217        geo_reference = Geo_reference(NetCDFObject=fid)
1218        points = geo_reference.get_absolute(map(None, x, y))
1219        points = ensure_numeric(points)
1220
1221        x = points[:,0]
1222        y = points[:,1]
1223
1224        for i, index in enumerate(permutation):
1225            # Check that STS points are stored in the correct order
1226           
1227            # Work out the UTM coordinates sts point i
1228            zone, e, n = redfearn(lat_long_points[index][0], 
1229                                  lat_long_points[index][1])             
1230
1231            #print i, [x[i],y[i]], [e,n]
1232            assert num.allclose([x[i],y[i]], [e,n])
1233           
1234                       
1235        # Check the time vector
1236        times = fid.variables['time'][:]
1237        assert num.allclose(ensure_numeric(times),
1238                            ensure_numeric(times_ref))
1239                       
1240
1241        # Check sts values
1242        stage = fid.variables['stage'][:].copy()
1243        xmomentum = fid.variables['xmomentum'][:].copy()
1244        ymomentum = fid.variables['ymomentum'][:].copy()
1245        elevation = fid.variables['elevation'][:].copy()
1246
1247       
1248        #print 'stage', stage
1249        #print 'elevation', elevation
1250       
1251        # The quantities stored in the .sts file should be the weighted sum of the
1252        # quantities written to the mux2 files subject to the permutation vector.
1253       
1254        ha_ref = (weights[0]*num.take(ha0, permutation, axis=0)
1255                  + weights[1]*num.take(ha1, permutation, axis=0))
1256        ua_ref = (weights[0]*num.take(ua0, permutation, axis=0)
1257                  + weights[1]*num.take(ua1, permutation, axis=0))
1258        va_ref = (weights[0]*num.take(va0, permutation, axis=0)
1259                  + weights[1]*num.take(va1, permutation, axis=0))
1260
1261        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
1262
1263
1264        #print
1265        #print stage
1266        #print transpose(ha_ref)+tide - stage
1267
1268        assert num.allclose(num.transpose(ha_ref)+tide, stage)  # Meters
1269
1270        #Check the momentums - ua
1271        #momentum = velocity*(stage-elevation)
1272        # elevation = - depth
1273        #momentum = velocity_ua *(stage+depth)
1274
1275        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
1276        for i in range(len(permutation)):
1277            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
1278
1279           
1280       
1281
1282        # The xmomentum stored in the .sts file should be the sum of the ua
1283        # in the two mux2 files multiplied by the depth.
1284        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum) 
1285
1286        #Check the momentums - va
1287        #momentum = velocity*(stage-elevation)
1288        # elevation = - depth
1289        #momentum = velocity_va *(stage+depth)
1290
1291        # The ymomentum stored in the .sts file should be the sum of the va
1292        # in the two mux2 files multiplied by the depth.
1293       
1294       
1295        #print transpose(va_ref*depth_ref)
1296        #print ymomentum
1297
1298        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum)
1299
1300        # check the elevation values.
1301        # -ve since urs measures depth, sww meshers height,
1302        assert num.allclose(-gauge_depth_ref, elevation)  #Meters
1303
1304        fid.close()
1305        self.delete_mux(filesI)
1306        self.delete_mux(filesII)
1307        os.remove(sts_file)
1308       
1309        #---------------
1310        # "Manually" add the timeseries up with weights and test
1311        # Tide is discounted from individual results and added back in       
1312        #
1313
1314        stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
1315        assert num.allclose(stage_man, stage)
1316               
1317       
1318    def test_file_boundary_stsI(self):
1319        """test_file_boundary_stsI(self):
1320        """
1321       
1322        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
1323        tide = 0.37
1324        time_step_count = 5
1325        time_step = 2
1326        lat_long_points =bounding_polygon[0:3]
1327        n=len(lat_long_points)
1328        first_tstep=num.ones(n,num.int)
1329        last_tstep=(time_step_count)*num.ones(n,num.int)
1330
1331        h = 20       
1332        w = 2
1333        u = 10
1334        v = -10
1335        gauge_depth=h*num.ones(n,num.float)
1336        ha=w*num.ones((n,time_step_count),num.float)
1337        ua=u*num.ones((n,time_step_count),num.float)
1338        va=v*num.ones((n,time_step_count),num.float)
1339        base_name, files = self.write_mux2(lat_long_points,
1340                                           time_step_count, time_step,
1341                                           first_tstep, last_tstep,
1342                                           depth=gauge_depth,
1343                                           ha=ha,
1344                                           ua=ua,
1345                                           va=va)
1346
1347        sts_file=base_name
1348        urs2sts(base_name,
1349                sts_file,
1350                mean_stage=tide,
1351                verbose=False)
1352        self.delete_mux(files)
1353
1354        #print 'start create mesh from regions'
1355        for i in range(len(bounding_polygon)):
1356            zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1])
1357        extent_res=1000000
1358        meshname = 'urs_test_mesh' + '.tsh'
1359        interior_regions=None
1360        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1361        create_mesh_from_regions(bounding_polygon,
1362                                 boundary_tags=boundary_tags,
1363                                 maximum_triangle_area=extent_res,
1364                                 filename=meshname,
1365                                 interior_regions=interior_regions,
1366                                 verbose=False)
1367       
1368        domain_fbound = Domain(meshname)
1369        domain_fbound.set_quantity('stage', tide)
1370        Bf = File_boundary(sts_file+'.sts', 
1371                           domain_fbound, 
1372                           boundary_polygon=bounding_polygon)
1373        Br = Reflective_boundary(domain_fbound)
1374        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])       
1375
1376        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1377       
1378        # Check boundary object evaluates as it should
1379        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
1380            if B is Bf:
1381           
1382                qf = B.evaluate(vol_id, edge_id)  # File boundary
1383                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
1384
1385                assert num.allclose(qf, qd) 
1386               
1387       
1388        # Evolve
1389        finaltime=time_step*(time_step_count-1)
1390        yieldstep=time_step
1391        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
1392
1393        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1394                                                   finaltime=finaltime, 
1395                                                   skip_initial_step=False)):
1396                                                   
1397            D = domain_fbound
1398            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
1399
1400            # Check that file boundary object has populated
1401            # boundary array correctly 
1402            # FIXME (Ole): Do this for the other tests too!
1403            for j, val in enumerate(D.get_quantity('stage').boundary_values):
1404           
1405                (vol_id, edge_id), B = D.boundary_objects[j]
1406                if isinstance(B, File_boundary):
1407                    #print j, val
1408                    assert num.allclose(val, w + tide)
1409
1410
1411       
1412        domain_drchlt = Domain(meshname)
1413        domain_drchlt.set_quantity('stage', tide)
1414        Br = Reflective_boundary(domain_drchlt)
1415
1416        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
1417        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
1418
1419        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
1420                                                   finaltime=finaltime, 
1421                                                   skip_initial_step=False)):
1422            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
1423
1424        #print domain_fbound.quantities['stage'].vertex_values
1425        #print domain_drchlt.quantities['stage'].vertex_values
1426                   
1427        assert num.allclose(temp_fbound,temp_drchlt)
1428       
1429        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
1430                            domain_drchlt.quantities['stage'].vertex_values)
1431                       
1432        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
1433                            domain_drchlt.quantities['xmomentum'].vertex_values) 
1434                       
1435        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
1436                            domain_drchlt.quantities['ymomentum'].vertex_values)
1437       
1438       
1439        os.remove(sts_file+'.sts')
1440        os.remove(meshname)
1441               
1442       
1443    def test_file_boundary_stsI_beyond_model_time(self):
1444        """test_file_boundary_stsI(self):
1445       
1446        Test that file_boundary can work when model time
1447        exceeds available data.
1448        This is optional and requires the user to specify a default
1449        boundary object.
1450        """
1451       
1452        # Don't do warnings in unit test
1453        import warnings
1454        warnings.simplefilter('ignore')
1455
1456        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
1457                          [6.02,97.02],[6.00,97.02]]
1458        tide = 0.37
1459        time_step_count = 5
1460        time_step = 2
1461        lat_long_points = bounding_polygon[0:3]
1462        n=len(lat_long_points)
1463        first_tstep=num.ones(n,num.int)
1464        last_tstep=(time_step_count)*num.ones(n,num.int)
1465
1466        h = 20       
1467        w = 2
1468        u = 10
1469        v = -10
1470        gauge_depth=h*num.ones(n,num.float)
1471        ha=w*num.ones((n,time_step_count),num.float)
1472        ua=u*num.ones((n,time_step_count),num.float)
1473        va=v*num.ones((n,time_step_count),num.float)
1474        base_name, files = self.write_mux2(lat_long_points,
1475                                           time_step_count, time_step,
1476                                           first_tstep, last_tstep,
1477                                           depth=gauge_depth,
1478                                           ha=ha,
1479                                           ua=ua,
1480                                           va=va)
1481
1482        sts_file=base_name
1483        urs2sts(base_name,
1484                sts_file,
1485                mean_stage=tide,
1486                verbose=False)
1487        self.delete_mux(files)
1488
1489        #print 'start create mesh from regions'
1490        for i in range(len(bounding_polygon)):
1491            zone,\
1492            bounding_polygon[i][0],\
1493            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
1494                                            bounding_polygon[i][1])
1495                                           
1496        extent_res=1000000
1497        meshname = 'urs_test_mesh' + '.tsh'
1498        interior_regions=None
1499        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1500        create_mesh_from_regions(bounding_polygon,
1501                                 boundary_tags=boundary_tags,
1502                                 maximum_triangle_area=extent_res,
1503                                 filename=meshname,
1504                                 interior_regions=interior_regions,
1505                                 verbose=False)
1506       
1507        domain_fbound = Domain(meshname)
1508        domain_fbound.set_quantity('stage', tide)
1509       
1510        Br = Reflective_boundary(domain_fbound)
1511        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])       
1512        Bf = File_boundary(sts_file+'.sts', 
1513                           domain_fbound, 
1514                           boundary_polygon=bounding_polygon,
1515                           default_boundary=Bd) # Condition to be used
1516                                                # if model time exceeds
1517                                                # available data
1518
1519        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1520       
1521        # Check boundary object evaluates as it should
1522        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
1523            if B is Bf:
1524           
1525                qf = B.evaluate(vol_id, edge_id)  # File boundary
1526                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
1527
1528                assert num.allclose(qf, qd) 
1529               
1530       
1531        # Evolve
1532        data_finaltime = time_step*(time_step_count-1)
1533        finaltime = data_finaltime + 10 # Let model time exceed available data
1534        yieldstep = time_step
1535        temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float)
1536
1537        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1538                                                   finaltime=finaltime, 
1539                                                   skip_initial_step=False)):
1540                                                   
1541            D = domain_fbound
1542            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
1543
1544            # Check that file boundary object has populated
1545            # boundary array correctly 
1546            # FIXME (Ole): Do this for the other tests too!
1547            for j, val in enumerate(D.get_quantity('stage').boundary_values):
1548           
1549                (vol_id, edge_id), B = D.boundary_objects[j]
1550                if isinstance(B, File_boundary):
1551                    #print j, val
1552                    assert num.allclose(val, w + tide)
1553
1554
1555        domain_drchlt = Domain(meshname)
1556        domain_drchlt.set_quantity('stage', tide)
1557        Br = Reflective_boundary(domain_drchlt)
1558
1559        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
1560        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
1561
1562        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
1563                                                   finaltime=finaltime, 
1564                                                   skip_initial_step=False)):
1565            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
1566
1567        #print domain_fbound.quantities['stage'].vertex_values
1568        #print domain_drchlt.quantities['stage'].vertex_values
1569                   
1570        assert num.allclose(temp_fbound,temp_drchlt)
1571       
1572        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
1573                            domain_drchlt.quantities['stage'].vertex_values)
1574                       
1575        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
1576                            domain_drchlt.quantities['xmomentum'].vertex_values) 
1577                       
1578        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
1579                            domain_drchlt.quantities['ymomentum'].vertex_values) 
1580       
1581        os.remove(sts_file+'.sts')
1582        os.remove(meshname)
1583               
1584       
1585    def test_field_boundary_stsI_beyond_model_time(self):
1586        """test_field_boundary(self):
1587       
1588        Test that field_boundary can work when model time
1589        exceeds available data whilst adjusting mean_stage.
1590       
1591        """
1592       
1593        # Don't do warnings in unit test
1594        import warnings
1595        warnings.simplefilter('ignore')
1596
1597        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
1598                          [6.02,97.02],[6.00,97.02]]
1599        tide = 0.37
1600        time_step_count = 5
1601        time_step = 2
1602        lat_long_points = bounding_polygon[0:3]
1603        n=len(lat_long_points)
1604        first_tstep=num.ones(n,num.int)
1605        last_tstep=(time_step_count)*num.ones(n,num.int)
1606
1607        h = 20       
1608        w = 2
1609        u = 10
1610        v = -10
1611        gauge_depth=h*num.ones(n,num.float)
1612        ha=w*num.ones((n,time_step_count),num.float)
1613        ua=u*num.ones((n,time_step_count),num.float)
1614        va=v*num.ones((n,time_step_count),num.float)
1615        base_name, files = self.write_mux2(lat_long_points,
1616                                           time_step_count, time_step,
1617                                           first_tstep, last_tstep,
1618                                           depth=gauge_depth,
1619                                           ha=ha,
1620                                           ua=ua,
1621                                           va=va)
1622
1623        sts_file=base_name
1624        urs2sts(base_name,
1625                sts_file,
1626                mean_stage=0.0, # Deliberately let Field_boundary do the adjustment
1627                verbose=False)
1628        self.delete_mux(files)
1629
1630        #print 'start create mesh from regions'
1631        for i in range(len(bounding_polygon)):
1632            zone,\
1633            bounding_polygon[i][0],\
1634            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
1635                                            bounding_polygon[i][1])
1636                                           
1637        extent_res=1000000
1638        meshname = 'urs_test_mesh' + '.tsh'
1639        interior_regions=None
1640        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1641        create_mesh_from_regions(bounding_polygon,
1642                                 boundary_tags=boundary_tags,
1643                                 maximum_triangle_area=extent_res,
1644                                 filename=meshname,
1645                                 interior_regions=interior_regions,
1646                                 verbose=False)
1647       
1648        domain_fbound = Domain(meshname)
1649        domain_fbound.set_quantity('stage', tide)
1650       
1651        Br = Reflective_boundary(domain_fbound)
1652        Bd = Dirichlet_boundary([w+tide, u*(w+h), v*(w+h)])
1653        Bdefault = Dirichlet_boundary([w, u*(w+h), v*(w+h)])       
1654               
1655        Bf = Field_boundary(sts_file+'.sts', 
1656                           domain_fbound, 
1657                           mean_stage=tide, # Field boundary to adjust for tide
1658                           boundary_polygon=bounding_polygon,
1659                           default_boundary=Bdefault) # Condition to be used
1660                                                      # if model time exceeds
1661                                                      # available data
1662
1663        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1664       
1665        # Check boundary object evaluates as it should
1666        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
1667            if B is Bf:
1668           
1669                qf = B.evaluate(vol_id, edge_id)  # Field boundary
1670                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
1671               
1672                msg = 'Got %s, should have been %s' %(qf, qd)
1673                assert num.allclose(qf, qd), msg
1674               
1675        # Evolve
1676        data_finaltime = time_step*(time_step_count-1)
1677        finaltime = data_finaltime + 10 # Let model time exceed available data
1678        yieldstep = time_step
1679        temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float)
1680
1681        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1682                                                   finaltime=finaltime, 
1683                                                   skip_initial_step=False)):
1684                                                   
1685            D = domain_fbound
1686            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
1687
1688            # Check that file boundary object has populated
1689            # boundary array correctly 
1690            # FIXME (Ole): Do this for the other tests too!
1691            for j, val in enumerate(D.get_quantity('stage').boundary_values):
1692           
1693                (vol_id, edge_id), B = D.boundary_objects[j]
1694                if isinstance(B, Field_boundary):
1695                    msg = 'Got %f should have been %f' %(val, w+tide)
1696                    assert num.allclose(val, w + tide), msg
1697
1698
1699    def test_file_boundary_stsII(self):
1700        """test_file_boundary_stsII(self):
1701       
1702         mux2 file has points not included in boundary
1703         mux2 gauges are not stored with the same order as they are
1704         found in bounding_polygon. This does not matter as long as bounding
1705         polygon passed to file_function contains the mux2 points needed (in
1706         the correct order).
1707         """
1708
1709        bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]]
1710        tide = -2.20 
1711        time_step_count = 20
1712        time_step = 2
1713        lat_long_points=bounding_polygon[0:2]
1714        lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1])
1715        lat_long_points.insert(0,[6.0,97.01])
1716        n=len(lat_long_points)
1717        first_tstep=num.ones(n,num.int)
1718        last_tstep=(time_step_count)*num.ones(n,num.int)
1719        gauge_depth=20*num.ones(n,num.float)
1720        ha=2*num.ones((n,time_step_count),num.float)
1721        ua=10*num.ones((n,time_step_count),num.float)
1722        va=-10*num.ones((n,time_step_count),num.float)
1723        base_name, files = self.write_mux2(lat_long_points,
1724                                           time_step_count,
1725                                           time_step,
1726                                           first_tstep,
1727                                           last_tstep,
1728                                           depth=gauge_depth,
1729                                           ha=ha,
1730                                           ua=ua,
1731                                           va=va)
1732
1733        sts_file=base_name
1734        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
1735        self.delete_mux(files)
1736
1737        #print 'start create mesh from regions'
1738        for i in range(len(bounding_polygon)):
1739            zone,\
1740            bounding_polygon[i][0],\
1741            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
1742                                            bounding_polygon[i][1])
1743           
1744        extent_res=1000000
1745        meshname = 'urs_test_mesh' + '.tsh'
1746        interior_regions=None
1747        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1748        # have to change boundary tags from last example because now bounding
1749        # polygon starts in different place.
1750        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
1751                         maximum_triangle_area=extent_res,filename=meshname,
1752                         interior_regions=interior_regions,verbose=False)
1753       
1754        domain_fbound = Domain(meshname)
1755        domain_fbound.set_quantity('stage', tide)
1756        Bf = File_boundary(sts_file+'.sts',
1757                           domain_fbound,
1758                           boundary_polygon=bounding_polygon)
1759        Br = Reflective_boundary(domain_fbound)
1760
1761        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1762        finaltime=time_step*(time_step_count-1)
1763        yieldstep=time_step
1764        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
1765       
1766        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1767                                                   finaltime=finaltime, 
1768                                                   skip_initial_step = False)):
1769            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
1770       
1771        domain_drchlt = Domain(meshname)
1772        domain_drchlt.set_quantity('stage', tide)
1773        Br = Reflective_boundary(domain_drchlt)
1774        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
1775        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
1776        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
1777
1778        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
1779                                                   finaltime=finaltime, 
1780                                                   skip_initial_step = False)):
1781            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
1782
1783
1784        assert num.allclose(temp_fbound,temp_drchlt)           
1785           
1786        #print domain_fbound.quantities['stage'].vertex_values
1787        #print domain_drchlt.quantities['stage'].vertex_values
1788                   
1789           
1790        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
1791                            domain_drchlt.quantities['stage'].vertex_values)
1792                       
1793        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
1794                            domain_drchlt.quantities['xmomentum'].vertex_values)
1795                       
1796        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
1797                            domain_drchlt.quantities['ymomentum'].vertex_values)
1798           
1799           
1800
1801        os.remove(sts_file+'.sts')
1802        os.remove(meshname)
1803
1804       
1805       
1806    def test_file_boundary_stsIII_ordering(self):
1807        """test_file_boundary_stsIII_ordering(self):
1808        Read correct points from ordering file and apply sts to boundary
1809        """
1810
1811        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
1812        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
1813                          [6.02,97.02],[6.00,97.02]]
1814        tide = 3.0
1815        time_step_count = 50
1816        time_step = 2
1817        n=len(lat_long_points)
1818        first_tstep=num.ones(n,num.int)
1819        last_tstep=(time_step_count)*num.ones(n,num.int)
1820        gauge_depth=20*num.ones(n,num.float)
1821        ha=2*num.ones((n,time_step_count),num.float)
1822        ua=10*num.ones((n,time_step_count),num.float)
1823        va=-10*num.ones((n,time_step_count),num.float)
1824        base_name, files = self.write_mux2(lat_long_points,
1825                                           time_step_count,
1826                                           time_step,
1827                                           first_tstep,
1828                                           last_tstep,
1829                                           depth=gauge_depth,
1830                                           ha=ha,
1831                                           ua=ua,
1832                                           va=va)
1833        # base name will not exist, but 3 other files are created
1834
1835        # Write order file
1836        file_handle, order_base_name = tempfile.mkstemp("")
1837        os.close(file_handle)
1838        os.remove(order_base_name)
1839        d=","
1840        order_file=order_base_name+'order.txt'
1841        fid=open(order_file,'w')
1842       
1843        # Write Header
1844        header='index, longitude, latitude\n'
1845        fid.write(header)
1846        indices=[3,0,1]
1847        for i in indices:
1848            line=str(i)+d+str(lat_long_points[i][1])+d+\
1849                str(lat_long_points[i][0])+"\n"
1850            fid.write(line)
1851        fid.close()
1852
1853        sts_file=base_name
1854        urs2sts(base_name,
1855                basename_out=sts_file,
1856                ordering_filename=order_file,
1857                mean_stage=tide,
1858                verbose=False)
1859        self.delete_mux(files)
1860
1861        assert(os.access(sts_file+'.sts', os.F_OK))
1862
1863        boundary_polygon = create_sts_boundary(base_name)
1864
1865        os.remove(order_file)
1866
1867        # Append the remaining part of the boundary polygon to be defined by
1868        # the user
1869        bounding_polygon_utm=[]
1870        for point in bounding_polygon:
1871            zone,easting,northing=redfearn(point[0],point[1])
1872            bounding_polygon_utm.append([easting,northing])
1873
1874        boundary_polygon.append(bounding_polygon_utm[3])
1875        boundary_polygon.append(bounding_polygon_utm[4])
1876
1877        plot=False
1878        if plot:
1879            from pylab import plot,show,axis
1880            boundary_polygon=ensure_numeric(boundary_polygon)
1881            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
1882            #lat_long_points=ensure_numeric(lat_long_points)
1883            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
1884            plot(boundary_polygon[:,0], boundary_polygon[:,1],'d')
1885            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1],'o')
1886            show()
1887
1888        assert num.allclose(bounding_polygon_utm,boundary_polygon)
1889
1890
1891        extent_res=1000000
1892        meshname = 'urs_test_mesh' + '.tsh'
1893        interior_regions=None
1894        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
1895       
1896        # have to change boundary tags from last example because now bounding
1897        # polygon starts in different place.
1898        create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags,
1899                         maximum_triangle_area=extent_res,filename=meshname,
1900                         interior_regions=interior_regions,verbose=False)
1901       
1902        domain_fbound = Domain(meshname)
1903        domain_fbound.set_quantity('stage', tide)
1904        Bf = File_boundary(sts_file+'.sts',
1905                           domain_fbound,
1906                           boundary_polygon=boundary_polygon)
1907        Br = Reflective_boundary(domain_fbound)
1908
1909        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
1910        finaltime=time_step*(time_step_count-1)
1911        yieldstep=time_step
1912        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
1913   
1914        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
1915                                                   finaltime=finaltime, 
1916                                                   skip_initial_step = False)):
1917            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
1918   
1919       
1920        domain_drchlt = Domain(meshname)
1921        domain_drchlt.set_quantity('stage', tide)
1922        Br = Reflective_boundary(domain_drchlt)
1923        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
1924        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
1925        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
1926       
1927        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
1928                                                   finaltime=finaltime, 
1929                                                   skip_initial_step = False)):
1930            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
1931
1932       
1933        #print domain_fbound.quantities['stage'].vertex_values
1934        #print domain_drchlt.quantities['stage'].vertex_values
1935                   
1936        assert num.allclose(temp_fbound,temp_drchlt)
1937
1938       
1939        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
1940                            domain_drchlt.quantities['stage'].vertex_values)
1941                       
1942        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
1943                            domain_drchlt.quantities['xmomentum'].vertex_values)                       
1944                       
1945        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
1946                            domain_drchlt.quantities['ymomentum'].vertex_values)
1947       
1948        # Use known Dirichlet condition (if sufficient timesteps have been taken)
1949
1950        #FIXME: What do these assertions test? Also do they assume tide =0
1951        #print domain_fbound.quantities['stage'].vertex_values
1952        #assert allclose(domain_drchlt.quantities['stage'].vertex_values[6], 2)       
1953        #assert allclose(domain_fbound.quantities['stage'].vertex_values[6], 2)
1954
1955        try:
1956            os.remove(sts_file+'.sts')
1957        except IOError:
1958            # Windoze can't remove this file for some reason
1959            pass
1960       
1961        os.remove(meshname)
1962       
1963
1964
1965#-------------------------------------------------------------
1966
1967if __name__ == "__main__":
1968    suite = unittest.makeSuite(Test_Urs2Sts,'test')
1969    runner = unittest.TextTestRunner()
1970    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.