Context Navigation

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

Last change on this file since 7772 was 7772, checked in by hudson, 14 years ago
Fixed a bunch of failing tests - only 3 still failing.
File size: 87.2 KB

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

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

Download in other formats: