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

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

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

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