Context Navigation

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, 14 years ago
New filename conventions for file conversion. Filenames must always be passed in with the correct extension.
File size: 79.0 KB

Line
1	import numpy as num
2	import unittest
3	import tempfile
4	import os
5	from Scientific.IO.NetCDF import NetCDFFile
6
7	from anuga.utilities.system_tools import get_pathname_from_package
8	from anuga.coordinate_transforms.geo_reference import Geo_reference
9	from anuga.coordinate_transforms.redfearn import redfearn
10	from anuga.utilities.numerical_tools import ensure_numeric
11	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
12	from anuga.file.sts import create_sts_boundary
13	from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array
14
15	from anuga.shallow_water.shallow_water_domain import Domain
16
17	# boundary functions
18	from anuga.shallow_water.boundaries import Reflective_boundary, \
19	Field_boundary, Transmissive_momentum_set_stage_boundary, \
20	Transmissive_stage_zero_momentum_boundary
21	from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
22	import Transmissive_boundary, Dirichlet_boundary, \
23	Time_boundary, File_boundary, AWI_boundary
24
25	from anuga.pmesh.mesh_interface import create_mesh_from_regions
26
27	from urs2sts import urs2sts
28
29	# Allow us to use helper methods from this test.
30	from anuga.file.test_mux import Test_Mux
31
32	class 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(2time_step_count,3time_step_count)
59	ha[3]=num.arange(3time_step_count,4time_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(2time_step_count,3time_step_count)
174	ha[3]=num.arange(3time_step_count,4time_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(2time_step_count,3time_step_count)
575	ha[3]=num.arange(3time_step_count,4time_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.0num.transpose(ua_permutationdepth_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.0num.transpose(va_permutationdepth_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(2time_step_count,3time_step_count)
762	ha[3]=num.arange(3time_step_count,4time_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(2time_step_count,3time_step_count)
879	ha0[3]=num.arange(3time_step_count,4time_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]=5num.sin(4times_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(3time_step_count,4time_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+10tide,-220-10tide])
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+10tide,-220-10tide])
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
1966	if __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.

Download in other formats: