Context Navigation

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

Last change on this file since 7780 was 7780, checked in by hudson, 14 years ago
Almost all failing tests fixed.
File size: 76.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: