Context Navigation

test_data_manager.py @ 5274

Last change on this file since 5274 was 5274, checked in by ole, 16 years ago
Work in progress regarding flow through a cross section, ticket:175
File size: 253.1 KB

Line
1	#!/usr/bin/env python
2	#
3
4
5	import unittest
6	import copy
7	from Numeric import zeros, array, allclose, Float
8	from anuga.utilities.numerical_tools import mean
9	import tempfile
10	import os
11	from Scientific.IO.NetCDF import NetCDFFile
12	from struct import pack
13	from sets import ImmutableSet
14
15	from anuga.shallow_water import *
16	from anuga.shallow_water.data_manager import *
17	from anuga.config import epsilon
18	from anuga.utilities.anuga_exceptions import ANUGAError
19	from anuga.utilities.numerical_tools import ensure_numeric
20	from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
21	from anuga.abstract_2d_finite_volumes.util import file_function
22
23	# This is needed to run the tests of local functions
24	import data_manager
25	from anuga.coordinate_transforms.redfearn import redfearn
26	from anuga.coordinate_transforms.geo_reference import Geo_reference, \
27	DEFAULT_ZONE
28	from anuga.geospatial_data.geospatial_data import Geospatial_data
29
30	class Test_Data_Manager(unittest.TestCase):
31	# Class variable
32	verbose = False
33
34	def set_verbose(self):
35	Test_Data_Manager.verbose = True
36
37	def setUp(self):
38	import time
39	from mesh_factory import rectangular
40
41
42	self.verbose = Test_Data_Manager.verbose
43	#Create basic mesh
44	points, vertices, boundary = rectangular(2, 2)
45
46	#Create shallow water domain
47	domain = Domain(points, vertices, boundary)
48	domain.default_order = 2
49
50	#Set some field values
51	domain.set_quantity('elevation', lambda x,y: -x)
52	domain.set_quantity('friction', 0.03)
53
54
55	######################
56	# Boundary conditions
57	B = Transmissive_boundary(domain)
58	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
59
60
61	######################
62	#Initial condition - with jumps
63	bed = domain.quantities['elevation'].vertex_values
64	stage = zeros(bed.shape, Float)
65
66	h = 0.3
67	for i in range(stage.shape[0]):
68	if i % 2 == 0:
69	stage[i,:] = bed[i,:] + h
70	else:
71	stage[i,:] = bed[i,:]
72
73	domain.set_quantity('stage', stage)
74
75
76	domain.distribute_to_vertices_and_edges()
77	self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
78
79
80
81	self.domain = domain
82
83	C = domain.get_vertex_coordinates()
84	self.X = C[:,0:6:2].copy()
85	self.Y = C[:,1:6:2].copy()
86
87	self.F = bed
88
89	#Write A testfile (not realistic. Values aren't realistic)
90	self.test_MOST_file = 'most_small'
91
92	longitudes = [150.66667, 150.83334, 151., 151.16667]
93	latitudes = [-34.5, -34.33333, -34.16667, -34]
94
95	long_name = 'LON'
96	lat_name = 'LAT'
97
98	nx = 4
99	ny = 4
100	six = 6
101
102
103	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
104	fid = NetCDFFile(self.test_MOST_file + ext, 'w')
105
106	fid.createDimension(long_name,nx)
107	fid.createVariable(long_name,'d',(long_name,))
108	fid.variables[long_name].point_spacing='uneven'
109	fid.variables[long_name].units='degrees_east'
110	fid.variables[long_name].assignValue(longitudes)
111
112	fid.createDimension(lat_name,ny)
113	fid.createVariable(lat_name,'d',(lat_name,))
114	fid.variables[lat_name].point_spacing='uneven'
115	fid.variables[lat_name].units='degrees_north'
116	fid.variables[lat_name].assignValue(latitudes)
117
118	fid.createDimension('TIME',six)
119	fid.createVariable('TIME','d',('TIME',))
120	fid.variables['TIME'].point_spacing='uneven'
121	fid.variables['TIME'].units='seconds'
122	fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
123
124
125	name = ext[1:3].upper()
126	if name == 'E.': name = 'ELEVATION'
127	fid.createVariable(name,'d',('TIME', lat_name, long_name))
128	fid.variables[name].units='CENTIMETERS'
129	fid.variables[name].missing_value=-1.e+034
130
131	fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
132	[-0.1214216, 0, 0, 0],
133	[0, 0, 0, 0],
134	[0, 0, 0, 0]],
135	[[0.3400644, 2.291054e-005, -23.33335, -6.50198],
136	[-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
137	[2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
138	[0, 2.291054e-005, 2.291054e-005, 0]],
139	[[0.3400644, 0.0001374632, -23.31503, -6.50198],
140	[-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
141	[0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
142	[0, 0.0001374632, 0.0001374632, 0]],
143	[[0.3400644, 0.0002520159, -23.29672, -6.50198],
144	[-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
145	[0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
146	[0, 0.0002520159, 0.0002520159, 0]],
147	[[0.3400644, 0.0003665686, -23.27842, -6.50198],
148	[-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
149	[0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
150	[0, 0.0003665686, 0.0003665686, 0]],
151	[[0.3400644, 0.0004811212, -23.26012, -6.50198],
152	[-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
153	[0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
154	[0, 0.0004811212, 0.0004811212, 0]]])
155
156
157	fid.close()
158
159
160
161
162	def tearDown(self):
163	import os
164	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
165	#print 'Trying to remove', self.test_MOST_file + ext
166	os.remove(self.test_MOST_file + ext)
167
168	def test_sww_constant(self):
169	"""Test that constant sww information can be written correctly
170	(non smooth)
171	"""
172	self.domain.set_name('datatest' + str(id(self)))
173	self.domain.format = 'sww' #Remove??
174	self.domain.smooth = False
175
176	sww = get_dataobject(self.domain)
177	sww.store_connectivity()
178
179	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
180
181	# Get the variables
182	x = fid.variables['x']
183	y = fid.variables['y']
184	z = fid.variables['elevation']
185	V = fid.variables['volumes']
186
187	assert allclose (x[:], self.X.flat)
188	assert allclose (y[:], self.Y.flat)
189	assert allclose (z[:], self.F.flat)
190
191	P = len(self.domain)
192	for k in range(P):
193	assert V[k, 0] == 3*k
194	assert V[k, 1] == 3*k+1
195	assert V[k, 2] == 3*k+2
196
197	fid.close()
198	os.remove(sww.filename)
199
200	def test_sww_header(self):
201	"""Test that constant sww information can be written correctly
202	(non smooth)
203	"""
204	self.domain.set_name('datatest' + str(id(self)))
205	self.domain.format = 'sww' #Remove??
206	self.domain.smooth = False
207
208	sww = get_dataobject(self.domain)
209	sww.store_connectivity()
210
211	#Check contents
212	#Get NetCDF
213	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
214
215	# Get the variables
216	sww_revision = fid.revision_number
217	try:
218	revision_number = get_revision_number()
219	except:
220	revision_number = None
221
222	assert str(revision_number) == sww_revision
223	fid.close()
224
225	#print "sww.filename", sww.filename
226	os.remove(sww.filename)
227
228	def test_sww_range(self):
229	"""Test that constant sww information can be written correctly
230	(non smooth)
231	"""
232	self.domain.set_name('datatest' + str(id(self)))
233	self.domain.format = 'sww'
234	self.domain.smooth = True
235
236	self.domain.tight_slope_limiters = 0 # Backwards compatibility
237
238	sww = get_dataobject(self.domain)
239
240	for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
241	pass
242
243	# Get NetCDF
244	fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
245
246	# Get the variables
247	range = fid.variables['stage_range'][:]
248	#print range
249	assert allclose(range,[-0.93519, 0.15]) or\
250	allclose(range,[-0.9352743, 0.15]) or\
251	allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
252
253	range = fid.variables['xmomentum_range'][:]
254	#print range
255	assert allclose(range,[0,0.4695096]) or \
256	allclose(range,[0,0.47790655]) or\
257	allclose(range,[0,0.46941957]) or\
258	allclose(range,[0,0.47769409])
259
260
261	range = fid.variables['ymomentum_range'][:]
262	#print range
263	assert allclose(range,[0,0.02174380]) or\
264	allclose(range,[0,0.02174439]) or\
265	allclose(range,[0,0.02283983]) or\
266	allclose(range,[0,0.0217342]) or\
267	allclose(range,[0,0.0227564]) # Old slope limiters
268
269	fid.close()
270	os.remove(sww.filename)
271
272	def test_sww_extrema(self):
273	"""Test that extrema of quantities can be retrieved at every vertex
274	Extrema are updated at every internal timestep
275	"""
276
277	domain = self.domain
278
279	domain.set_name('extrema_test' + str(id(self)))
280	domain.format = 'sww'
281	domain.smooth = True
282
283	assert domain.quantities_to_be_monitored is None
284	assert domain.monitor_polygon is None
285	assert domain.monitor_time_interval is None
286
287	domain.set_quantities_to_be_monitored(['xmomentum',
288	'ymomentum',
289	'stage-elevation'])
290
291	assert domain.monitor_polygon is None
292	assert domain.monitor_time_interval is None
293
294
295	domain.set_quantities_to_be_monitored(['xmomentum',
296	'ymomentum',
297	'stage-elevation'],
298	polygon=domain.get_boundary_polygon(),
299	time_interval=[0,1])
300
301
302	assert len(domain.quantities_to_be_monitored) == 3
303	assert domain.quantities_to_be_monitored.has_key('stage-elevation')
304	assert domain.quantities_to_be_monitored.has_key('xmomentum')
305	assert domain.quantities_to_be_monitored.has_key('ymomentum')
306
307
308	#domain.protect_against_isolated_degenerate_timesteps = True
309	#domain.tight_slope_limiters = 1
310	domain.tight_slope_limiters = 0 # Backwards compatibility
311
312	sww = get_dataobject(domain)
313
314	for t in domain.evolve(yieldstep = 1, finaltime = 1):
315	pass
316	#print domain.timestepping_statistics()
317	domain.quantity_statistics(precision = '%.8f') # Silent
318
319
320	# Get NetCDF
321	fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
322
323	# Get the variables
324	extrema = fid.variables['stage-elevation.extrema'][:]
325	assert allclose(extrema, [0.00, 0.30])
326
327	loc = fid.variables['stage-elevation.min_location'][:]
328	assert allclose(loc, [0.16666667, 0.33333333])
329
330	loc = fid.variables['stage-elevation.max_location'][:]
331	assert allclose(loc, [0.8333333, 0.16666667])
332
333	time = fid.variables['stage-elevation.max_time'][:]
334	assert allclose(time, 0.0)
335
336	extrema = fid.variables['xmomentum.extrema'][:]
337	assert allclose(extrema,[-0.06062178, 0.47873023]) or allclose(extrema, [-0.06062178, 0.47847986])
338
339	extrema = fid.variables['ymomentum.extrema'][:]
340	assert allclose(extrema,[0.00, 0.0625786]) or allclose(extrema,[0.00, 0.06062178])
341
342	time_interval = fid.variables['extrema.time_interval'][:]
343	assert allclose(time_interval, [0,1])
344
345	polygon = fid.variables['extrema.polygon'][:]
346	assert allclose(polygon, domain.get_boundary_polygon())
347
348	fid.close()
349	#print "sww.filename", sww.filename
350	os.remove(sww.filename)
351
352
353
354	def test_sww_constant_smooth(self):
355	"""Test that constant sww information can be written correctly
356	(non smooth)
357	"""
358	self.domain.set_name('datatest' + str(id(self)))
359	self.domain.format = 'sww'
360	self.domain.smooth = True
361
362	sww = get_dataobject(self.domain)
363	sww.store_connectivity()
364
365	#Check contents
366	#Get NetCDF
367	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
368
369	# Get the variables
370	X = fid.variables['x'][:]
371	Y = fid.variables['y'][:]
372	Z = fid.variables['elevation'][:]
373	V = fid.variables['volumes']
374
375	assert allclose([X[0], Y[0]], array([0.0, 0.0]))
376	assert allclose([X[1], Y[1]], array([0.0, 0.5]))
377	assert allclose([X[2], Y[2]], array([0.0, 1.0]))
378	assert allclose([X[4], Y[4]], array([0.5, 0.5]))
379	assert allclose([X[7], Y[7]], array([1.0, 0.5]))
380
381	assert Z[4] == -0.5
382
383	assert V[2,0] == 4
384	assert V[2,1] == 5
385	assert V[2,2] == 1
386	assert V[4,0] == 6
387	assert V[4,1] == 7
388	assert V[4,2] == 3
389
390	fid.close()
391	os.remove(sww.filename)
392
393
394	def test_sww_variable(self):
395	"""Test that sww information can be written correctly
396	"""
397	self.domain.set_name('datatest' + str(id(self)))
398	self.domain.format = 'sww'
399	self.domain.smooth = True
400	self.domain.reduction = mean
401
402	sww = get_dataobject(self.domain)
403	sww.store_connectivity()
404	sww.store_timestep()
405
406	#Check contents
407	#Get NetCDF
408	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
409
410
411	# Get the variables
412	time = fid.variables['time']
413	stage = fid.variables['stage']
414
415	Q = self.domain.quantities['stage']
416	Q0 = Q.vertex_values[:,0]
417	Q1 = Q.vertex_values[:,1]
418	Q2 = Q.vertex_values[:,2]
419
420	A = stage[0,:]
421	#print A[0], (Q2[0,0] + Q1[1,0])/2
422	assert allclose(A[0], (Q2[0] + Q1[1])/2)
423	assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
424	assert allclose(A[2], Q0[3])
425	assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
426
427	#Center point
428	assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
429	Q0[5] + Q2[6] + Q1[7])/6)
430
431	fid.close()
432	os.remove(sww.filename)
433
434
435	def test_sww_variable2(self):
436	"""Test that sww information can be written correctly
437	multiple timesteps. Use average as reduction operator
438	"""
439
440	import time, os
441	from Numeric import array, zeros, allclose, Float, concatenate
442	from Scientific.IO.NetCDF import NetCDFFile
443
444	self.domain.set_name('datatest' + str(id(self)))
445	self.domain.format = 'sww'
446	self.domain.smooth = True
447
448	self.domain.reduction = mean
449
450	sww = get_dataobject(self.domain)
451	sww.store_connectivity()
452	sww.store_timestep()
453	#self.domain.tight_slope_limiters = 1
454	self.domain.evolve_to_end(finaltime = 0.01)
455	sww.store_timestep()
456
457
458	#Check contents
459	#Get NetCDF
460	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
461
462	# Get the variables
463	x = fid.variables['x']
464	y = fid.variables['y']
465	z = fid.variables['elevation']
466	time = fid.variables['time']
467	stage = fid.variables['stage']
468
469	#Check values
470	Q = self.domain.quantities['stage']
471	Q0 = Q.vertex_values[:,0]
472	Q1 = Q.vertex_values[:,1]
473	Q2 = Q.vertex_values[:,2]
474
475	A = stage[1,:]
476	assert allclose(A[0], (Q2[0] + Q1[1])/2)
477	assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
478	assert allclose(A[2], Q0[3])
479	assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
480
481	#Center point
482	assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
483	Q0[5] + Q2[6] + Q1[7])/6)
484
485
486	fid.close()
487
488	#Cleanup
489	os.remove(sww.filename)
490
491	def no_test_sww_variable3(self):
492	"""Test that sww information can be written correctly
493	multiple timesteps using a different reduction operator (min)
494	"""
495
496	# Different reduction in sww files has been made obsolete.
497
498	import time, os
499	from Numeric import array, zeros, allclose, Float, concatenate
500	from Scientific.IO.NetCDF import NetCDFFile
501
502	self.domain.set_name('datatest' + str(id(self)))
503	self.domain.format = 'sww'
504	self.domain.smooth = True
505	self.domain.reduction = min
506
507	sww = get_dataobject(self.domain)
508	sww.store_connectivity()
509	sww.store_timestep()
510	#self.domain.tight_slope_limiters = 1
511	self.domain.evolve_to_end(finaltime = 0.01)
512	sww.store_timestep()
513
514
515	#Check contents
516	#Get NetCDF
517	fid = NetCDFFile(sww.filename, 'r')
518
519	# Get the variables
520	x = fid.variables['x']
521	y = fid.variables['y']
522	z = fid.variables['elevation']
523	time = fid.variables['time']
524	stage = fid.variables['stage']
525
526	#Check values
527	Q = self.domain.quantities['stage']
528	Q0 = Q.vertex_values[:,0]
529	Q1 = Q.vertex_values[:,1]
530	Q2 = Q.vertex_values[:,2]
531
532	A = stage[1,:]
533	assert allclose(A[0], min(Q2[0], Q1[1]))
534	assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
535	assert allclose(A[2], Q0[3])
536	assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
537
538	#Center point
539	assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
540	Q0[5], Q2[6], Q1[7]))
541
542
543	fid.close()
544
545	#Cleanup
546	os.remove(sww.filename)
547
548
549	def test_sync(self):
550	"""test_sync - Test info stored at each timestep is as expected (incl initial condition)
551	"""
552
553	import time, os, config
554	from Numeric import array, zeros, allclose, Float, concatenate
555	from Scientific.IO.NetCDF import NetCDFFile
556
557	self.domain.set_name('synctest')
558	self.domain.format = 'sww'
559	self.domain.smooth = False
560	self.domain.store = True
561	self.domain.beta_h = 0
562
563	# In this case tight_slope_limiters as default
564	# in conjunction with protection
565	# against isolated degenerate timesteps works.
566	#self.domain.tight_slope_limiters = 1
567	#self.domain.protect_against_isolated_degenerate_timesteps = True
568
569	#print 'tight_sl', self.domain.tight_slope_limiters
570
571
572	#Evolution
573	for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
574
575	#########self.domain.write_time(track_speeds=True)
576	stage = self.domain.quantities['stage'].vertex_values
577
578	#Get NetCDF
579	fid = NetCDFFile(self.domain.writer.filename, 'r')
580	stage_file = fid.variables['stage']
581
582	if t == 0.0:
583	assert allclose(stage, self.initial_stage)
584	assert allclose(stage_file[:], stage.flat)
585	else:
586	assert not allclose(stage, self.initial_stage)
587	assert not allclose(stage_file[:], stage.flat)
588
589	fid.close()
590
591	os.remove(self.domain.writer.filename)
592
593
594	def test_sww_minimum_storable_height(self):
595	"""Test that sww information can be written correctly
596	multiple timesteps using a different reduction operator (min)
597	"""
598
599	import time, os
600	from Numeric import array, zeros, allclose, Float, concatenate
601	from Scientific.IO.NetCDF import NetCDFFile
602
603	self.domain.set_name('datatest' + str(id(self)))
604	self.domain.format = 'sww'
605	self.domain.smooth = True
606	self.domain.reduction = min
607	self.domain.minimum_storable_height = 100
608
609	sww = get_dataobject(self.domain)
610	sww.store_connectivity()
611	sww.store_timestep()
612
613	#self.domain.tight_slope_limiters = 1
614	self.domain.evolve_to_end(finaltime = 0.01)
615	sww.store_timestep()
616
617
618	#Check contents
619	#Get NetCDF
620	fid = NetCDFFile(sww.filename, 'r')
621
622
623	# Get the variables
624	x = fid.variables['x']
625	y = fid.variables['y']
626	z = fid.variables['elevation']
627	time = fid.variables['time']
628	stage = fid.variables['stage']
629	xmomentum = fid.variables['xmomentum']
630	ymomentum = fid.variables['ymomentum']
631
632	#Check values
633	Q = self.domain.quantities['stage']
634	Q0 = Q.vertex_values[:,0]
635	Q1 = Q.vertex_values[:,1]
636	Q2 = Q.vertex_values[:,2]
637
638	A = stage[1,:]
639	assert allclose(stage[1,:], z[:])
640
641
642	assert allclose(xmomentum, 0.0)
643	assert allclose(ymomentum, 0.0)
644
645	fid.close()
646
647	#Cleanup
648	os.remove(sww.filename)
649
650
651	def Not_a_test_sww_DSG(self):
652	"""Not a test, rather a look at the sww format
653	"""
654
655	import time, os
656	from Numeric import array, zeros, allclose, Float, concatenate
657	from Scientific.IO.NetCDF import NetCDFFile
658
659	self.domain.set_name('datatest' + str(id(self)))
660	self.domain.format = 'sww'
661	self.domain.smooth = True
662	self.domain.reduction = mean
663
664	sww = get_dataobject(self.domain)
665	sww.store_connectivity()
666	sww.store_timestep()
667
668	#Check contents
669	#Get NetCDF
670	fid = NetCDFFile(sww.filename, 'r')
671
672	# Get the variables
673	x = fid.variables['x']
674	y = fid.variables['y']
675	z = fid.variables['elevation']
676
677	volumes = fid.variables['volumes']
678	time = fid.variables['time']
679
680	# 2D
681	stage = fid.variables['stage']
682
683	X = x[:]
684	Y = y[:]
685	Z = z[:]
686	V = volumes[:]
687	T = time[:]
688	S = stage[:,:]
689
690	# print "****************************"
691	# print "X ",X
692	# print "****************************"
693	# print "Y ",Y
694	# print "****************************"
695	# print "Z ",Z
696	# print "****************************"
697	# print "V ",V
698	# print "****************************"
699	# print "Time ",T
700	# print "****************************"
701	# print "Stage ",S
702	# print "****************************"
703
704
705	fid.close()
706
707	#Cleanup
708	os.remove(sww.filename)
709
710
711
712	def test_dem2pts_bounding_box_v2(self):
713	"""Test conversion from dem in ascii format to native NetCDF format
714	"""
715
716	import time, os
717	from Numeric import array, zeros, allclose, Float, concatenate, ones
718	from Scientific.IO.NetCDF import NetCDFFile
719
720	#Write test asc file
721	root = 'demtest'
722
723	filename = root+'.asc'
724	fid = open(filename, 'w')
725	fid.write("""ncols 10
726	nrows 10
727	xllcorner 2000
728	yllcorner 3000
729	cellsize 1
730	NODATA_value -9999
731	""")
732	#Create linear function
733	ref_points = []
734	ref_elevation = []
735	x0 = 2000
736	y = 3010
737	yvec = range(10)
738	xvec = range(10)
739	z = -1
740	for i in range(10):
741	y = y - 1
742	for j in range(10):
743	x = x0 + xvec[j]
744	z += 1
745	ref_points.append ([x,y])
746	ref_elevation.append(z)
747	fid.write('%f ' %z)
748	fid.write('\n')
749
750	fid.close()
751
752	#print 'sending pts', ref_points
753	#print 'sending elev', ref_elevation
754
755	#Write prj file with metadata
756	metafilename = root+'.prj'
757	fid = open(metafilename, 'w')
758
759
760	fid.write("""Projection UTM
761	Zone 56
762	Datum WGS84
763	Zunits NO
764	Units METERS
765	Spheroid WGS84
766	Xshift 0.0000000000
767	Yshift 10000000.0000000000
768	Parameters
769	""")
770	fid.close()
771
772	#Convert to NetCDF pts
773	convert_dem_from_ascii2netcdf(root)
774	dem2pts(root, easting_min=2002.0, easting_max=2007.0,
775	northing_min=3003.0, northing_max=3006.0,
776	verbose=self.verbose)
777
778	#Check contents
779	#Get NetCDF
780	fid = NetCDFFile(root+'.pts', 'r')
781
782	# Get the variables
783	#print fid.variables.keys()
784	points = fid.variables['points']
785	elevation = fid.variables['elevation']
786
787	#Check values
788	assert fid.xllcorner[0] == 2002.0
789	assert fid.yllcorner[0] == 3003.0
790
791	#create new reference points
792	newz = []
793	newz[0:5] = ref_elevation[32:38]
794	newz[6:11] = ref_elevation[42:48]
795	newz[12:17] = ref_elevation[52:58]
796	newz[18:23] = ref_elevation[62:68]
797	ref_elevation = []
798	ref_elevation = newz
799	ref_points = []
800	x0 = 2002
801	y = 3007
802	yvec = range(4)
803	xvec = range(6)
804	for i in range(4):
805	y = y - 1
806	ynew = y - 3003.0
807	for j in range(6):
808	x = x0 + xvec[j]
809	xnew = x - 2002.0
810	ref_points.append ([xnew,ynew]) #Relative point values
811
812	assert allclose(points, ref_points)
813
814	assert allclose(elevation, ref_elevation)
815
816	#Cleanup
817	fid.close()
818
819
820	os.remove(root + '.pts')
821	os.remove(root + '.dem')
822	os.remove(root + '.asc')
823	os.remove(root + '.prj')
824
825
826	def test_dem2pts_bounding_box_removeNullvalues_v2(self):
827	"""Test conversion from dem in ascii format to native NetCDF format
828	"""
829
830	import time, os
831	from Numeric import array, zeros, allclose, Float, concatenate, ones
832	from Scientific.IO.NetCDF import NetCDFFile
833
834	#Write test asc file
835	root = 'demtest'
836
837	filename = root+'.asc'
838	fid = open(filename, 'w')
839	fid.write("""ncols 10
840	nrows 10
841	xllcorner 2000
842	yllcorner 3000
843	cellsize 1
844	NODATA_value -9999
845	""")
846	#Create linear function
847	ref_points = []
848	ref_elevation = []
849	x0 = 2000
850	y = 3010
851	yvec = range(10)
852	xvec = range(10)
853	#z = range(100)
854	z = zeros(100)
855	NODATA_value = -9999
856	count = -1
857	for i in range(10):
858	y = y - 1
859	for j in range(10):
860	x = x0 + xvec[j]
861	ref_points.append ([x,y])
862	count += 1
863	z[count] = (4i - 3j)%13
864	if j == 4: z[count] = NODATA_value #column inside clipping region
865	if j == 8: z[count] = NODATA_value #column outside clipping region
866	if i == 9: z[count] = NODATA_value #row outside clipping region
867	if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
868	ref_elevation.append( z[count] )
869	fid.write('%f ' %z[count])
870	fid.write('\n')
871
872	fid.close()
873
874	#print 'sending elev', ref_elevation
875
876	#Write prj file with metadata
877	metafilename = root+'.prj'
878	fid = open(metafilename, 'w')
879
880
881	fid.write("""Projection UTM
882	Zone 56
883	Datum WGS84
884	Zunits NO
885	Units METERS
886	Spheroid WGS84
887	Xshift 0.0000000000
888	Yshift 10000000.0000000000
889	Parameters
890	""")
891	fid.close()
892
893	#Convert to NetCDF pts
894	convert_dem_from_ascii2netcdf(root)
895	dem2pts(root, easting_min=2002.0, easting_max=2007.0,
896	northing_min=3003.0, northing_max=3006.0,
897	verbose=self.verbose)
898
899	#Check contents
900	#Get NetCDF
901	fid = NetCDFFile(root+'.pts', 'r')
902
903	# Get the variables
904	#print fid.variables.keys()
905	points = fid.variables['points']
906	elevation = fid.variables['elevation']
907
908	#Check values
909	assert fid.xllcorner[0] == 2002.0
910	assert fid.yllcorner[0] == 3003.0
911
912	#create new reference points
913	newz = zeros(19)
914	newz[0:2] = ref_elevation[32:34]
915	newz[2:5] = ref_elevation[35:38]
916	newz[5:7] = ref_elevation[42:44]
917	newz[7] = ref_elevation[45]
918	newz[8] = ref_elevation[47]
919	newz[9:11] = ref_elevation[52:54]
920	newz[11:14] = ref_elevation[55:58]
921	newz[14:16] = ref_elevation[62:64]
922	newz[16:19] = ref_elevation[65:68]
923
924
925	ref_elevation = newz
926	ref_points = []
927	new_ref_points = []
928	x0 = 2002
929	y = 3007
930	yvec = range(4)
931	xvec = range(6)
932	for i in range(4):
933	y = y - 1
934	ynew = y - 3003.0
935	for j in range(6):
936	x = x0 + xvec[j]
937	xnew = x - 2002.0
938	if j <> 2 and (i<>1 or j<>4):
939	ref_points.append([x,y])
940	new_ref_points.append ([xnew,ynew])
941
942
943	assert allclose(points, new_ref_points)
944	assert allclose(elevation, ref_elevation)
945
946	#Cleanup
947	fid.close()
948
949
950	os.remove(root + '.pts')
951	os.remove(root + '.dem')
952	os.remove(root + '.asc')
953	os.remove(root + '.prj')
954
955
956	def test_dem2pts_bounding_box_removeNullvalues_v3(self):
957	"""Test conversion from dem in ascii format to native NetCDF format
958	Check missing values on clipping boundary
959	"""
960
961	import time, os
962	from Numeric import array, zeros, allclose, Float, concatenate, ones
963	from Scientific.IO.NetCDF import NetCDFFile
964
965	#Write test asc file
966	root = 'demtest'
967
968	filename = root+'.asc'
969	fid = open(filename, 'w')
970	fid.write("""ncols 10
971	nrows 10
972	xllcorner 2000
973	yllcorner 3000
974	cellsize 1
975	NODATA_value -9999
976	""")
977	#Create linear function
978	ref_points = []
979	ref_elevation = []
980	x0 = 2000
981	y = 3010
982	yvec = range(10)
983	xvec = range(10)
984	#z = range(100)
985	z = zeros(100)
986	NODATA_value = -9999
987	count = -1
988	for i in range(10):
989	y = y - 1
990	for j in range(10):
991	x = x0 + xvec[j]
992	ref_points.append ([x,y])
993	count += 1
994	z[count] = (4i - 3j)%13
995	if j == 4: z[count] = NODATA_value #column inside clipping region
996	if j == 8: z[count] = NODATA_value #column outside clipping region
997	if i == 6: z[count] = NODATA_value #row on clipping boundary
998	if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
999	ref_elevation.append( z[count] )
1000	fid.write('%f ' %z[count])
1001	fid.write('\n')
1002
1003	fid.close()
1004
1005	#print 'sending elev', ref_elevation
1006
1007	#Write prj file with metadata
1008	metafilename = root+'.prj'
1009	fid = open(metafilename, 'w')
1010
1011
1012	fid.write("""Projection UTM
1013	Zone 56
1014	Datum WGS84
1015	Zunits NO
1016	Units METERS
1017	Spheroid WGS84
1018	Xshift 0.0000000000
1019	Yshift 10000000.0000000000
1020	Parameters
1021	""")
1022	fid.close()
1023
1024	#Convert to NetCDF pts
1025	convert_dem_from_ascii2netcdf(root)
1026	dem2pts(root, easting_min=2002.0, easting_max=2007.0,
1027	northing_min=3003.0, northing_max=3006.0,
1028	verbose=self.verbose)
1029
1030	#Check contents
1031	#Get NetCDF
1032	fid = NetCDFFile(root+'.pts', 'r')
1033
1034	# Get the variables
1035	#print fid.variables.keys()
1036	points = fid.variables['points']
1037	elevation = fid.variables['elevation']
1038
1039	#Check values
1040	assert fid.xllcorner[0] == 2002.0
1041	assert fid.yllcorner[0] == 3003.0
1042
1043	#create new reference points
1044	newz = zeros(14)
1045	newz[0:2] = ref_elevation[32:34]
1046	newz[2:5] = ref_elevation[35:38]
1047	newz[5:7] = ref_elevation[42:44]
1048	newz[7] = ref_elevation[45]
1049	newz[8] = ref_elevation[47]
1050	newz[9:11] = ref_elevation[52:54]
1051	newz[11:14] = ref_elevation[55:58]
1052
1053
1054
1055	ref_elevation = newz
1056	ref_points = []
1057	new_ref_points = []
1058	x0 = 2002
1059	y = 3007
1060	yvec = range(4)
1061	xvec = range(6)
1062	for i in range(4):
1063	y = y - 1
1064	ynew = y - 3003.0
1065	for j in range(6):
1066	x = x0 + xvec[j]
1067	xnew = x - 2002.0
1068	if j <> 2 and (i<>1 or j<>4) and i<>3:
1069	ref_points.append([x,y])
1070	new_ref_points.append ([xnew,ynew])
1071
1072
1073	#print points[:],points[:].shape
1074	#print new_ref_points, len(new_ref_points)
1075
1076	assert allclose(elevation, ref_elevation)
1077	assert allclose(points, new_ref_points)
1078
1079
1080	#Cleanup
1081	fid.close()
1082
1083
1084	os.remove(root + '.pts')
1085	os.remove(root + '.dem')
1086	os.remove(root + '.asc')
1087	os.remove(root + '.prj')
1088
1089
1090	def test_hecras_cross_sections2pts(self):
1091	"""Test conversion from HECRAS cross sections in ascii format
1092	to native NetCDF pts format
1093	"""
1094
1095	import time, os
1096	from Numeric import array, zeros, allclose, Float, concatenate
1097	from Scientific.IO.NetCDF import NetCDFFile
1098
1099	#Write test asc file
1100	root = 'hecrastest'
1101
1102	filename = root+'.sdf'
1103	fid = open(filename, 'w')
1104	fid.write("""
1105	# RAS export file created on Mon 15Aug2005 11:42
1106	# by HEC-RAS Version 3.1.1
1107
1108	BEGIN HEADER:
1109	UNITS: METRIC
1110	DTM TYPE: TIN
1111	DTM: v:\1\cit\perth_topo\river_tin
1112	STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
1113	CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
1114	MAP PROJECTION: UTM
1115	PROJECTION ZONE: 50
1116	DATUM: AGD66
1117	VERTICAL DATUM:
1118	NUMBER OF REACHES: 19
1119	NUMBER OF CROSS-SECTIONS: 2
1120	END HEADER:
1121
1122
1123	BEGIN CROSS-SECTIONS:
1124
1125	CROSS-SECTION:
1126	STREAM ID:Southern-Wungong
1127	REACH ID:Southern-Wungong
1128	STATION:21410
1129	CUT LINE:
1130	407546.08 , 6437277.542
1131	407329.32 , 6437489.482
1132	407283.11 , 6437541.232
1133	SURFACE LINE:
1134	407546.08, 6437277.54, 52.14
1135	407538.88, 6437284.58, 51.07
1136	407531.68, 6437291.62, 50.56
1137	407524.48, 6437298.66, 49.58
1138	407517.28, 6437305.70, 49.09
1139	407510.08, 6437312.74, 48.76
1140	END:
1141
1142	CROSS-SECTION:
1143	STREAM ID:Swan River
1144	REACH ID:Swan Mouth
1145	STATION:840.*
1146	CUT LINE:
1147	381178.0855 , 6452559.0685
1148	380485.4755 , 6453169.272
1149	SURFACE LINE:
1150	381178.09, 6452559.07, 4.17
1151	381169.49, 6452566.64, 4.26
1152	381157.78, 6452576.96, 4.34
1153	381155.97, 6452578.56, 4.35
1154	381143.72, 6452589.35, 4.43
1155	381136.69, 6452595.54, 4.58
1156	381114.74, 6452614.88, 4.41
1157	381075.53, 6452649.43, 4.17
1158	381071.47, 6452653.00, 3.99
1159	381063.46, 6452660.06, 3.67
1160	381054.41, 6452668.03, 3.67
1161	END:
1162	END CROSS-SECTIONS:
1163	""")
1164
1165	fid.close()
1166
1167
1168	#Convert to NetCDF pts
1169	hecras_cross_sections2pts(root)
1170
1171	#Check contents
1172	#Get NetCDF
1173	fid = NetCDFFile(root+'.pts', 'r')
1174
1175	# Get the variables
1176	#print fid.variables.keys()
1177	points = fid.variables['points']
1178	elevation = fid.variables['elevation']
1179
1180	#Check values
1181	ref_points = [[407546.08, 6437277.54],
1182	[407538.88, 6437284.58],
1183	[407531.68, 6437291.62],
1184	[407524.48, 6437298.66],
1185	[407517.28, 6437305.70],
1186	[407510.08, 6437312.74]]
1187
1188	ref_points += [[381178.09, 6452559.07],
1189	[381169.49, 6452566.64],
1190	[381157.78, 6452576.96],
1191	[381155.97, 6452578.56],
1192	[381143.72, 6452589.35],
1193	[381136.69, 6452595.54],
1194	[381114.74, 6452614.88],
1195	[381075.53, 6452649.43],
1196	[381071.47, 6452653.00],
1197	[381063.46, 6452660.06],
1198	[381054.41, 6452668.03]]
1199
1200
1201	ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
1202	ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
1203
1204	#print points[:]
1205	#print ref_points
1206	assert allclose(points, ref_points)
1207
1208	#print attributes[:]
1209	#print ref_elevation
1210	assert allclose(elevation, ref_elevation)
1211
1212	#Cleanup
1213	fid.close()
1214
1215
1216	os.remove(root + '.sdf')
1217	os.remove(root + '.pts')
1218
1219
1220	def test_sww2dem_asc_elevation_depth(self):
1221	"""
1222	test_sww2dem_asc_elevation_depth(self):
1223	Test that sww information can be converted correctly to asc/prj
1224	format readable by e.g. ArcView
1225	"""
1226
1227	import time, os
1228	from Numeric import array, zeros, allclose, Float, concatenate
1229	from Scientific.IO.NetCDF import NetCDFFile
1230
1231	#Setup
1232	self.domain.set_name('datatest')
1233
1234	prjfile = self.domain.get_name() + '_elevation.prj'
1235	ascfile = self.domain.get_name() + '_elevation.asc'
1236	swwfile = self.domain.get_name() + '.sww'
1237
1238	self.domain.set_datadir('.')
1239	self.domain.format = 'sww'
1240	self.domain.smooth = True
1241	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1242	self.domain.set_quantity('stage', 1.0)
1243
1244	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1245
1246	sww = get_dataobject(self.domain)
1247	sww.store_connectivity()
1248	sww.store_timestep()
1249
1250	#self.domain.tight_slope_limiters = 1
1251
1252	self.domain.evolve_to_end(finaltime = 0.01)
1253	sww.store_timestep()
1254
1255	cellsize = 0.25
1256	#Check contents
1257	#Get NetCDF
1258
1259	fid = NetCDFFile(sww.filename, 'r')
1260
1261	# Get the variables
1262	x = fid.variables['x'][:]
1263	y = fid.variables['y'][:]
1264	z = fid.variables['elevation'][:]
1265	time = fid.variables['time'][:]
1266	stage = fid.variables['stage'][:]
1267
1268	fid.close()
1269
1270	#Export to ascii/prj files
1271	sww2dem(self.domain.get_name(),
1272	quantity = 'elevation',
1273	cellsize = cellsize,
1274	verbose = self.verbose,
1275	format = 'asc')
1276
1277	#Check prj (meta data)
1278	prjid = open(prjfile)
1279	lines = prjid.readlines()
1280	prjid.close()
1281
1282	L = lines[0].strip().split()
1283	assert L[0].strip().lower() == 'projection'
1284	assert L[1].strip().lower() == 'utm'
1285
1286	L = lines[1].strip().split()
1287	assert L[0].strip().lower() == 'zone'
1288	assert L[1].strip().lower() == '56'
1289
1290	L = lines[2].strip().split()
1291	assert L[0].strip().lower() == 'datum'
1292	assert L[1].strip().lower() == 'wgs84'
1293
1294	L = lines[3].strip().split()
1295	assert L[0].strip().lower() == 'zunits'
1296	assert L[1].strip().lower() == 'no'
1297
1298	L = lines[4].strip().split()
1299	assert L[0].strip().lower() == 'units'
1300	assert L[1].strip().lower() == 'meters'
1301
1302	L = lines[5].strip().split()
1303	assert L[0].strip().lower() == 'spheroid'
1304	assert L[1].strip().lower() == 'wgs84'
1305
1306	L = lines[6].strip().split()
1307	assert L[0].strip().lower() == 'xshift'
1308	assert L[1].strip().lower() == '500000'
1309
1310	L = lines[7].strip().split()
1311	assert L[0].strip().lower() == 'yshift'
1312	assert L[1].strip().lower() == '10000000'
1313
1314	L = lines[8].strip().split()
1315	assert L[0].strip().lower() == 'parameters'
1316
1317
1318	#Check asc file
1319	ascid = open(ascfile)
1320	lines = ascid.readlines()
1321	ascid.close()
1322
1323	L = lines[0].strip().split()
1324	assert L[0].strip().lower() == 'ncols'
1325	assert L[1].strip().lower() == '5'
1326
1327	L = lines[1].strip().split()
1328	assert L[0].strip().lower() == 'nrows'
1329	assert L[1].strip().lower() == '5'
1330
1331	L = lines[2].strip().split()
1332	assert L[0].strip().lower() == 'xllcorner'
1333	assert allclose(float(L[1].strip().lower()), 308500)
1334
1335	L = lines[3].strip().split()
1336	assert L[0].strip().lower() == 'yllcorner'
1337	assert allclose(float(L[1].strip().lower()), 6189000)
1338
1339	L = lines[4].strip().split()
1340	assert L[0].strip().lower() == 'cellsize'
1341	assert allclose(float(L[1].strip().lower()), cellsize)
1342
1343	L = lines[5].strip().split()
1344	assert L[0].strip() == 'NODATA_value'
1345	assert L[1].strip().lower() == '-9999'
1346
1347	#Check grid values
1348	for j in range(5):
1349	L = lines[6+j].strip().split()
1350	y = (4-j) * cellsize
1351	for i in range(5):
1352	assert allclose(float(L[i]), -i*cellsize - y)
1353
1354	#Cleanup
1355	os.remove(prjfile)
1356	os.remove(ascfile)
1357
1358	#Export to ascii/prj files
1359	sww2dem(self.domain.get_name(),
1360	quantity = 'depth',
1361	cellsize = cellsize,
1362	verbose = self.verbose,
1363	format = 'asc')
1364
1365	#Check asc file
1366	ascfile = self.domain.get_name() + '_depth.asc'
1367	prjfile = self.domain.get_name() + '_depth.prj'
1368	ascid = open(ascfile)
1369	lines = ascid.readlines()
1370	ascid.close()
1371
1372	L = lines[0].strip().split()
1373	assert L[0].strip().lower() == 'ncols'
1374	assert L[1].strip().lower() == '5'
1375
1376	L = lines[1].strip().split()
1377	assert L[0].strip().lower() == 'nrows'
1378	assert L[1].strip().lower() == '5'
1379
1380	L = lines[2].strip().split()
1381	assert L[0].strip().lower() == 'xllcorner'
1382	assert allclose(float(L[1].strip().lower()), 308500)
1383
1384	L = lines[3].strip().split()
1385	assert L[0].strip().lower() == 'yllcorner'
1386	assert allclose(float(L[1].strip().lower()), 6189000)
1387
1388	L = lines[4].strip().split()
1389	assert L[0].strip().lower() == 'cellsize'
1390	assert allclose(float(L[1].strip().lower()), cellsize)
1391
1392	L = lines[5].strip().split()
1393	assert L[0].strip() == 'NODATA_value'
1394	assert L[1].strip().lower() == '-9999'
1395
1396	#Check grid values
1397	for j in range(5):
1398	L = lines[6+j].strip().split()
1399	y = (4-j) * cellsize
1400	for i in range(5):
1401	assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1402
1403
1404	#Cleanup
1405	os.remove(prjfile)
1406	os.remove(ascfile)
1407	os.remove(swwfile)
1408
1409
1410	def test_export_grid(self):
1411	"""
1412	test_export_grid(self):
1413	Test that sww information can be converted correctly to asc/prj
1414	format readable by e.g. ArcView
1415	"""
1416
1417	import time, os
1418	from Numeric import array, zeros, allclose, Float, concatenate
1419	from Scientific.IO.NetCDF import NetCDFFile
1420
1421	try:
1422	os.remove('teg*.sww')
1423	except:
1424	pass
1425
1426	#Setup
1427	self.domain.set_name('teg')
1428
1429	prjfile = self.domain.get_name() + '_elevation.prj'
1430	ascfile = self.domain.get_name() + '_elevation.asc'
1431	swwfile = self.domain.get_name() + '.sww'
1432
1433	self.domain.set_datadir('.')
1434	self.domain.format = 'sww'
1435	self.domain.smooth = True
1436	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1437	self.domain.set_quantity('stage', 1.0)
1438
1439	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1440
1441	sww = get_dataobject(self.domain)
1442	sww.store_connectivity()
1443	sww.store_timestep()
1444	self.domain.evolve_to_end(finaltime = 0.01)
1445	sww.store_timestep()
1446
1447	cellsize = 0.25
1448	#Check contents
1449	#Get NetCDF
1450
1451	fid = NetCDFFile(sww.filename, 'r')
1452
1453	# Get the variables
1454	x = fid.variables['x'][:]
1455	y = fid.variables['y'][:]
1456	z = fid.variables['elevation'][:]
1457	time = fid.variables['time'][:]
1458	stage = fid.variables['stage'][:]
1459
1460	fid.close()
1461
1462	#Export to ascii/prj files
1463	export_grid(self.domain.get_name(),
1464	quantities = 'elevation',
1465	cellsize = cellsize,
1466	verbose = self.verbose,
1467	format = 'asc')
1468
1469	#Check asc file
1470	ascid = open(ascfile)
1471	lines = ascid.readlines()
1472	ascid.close()
1473
1474	L = lines[2].strip().split()
1475	assert L[0].strip().lower() == 'xllcorner'
1476	assert allclose(float(L[1].strip().lower()), 308500)
1477
1478	L = lines[3].strip().split()
1479	assert L[0].strip().lower() == 'yllcorner'
1480	assert allclose(float(L[1].strip().lower()), 6189000)
1481
1482	#Check grid values
1483	for j in range(5):
1484	L = lines[6+j].strip().split()
1485	y = (4-j) * cellsize
1486	for i in range(5):
1487	assert allclose(float(L[i]), -i*cellsize - y)
1488
1489	#Cleanup
1490	os.remove(prjfile)
1491	os.remove(ascfile)
1492	os.remove(swwfile)
1493
1494	def test_export_gridII(self):
1495	"""
1496	test_export_gridII(self):
1497	Test that sww information can be converted correctly to asc/prj
1498	format readable by e.g. ArcView
1499	"""
1500
1501	import time, os
1502	from Numeric import array, zeros, allclose, Float, concatenate
1503	from Scientific.IO.NetCDF import NetCDFFile
1504
1505	try:
1506	os.remove('teg*.sww')
1507	except:
1508	pass
1509
1510	#Setup
1511	self.domain.set_name('tegII')
1512
1513	swwfile = self.domain.get_name() + '.sww'
1514
1515	self.domain.set_datadir('.')
1516	self.domain.format = 'sww'
1517	self.domain.smooth = True
1518	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1519	self.domain.set_quantity('stage', 1.0)
1520
1521	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1522
1523	sww = get_dataobject(self.domain)
1524	sww.store_connectivity()
1525	sww.store_timestep() #'stage')
1526	self.domain.evolve_to_end(finaltime = 0.01)
1527	sww.store_timestep() #'stage')
1528
1529	cellsize = 0.25
1530	#Check contents
1531	#Get NetCDF
1532
1533	fid = NetCDFFile(sww.filename, 'r')
1534
1535	# Get the variables
1536	x = fid.variables['x'][:]
1537	y = fid.variables['y'][:]
1538	z = fid.variables['elevation'][:]
1539	time = fid.variables['time'][:]
1540	stage = fid.variables['stage'][:]
1541	xmomentum = fid.variables['xmomentum'][:]
1542	ymomentum = fid.variables['ymomentum'][:]
1543
1544	#print 'stage', stage
1545	#print 'xmom', xmomentum
1546	#print 'ymom', ymomentum
1547
1548	fid.close()
1549
1550	#Export to ascii/prj files
1551	if True:
1552	export_grid(self.domain.get_name(),
1553	quantities = ['elevation', 'depth'],
1554	cellsize = cellsize,
1555	verbose = self.verbose,
1556	format = 'asc')
1557
1558	else:
1559	export_grid(self.domain.get_name(),
1560	quantities = ['depth'],
1561	cellsize = cellsize,
1562	verbose = self.verbose,
1563	format = 'asc')
1564
1565
1566	export_grid(self.domain.get_name(),
1567	quantities = ['elevation'],
1568	cellsize = cellsize,
1569	verbose = self.verbose,
1570	format = 'asc')
1571
1572	prjfile = self.domain.get_name() + '_elevation.prj'
1573	ascfile = self.domain.get_name() + '_elevation.asc'
1574
1575	#Check asc file
1576	ascid = open(ascfile)
1577	lines = ascid.readlines()
1578	ascid.close()
1579
1580	L = lines[2].strip().split()
1581	assert L[0].strip().lower() == 'xllcorner'
1582	assert allclose(float(L[1].strip().lower()), 308500)
1583
1584	L = lines[3].strip().split()
1585	assert L[0].strip().lower() == 'yllcorner'
1586	assert allclose(float(L[1].strip().lower()), 6189000)
1587
1588	#print "ascfile", ascfile
1589	#Check grid values
1590	for j in range(5):
1591	L = lines[6+j].strip().split()
1592	y = (4-j) * cellsize
1593	for i in range(5):
1594	#print " -icellsize - y", -icellsize - y
1595	#print "float(L[i])", float(L[i])
1596	assert allclose(float(L[i]), -i*cellsize - y)
1597
1598	#Cleanup
1599	os.remove(prjfile)
1600	os.remove(ascfile)
1601
1602	#Check asc file
1603	ascfile = self.domain.get_name() + '_depth.asc'
1604	prjfile = self.domain.get_name() + '_depth.prj'
1605	ascid = open(ascfile)
1606	lines = ascid.readlines()
1607	ascid.close()
1608
1609	L = lines[2].strip().split()
1610	assert L[0].strip().lower() == 'xllcorner'
1611	assert allclose(float(L[1].strip().lower()), 308500)
1612
1613	L = lines[3].strip().split()
1614	assert L[0].strip().lower() == 'yllcorner'
1615	assert allclose(float(L[1].strip().lower()), 6189000)
1616
1617	#Check grid values
1618	for j in range(5):
1619	L = lines[6+j].strip().split()
1620	y = (4-j) * cellsize
1621	for i in range(5):
1622	#print " -icellsize - y", -icellsize - y
1623	#print "float(L[i])", float(L[i])
1624	assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1625
1626	#Cleanup
1627	os.remove(prjfile)
1628	os.remove(ascfile)
1629	os.remove(swwfile)
1630
1631
1632	def test_export_gridIII(self):
1633	"""
1634	test_export_gridIII
1635	Test that sww information can be converted correctly to asc/prj
1636	format readable by e.g. ArcView
1637	"""
1638
1639	import time, os
1640	from Numeric import array, zeros, allclose, Float, concatenate
1641	from Scientific.IO.NetCDF import NetCDFFile
1642
1643	try:
1644	os.remove('teg*.sww')
1645	except:
1646	pass
1647
1648	#Setup
1649
1650	self.domain.set_name('tegIII')
1651
1652	swwfile = self.domain.get_name() + '.sww'
1653
1654	self.domain.set_datadir('.')
1655	self.domain.format = 'sww'
1656	self.domain.smooth = True
1657	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1658	self.domain.set_quantity('stage', 1.0)
1659
1660	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1661
1662	sww = get_dataobject(self.domain)
1663	sww.store_connectivity()
1664	sww.store_timestep() #'stage')
1665	self.domain.evolve_to_end(finaltime = 0.01)
1666	sww.store_timestep() #'stage')
1667
1668	cellsize = 0.25
1669	#Check contents
1670	#Get NetCDF
1671
1672	fid = NetCDFFile(sww.filename, 'r')
1673
1674	# Get the variables
1675	x = fid.variables['x'][:]
1676	y = fid.variables['y'][:]
1677	z = fid.variables['elevation'][:]
1678	time = fid.variables['time'][:]
1679	stage = fid.variables['stage'][:]
1680
1681	fid.close()
1682
1683	#Export to ascii/prj files
1684	extra_name_out = 'yeah'
1685	if True:
1686	export_grid(self.domain.get_name(),
1687	quantities = ['elevation', 'depth'],
1688	extra_name_out = extra_name_out,
1689	cellsize = cellsize,
1690	verbose = self.verbose,
1691	format = 'asc')
1692
1693	else:
1694	export_grid(self.domain.get_name(),
1695	quantities = ['depth'],
1696	cellsize = cellsize,
1697	verbose = self.verbose,
1698	format = 'asc')
1699
1700
1701	export_grid(self.domain.get_name(),
1702	quantities = ['elevation'],
1703	cellsize = cellsize,
1704	verbose = self.verbose,
1705	format = 'asc')
1706
1707	prjfile = self.domain.get_name() + '_elevation_yeah.prj'
1708	ascfile = self.domain.get_name() + '_elevation_yeah.asc'
1709
1710	#Check asc file
1711	ascid = open(ascfile)
1712	lines = ascid.readlines()
1713	ascid.close()
1714
1715	L = lines[2].strip().split()
1716	assert L[0].strip().lower() == 'xllcorner'
1717	assert allclose(float(L[1].strip().lower()), 308500)
1718
1719	L = lines[3].strip().split()
1720	assert L[0].strip().lower() == 'yllcorner'
1721	assert allclose(float(L[1].strip().lower()), 6189000)
1722
1723	#print "ascfile", ascfile
1724	#Check grid values
1725	for j in range(5):
1726	L = lines[6+j].strip().split()
1727	y = (4-j) * cellsize
1728	for i in range(5):
1729	#print " -icellsize - y", -icellsize - y
1730	#print "float(L[i])", float(L[i])
1731	assert allclose(float(L[i]), -i*cellsize - y)
1732
1733	#Cleanup
1734	os.remove(prjfile)
1735	os.remove(ascfile)
1736
1737	#Check asc file
1738	ascfile = self.domain.get_name() + '_depth_yeah.asc'
1739	prjfile = self.domain.get_name() + '_depth_yeah.prj'
1740	ascid = open(ascfile)
1741	lines = ascid.readlines()
1742	ascid.close()
1743
1744	L = lines[2].strip().split()
1745	assert L[0].strip().lower() == 'xllcorner'
1746	assert allclose(float(L[1].strip().lower()), 308500)
1747
1748	L = lines[3].strip().split()
1749	assert L[0].strip().lower() == 'yllcorner'
1750	assert allclose(float(L[1].strip().lower()), 6189000)
1751
1752	#Check grid values
1753	for j in range(5):
1754	L = lines[6+j].strip().split()
1755	y = (4-j) * cellsize
1756	for i in range(5):
1757	assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1758
1759	#Cleanup
1760	os.remove(prjfile)
1761	os.remove(ascfile)
1762	os.remove(swwfile)
1763
1764	def test_export_grid_bad(self):
1765	"""Test that sww information can be converted correctly to asc/prj
1766	format readable by e.g. ArcView
1767	"""
1768
1769	try:
1770	export_grid('a_small_round-egg',
1771	quantities = ['elevation', 'depth'],
1772	cellsize = 99,
1773	verbose = self.verbose,
1774	format = 'asc')
1775	except IOError:
1776	pass
1777	else:
1778	self.failUnless(0 ==1, 'Bad input did not throw exception error!')
1779
1780	def test_export_grid_parallel(self):
1781	"""Test that sww information can be converted correctly to asc/prj
1782	format readable by e.g. ArcView
1783	"""
1784
1785	import time, os
1786	from Numeric import array, zeros, allclose, Float, concatenate
1787	from Scientific.IO.NetCDF import NetCDFFile
1788
1789	base_name = 'tegp'
1790	#Setup
1791	self.domain.set_name(base_name+'_P0_8')
1792	swwfile = self.domain.get_name() + '.sww'
1793
1794	self.domain.set_datadir('.')
1795	self.domain.format = 'sww'
1796	self.domain.smooth = True
1797	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1798	self.domain.set_quantity('stage', 1.0)
1799
1800	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1801
1802	sww = get_dataobject(self.domain)
1803	sww.store_connectivity()
1804	sww.store_timestep()
1805	self.domain.evolve_to_end(finaltime = 0.0001)
1806	#Setup
1807	self.domain.set_name(base_name+'_P1_8')
1808	swwfile2 = self.domain.get_name() + '.sww'
1809	sww = get_dataobject(self.domain)
1810	sww.store_connectivity()
1811	sww.store_timestep()
1812	self.domain.evolve_to_end(finaltime = 0.0002)
1813	sww.store_timestep()
1814
1815	cellsize = 0.25
1816	#Check contents
1817	#Get NetCDF
1818
1819	fid = NetCDFFile(sww.filename, 'r')
1820
1821	# Get the variables
1822	x = fid.variables['x'][:]
1823	y = fid.variables['y'][:]
1824	z = fid.variables['elevation'][:]
1825	time = fid.variables['time'][:]
1826	stage = fid.variables['stage'][:]
1827
1828	fid.close()
1829
1830	#Export to ascii/prj files
1831	extra_name_out = 'yeah'
1832	export_grid(base_name,
1833	quantities = ['elevation', 'depth'],
1834	extra_name_out = extra_name_out,
1835	cellsize = cellsize,
1836	verbose = self.verbose,
1837	format = 'asc')
1838
1839	prjfile = base_name + '_P0_8_elevation_yeah.prj'
1840	ascfile = base_name + '_P0_8_elevation_yeah.asc'
1841	#Check asc file
1842	ascid = open(ascfile)
1843	lines = ascid.readlines()
1844	ascid.close()
1845	#Check grid values
1846	for j in range(5):
1847	L = lines[6+j].strip().split()
1848	y = (4-j) * cellsize
1849	for i in range(5):
1850	#print " -icellsize - y", -icellsize - y
1851	#print "float(L[i])", float(L[i])
1852	assert allclose(float(L[i]), -i*cellsize - y)
1853	#Cleanup
1854	os.remove(prjfile)
1855	os.remove(ascfile)
1856
1857	prjfile = base_name + '_P1_8_elevation_yeah.prj'
1858	ascfile = base_name + '_P1_8_elevation_yeah.asc'
1859	#Check asc file
1860	ascid = open(ascfile)
1861	lines = ascid.readlines()
1862	ascid.close()
1863	#Check grid values
1864	for j in range(5):
1865	L = lines[6+j].strip().split()
1866	y = (4-j) * cellsize
1867	for i in range(5):
1868	#print " -icellsize - y", -icellsize - y
1869	#print "float(L[i])", float(L[i])
1870	assert allclose(float(L[i]), -i*cellsize - y)
1871	#Cleanup
1872	os.remove(prjfile)
1873	os.remove(ascfile)
1874	os.remove(swwfile)
1875
1876	#Check asc file
1877	ascfile = base_name + '_P0_8_depth_yeah.asc'
1878	prjfile = base_name + '_P0_8_depth_yeah.prj'
1879	ascid = open(ascfile)
1880	lines = ascid.readlines()
1881	ascid.close()
1882	#Check grid values
1883	for j in range(5):
1884	L = lines[6+j].strip().split()
1885	y = (4-j) * cellsize
1886	for i in range(5):
1887	assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1888	#Cleanup
1889	os.remove(prjfile)
1890	os.remove(ascfile)
1891
1892	#Check asc file
1893	ascfile = base_name + '_P1_8_depth_yeah.asc'
1894	prjfile = base_name + '_P1_8_depth_yeah.prj'
1895	ascid = open(ascfile)
1896	lines = ascid.readlines()
1897	ascid.close()
1898	#Check grid values
1899	for j in range(5):
1900	L = lines[6+j].strip().split()
1901	y = (4-j) * cellsize
1902	for i in range(5):
1903	assert allclose(float(L[i]), 1 - (-i*cellsize - y))
1904	#Cleanup
1905	os.remove(prjfile)
1906	os.remove(ascfile)
1907	os.remove(swwfile2)
1908
1909	def test_sww2dem_larger(self):
1910	"""Test that sww information can be converted correctly to asc/prj
1911	format readable by e.g. ArcView. Here:
1912
1913	ncols 11
1914	nrows 11
1915	xllcorner 308500
1916	yllcorner 6189000
1917	cellsize 10.000000
1918	NODATA_value -9999
1919	-100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1920	-90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1921	-80 -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1922	-70 -80 -90 -100 -110 -120 -130 -140 -150 -160 -170
1923	-60 -70 -80 -90 -100 -110 -120 -130 -140 -150 -160
1924	-50 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150
1925	-40 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140
1926	-30 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130
1927	-20 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120
1928	-10 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110
1929	0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100
1930
1931	"""
1932
1933	import time, os
1934	from Numeric import array, zeros, allclose, Float, concatenate
1935	from Scientific.IO.NetCDF import NetCDFFile
1936
1937	#Setup
1938
1939	from mesh_factory import rectangular
1940
1941	#Create basic mesh (100m x 100m)
1942	points, vertices, boundary = rectangular(2, 2, 100, 100)
1943
1944	#Create shallow water domain
1945	domain = Domain(points, vertices, boundary)
1946	domain.default_order = 2
1947
1948	domain.set_name('datatest')
1949
1950	prjfile = domain.get_name() + '_elevation.prj'
1951	ascfile = domain.get_name() + '_elevation.asc'
1952	swwfile = domain.get_name() + '.sww'
1953
1954	domain.set_datadir('.')
1955	domain.format = 'sww'
1956	domain.smooth = True
1957	domain.geo_reference = Geo_reference(56, 308500, 6189000)
1958
1959	#
1960	domain.set_quantity('elevation', lambda x,y: -x-y)
1961	domain.set_quantity('stage', 0)
1962
1963	B = Transmissive_boundary(domain)
1964	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1965
1966
1967	#
1968	sww = get_dataobject(domain)
1969	sww.store_connectivity()
1970	sww.store_timestep()
1971
1972	domain.tight_slope_limiters = 1
1973	domain.evolve_to_end(finaltime = 0.01)
1974	sww.store_timestep()
1975
1976	cellsize = 10 #10m grid
1977
1978
1979	#Check contents
1980	#Get NetCDF
1981
1982	fid = NetCDFFile(sww.filename, 'r')
1983
1984	# Get the variables
1985	x = fid.variables['x'][:]
1986	y = fid.variables['y'][:]
1987	z = fid.variables['elevation'][:]
1988	time = fid.variables['time'][:]
1989	stage = fid.variables['stage'][:]
1990
1991
1992	#Export to ascii/prj files
1993	sww2dem(domain.get_name(),
1994	quantity = 'elevation',
1995	cellsize = cellsize,
1996	verbose = self.verbose,
1997	format = 'asc')
1998
1999
2000	#Check prj (meta data)
2001	prjid = open(prjfile)
2002	lines = prjid.readlines()
2003	prjid.close()
2004
2005	L = lines[0].strip().split()
2006	assert L[0].strip().lower() == 'projection'
2007	assert L[1].strip().lower() == 'utm'
2008
2009	L = lines[1].strip().split()
2010	assert L[0].strip().lower() == 'zone'
2011	assert L[1].strip().lower() == '56'
2012
2013	L = lines[2].strip().split()
2014	assert L[0].strip().lower() == 'datum'
2015	assert L[1].strip().lower() == 'wgs84'
2016
2017	L = lines[3].strip().split()
2018	assert L[0].strip().lower() == 'zunits'
2019	assert L[1].strip().lower() == 'no'
2020
2021	L = lines[4].strip().split()
2022	assert L[0].strip().lower() == 'units'
2023	assert L[1].strip().lower() == 'meters'
2024
2025	L = lines[5].strip().split()
2026	assert L[0].strip().lower() == 'spheroid'
2027	assert L[1].strip().lower() == 'wgs84'
2028
2029	L = lines[6].strip().split()
2030	assert L[0].strip().lower() == 'xshift'
2031	assert L[1].strip().lower() == '500000'
2032
2033	L = lines[7].strip().split()
2034	assert L[0].strip().lower() == 'yshift'
2035	assert L[1].strip().lower() == '10000000'
2036
2037	L = lines[8].strip().split()
2038	assert L[0].strip().lower() == 'parameters'
2039
2040
2041	#Check asc file
2042	ascid = open(ascfile)
2043	lines = ascid.readlines()
2044	ascid.close()
2045
2046	L = lines[0].strip().split()
2047	assert L[0].strip().lower() == 'ncols'
2048	assert L[1].strip().lower() == '11'
2049
2050	L = lines[1].strip().split()
2051	assert L[0].strip().lower() == 'nrows'
2052	assert L[1].strip().lower() == '11'
2053
2054	L = lines[2].strip().split()
2055	assert L[0].strip().lower() == 'xllcorner'
2056	assert allclose(float(L[1].strip().lower()), 308500)
2057
2058	L = lines[3].strip().split()
2059	assert L[0].strip().lower() == 'yllcorner'
2060	assert allclose(float(L[1].strip().lower()), 6189000)
2061
2062	L = lines[4].strip().split()
2063	assert L[0].strip().lower() == 'cellsize'
2064	assert allclose(float(L[1].strip().lower()), cellsize)
2065
2066	L = lines[5].strip().split()
2067	assert L[0].strip() == 'NODATA_value'
2068	assert L[1].strip().lower() == '-9999'
2069
2070	#Check grid values (FIXME: Use same strategy for other sww2dem tests)
2071	for i, line in enumerate(lines[6:]):
2072	for j, value in enumerate( line.split() ):
2073	#assert allclose(float(value), -(10-i+j)*cellsize)
2074	assert float(value) == -(10-i+j)*cellsize
2075
2076
2077	fid.close()
2078
2079	#Cleanup
2080	os.remove(prjfile)
2081	os.remove(ascfile)
2082	os.remove(swwfile)
2083
2084
2085
2086
2087	def test_sww2dem_boundingbox(self):
2088	"""Test that sww information can be converted correctly to asc/prj
2089	format readable by e.g. ArcView.
2090	This will test that mesh can be restricted by bounding box
2091
2092	Original extent is 100m x 100m:
2093
2094	Eastings: 308500 - 308600
2095	Northings: 6189000 - 6189100
2096
2097	Bounding box changes this to the 50m x 50m square defined by
2098
2099	Eastings: 308530 - 308570
2100	Northings: 6189050 - 6189100
2101
2102	The cropped values should be
2103
2104	-130 -140 -150 -160 -170
2105	-120 -130 -140 -150 -160
2106	-110 -120 -130 -140 -150
2107	-100 -110 -120 -130 -140
2108	-90 -100 -110 -120 -130
2109	-80 -90 -100 -110 -120
2110
2111	and the new lower reference point should be
2112	Eastings: 308530
2113	Northings: 6189050
2114
2115	Original dataset is the same as in test_sww2dem_larger()
2116
2117	"""
2118
2119	import time, os
2120	from Numeric import array, zeros, allclose, Float, concatenate
2121	from Scientific.IO.NetCDF import NetCDFFile
2122
2123	#Setup
2124
2125	from mesh_factory import rectangular
2126
2127	#Create basic mesh (100m x 100m)
2128	points, vertices, boundary = rectangular(2, 2, 100, 100)
2129
2130	#Create shallow water domain
2131	domain = Domain(points, vertices, boundary)
2132	domain.default_order = 2
2133
2134	domain.set_name('datatest')
2135
2136	prjfile = domain.get_name() + '_elevation.prj'
2137	ascfile = domain.get_name() + '_elevation.asc'
2138	swwfile = domain.get_name() + '.sww'
2139
2140	domain.set_datadir('.')
2141	domain.format = 'sww'
2142	domain.smooth = True
2143	domain.geo_reference = Geo_reference(56, 308500, 6189000)
2144
2145	#
2146	domain.set_quantity('elevation', lambda x,y: -x-y)
2147	domain.set_quantity('stage', 0)
2148
2149	B = Transmissive_boundary(domain)
2150	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
2151
2152
2153	#
2154	sww = get_dataobject(domain)
2155	sww.store_connectivity()
2156	sww.store_timestep()
2157
2158	#domain.tight_slope_limiters = 1
2159	domain.evolve_to_end(finaltime = 0.01)
2160	sww.store_timestep()
2161
2162	cellsize = 10 #10m grid
2163
2164
2165	#Check contents
2166	#Get NetCDF
2167
2168	fid = NetCDFFile(sww.filename, 'r')
2169
2170	# Get the variables
2171	x = fid.variables['x'][:]
2172	y = fid.variables['y'][:]
2173	z = fid.variables['elevation'][:]
2174	time = fid.variables['time'][:]
2175	stage = fid.variables['stage'][:]
2176
2177
2178	#Export to ascii/prj files
2179	sww2dem(domain.get_name(),
2180	quantity = 'elevation',
2181	cellsize = cellsize,
2182	easting_min = 308530,
2183	easting_max = 308570,
2184	northing_min = 6189050,
2185	northing_max = 6189100,
2186	verbose = self.verbose,
2187	format = 'asc')
2188
2189	fid.close()
2190
2191
2192	#Check prj (meta data)
2193	prjid = open(prjfile)
2194	lines = prjid.readlines()
2195	prjid.close()
2196
2197	L = lines[0].strip().split()
2198	assert L[0].strip().lower() == 'projection'
2199	assert L[1].strip().lower() == 'utm'
2200
2201	L = lines[1].strip().split()
2202	assert L[0].strip().lower() == 'zone'
2203	assert L[1].strip().lower() == '56'
2204
2205	L = lines[2].strip().split()
2206	assert L[0].strip().lower() == 'datum'
2207	assert L[1].strip().lower() == 'wgs84'
2208
2209	L = lines[3].strip().split()
2210	assert L[0].strip().lower() == 'zunits'
2211	assert L[1].strip().lower() == 'no'
2212
2213	L = lines[4].strip().split()
2214	assert L[0].strip().lower() == 'units'
2215	assert L[1].strip().lower() == 'meters'
2216
2217	L = lines[5].strip().split()
2218	assert L[0].strip().lower() == 'spheroid'
2219	assert L[1].strip().lower() == 'wgs84'
2220
2221	L = lines[6].strip().split()
2222	assert L[0].strip().lower() == 'xshift'
2223	assert L[1].strip().lower() == '500000'
2224
2225	L = lines[7].strip().split()
2226	assert L[0].strip().lower() == 'yshift'
2227	assert L[1].strip().lower() == '10000000'
2228
2229	L = lines[8].strip().split()
2230	assert L[0].strip().lower() == 'parameters'
2231
2232
2233	#Check asc file
2234	ascid = open(ascfile)
2235	lines = ascid.readlines()
2236	ascid.close()
2237
2238	L = lines[0].strip().split()
2239	assert L[0].strip().lower() == 'ncols'
2240	assert L[1].strip().lower() == '5'
2241
2242	L = lines[1].strip().split()
2243	assert L[0].strip().lower() == 'nrows'
2244	assert L[1].strip().lower() == '6'
2245
2246	L = lines[2].strip().split()
2247	assert L[0].strip().lower() == 'xllcorner'
2248	assert allclose(float(L[1].strip().lower()), 308530)
2249
2250	L = lines[3].strip().split()
2251	assert L[0].strip().lower() == 'yllcorner'
2252	assert allclose(float(L[1].strip().lower()), 6189050)
2253
2254	L = lines[4].strip().split()
2255	assert L[0].strip().lower() == 'cellsize'
2256	assert allclose(float(L[1].strip().lower()), cellsize)
2257
2258	L = lines[5].strip().split()
2259	assert L[0].strip() == 'NODATA_value'
2260	assert L[1].strip().lower() == '-9999'
2261
2262	#Check grid values
2263	for i, line in enumerate(lines[6:]):
2264	for j, value in enumerate( line.split() ):
2265	#assert float(value) == -(10-i+j)*cellsize
2266	assert float(value) == -(10-i+j+3)*cellsize
2267
2268
2269
2270	#Cleanup
2271	os.remove(prjfile)
2272	os.remove(ascfile)
2273	os.remove(swwfile)
2274
2275
2276
2277	def test_sww2dem_asc_stage_reduction(self):
2278	"""Test that sww information can be converted correctly to asc/prj
2279	format readable by e.g. ArcView
2280
2281	This tests the reduction of quantity stage using min
2282	"""
2283
2284	import time, os
2285	from Numeric import array, zeros, allclose, Float, concatenate
2286	from Scientific.IO.NetCDF import NetCDFFile
2287
2288	#Setup
2289	self.domain.set_name('datatest')
2290
2291	prjfile = self.domain.get_name() + '_stage.prj'
2292	ascfile = self.domain.get_name() + '_stage.asc'
2293	swwfile = self.domain.get_name() + '.sww'
2294
2295	self.domain.set_datadir('.')
2296	self.domain.format = 'sww'
2297	self.domain.smooth = True
2298	self.domain.set_quantity('elevation', lambda x,y: -x-y)
2299
2300	self.domain.geo_reference = Geo_reference(56,308500,6189000)
2301
2302
2303	sww = get_dataobject(self.domain)
2304	sww.store_connectivity()
2305	sww.store_timestep()
2306
2307	#self.domain.tight_slope_limiters = 1
2308	self.domain.evolve_to_end(finaltime = 0.01)
2309	sww.store_timestep()
2310
2311	cellsize = 0.25
2312	#Check contents
2313	#Get NetCDF
2314
2315	fid = NetCDFFile(sww.filename, 'r')
2316
2317	# Get the variables
2318	x = fid.variables['x'][:]
2319	y = fid.variables['y'][:]
2320	z = fid.variables['elevation'][:]
2321	time = fid.variables['time'][:]
2322	stage = fid.variables['stage'][:]
2323
2324
2325	#Export to ascii/prj files
2326	sww2dem(self.domain.get_name(),
2327	quantity = 'stage',
2328	cellsize = cellsize,
2329	reduction = min,
2330	format = 'asc',
2331	verbose=self.verbose)
2332
2333
2334	#Check asc file
2335	ascid = open(ascfile)
2336	lines = ascid.readlines()
2337	ascid.close()
2338
2339	L = lines[0].strip().split()
2340	assert L[0].strip().lower() == 'ncols'
2341	assert L[1].strip().lower() == '5'
2342
2343	L = lines[1].strip().split()
2344	assert L[0].strip().lower() == 'nrows'
2345	assert L[1].strip().lower() == '5'
2346
2347	L = lines[2].strip().split()
2348	assert L[0].strip().lower() == 'xllcorner'
2349	assert allclose(float(L[1].strip().lower()), 308500)
2350
2351	L = lines[3].strip().split()
2352	assert L[0].strip().lower() == 'yllcorner'
2353	assert allclose(float(L[1].strip().lower()), 6189000)
2354
2355	L = lines[4].strip().split()
2356	assert L[0].strip().lower() == 'cellsize'
2357	assert allclose(float(L[1].strip().lower()), cellsize)
2358
2359	L = lines[5].strip().split()
2360	assert L[0].strip() == 'NODATA_value'
2361	assert L[1].strip().lower() == '-9999'
2362
2363
2364	#Check grid values (where applicable)
2365	for j in range(5):
2366	if j%2 == 0:
2367	L = lines[6+j].strip().split()
2368	jj = 4-j
2369	for i in range(5):
2370	if i%2 == 0:
2371	index = jj/2 + i/2*3
2372	val0 = stage[0,index]
2373	val1 = stage[1,index]
2374
2375	#print i, j, index, ':', L[i], val0, val1
2376	assert allclose(float(L[i]), min(val0, val1))
2377
2378
2379	fid.close()
2380
2381	#Cleanup
2382	os.remove(prjfile)
2383	os.remove(ascfile)
2384	os.remove(swwfile)
2385
2386
2387
2388	def test_sww2dem_asc_derived_quantity(self):
2389	"""Test that sww information can be converted correctly to asc/prj
2390	format readable by e.g. ArcView
2391
2392	This tests the use of derived quantities
2393	"""
2394
2395	import time, os
2396	from Numeric import array, zeros, allclose, Float, concatenate
2397	from Scientific.IO.NetCDF import NetCDFFile
2398
2399	#Setup
2400	self.domain.set_name('datatest')
2401
2402	prjfile = self.domain.get_name() + '_depth.prj'
2403	ascfile = self.domain.get_name() + '_depth.asc'
2404	swwfile = self.domain.get_name() + '.sww'
2405
2406	self.domain.set_datadir('.')
2407	self.domain.format = 'sww'
2408	self.domain.smooth = True
2409	self.domain.set_quantity('elevation', lambda x,y: -x-y)
2410	self.domain.set_quantity('stage', 0.0)
2411
2412	self.domain.geo_reference = Geo_reference(56,308500,6189000)
2413
2414
2415	sww = get_dataobject(self.domain)
2416	sww.store_connectivity()
2417	sww.store_timestep()
2418
2419	#self.domain.tight_slope_limiters = 1
2420	self.domain.evolve_to_end(finaltime = 0.01)
2421	sww.store_timestep()
2422
2423	cellsize = 0.25
2424	#Check contents
2425	#Get NetCDF
2426
2427	fid = NetCDFFile(sww.filename, 'r')
2428
2429	# Get the variables
2430	x = fid.variables['x'][:]
2431	y = fid.variables['y'][:]
2432	z = fid.variables['elevation'][:]
2433	time = fid.variables['time'][:]
2434	stage = fid.variables['stage'][:]
2435
2436
2437	#Export to ascii/prj files
2438	sww2dem(self.domain.get_name(),
2439	basename_out = 'datatest_depth',
2440	quantity = 'stage - elevation',
2441	cellsize = cellsize,
2442	reduction = min,
2443	format = 'asc',
2444	verbose = self.verbose)
2445
2446
2447	#Check asc file
2448	ascid = open(ascfile)
2449	lines = ascid.readlines()
2450	ascid.close()
2451
2452	L = lines[0].strip().split()
2453	assert L[0].strip().lower() == 'ncols'
2454	assert L[1].strip().lower() == '5'
2455
2456	L = lines[1].strip().split()
2457	assert L[0].strip().lower() == 'nrows'
2458	assert L[1].strip().lower() == '5'
2459
2460	L = lines[2].strip().split()
2461	assert L[0].strip().lower() == 'xllcorner'
2462	assert allclose(float(L[1].strip().lower()), 308500)
2463
2464	L = lines[3].strip().split()
2465	assert L[0].strip().lower() == 'yllcorner'
2466	assert allclose(float(L[1].strip().lower()), 6189000)
2467
2468	L = lines[4].strip().split()
2469	assert L[0].strip().lower() == 'cellsize'
2470	assert allclose(float(L[1].strip().lower()), cellsize)
2471
2472	L = lines[5].strip().split()
2473	assert L[0].strip() == 'NODATA_value'
2474	assert L[1].strip().lower() == '-9999'
2475
2476
2477	#Check grid values (where applicable)
2478	for j in range(5):
2479	if j%2 == 0:
2480	L = lines[6+j].strip().split()
2481	jj = 4-j
2482	for i in range(5):
2483	if i%2 == 0:
2484	index = jj/2 + i/2*3
2485	val0 = stage[0,index] - z[index]
2486	val1 = stage[1,index] - z[index]
2487
2488	#print i, j, index, ':', L[i], val0, val1
2489	assert allclose(float(L[i]), min(val0, val1))
2490
2491
2492	fid.close()
2493
2494	#Cleanup
2495	os.remove(prjfile)
2496	os.remove(ascfile)
2497	os.remove(swwfile)
2498
2499
2500
2501
2502
2503	def test_sww2dem_asc_missing_points(self):
2504	"""Test that sww information can be converted correctly to asc/prj
2505	format readable by e.g. ArcView
2506
2507	This test includes the writing of missing values
2508	"""
2509
2510	import time, os
2511	from Numeric import array, zeros, allclose, Float, concatenate
2512	from Scientific.IO.NetCDF import NetCDFFile
2513
2514	#Setup mesh not coinciding with rectangle.
2515	#This will cause missing values to occur in gridded data
2516
2517
2518	points = [ [1.0, 1.0],
2519	[0.5, 0.5], [1.0, 0.5],
2520	[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
2521
2522	vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
2523
2524	#Create shallow water domain
2525	domain = Domain(points, vertices)
2526	domain.default_order=2
2527
2528
2529	#Set some field values
2530	domain.set_quantity('elevation', lambda x,y: -x-y)
2531	domain.set_quantity('friction', 0.03)
2532
2533
2534	######################
2535	# Boundary conditions
2536	B = Transmissive_boundary(domain)
2537	domain.set_boundary( {'exterior': B} )
2538
2539
2540	######################
2541	#Initial condition - with jumps
2542
2543	bed = domain.quantities['elevation'].vertex_values
2544	stage = zeros(bed.shape, Float)
2545
2546	h = 0.3
2547	for i in range(stage.shape[0]):
2548	if i % 2 == 0:
2549	stage[i,:] = bed[i,:] + h
2550	else:
2551	stage[i,:] = bed[i,:]
2552
2553	domain.set_quantity('stage', stage)
2554	domain.distribute_to_vertices_and_edges()
2555
2556	domain.set_name('datatest')
2557
2558	prjfile = domain.get_name() + '_elevation.prj'
2559	ascfile = domain.get_name() + '_elevation.asc'
2560	swwfile = domain.get_name() + '.sww'
2561
2562	domain.set_datadir('.')
2563	domain.format = 'sww'
2564	domain.smooth = True
2565
2566	domain.geo_reference = Geo_reference(56,308500,6189000)
2567
2568	sww = get_dataobject(domain)
2569	sww.store_connectivity()
2570	sww.store_timestep()
2571
2572	cellsize = 0.25
2573	#Check contents
2574	#Get NetCDF
2575
2576	fid = NetCDFFile(swwfile, 'r')
2577
2578	# Get the variables
2579	x = fid.variables['x'][:]
2580	y = fid.variables['y'][:]
2581	z = fid.variables['elevation'][:]
2582	time = fid.variables['time'][:]
2583
2584	try:
2585	geo_reference = Geo_reference(NetCDFObject=fid)
2586	except AttributeError, e:
2587	geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
2588
2589	#Export to ascii/prj files
2590	sww2dem(domain.get_name(),
2591	quantity = 'elevation',
2592	cellsize = cellsize,
2593	verbose = self.verbose,
2594	format = 'asc')
2595
2596
2597	#Check asc file
2598	ascid = open(ascfile)
2599	lines = ascid.readlines()
2600	ascid.close()
2601
2602	L = lines[0].strip().split()
2603	assert L[0].strip().lower() == 'ncols'
2604	assert L[1].strip().lower() == '5'
2605
2606	L = lines[1].strip().split()
2607	assert L[0].strip().lower() == 'nrows'
2608	assert L[1].strip().lower() == '5'
2609
2610	L = lines[2].strip().split()
2611	assert L[0].strip().lower() == 'xllcorner'
2612	assert allclose(float(L[1].strip().lower()), 308500)
2613
2614	L = lines[3].strip().split()
2615	assert L[0].strip().lower() == 'yllcorner'
2616	assert allclose(float(L[1].strip().lower()), 6189000)
2617
2618	L = lines[4].strip().split()
2619	assert L[0].strip().lower() == 'cellsize'
2620	assert allclose(float(L[1].strip().lower()), cellsize)
2621
2622	L = lines[5].strip().split()
2623	assert L[0].strip() == 'NODATA_value'
2624	assert L[1].strip().lower() == '-9999'
2625
2626	#Check grid values
2627	for j in range(5):
2628	L = lines[6+j].strip().split()
2629	assert len(L) == 5
2630	y = (4-j) * cellsize
2631
2632	for i in range(5):
2633	#print i
2634	if i+j >= 4:
2635	assert allclose(float(L[i]), -i*cellsize - y)
2636	else:
2637	#Missing values
2638	assert allclose(float(L[i]), -9999)
2639
2640
2641
2642	fid.close()
2643
2644	#Cleanup
2645	os.remove(prjfile)
2646	os.remove(ascfile)
2647	os.remove(swwfile)
2648
2649	def test_sww2ers_simple(self):
2650	"""Test that sww information can be converted correctly to asc/prj
2651	format readable by e.g. ArcView
2652	"""
2653
2654	import time, os
2655	from Numeric import array, zeros, allclose, Float, concatenate
2656	from Scientific.IO.NetCDF import NetCDFFile
2657
2658
2659	NODATA_value = 1758323
2660
2661	#Setup
2662	self.domain.set_name('datatest')
2663
2664	headerfile = self.domain.get_name() + '.ers'
2665	swwfile = self.domain.get_name() + '.sww'
2666
2667	self.domain.set_datadir('.')
2668	self.domain.format = 'sww'
2669	self.domain.smooth = True
2670	self.domain.set_quantity('elevation', lambda x,y: -x-y)
2671
2672	self.domain.geo_reference = Geo_reference(56,308500,6189000)
2673
2674	sww = get_dataobject(self.domain)
2675	sww.store_connectivity()
2676	sww.store_timestep()
2677
2678	#self.domain.tight_slope_limiters = 1
2679	self.domain.evolve_to_end(finaltime = 0.01)
2680	sww.store_timestep()
2681
2682	cellsize = 0.25
2683	#Check contents
2684	#Get NetCDF
2685
2686	fid = NetCDFFile(sww.filename, 'r')
2687
2688	# Get the variables
2689	x = fid.variables['x'][:]
2690	y = fid.variables['y'][:]
2691	z = fid.variables['elevation'][:]
2692	time = fid.variables['time'][:]
2693	stage = fid.variables['stage'][:]
2694
2695
2696	#Export to ers files
2697	sww2dem(self.domain.get_name(),
2698	quantity = 'elevation',
2699	cellsize = cellsize,
2700	NODATA_value = NODATA_value,
2701	verbose = self.verbose,
2702	format = 'ers')
2703
2704	#Check header data
2705	from ermapper_grids import read_ermapper_header, read_ermapper_data
2706
2707	header = read_ermapper_header(self.domain.get_name() + '_elevation.ers')
2708	#print header
2709	assert header['projection'].lower() == '"utm-56"'
2710	assert header['datum'].lower() == '"wgs84"'
2711	assert header['units'].lower() == '"meters"'
2712	assert header['value'].lower() == '"elevation"'
2713	assert header['xdimension'] == '0.25'
2714	assert header['ydimension'] == '0.25'
2715	assert float(header['eastings']) == 308500.0 #xllcorner
2716	assert float(header['northings']) == 6189000.0 #yllcorner
2717	assert int(header['nroflines']) == 5
2718	assert int(header['nrofcellsperline']) == 5
2719	assert int(header['nullcellvalue']) == NODATA_value
2720	#FIXME - there is more in the header
2721
2722
2723	#Check grid data
2724	grid = read_ermapper_data(self.domain.get_name() + '_elevation')
2725
2726	#FIXME (Ole): Why is this the desired reference grid for -x-y?
2727	ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
2728	-1, -1.25, -1.5, -1.75, -2.0,
2729	-0.75, -1.0, -1.25, -1.5, -1.75,
2730	-0.5, -0.75, -1.0, -1.25, -1.5,
2731	-0.25, -0.5, -0.75, -1.0, -1.25]
2732
2733
2734	#print grid
2735	assert allclose(grid, ref_grid)
2736
2737	fid.close()
2738
2739	#Cleanup
2740	#FIXME the file clean-up doesn't work (eg Permission Denied Error)
2741	#Done (Ole) - it was because sww2ers didn't close it's sww file
2742	os.remove(sww.filename)
2743	os.remove(self.domain.get_name() + '_elevation')
2744	os.remove(self.domain.get_name() + '_elevation.ers')
2745
2746
2747
2748	def test_sww2pts_centroids(self):
2749	"""Test that sww information can be converted correctly to pts data at specified coordinates
2750	- in this case, the centroids.
2751	"""
2752
2753	import time, os
2754	from Numeric import array, zeros, allclose, Float, concatenate, NewAxis
2755	from Scientific.IO.NetCDF import NetCDFFile
2756	# Used for points that lie outside mesh
2757	NODATA_value = 1758323
2758
2759	# Setup
2760	self.domain.set_name('datatest')
2761
2762	ptsfile = self.domain.get_name() + '_elevation.pts'
2763	swwfile = self.domain.get_name() + '.sww'
2764
2765	self.domain.set_datadir('.')
2766	self.domain.format = 'sww'
2767	self.smooth = True #self.set_store_vertices_uniquely(False)
2768	self.domain.set_quantity('elevation', lambda x,y: -x-y)
2769
2770	self.domain.geo_reference = Geo_reference(56,308500,6189000)
2771
2772	sww = get_dataobject(self.domain)
2773	sww.store_connectivity()
2774	sww.store_timestep()
2775
2776	#self.domain.tight_slope_limiters = 1
2777	self.domain.evolve_to_end(finaltime = 0.01)
2778	sww.store_timestep()
2779
2780	# Check contents in NetCDF
2781	fid = NetCDFFile(sww.filename, 'r')
2782
2783	# Get the variables
2784	x = fid.variables['x'][:]
2785	y = fid.variables['y'][:]
2786	elevation = fid.variables['elevation'][:]
2787	time = fid.variables['time'][:]
2788	stage = fid.variables['stage'][:]
2789
2790	volumes = fid.variables['volumes'][:]
2791
2792
2793	# Invoke interpolation for vertex points
2794	points = concatenate( (x[:,NewAxis],y[:,NewAxis]), axis=1 )
2795	sww2pts(self.domain.get_name(),
2796	quantity = 'elevation',
2797	data_points = points,
2798	NODATA_value = NODATA_value,
2799	verbose = self.verbose)
2800	ref_point_values = elevation
2801	point_values = Geospatial_data(ptsfile).get_attributes()
2802	#print 'P', point_values
2803	#print 'Ref', ref_point_values
2804	assert allclose(point_values, ref_point_values)
2805
2806
2807
2808	# Invoke interpolation for centroids
2809	points = self.domain.get_centroid_coordinates()
2810	#print points
2811	sww2pts(self.domain.get_name(),
2812	quantity = 'elevation',
2813	data_points = points,
2814	NODATA_value = NODATA_value,
2815	verbose = self.verbose)
2816	ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5] #At centroids
2817
2818
2819	point_values = Geospatial_data(ptsfile).get_attributes()
2820	#print 'P', point_values
2821	#print 'Ref', ref_point_values
2822	assert allclose(point_values, ref_point_values)
2823
2824
2825
2826	fid.close()
2827
2828	#Cleanup
2829	os.remove(sww.filename)
2830	os.remove(ptsfile)
2831
2832
2833
2834
2835	def test_ferret2sww1(self):
2836	"""Test that georeferencing etc works when converting from
2837	ferret format (lat/lon) to sww format (UTM)
2838	"""
2839	from Scientific.IO.NetCDF import NetCDFFile
2840	import os, sys
2841
2842	#The test file has
2843	# LON = 150.66667, 150.83334, 151, 151.16667
2844	# LAT = -34.5, -34.33333, -34.16667, -34 ;
2845	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2846	#
2847	# First value (index=0) in small_ha.nc is 0.3400644 cm,
2848	# Fourth value (index==3) is -6.50198 cm
2849
2850
2851
2852	#Read
2853	from anuga.coordinate_transforms.redfearn import redfearn
2854	#fid = NetCDFFile(self.test_MOST_file)
2855	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2856	first_value = fid.variables['HA'][:][0,0,0]
2857	fourth_value = fid.variables['HA'][:][0,0,3]
2858	fid.close()
2859
2860
2861	#Call conversion (with zero origin)
2862	#ferret2sww('small', verbose=False,
2863	# origin = (56, 0, 0))
2864	ferret2sww(self.test_MOST_file, verbose=self.verbose,
2865	origin = (56, 0, 0))
2866
2867	#Work out the UTM coordinates for first point
2868	zone, e, n = redfearn(-34.5, 150.66667)
2869	#print zone, e, n
2870
2871	#Read output file 'small.sww'
2872	#fid = NetCDFFile('small.sww')
2873	fid = NetCDFFile(self.test_MOST_file + '.sww')
2874
2875	x = fid.variables['x'][:]
2876	y = fid.variables['y'][:]
2877
2878	#Check that first coordinate is correctly represented
2879	assert allclose(x[0], e)
2880	assert allclose(y[0], n)
2881
2882	#Check first value
2883	stage = fid.variables['stage'][:]
2884	xmomentum = fid.variables['xmomentum'][:]
2885	ymomentum = fid.variables['ymomentum'][:]
2886
2887	#print ymomentum
2888
2889	assert allclose(stage[0,0], first_value/100) #Meters
2890
2891	#Check fourth value
2892	assert allclose(stage[0,3], fourth_value/100) #Meters
2893
2894	fid.close()
2895
2896	#Cleanup
2897	import os
2898	os.remove(self.test_MOST_file + '.sww')
2899
2900
2901
2902	def test_ferret2sww_zscale(self):
2903	"""Test that zscale workse
2904	"""
2905	from Scientific.IO.NetCDF import NetCDFFile
2906	import os, sys
2907
2908	#The test file has
2909	# LON = 150.66667, 150.83334, 151, 151.16667
2910	# LAT = -34.5, -34.33333, -34.16667, -34 ;
2911	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2912	#
2913	# First value (index=0) in small_ha.nc is 0.3400644 cm,
2914	# Fourth value (index==3) is -6.50198 cm
2915
2916
2917	#Read
2918	from anuga.coordinate_transforms.redfearn import redfearn
2919	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2920	first_value = fid.variables['HA'][:][0,0,0]
2921	fourth_value = fid.variables['HA'][:][0,0,3]
2922	fid.close()
2923
2924	#Call conversion (with no scaling)
2925	ferret2sww(self.test_MOST_file, verbose=self.verbose,
2926	origin = (56, 0, 0))
2927
2928	#Work out the UTM coordinates for first point
2929	fid = NetCDFFile(self.test_MOST_file + '.sww')
2930
2931	#Check values
2932	stage_1 = fid.variables['stage'][:]
2933	xmomentum_1 = fid.variables['xmomentum'][:]
2934	ymomentum_1 = fid.variables['ymomentum'][:]
2935
2936	assert allclose(stage_1[0,0], first_value/100) #Meters
2937	assert allclose(stage_1[0,3], fourth_value/100) #Meters
2938
2939	fid.close()
2940
2941	#Call conversion (with scaling)
2942	ferret2sww(self.test_MOST_file,
2943	zscale = 5,
2944	verbose=self.verbose,
2945	origin = (56, 0, 0))
2946
2947	#Work out the UTM coordinates for first point
2948	fid = NetCDFFile(self.test_MOST_file + '.sww')
2949
2950	#Check values
2951	stage_5 = fid.variables['stage'][:]
2952	xmomentum_5 = fid.variables['xmomentum'][:]
2953	ymomentum_5 = fid.variables['ymomentum'][:]
2954	elevation = fid.variables['elevation'][:]
2955
2956	assert allclose(stage_5[0,0], 5*first_value/100) #Meters
2957	assert allclose(stage_5[0,3], 5*fourth_value/100) #Meters
2958
2959	assert allclose(5*stage_1, stage_5)
2960
2961	# Momentum will also be changed due to new depth
2962
2963	depth_1 = stage_1-elevation
2964	depth_5 = stage_5-elevation
2965
2966
2967	for i in range(stage_1.shape[0]):
2968	for j in range(stage_1.shape[1]):
2969	if depth_1[i,j] > epsilon:
2970
2971	scale = depth_5[i,j]/depth_1[i,j]
2972	ref_xmomentum = xmomentum_1[i,j] * scale
2973	ref_ymomentum = ymomentum_1[i,j] * scale
2974
2975	#print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
2976
2977	assert allclose(xmomentum_5[i,j], ref_xmomentum)
2978	assert allclose(ymomentum_5[i,j], ref_ymomentum)
2979
2980
2981
2982	fid.close()
2983
2984
2985	#Cleanup
2986	import os
2987	os.remove(self.test_MOST_file + '.sww')
2988
2989
2990
2991	def test_ferret2sww_2(self):
2992	"""Test that georeferencing etc works when converting from
2993	ferret format (lat/lon) to sww format (UTM)
2994	"""
2995	from Scientific.IO.NetCDF import NetCDFFile
2996
2997	#The test file has
2998	# LON = 150.66667, 150.83334, 151, 151.16667
2999	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3000	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
3001	#
3002	# First value (index=0) in small_ha.nc is 0.3400644 cm,
3003	# Fourth value (index==3) is -6.50198 cm
3004
3005
3006	from anuga.coordinate_transforms.redfearn import redfearn
3007
3008	#fid = NetCDFFile('small_ha.nc')
3009	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3010
3011	#Pick a coordinate and a value
3012
3013	time_index = 1
3014	lat_index = 0
3015	lon_index = 2
3016
3017	test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
3018	test_time = fid.variables['TIME'][:][time_index]
3019	test_lat = fid.variables['LAT'][:][lat_index]
3020	test_lon = fid.variables['LON'][:][lon_index]
3021
3022	linear_point_index = lat_index*4 + lon_index
3023	fid.close()
3024
3025	#Call conversion (with zero origin)
3026	ferret2sww(self.test_MOST_file, verbose=self.verbose,
3027	origin = (56, 0, 0))
3028
3029
3030	#Work out the UTM coordinates for test point
3031	zone, e, n = redfearn(test_lat, test_lon)
3032
3033	#Read output file 'small.sww'
3034	fid = NetCDFFile(self.test_MOST_file + '.sww')
3035
3036	x = fid.variables['x'][:]
3037	y = fid.variables['y'][:]
3038
3039	#Check that test coordinate is correctly represented
3040	assert allclose(x[linear_point_index], e)
3041	assert allclose(y[linear_point_index], n)
3042
3043	#Check test value
3044	stage = fid.variables['stage'][:]
3045
3046	assert allclose(stage[time_index, linear_point_index], test_value/100)
3047
3048	fid.close()
3049
3050	#Cleanup
3051	import os
3052	os.remove(self.test_MOST_file + '.sww')
3053
3054
3055	def test_ferret2sww_lat_long(self):
3056	# Test that min lat long works
3057
3058	#The test file has
3059	# LON = 150.66667, 150.83334, 151, 151.16667
3060	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3061
3062	#Read
3063	from anuga.coordinate_transforms.redfearn import redfearn
3064	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3065	first_value = fid.variables['HA'][:][0,0,0]
3066	fourth_value = fid.variables['HA'][:][0,0,3]
3067	fid.close()
3068
3069
3070	#Call conversion (with zero origin)
3071	#ferret2sww('small', verbose=self.verbose,
3072	# origin = (56, 0, 0))
3073	ferret2sww(self.test_MOST_file, verbose=self.verbose,
3074	origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
3075
3076	#Work out the UTM coordinates for first point
3077	zone, e, n = redfearn(-34.5, 150.66667)
3078	#print zone, e, n
3079
3080	#Read output file 'small.sww'
3081	#fid = NetCDFFile('small.sww')
3082	fid = NetCDFFile(self.test_MOST_file + '.sww')
3083
3084	x = fid.variables['x'][:]
3085	y = fid.variables['y'][:]
3086	#Check that first coordinate is correctly represented
3087	assert 16 == len(x)
3088
3089	fid.close()
3090
3091	#Cleanup
3092	import os
3093	os.remove(self.test_MOST_file + '.sww')
3094
3095
3096	def test_ferret2sww_lat_longII(self):
3097	# Test that min lat long works
3098
3099	#The test file has
3100	# LON = 150.66667, 150.83334, 151, 151.16667
3101	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3102
3103	#Read
3104	from anuga.coordinate_transforms.redfearn import redfearn
3105	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3106	first_value = fid.variables['HA'][:][0,0,0]
3107	fourth_value = fid.variables['HA'][:][0,0,3]
3108	fid.close()
3109
3110
3111	#Call conversion (with zero origin)
3112	#ferret2sww('small', verbose=False,
3113	# origin = (56, 0, 0))
3114	ferret2sww(self.test_MOST_file, verbose=False,
3115	origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
3116
3117	#Work out the UTM coordinates for first point
3118	zone, e, n = redfearn(-34.5, 150.66667)
3119	#print zone, e, n
3120
3121	#Read output file 'small.sww'
3122	#fid = NetCDFFile('small.sww')
3123	fid = NetCDFFile(self.test_MOST_file + '.sww')
3124
3125	x = fid.variables['x'][:]
3126	y = fid.variables['y'][:]
3127	#Check that first coordinate is correctly represented
3128	assert 12 == len(x)
3129
3130	fid.close()
3131
3132	#Cleanup
3133	import os
3134	os.remove(self.test_MOST_file + '.sww')
3135
3136	def test_ferret2sww3(self):
3137	"""Elevation included
3138	"""
3139	from Scientific.IO.NetCDF import NetCDFFile
3140
3141	#The test file has
3142	# LON = 150.66667, 150.83334, 151, 151.16667
3143	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3144	# ELEVATION = [-1 -2 -3 -4
3145	# -5 -6 -7 -8
3146	# ...
3147	# ... -16]
3148	# where the top left corner is -1m,
3149	# and the ll corner is -13.0m
3150	#
3151	# First value (index=0) in small_ha.nc is 0.3400644 cm,
3152	# Fourth value (index==3) is -6.50198 cm
3153
3154	from anuga.coordinate_transforms.redfearn import redfearn
3155	import os
3156	fid1 = NetCDFFile('test_ha.nc','w')
3157	fid2 = NetCDFFile('test_ua.nc','w')
3158	fid3 = NetCDFFile('test_va.nc','w')
3159	fid4 = NetCDFFile('test_e.nc','w')
3160
3161	h1_list = [150.66667,150.83334,151.]
3162	h2_list = [-34.5,-34.33333]
3163
3164	long_name = 'LON'
3165	lat_name = 'LAT'
3166	time_name = 'TIME'
3167
3168	nx = 3
3169	ny = 2
3170
3171	for fid in [fid1,fid2,fid3]:
3172	fid.createDimension(long_name,nx)
3173	fid.createVariable(long_name,'d',(long_name,))
3174	fid.variables[long_name].point_spacing='uneven'
3175	fid.variables[long_name].units='degrees_east'
3176	fid.variables[long_name].assignValue(h1_list)
3177
3178	fid.createDimension(lat_name,ny)
3179	fid.createVariable(lat_name,'d',(lat_name,))
3180	fid.variables[lat_name].point_spacing='uneven'
3181	fid.variables[lat_name].units='degrees_north'
3182	fid.variables[lat_name].assignValue(h2_list)
3183
3184	fid.createDimension(time_name,2)
3185	fid.createVariable(time_name,'d',(time_name,))
3186	fid.variables[time_name].point_spacing='uneven'
3187	fid.variables[time_name].units='seconds'
3188	fid.variables[time_name].assignValue([0.,1.])
3189	if fid == fid3: break
3190
3191
3192	for fid in [fid4]:
3193	fid.createDimension(long_name,nx)
3194	fid.createVariable(long_name,'d',(long_name,))
3195	fid.variables[long_name].point_spacing='uneven'
3196	fid.variables[long_name].units='degrees_east'
3197	fid.variables[long_name].assignValue(h1_list)
3198
3199	fid.createDimension(lat_name,ny)
3200	fid.createVariable(lat_name,'d',(lat_name,))
3201	fid.variables[lat_name].point_spacing='uneven'
3202	fid.variables[lat_name].units='degrees_north'
3203	fid.variables[lat_name].assignValue(h2_list)
3204
3205	name = {}
3206	name[fid1]='HA'
3207	name[fid2]='UA'
3208	name[fid3]='VA'
3209	name[fid4]='ELEVATION'
3210
3211	units = {}
3212	units[fid1]='cm'
3213	units[fid2]='cm/s'
3214	units[fid3]='cm/s'
3215	units[fid4]='m'
3216
3217	values = {}
3218	values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
3219	values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
3220	values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
3221	values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
3222
3223	for fid in [fid1,fid2,fid3]:
3224	fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
3225	fid.variables[name[fid]].point_spacing='uneven'
3226	fid.variables[name[fid]].units=units[fid]
3227	fid.variables[name[fid]].assignValue(values[fid])
3228	fid.variables[name[fid]].missing_value = -99999999.
3229	if fid == fid3: break
3230
3231	for fid in [fid4]:
3232	fid.createVariable(name[fid],'d',(lat_name,long_name))
3233	fid.variables[name[fid]].point_spacing='uneven'
3234	fid.variables[name[fid]].units=units[fid]
3235	fid.variables[name[fid]].assignValue(values[fid])
3236	fid.variables[name[fid]].missing_value = -99999999.
3237
3238
3239	fid1.sync(); fid1.close()
3240	fid2.sync(); fid2.close()
3241	fid3.sync(); fid3.close()
3242	fid4.sync(); fid4.close()
3243
3244	fid1 = NetCDFFile('test_ha.nc','r')
3245	fid2 = NetCDFFile('test_e.nc','r')
3246	fid3 = NetCDFFile('test_va.nc','r')
3247
3248
3249	first_amp = fid1.variables['HA'][:][0,0,0]
3250	third_amp = fid1.variables['HA'][:][0,0,2]
3251	first_elevation = fid2.variables['ELEVATION'][0,0]
3252	third_elevation= fid2.variables['ELEVATION'][:][0,2]
3253	first_speed = fid3.variables['VA'][0,0,0]
3254	third_speed = fid3.variables['VA'][:][0,0,2]
3255
3256	fid1.close()
3257	fid2.close()
3258	fid3.close()
3259
3260	#Call conversion (with zero origin)
3261	ferret2sww('test', verbose=self.verbose,
3262	origin = (56, 0, 0), inverted_bathymetry=False)
3263
3264	os.remove('test_va.nc')
3265	os.remove('test_ua.nc')
3266	os.remove('test_ha.nc')
3267	os.remove('test_e.nc')
3268
3269	#Read output file 'test.sww'
3270	fid = NetCDFFile('test.sww')
3271
3272
3273	#Check first value
3274	elevation = fid.variables['elevation'][:]
3275	stage = fid.variables['stage'][:]
3276	xmomentum = fid.variables['xmomentum'][:]
3277	ymomentum = fid.variables['ymomentum'][:]
3278
3279	#print ymomentum
3280	first_height = first_amp/100 - first_elevation
3281	third_height = third_amp/100 - third_elevation
3282	first_momentum=first_speed*first_height/100
3283	third_momentum=third_speed*third_height/100
3284
3285	assert allclose(ymomentum[0][0],first_momentum) #Meters
3286	assert allclose(ymomentum[0][2],third_momentum) #Meters
3287
3288	fid.close()
3289
3290	#Cleanup
3291	os.remove('test.sww')
3292
3293
3294
3295	def test_ferret2sww4(self):
3296	"""Like previous but with augmented variable names as
3297	in files produced by ferret as opposed to MOST
3298	"""
3299	from Scientific.IO.NetCDF import NetCDFFile
3300
3301	#The test file has
3302	# LON = 150.66667, 150.83334, 151, 151.16667
3303	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3304	# ELEVATION = [-1 -2 -3 -4
3305	# -5 -6 -7 -8
3306	# ...
3307	# ... -16]
3308	# where the top left corner is -1m,
3309	# and the ll corner is -13.0m
3310	#
3311	# First value (index=0) in small_ha.nc is 0.3400644 cm,
3312	# Fourth value (index==3) is -6.50198 cm
3313
3314	from anuga.coordinate_transforms.redfearn import redfearn
3315	import os
3316	fid1 = NetCDFFile('test_ha.nc','w')
3317	fid2 = NetCDFFile('test_ua.nc','w')
3318	fid3 = NetCDFFile('test_va.nc','w')
3319	fid4 = NetCDFFile('test_e.nc','w')
3320
3321	h1_list = [150.66667,150.83334,151.]
3322	h2_list = [-34.5,-34.33333]
3323
3324	# long_name = 'LON961_1261'
3325	# lat_name = 'LAT481_841'
3326	# time_name = 'TIME1'
3327
3328	long_name = 'LON'
3329	lat_name = 'LAT'
3330	time_name = 'TIME'
3331
3332	nx = 3
3333	ny = 2
3334
3335	for fid in [fid1,fid2,fid3]:
3336	fid.createDimension(long_name,nx)
3337	fid.createVariable(long_name,'d',(long_name,))
3338	fid.variables[long_name].point_spacing='uneven'
3339	fid.variables[long_name].units='degrees_east'
3340	fid.variables[long_name].assignValue(h1_list)
3341
3342	fid.createDimension(lat_name,ny)
3343	fid.createVariable(lat_name,'d',(lat_name,))
3344	fid.variables[lat_name].point_spacing='uneven'
3345	fid.variables[lat_name].units='degrees_north'
3346	fid.variables[lat_name].assignValue(h2_list)
3347
3348	fid.createDimension(time_name,2)
3349	fid.createVariable(time_name,'d',(time_name,))
3350	fid.variables[time_name].point_spacing='uneven'
3351	fid.variables[time_name].units='seconds'
3352	fid.variables[time_name].assignValue([0.,1.])
3353	if fid == fid3: break
3354
3355
3356	for fid in [fid4]:
3357	fid.createDimension(long_name,nx)
3358	fid.createVariable(long_name,'d',(long_name,))
3359	fid.variables[long_name].point_spacing='uneven'
3360	fid.variables[long_name].units='degrees_east'
3361	fid.variables[long_name].assignValue(h1_list)
3362
3363	fid.createDimension(lat_name,ny)
3364	fid.createVariable(lat_name,'d',(lat_name,))
3365	fid.variables[lat_name].point_spacing='uneven'
3366	fid.variables[lat_name].units='degrees_north'
3367	fid.variables[lat_name].assignValue(h2_list)
3368
3369	name = {}
3370	name[fid1]='HA'
3371	name[fid2]='UA'
3372	name[fid3]='VA'
3373	name[fid4]='ELEVATION'
3374
3375	units = {}
3376	units[fid1]='cm'
3377	units[fid2]='cm/s'
3378	units[fid3]='cm/s'
3379	units[fid4]='m'
3380
3381	values = {}
3382	values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
3383	values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
3384	values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
3385	values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
3386
3387	for fid in [fid1,fid2,fid3]:
3388	fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
3389	fid.variables[name[fid]].point_spacing='uneven'
3390	fid.variables[name[fid]].units=units[fid]
3391	fid.variables[name[fid]].assignValue(values[fid])
3392	fid.variables[name[fid]].missing_value = -99999999.
3393	if fid == fid3: break
3394
3395	for fid in [fid4]:
3396	fid.createVariable(name[fid],'d',(lat_name,long_name))
3397	fid.variables[name[fid]].point_spacing='uneven'
3398	fid.variables[name[fid]].units=units[fid]
3399	fid.variables[name[fid]].assignValue(values[fid])
3400	fid.variables[name[fid]].missing_value = -99999999.
3401
3402
3403	fid1.sync(); fid1.close()
3404	fid2.sync(); fid2.close()
3405	fid3.sync(); fid3.close()
3406	fid4.sync(); fid4.close()
3407
3408	fid1 = NetCDFFile('test_ha.nc','r')
3409	fid2 = NetCDFFile('test_e.nc','r')
3410	fid3 = NetCDFFile('test_va.nc','r')
3411
3412
3413	first_amp = fid1.variables['HA'][:][0,0,0]
3414	third_amp = fid1.variables['HA'][:][0,0,2]
3415	first_elevation = fid2.variables['ELEVATION'][0,0]
3416	third_elevation= fid2.variables['ELEVATION'][:][0,2]
3417	first_speed = fid3.variables['VA'][0,0,0]
3418	third_speed = fid3.variables['VA'][:][0,0,2]
3419
3420	fid1.close()
3421	fid2.close()
3422	fid3.close()
3423
3424	#Call conversion (with zero origin)
3425	ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
3426	, inverted_bathymetry=False)
3427
3428	os.remove('test_va.nc')
3429	os.remove('test_ua.nc')
3430	os.remove('test_ha.nc')
3431	os.remove('test_e.nc')
3432
3433	#Read output file 'test.sww'
3434	fid = NetCDFFile('test.sww')
3435
3436
3437	#Check first value
3438	elevation = fid.variables['elevation'][:]
3439	stage = fid.variables['stage'][:]
3440	xmomentum = fid.variables['xmomentum'][:]
3441	ymomentum = fid.variables['ymomentum'][:]
3442
3443	#print ymomentum
3444	first_height = first_amp/100 - first_elevation
3445	third_height = third_amp/100 - third_elevation
3446	first_momentum=first_speed*first_height/100
3447	third_momentum=third_speed*third_height/100
3448
3449	assert allclose(ymomentum[0][0],first_momentum) #Meters
3450	assert allclose(ymomentum[0][2],third_momentum) #Meters
3451
3452	fid.close()
3453
3454	#Cleanup
3455	os.remove('test.sww')
3456
3457
3458
3459
3460	def test_ferret2sww_nz_origin(self):
3461	from Scientific.IO.NetCDF import NetCDFFile
3462	from anuga.coordinate_transforms.redfearn import redfearn
3463
3464	#Call conversion (with nonzero origin)
3465	ferret2sww(self.test_MOST_file, verbose=self.verbose,
3466	origin = (56, 100000, 200000))
3467
3468
3469	#Work out the UTM coordinates for first point
3470	zone, e, n = redfearn(-34.5, 150.66667)
3471
3472	#Read output file 'small.sww'
3473	#fid = NetCDFFile('small.sww', 'r')
3474	fid = NetCDFFile(self.test_MOST_file + '.sww')
3475
3476	x = fid.variables['x'][:]
3477	y = fid.variables['y'][:]
3478
3479	#Check that first coordinate is correctly represented
3480	assert allclose(x[0], e-100000)
3481	assert allclose(y[0], n-200000)
3482
3483	fid.close()
3484
3485	#Cleanup
3486	os.remove(self.test_MOST_file + '.sww')
3487
3488
3489	def test_ferret2sww_lat_longII(self):
3490	# Test that min lat long works
3491
3492	#The test file has
3493	# LON = 150.66667, 150.83334, 151, 151.16667
3494	# LAT = -34.5, -34.33333, -34.16667, -34 ;
3495
3496	#Read
3497	from anuga.coordinate_transforms.redfearn import redfearn
3498	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
3499	first_value = fid.variables['HA'][:][0,0,0]
3500	fourth_value = fid.variables['HA'][:][0,0,3]
3501	fid.close()
3502
3503
3504	#Call conversion (with zero origin)
3505	#ferret2sww('small', verbose=self.verbose,
3506	# origin = (56, 0, 0))
3507	try:
3508	ferret2sww(self.test_MOST_file, verbose=self.verbose,
3509	origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
3510	except AssertionError:
3511	pass
3512	else:
3513	self.failUnless(0 ==1, 'Bad input did not throw exception error!')
3514
3515	def test_sww_extent(self):
3516	"""Not a test, rather a look at the sww format
3517	"""
3518
3519	import time, os
3520	from Numeric import array, zeros, allclose, Float, concatenate
3521	from Scientific.IO.NetCDF import NetCDFFile
3522
3523	self.domain.set_name('datatest' + str(id(self)))
3524	self.domain.format = 'sww'
3525	self.domain.smooth = True
3526	self.domain.reduction = mean
3527	self.domain.set_datadir('.')
3528	#self.domain.tight_slope_limiters = 1
3529
3530
3531	sww = get_dataobject(self.domain)
3532	sww.store_connectivity()
3533	sww.store_timestep()
3534	self.domain.time = 2.
3535
3536	#Modify stage at second timestep
3537	stage = self.domain.quantities['stage'].vertex_values
3538	self.domain.set_quantity('stage', stage/2)
3539
3540	sww.store_timestep()
3541
3542	file_and_extension_name = self.domain.get_name() + ".sww"
3543	#print "file_and_extension_name",file_and_extension_name
3544	[xmin, xmax, ymin, ymax, stagemin, stagemax] = \
3545	extent_sww(file_and_extension_name )
3546
3547	assert allclose(xmin, 0.0)
3548	assert allclose(xmax, 1.0)
3549	assert allclose(ymin, 0.0)
3550	assert allclose(ymax, 1.0)
3551
3552	# FIXME (Ole): Revisit these numbers
3553	#assert allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
3554	#assert allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
3555
3556
3557	#Cleanup
3558	os.remove(sww.filename)
3559
3560
3561
3562	def test_sww2domain1(self):
3563	################################################
3564	#Create a test domain, and evolve and save it.
3565	################################################
3566	from mesh_factory import rectangular
3567	from Numeric import array
3568
3569	#Create basic mesh
3570
3571	yiel=0.01
3572	points, vertices, boundary = rectangular(10,10)
3573
3574	#Create shallow water domain
3575	domain = Domain(points, vertices, boundary)
3576	domain.geo_reference = Geo_reference(56,11,11)
3577	domain.smooth = False
3578	domain.store = True
3579	domain.set_name('bedslope')
3580	domain.default_order=2
3581	#Bed-slope and friction
3582	domain.set_quantity('elevation', lambda x,y: -x/3)
3583	domain.set_quantity('friction', 0.1)
3584	# Boundary conditions
3585	from math import sin, pi
3586	Br = Reflective_boundary(domain)
3587	Bt = Transmissive_boundary(domain)
3588	Bd = Dirichlet_boundary([0.2,0.,0.])
3589	Bw = Time_boundary(domain=domain,f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
3590
3591	#domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3592	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3593
3594	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
3595	#Initial condition
3596	h = 0.05
3597	elevation = domain.quantities['elevation'].vertex_values
3598	domain.set_quantity('stage', elevation + h)
3599
3600	domain.check_integrity()
3601	#Evolution
3602	#domain.tight_slope_limiters = 1
3603	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3604	#domain.write_time()
3605	pass
3606
3607
3608	##########################################
3609	#Import the example's file as a new domain
3610	##########################################
3611	from data_manager import sww2domain
3612	from Numeric import allclose
3613	import os
3614
3615	filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3616	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
3617	#points, vertices, boundary = rectangular(15,15)
3618	#domain2.boundary = boundary
3619	###################
3620	##NOW TEST IT!!!
3621	###################
3622
3623	os.remove(filename)
3624
3625	bits = ['vertex_coordinates']
3626	for quantity in ['elevation']+domain.quantities_to_be_stored:
3627	bits.append('get_quantity("%s").get_integral()' %quantity)
3628	bits.append('get_quantity("%s").get_values()' %quantity)
3629
3630	for bit in bits:
3631	#print 'testing that domain.'+bit+' has been restored'
3632	#print bit
3633	#print 'done'
3634	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3635
3636	######################################
3637	#Now evolve them both, just to be sure
3638	######################################x
3639	domain.time = 0.
3640	from time import sleep
3641
3642	final = .1
3643	domain.set_quantity('friction', 0.1)
3644	domain.store = False
3645	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3646
3647
3648	for t in domain.evolve(yieldstep = yiel, finaltime = final):
3649	#domain.write_time()
3650	pass
3651
3652	final = final - (domain2.starttime-domain.starttime)
3653	#BUT since domain1 gets time hacked back to 0:
3654	final = final + (domain2.starttime-domain.starttime)
3655
3656	domain2.smooth = False
3657	domain2.store = False
3658	domain2.default_order=2
3659	domain2.set_quantity('friction', 0.1)
3660	#Bed-slope and friction
3661	# Boundary conditions
3662	Bd2=Dirichlet_boundary([0.2,0.,0.])
3663	domain2.boundary = domain.boundary
3664	#print 'domain2.boundary'
3665	#print domain2.boundary
3666	domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3667	#domain2.set_boundary({'exterior': Bd})
3668
3669	domain2.check_integrity()
3670
3671	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3672	#domain2.write_time()
3673	pass
3674
3675	###################
3676	##NOW TEST IT!!!
3677	##################
3678
3679	bits = ['vertex_coordinates']
3680
3681	for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
3682	bits.append('get_quantity("%s").get_integral()' %quantity)
3683	bits.append('get_quantity("%s").get_values()' %quantity)
3684
3685	#print bits
3686	for bit in bits:
3687	#print bit
3688	#print eval('domain.'+bit)
3689	#print eval('domain2.'+bit)
3690
3691	#print eval('domain.'+bit+'-domain2.'+bit)
3692	msg = 'Values in the two domains are different for ' + bit
3693	assert allclose(eval('domain.'+bit),eval('domain2.'+bit),
3694	rtol=1.e-5, atol=3.e-8), msg
3695
3696
3697	def DISABLEDtest_sww2domain2(self):
3698	##################################################################
3699	#Same as previous test, but this checks how NaNs are handled.
3700	##################################################################
3701
3702
3703	from mesh_factory import rectangular
3704	from Numeric import array
3705
3706	#Create basic mesh
3707	points, vertices, boundary = rectangular(2,2)
3708
3709	#Create shallow water domain
3710	domain = Domain(points, vertices, boundary)
3711	domain.smooth = False
3712	domain.store = True
3713	domain.set_name('test_file')
3714	domain.set_datadir('.')
3715	domain.default_order=2
3716	#domain.quantities_to_be_stored=['stage']
3717
3718	domain.set_quantity('elevation', lambda x,y: -x/3)
3719	domain.set_quantity('friction', 0.1)
3720
3721	from math import sin, pi
3722	Br = Reflective_boundary(domain)
3723	Bt = Transmissive_boundary(domain)
3724	Bd = Dirichlet_boundary([0.2,0.,0.])
3725	Bw = Time_boundary(domain=domain,
3726	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
3727
3728	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3729
3730	h = 0.05
3731	elevation = domain.quantities['elevation'].vertex_values
3732	domain.set_quantity('stage', elevation + h)
3733
3734	domain.check_integrity()
3735
3736	for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
3737	pass
3738	#domain.write_time()
3739
3740
3741
3742	##################################
3743	#Import the file as a new domain
3744	##################################
3745	from data_manager import sww2domain
3746	from Numeric import allclose
3747	import os
3748
3749	filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3750
3751	#Fail because NaNs are present
3752	try:
3753	domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=self.verbose)
3754	except:
3755	#Now import it, filling NaNs to be 0
3756	filler = 0
3757	domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=self.verbose)
3758
3759	#Clean up
3760	os.remove(filename)
3761
3762
3763	bits = ['geo_reference.get_xllcorner()',
3764	'geo_reference.get_yllcorner()',
3765	'vertex_coordinates']
3766
3767	for quantity in ['elevation']+domain.quantities_to_be_stored:
3768	bits.append('get_quantity("%s").get_integral()' %quantity)
3769	bits.append('get_quantity("%s").get_values()' %quantity)
3770
3771	for bit in bits:
3772	# print 'testing that domain.'+bit+' has been restored'
3773	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3774
3775	#print filler
3776	#print max(max(domain2.get_quantity('xmomentum').get_values()))
3777
3778	assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
3779	assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
3780	assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
3781	assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
3782
3783
3784
3785	def test_weed(self):
3786	from data_manager import weed
3787
3788	coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
3789	volumes1 = [[0,1,2],[3,4,5]]
3790	boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
3791	coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
3792
3793	points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
3794
3795	assert len(points2)==len(coordinates2)
3796	for i in range(len(coordinates2)):
3797	coordinate = tuple(coordinates2[i])
3798	assert points2.has_key(coordinate)
3799	points2[coordinate]=i
3800
3801	for triangle in volumes1:
3802	for coordinate in triangle:
3803	assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
3804	assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
3805
3806
3807	#FIXME This fails - smooth makes the comparism too hard for allclose
3808	def ztest_sww2domain3(self):
3809	################################################
3810	#DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
3811	################################################
3812	from mesh_factory import rectangular
3813	from Numeric import array
3814	#Create basic mesh
3815
3816	yiel=0.01
3817	points, vertices, boundary = rectangular(10,10)
3818
3819	#Create shallow water domain
3820	domain = Domain(points, vertices, boundary)
3821	domain.geo_reference = Geo_reference(56,11,11)
3822	domain.smooth = True
3823	domain.store = True
3824	domain.set_name('bedslope')
3825	domain.default_order=2
3826	#Bed-slope and friction
3827	domain.set_quantity('elevation', lambda x,y: -x/3)
3828	domain.set_quantity('friction', 0.1)
3829	# Boundary conditions
3830	from math import sin, pi
3831	Br = Reflective_boundary(domain)
3832	Bt = Transmissive_boundary(domain)
3833	Bd = Dirichlet_boundary([0.2,0.,0.])
3834	Bw = Time_boundary(domain=domain,
3835	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
3836
3837	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3838
3839	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
3840	#Initial condition
3841	h = 0.05
3842	elevation = domain.quantities['elevation'].vertex_values
3843	domain.set_quantity('stage', elevation + h)
3844
3845
3846	domain.check_integrity()
3847	#Evolution
3848	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3849	# domain.write_time()
3850	pass
3851
3852
3853	##########################################
3854	#Import the example's file as a new domain
3855	##########################################
3856	from data_manager import sww2domain
3857	from Numeric import allclose
3858	import os
3859
3860	filename = domain.datadir + os.sep + domain.get_name() + '.sww'
3861	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
3862	#points, vertices, boundary = rectangular(15,15)
3863	#domain2.boundary = boundary
3864	###################
3865	##NOW TEST IT!!!
3866	###################
3867
3868	os.remove(domain.get_name() + '.sww')
3869
3870	#FIXME smooth domain so that they can be compared
3871
3872
3873	bits = []#'vertex_coordinates']
3874	for quantity in ['elevation']+domain.quantities_to_be_stored:
3875	bits.append('quantities["%s"].get_integral()'%quantity)
3876
3877
3878	for bit in bits:
3879	#print 'testing that domain.'+bit+' has been restored'
3880	#print bit
3881	#print 'done'
3882	#print ('domain.'+bit), eval('domain.'+bit)
3883	#print ('domain2.'+bit), eval('domain2.'+bit)
3884	assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
3885	pass
3886
3887	######################################
3888	#Now evolve them both, just to be sure
3889	######################################x
3890	domain.time = 0.
3891	from time import sleep
3892
3893	final = .5
3894	domain.set_quantity('friction', 0.1)
3895	domain.store = False
3896	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
3897
3898	for t in domain.evolve(yieldstep = yiel, finaltime = final):
3899	#domain.write_time()
3900	pass
3901
3902	domain2.smooth = True
3903	domain2.store = False
3904	domain2.default_order=2
3905	domain2.set_quantity('friction', 0.1)
3906	#Bed-slope and friction
3907	# Boundary conditions
3908	Bd2=Dirichlet_boundary([0.2,0.,0.])
3909	Br2 = Reflective_boundary(domain2)
3910	domain2.boundary = domain.boundary
3911	#print 'domain2.boundary'
3912	#print domain2.boundary
3913	domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
3914	#domain2.boundary = domain.boundary
3915	#domain2.set_boundary({'exterior': Bd})
3916
3917	domain2.check_integrity()
3918
3919	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3920	#domain2.write_time()
3921	pass
3922
3923	###################
3924	##NOW TEST IT!!!
3925	##################
3926
3927	print '><><><><>>'
3928	bits = [ 'vertex_coordinates']
3929
3930	for quantity in ['elevation','xmomentum','ymomentum']:
3931	#bits.append('quantities["%s"].get_integral()'%quantity)
3932	bits.append('get_quantity("%s").get_values()' %quantity)
3933
3934	for bit in bits:
3935	print bit
3936	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3937
3938
3939	def test_decimate_dem(self):
3940	"""Test decimation of dem file
3941	"""
3942
3943	import os
3944	from Numeric import ones, allclose, Float, arange
3945	from Scientific.IO.NetCDF import NetCDFFile
3946
3947	#Write test dem file
3948	root = 'decdemtest'
3949
3950	filename = root + '.dem'
3951	fid = NetCDFFile(filename, 'w')
3952
3953	fid.institution = 'Geoscience Australia'
3954	fid.description = 'NetCDF DEM format for compact and portable ' +\
3955	'storage of spatial point data'
3956
3957	nrows = 15
3958	ncols = 18
3959
3960	fid.ncols = ncols
3961	fid.nrows = nrows
3962	fid.xllcorner = 2000.5
3963	fid.yllcorner = 3000.5
3964	fid.cellsize = 25
3965	fid.NODATA_value = -9999
3966
3967	fid.zone = 56
3968	fid.false_easting = 0.0
3969	fid.false_northing = 0.0
3970	fid.projection = 'UTM'
3971	fid.datum = 'WGS84'
3972	fid.units = 'METERS'
3973
3974	fid.createDimension('number_of_points', nrows*ncols)
3975
3976	fid.createVariable('elevation', Float, ('number_of_points',))
3977
3978	elevation = fid.variables['elevation']
3979
3980	elevation[:] = (arange(nrows*ncols))
3981
3982	fid.close()
3983
3984	#generate the elevation values expected in the decimated file
3985	ref_elevation = [( 0+ 1+ 2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3986	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3987	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3988	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3989	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3990	( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3991	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3992	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3993	(144+145+146+162+163+164+180+181+182) / 9.0,
3994	(148+149+150+166+167+168+184+185+186) / 9.0,
3995	(152+153+154+170+171+172+188+189+190) / 9.0,
3996	(156+157+158+174+175+176+192+193+194) / 9.0,
3997	(216+217+218+234+235+236+252+253+254) / 9.0,
3998	(220+221+222+238+239+240+256+257+258) / 9.0,
3999	(224+225+226+242+243+244+260+261+262) / 9.0,
4000	(228+229+230+246+247+248+264+265+266) / 9.0]
4001
4002	#generate a stencil for computing the decimated values
4003	stencil = ones((3,3), Float) / 9.0
4004
4005	decimate_dem(root, stencil=stencil, cellsize_new=100)
4006
4007	#Open decimated NetCDF file
4008	fid = NetCDFFile(root + '_100.dem', 'r')
4009
4010	# Get decimated elevation
4011	elevation = fid.variables['elevation']
4012
4013	#Check values
4014	assert allclose(elevation, ref_elevation)
4015
4016	#Cleanup
4017	fid.close()
4018
4019	os.remove(root + '.dem')
4020	os.remove(root + '_100.dem')
4021
4022	def test_decimate_dem_NODATA(self):
4023	"""Test decimation of dem file that includes NODATA values
4024	"""
4025
4026	import os
4027	from Numeric import ones, allclose, Float, arange, reshape
4028	from Scientific.IO.NetCDF import NetCDFFile
4029
4030	#Write test dem file
4031	root = 'decdemtest'
4032
4033	filename = root + '.dem'
4034	fid = NetCDFFile(filename, 'w')
4035
4036	fid.institution = 'Geoscience Australia'
4037	fid.description = 'NetCDF DEM format for compact and portable ' +\
4038	'storage of spatial point data'
4039
4040	nrows = 15
4041	ncols = 18
4042	NODATA_value = -9999
4043
4044	fid.ncols = ncols
4045	fid.nrows = nrows
4046	fid.xllcorner = 2000.5
4047	fid.yllcorner = 3000.5
4048	fid.cellsize = 25
4049	fid.NODATA_value = NODATA_value
4050
4051	fid.zone = 56
4052	fid.false_easting = 0.0
4053	fid.false_northing = 0.0
4054	fid.projection = 'UTM'
4055	fid.datum = 'WGS84'
4056	fid.units = 'METERS'
4057
4058	fid.createDimension('number_of_points', nrows*ncols)
4059
4060	fid.createVariable('elevation', Float, ('number_of_points',))
4061
4062	elevation = fid.variables['elevation']
4063
4064	#generate initial elevation values
4065	elevation_tmp = (arange(nrows*ncols))
4066	#add some NODATA values
4067	elevation_tmp[0] = NODATA_value
4068	elevation_tmp[95] = NODATA_value
4069	elevation_tmp[188] = NODATA_value
4070	elevation_tmp[189] = NODATA_value
4071	elevation_tmp[190] = NODATA_value
4072	elevation_tmp[209] = NODATA_value
4073	elevation_tmp[252] = NODATA_value
4074
4075	elevation[:] = elevation_tmp
4076
4077	fid.close()
4078
4079	#generate the elevation values expected in the decimated file
4080	ref_elevation = [NODATA_value,
4081	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
4082	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
4083	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
4084	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
4085	NODATA_value,
4086	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
4087	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
4088	(144+145+146+162+163+164+180+181+182) / 9.0,
4089	(148+149+150+166+167+168+184+185+186) / 9.0,
4090	NODATA_value,
4091	(156+157+158+174+175+176+192+193+194) / 9.0,
4092	NODATA_value,
4093	(220+221+222+238+239+240+256+257+258) / 9.0,
4094	(224+225+226+242+243+244+260+261+262) / 9.0,
4095	(228+229+230+246+247+248+264+265+266) / 9.0]
4096
4097	#generate a stencil for computing the decimated values
4098	stencil = ones((3,3), Float) / 9.0
4099
4100	decimate_dem(root, stencil=stencil, cellsize_new=100)
4101
4102	#Open decimated NetCDF file
4103	fid = NetCDFFile(root + '_100.dem', 'r')
4104
4105	# Get decimated elevation
4106	elevation = fid.variables['elevation']
4107
4108	#Check values
4109	assert allclose(elevation, ref_elevation)
4110
4111	#Cleanup
4112	fid.close()
4113
4114	os.remove(root + '.dem')
4115	os.remove(root + '_100.dem')
4116
4117	def xxxtestz_sww2ers_real(self):
4118	"""Test that sww information can be converted correctly to asc/prj
4119	format readable by e.g. ArcView
4120	"""
4121
4122	import time, os
4123	from Numeric import array, zeros, allclose, Float, concatenate
4124	from Scientific.IO.NetCDF import NetCDFFile
4125
4126	# the memory optimised least squares
4127	# cellsize = 20, # this one seems to hang
4128	# cellsize = 200000, # Ran 1 test in 269.703s
4129	#Ran 1 test in 267.344s
4130	# cellsize = 20000, # Ran 1 test in 460.922s
4131	# cellsize = 2000 #Ran 1 test in 5340.250s
4132	# cellsize = 200 #this one seems to hang, building matirx A
4133
4134	# not optimised
4135	# seems to hang
4136	# cellsize = 2000 # Ran 1 test in 5334.563s
4137	#Export to ascii/prj files
4138	sww2dem('karratha_100m',
4139	quantity = 'depth',
4140	cellsize = 200000,
4141	verbose = True)
4142
4143	def test_read_asc(self):
4144	"""Test conversion from dem in ascii format to native NetCDF format
4145	"""
4146
4147	import time, os
4148	from Numeric import array, zeros, allclose, Float, concatenate
4149	from Scientific.IO.NetCDF import NetCDFFile
4150
4151	from data_manager import _read_asc
4152	#Write test asc file
4153	filename = tempfile.mktemp(".000")
4154	fid = open(filename, 'w')
4155	fid.write("""ncols 7
4156	nrows 4
4157	xllcorner 2000.5
4158	yllcorner 3000.5
4159	cellsize 25
4160	NODATA_value -9999
4161	97.921 99.285 125.588 180.830 258.645 342.872 415.836
4162	473.157 514.391 553.893 607.120 678.125 777.283 883.038
4163	984.494 1040.349 1008.161 900.738 730.882 581.430 514.980
4164	502.645 516.230 504.739 450.604 388.500 338.097 514.980
4165	""")
4166	fid.close()
4167	bath_metadata, grid = _read_asc(filename, verbose=self.verbose)
4168	self.failUnless(bath_metadata['xllcorner'] == 2000.5, 'Failed')
4169	self.failUnless(bath_metadata['yllcorner'] == 3000.5, 'Failed')
4170	self.failUnless(bath_metadata['cellsize'] == 25, 'Failed')
4171	self.failUnless(bath_metadata['NODATA_value'] == -9999, 'Failed')
4172	self.failUnless(grid[0][0] == 97.921, 'Failed')
4173	self.failUnless(grid[3][6] == 514.980, 'Failed')
4174
4175	os.remove(filename)
4176
4177	def test_asc_csiro2sww(self):
4178	import tempfile
4179
4180	bath_dir = tempfile.mkdtemp()
4181	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4182	#bath_dir = 'bath_data_manager_test'
4183	#print "os.getcwd( )",os.getcwd( )
4184	elevation_dir = tempfile.mkdtemp()
4185	#elevation_dir = 'elev_expanded'
4186	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4187	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4188
4189	fid = open(bath_dir_filename, 'w')
4190	fid.write(""" ncols 3
4191	nrows 2
4192	xllcorner 148.00000
4193	yllcorner -38.00000
4194	cellsize 0.25
4195	nodata_value -9999.0
4196	9000.000 -1000.000 3000.0
4197	-1000.000 9000.000 -1000.000
4198	""")
4199	fid.close()
4200
4201	fid = open(elevation_dir_filename1, 'w')
4202	fid.write(""" ncols 3
4203	nrows 2
4204	xllcorner 148.00000
4205	yllcorner -38.00000
4206	cellsize 0.25
4207	nodata_value -9999.0
4208	9000.000 0.000 3000.0
4209	0.000 9000.000 0.000
4210	""")
4211	fid.close()
4212
4213	fid = open(elevation_dir_filename2, 'w')
4214	fid.write(""" ncols 3
4215	nrows 2
4216	xllcorner 148.00000
4217	yllcorner -38.00000
4218	cellsize 0.25
4219	nodata_value -9999.0
4220	9000.000 4000.000 4000.0
4221	4000.000 9000.000 4000.000
4222	""")
4223	fid.close()
4224
4225	ucur_dir = tempfile.mkdtemp()
4226	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4227	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4228
4229	fid = open(ucur_dir_filename1, 'w')
4230	fid.write(""" ncols 3
4231	nrows 2
4232	xllcorner 148.00000
4233	yllcorner -38.00000
4234	cellsize 0.25
4235	nodata_value -9999.0
4236	90.000 60.000 30.0
4237	10.000 10.000 10.000
4238	""")
4239	fid.close()
4240	fid = open(ucur_dir_filename2, 'w')
4241	fid.write(""" ncols 3
4242	nrows 2
4243	xllcorner 148.00000
4244	yllcorner -38.00000
4245	cellsize 0.25
4246	nodata_value -9999.0
4247	90.000 60.000 30.0
4248	10.000 10.000 10.000
4249	""")
4250	fid.close()
4251
4252	vcur_dir = tempfile.mkdtemp()
4253	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4254	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4255
4256	fid = open(vcur_dir_filename1, 'w')
4257	fid.write(""" ncols 3
4258	nrows 2
4259	xllcorner 148.00000
4260	yllcorner -38.00000
4261	cellsize 0.25
4262	nodata_value -9999.0
4263	90.000 60.000 30.0
4264	10.000 10.000 10.000
4265	""")
4266	fid.close()
4267	fid = open(vcur_dir_filename2, 'w')
4268	fid.write(""" ncols 3
4269	nrows 2
4270	xllcorner 148.00000
4271	yllcorner -38.00000
4272	cellsize 0.25
4273	nodata_value -9999.0
4274	90.000 60.000 30.0
4275	10.000 10.000 10.000
4276	""")
4277	fid.close()
4278
4279	sww_file = 'a_test.sww'
4280	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
4281	verbose=self.verbose)
4282
4283	# check the sww file
4284
4285	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
4286	x = fid.variables['x'][:]
4287	y = fid.variables['y'][:]
4288	z = fid.variables['z'][:]
4289	stage = fid.variables['stage'][:]
4290	xmomentum = fid.variables['xmomentum'][:]
4291	geo_ref = Geo_reference(NetCDFObject=fid)
4292	#print "geo_ref",geo_ref
4293	x_ref = geo_ref.get_xllcorner()
4294	y_ref = geo_ref.get_yllcorner()
4295	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
4296	assert allclose(x_ref, 587798.418) # (-38, 148)
4297	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
4298
4299	#Zone: 55
4300	#Easting: 588095.674 Northing: 5821451.722
4301	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
4302	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
4303
4304	#Zone: 55
4305	#Easting: 632145.632 Northing: 5820863.269
4306	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
4307	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
4308
4309	#Zone: 55
4310	#Easting: 609748.788 Northing: 5793447.860
4311	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
4312	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
4313
4314	assert allclose(z[0],9000.0 )
4315	assert allclose(stage[0][1],0.0 )
4316
4317	#(4000+1000)*60
4318	assert allclose(xmomentum[1][1],300000.0 )
4319
4320
4321	fid.close()
4322
4323	#tidy up
4324	os.remove(bath_dir_filename)
4325	os.rmdir(bath_dir)
4326
4327	os.remove(elevation_dir_filename1)
4328	os.remove(elevation_dir_filename2)
4329	os.rmdir(elevation_dir)
4330
4331	os.remove(ucur_dir_filename1)
4332	os.remove(ucur_dir_filename2)
4333	os.rmdir(ucur_dir)
4334
4335	os.remove(vcur_dir_filename1)
4336	os.remove(vcur_dir_filename2)
4337	os.rmdir(vcur_dir)
4338
4339
4340	# remove sww file
4341	os.remove(sww_file)
4342
4343	def test_asc_csiro2sww2(self):
4344	import tempfile
4345
4346	bath_dir = tempfile.mkdtemp()
4347	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4348	#bath_dir = 'bath_data_manager_test'
4349	#print "os.getcwd( )",os.getcwd( )
4350	elevation_dir = tempfile.mkdtemp()
4351	#elevation_dir = 'elev_expanded'
4352	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4353	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4354
4355	fid = open(bath_dir_filename, 'w')
4356	fid.write(""" ncols 3
4357	nrows 2
4358	xllcorner 148.00000
4359	yllcorner -38.00000
4360	cellsize 0.25
4361	nodata_value -9999.0
4362	9000.000 -1000.000 3000.0
4363	-1000.000 9000.000 -1000.000
4364	""")
4365	fid.close()
4366
4367	fid = open(elevation_dir_filename1, 'w')
4368	fid.write(""" ncols 3
4369	nrows 2
4370	xllcorner 148.00000
4371	yllcorner -38.00000
4372	cellsize 0.25
4373	nodata_value -9999.0
4374	9000.000 0.000 3000.0
4375	0.000 -9999.000 -9999.000
4376	""")
4377	fid.close()
4378
4379	fid = open(elevation_dir_filename2, 'w')
4380	fid.write(""" ncols 3
4381	nrows 2
4382	xllcorner 148.00000
4383	yllcorner -38.00000
4384	cellsize 0.25
4385	nodata_value -9999.0
4386	9000.000 4000.000 4000.0
4387	4000.000 9000.000 4000.000
4388	""")
4389	fid.close()
4390
4391	ucur_dir = tempfile.mkdtemp()
4392	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4393	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4394
4395	fid = open(ucur_dir_filename1, 'w')
4396	fid.write(""" ncols 3
4397	nrows 2
4398	xllcorner 148.00000
4399	yllcorner -38.00000
4400	cellsize 0.25
4401	nodata_value -9999.0
4402	90.000 60.000 30.0
4403	10.000 10.000 10.000
4404	""")
4405	fid.close()
4406	fid = open(ucur_dir_filename2, 'w')
4407	fid.write(""" ncols 3
4408	nrows 2
4409	xllcorner 148.00000
4410	yllcorner -38.00000
4411	cellsize 0.25
4412	nodata_value -9999.0
4413	90.000 60.000 30.0
4414	10.000 10.000 10.000
4415	""")
4416	fid.close()
4417
4418	vcur_dir = tempfile.mkdtemp()
4419	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4420	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4421
4422	fid = open(vcur_dir_filename1, 'w')
4423	fid.write(""" ncols 3
4424	nrows 2
4425	xllcorner 148.00000
4426	yllcorner -38.00000
4427	cellsize 0.25
4428	nodata_value -9999.0
4429	90.000 60.000 30.0
4430	10.000 10.000 10.000
4431	""")
4432	fid.close()
4433	fid = open(vcur_dir_filename2, 'w')
4434	fid.write(""" ncols 3
4435	nrows 2
4436	xllcorner 148.00000
4437	yllcorner -38.00000
4438	cellsize 0.25
4439	nodata_value -9999.0
4440	90.000 60.000 30.0
4441	10.000 10.000 10.000
4442	""")
4443	fid.close()
4444
4445	try:
4446	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
4447	vcur_dir, sww_file,
4448	verbose=self.verbose)
4449	except:
4450	#tidy up
4451	os.remove(bath_dir_filename)
4452	os.rmdir(bath_dir)
4453
4454	os.remove(elevation_dir_filename1)
4455	os.remove(elevation_dir_filename2)
4456	os.rmdir(elevation_dir)
4457
4458	os.remove(ucur_dir_filename1)
4459	os.remove(ucur_dir_filename2)
4460	os.rmdir(ucur_dir)
4461
4462	os.remove(vcur_dir_filename1)
4463	os.remove(vcur_dir_filename2)
4464	os.rmdir(vcur_dir)
4465	else:
4466	#tidy up
4467	os.remove(bath_dir_filename)
4468	os.rmdir(bath_dir)
4469
4470	os.remove(elevation_dir_filename1)
4471	os.remove(elevation_dir_filename2)
4472	os.rmdir(elevation_dir)
4473	raise 'Should raise exception'
4474
4475	os.remove(ucur_dir_filename1)
4476	os.remove(ucur_dir_filename2)
4477	os.rmdir(ucur_dir)
4478
4479	os.remove(vcur_dir_filename1)
4480	os.remove(vcur_dir_filename2)
4481	os.rmdir(vcur_dir)
4482
4483
4484
4485	def test_asc_csiro2sww3(self):
4486	import tempfile
4487
4488	bath_dir = tempfile.mkdtemp()
4489	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4490	#bath_dir = 'bath_data_manager_test'
4491	#print "os.getcwd( )",os.getcwd( )
4492	elevation_dir = tempfile.mkdtemp()
4493	#elevation_dir = 'elev_expanded'
4494	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4495	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4496
4497	fid = open(bath_dir_filename, 'w')
4498	fid.write(""" ncols 3
4499	nrows 2
4500	xllcorner 148.00000
4501	yllcorner -38.00000
4502	cellsize 0.25
4503	nodata_value -9999.0
4504	9000.000 -1000.000 3000.0
4505	-1000.000 9000.000 -1000.000
4506	""")
4507	fid.close()
4508
4509	fid = open(elevation_dir_filename1, 'w')
4510	fid.write(""" ncols 3
4511	nrows 2
4512	xllcorner 148.00000
4513	yllcorner -38.00000
4514	cellsize 0.25
4515	nodata_value -9999.0
4516	9000.000 0.000 3000.0
4517	0.000 -9999.000 -9999.000
4518	""")
4519	fid.close()
4520
4521	fid = open(elevation_dir_filename2, 'w')
4522	fid.write(""" ncols 3
4523	nrows 2
4524	xllcorner 148.00000
4525	yllcorner -38.00000
4526	cellsize 0.25
4527	nodata_value -9999.0
4528	9000.000 4000.000 4000.0
4529	4000.000 9000.000 4000.000
4530	""")
4531	fid.close()
4532
4533	ucur_dir = tempfile.mkdtemp()
4534	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4535	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4536
4537	fid = open(ucur_dir_filename1, 'w')
4538	fid.write(""" ncols 3
4539	nrows 2
4540	xllcorner 148.00000
4541	yllcorner -38.00000
4542	cellsize 0.25
4543	nodata_value -9999.0
4544	90.000 60.000 30.0
4545	10.000 10.000 10.000
4546	""")
4547	fid.close()
4548	fid = open(ucur_dir_filename2, 'w')
4549	fid.write(""" ncols 3
4550	nrows 2
4551	xllcorner 148.00000
4552	yllcorner -38.00000
4553	cellsize 0.25
4554	nodata_value -9999.0
4555	90.000 60.000 30.0
4556	10.000 10.000 10.000
4557	""")
4558	fid.close()
4559
4560	vcur_dir = tempfile.mkdtemp()
4561	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4562	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4563
4564	fid = open(vcur_dir_filename1, 'w')
4565	fid.write(""" ncols 3
4566	nrows 2
4567	xllcorner 148.00000
4568	yllcorner -38.00000
4569	cellsize 0.25
4570	nodata_value -9999.0
4571	90.000 60.000 30.0
4572	10.000 10.000 10.000
4573	""")
4574	fid.close()
4575	fid = open(vcur_dir_filename2, 'w')
4576	fid.write(""" ncols 3
4577	nrows 2
4578	xllcorner 148.00000
4579	yllcorner -38.00000
4580	cellsize 0.25
4581	nodata_value -9999.0
4582	90.000 60.000 30.0
4583	10.000 10.000 10.000
4584	""")
4585	fid.close()
4586
4587	sww_file = 'a_test.sww'
4588	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
4589	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
4590	mean_stage = 100,
4591	verbose=self.verbose)
4592
4593	# check the sww file
4594
4595	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
4596	x = fid.variables['x'][:]
4597	y = fid.variables['y'][:]
4598	z = fid.variables['z'][:]
4599	stage = fid.variables['stage'][:]
4600	xmomentum = fid.variables['xmomentum'][:]
4601	geo_ref = Geo_reference(NetCDFObject=fid)
4602	#print "geo_ref",geo_ref
4603	x_ref = geo_ref.get_xllcorner()
4604	y_ref = geo_ref.get_yllcorner()
4605	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
4606	assert allclose(x_ref, 587798.418) # (-38, 148)
4607	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
4608
4609	#Zone: 55
4610	#Easting: 588095.674 Northing: 5821451.722
4611	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
4612	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
4613
4614	#Zone: 55
4615	#Easting: 632145.632 Northing: 5820863.269
4616	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
4617	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
4618
4619	#Zone: 55
4620	#Easting: 609748.788 Northing: 5793447.860
4621	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
4622	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
4623
4624	assert allclose(z[0],9000.0 )
4625	assert allclose(stage[0][4],100.0 )
4626	assert allclose(stage[0][5],100.0 )
4627
4628	#(100.0 - 9000)*10
4629	assert allclose(xmomentum[0][4], -89000.0 )
4630
4631	#(100.0 - -1000.000)*10
4632	assert allclose(xmomentum[0][5], 11000.0 )
4633
4634	fid.close()
4635
4636	#tidy up
4637	os.remove(bath_dir_filename)
4638	os.rmdir(bath_dir)
4639
4640	os.remove(elevation_dir_filename1)
4641	os.remove(elevation_dir_filename2)
4642	os.rmdir(elevation_dir)
4643
4644	os.remove(ucur_dir_filename1)
4645	os.remove(ucur_dir_filename2)
4646	os.rmdir(ucur_dir)
4647
4648	os.remove(vcur_dir_filename1)
4649	os.remove(vcur_dir_filename2)
4650	os.rmdir(vcur_dir)
4651
4652	# remove sww file
4653	os.remove(sww_file)
4654
4655
4656	def test_asc_csiro2sww4(self):
4657	"""
4658	Test specifying the extent
4659	"""
4660
4661	import tempfile
4662
4663	bath_dir = tempfile.mkdtemp()
4664	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
4665	#bath_dir = 'bath_data_manager_test'
4666	#print "os.getcwd( )",os.getcwd( )
4667	elevation_dir = tempfile.mkdtemp()
4668	#elevation_dir = 'elev_expanded'
4669	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
4670	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
4671
4672	fid = open(bath_dir_filename, 'w')
4673	fid.write(""" ncols 4
4674	nrows 4
4675	xllcorner 148.00000
4676	yllcorner -38.00000
4677	cellsize 0.25
4678	nodata_value -9999.0
4679	-9000.000 -1000.000 -3000.0 -2000.000
4680	-1000.000 9000.000 -1000.000 -3000.000
4681	-4000.000 6000.000 2000.000 -5000.000
4682	-9000.000 -1000.000 -3000.0 -2000.000
4683	""")
4684	fid.close()
4685
4686	fid = open(elevation_dir_filename1, 'w')
4687	fid.write(""" ncols 4
4688	nrows 4
4689	xllcorner 148.00000
4690	yllcorner -38.00000
4691	cellsize 0.25
4692	nodata_value -9999.0
4693	-900.000 -100.000 -300.0 -200.000
4694	-100.000 900.000 -100.000 -300.000
4695	-400.000 600.000 200.000 -500.000
4696	-900.000 -100.000 -300.0 -200.000
4697	""")
4698	fid.close()
4699
4700	fid = open(elevation_dir_filename2, 'w')
4701	fid.write(""" ncols 4
4702	nrows 4
4703	xllcorner 148.00000
4704	yllcorner -38.00000
4705	cellsize 0.25
4706	nodata_value -9999.0
4707	-990.000 -110.000 -330.0 -220.000
4708	-110.000 990.000 -110.000 -330.000
4709	-440.000 660.000 220.000 -550.000
4710	-990.000 -110.000 -330.0 -220.000
4711	""")
4712	fid.close()
4713
4714	ucur_dir = tempfile.mkdtemp()
4715	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
4716	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
4717
4718	fid = open(ucur_dir_filename1, 'w')
4719	fid.write(""" ncols 4
4720	nrows 4
4721	xllcorner 148.00000
4722	yllcorner -38.00000
4723	cellsize 0.25
4724	nodata_value -9999.0
4725	-90.000 -10.000 -30.0 -20.000
4726	-10.000 90.000 -10.000 -30.000
4727	-40.000 60.000 20.000 -50.000
4728	-90.000 -10.000 -30.0 -20.000
4729	""")
4730	fid.close()
4731	fid = open(ucur_dir_filename2, 'w')
4732	fid.write(""" ncols 4
4733	nrows 4
4734	xllcorner 148.00000
4735	yllcorner -38.00000
4736	cellsize 0.25
4737	nodata_value -9999.0
4738	-90.000 -10.000 -30.0 -20.000
4739	-10.000 99.000 -11.000 -30.000
4740	-40.000 66.000 22.000 -50.000
4741	-90.000 -10.000 -30.0 -20.000
4742	""")
4743	fid.close()
4744
4745	vcur_dir = tempfile.mkdtemp()
4746	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
4747	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
4748
4749	fid = open(vcur_dir_filename1, 'w')
4750	fid.write(""" ncols 4
4751	nrows 4
4752	xllcorner 148.00000
4753	yllcorner -38.00000
4754	cellsize 0.25
4755	nodata_value -9999.0
4756	-90.000 -10.000 -30.0 -20.000
4757	-10.000 80.000 -20.000 -30.000
4758	-40.000 50.000 10.000 -50.000
4759	-90.000 -10.000 -30.0 -20.000
4760	""")
4761	fid.close()
4762	fid = open(vcur_dir_filename2, 'w')
4763	fid.write(""" ncols 4
4764	nrows 4
4765	xllcorner 148.00000
4766	yllcorner -38.00000
4767	cellsize 0.25
4768	nodata_value -9999.0
4769	-90.000 -10.000 -30.0 -20.000
4770	-10.000 88.000 -22.000 -30.000
4771	-40.000 55.000 11.000 -50.000
4772	-90.000 -10.000 -30.0 -20.000
4773	""")
4774	fid.close()
4775
4776	sww_file = tempfile.mktemp(".sww")
4777	#sww_file = 'a_test.sww'
4778	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
4779	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
4780	mean_stage = 100,
4781	minlat = -37.6, maxlat = -37.6,
4782	minlon = 148.3, maxlon = 148.3,
4783	verbose=self.verbose
4784	#,verbose = True
4785	)
4786
4787	# check the sww file
4788
4789	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
4790	x = fid.variables['x'][:]
4791	y = fid.variables['y'][:]
4792	z = fid.variables['z'][:]
4793	stage = fid.variables['stage'][:]
4794	xmomentum = fid.variables['xmomentum'][:]
4795	ymomentum = fid.variables['ymomentum'][:]
4796	geo_ref = Geo_reference(NetCDFObject=fid)
4797	#print "geo_ref",geo_ref
4798	x_ref = geo_ref.get_xllcorner()
4799	y_ref = geo_ref.get_yllcorner()
4800	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
4801
4802	assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
4803	assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
4804	assert allclose(y_ref, 5820863.269 )# (-37.45, 148.5)
4805
4806	#Easting: 632145.632 Northing: 5820863.269
4807	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
4808
4809	#print "x",x
4810	#print "y",y
4811	self.failUnless(len(x) == 4,'failed') # 2*2
4812	self.failUnless(len(x) == 4,'failed') # 2*2
4813
4814	#Zone: 55
4815	#Easting: 632145.632 Northing: 5820863.269
4816	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
4817	# magic number - y is close enough for me.
4818	assert allclose(x[3], 632145.63 - x_ref)
4819	assert allclose(y[3], 5820863.269 - y_ref + 5.22155314684e-005)
4820
4821	assert allclose(z[0],9000.0 ) #z is elevation info
4822	#print "z",z
4823	# 2 time steps, 4 points
4824	self.failUnless(xmomentum.shape == (2,4), 'failed')
4825	self.failUnless(ymomentum.shape == (2,4), 'failed')
4826
4827	#(100.0 - -1000.000)*10
4828	#assert allclose(xmomentum[0][5], 11000.0 )
4829
4830	fid.close()
4831
4832	# is the sww file readable?
4833	#Lets see if we can convert it to a dem!
4834	# if you uncomment, remember to delete the file
4835	#print "sww_file",sww_file
4836	#dem_file = tempfile.mktemp(".dem")
4837	domain = sww2domain(sww_file) ###, dem_file)
4838	domain.check_integrity()
4839
4840	#tidy up
4841	os.remove(bath_dir_filename)
4842	os.rmdir(bath_dir)
4843
4844	os.remove(elevation_dir_filename1)
4845	os.remove(elevation_dir_filename2)
4846	os.rmdir(elevation_dir)
4847
4848	os.remove(ucur_dir_filename1)
4849	os.remove(ucur_dir_filename2)
4850	os.rmdir(ucur_dir)
4851
4852	os.remove(vcur_dir_filename1)
4853	os.remove(vcur_dir_filename2)
4854	os.rmdir(vcur_dir)
4855
4856
4857
4858
4859	# remove sww file
4860	os.remove(sww_file)
4861
4862
4863	def test_get_min_max_indexes(self):
4864	latitudes = [3,2,1,0]
4865	longitudes = [0,10,20,30]
4866
4867	# k - lat
4868	# l - lon
4869	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4870	latitudes,longitudes,
4871	-10,4,-10,31)
4872
4873	#print "kmin",kmin;print "kmax",kmax
4874	#print "lmin",lmin;print "lmax",lmax
4875	latitudes_new = latitudes[kmin:kmax]
4876	longitudes_news = longitudes[lmin:lmax]
4877	#print "latitudes_new", latitudes_new
4878	#print "longitudes_news",longitudes_news
4879	self.failUnless(latitudes == latitudes_new and \
4880	longitudes == longitudes_news,
4881	'failed')
4882
4883	## 2nd test
4884	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4885	latitudes,longitudes,
4886	0.5,2.5,5,25)
4887	#print "kmin",kmin;print "kmax",kmax
4888	#print "lmin",lmin;print "lmax",lmax
4889	latitudes_new = latitudes[kmin:kmax]
4890	longitudes_news = longitudes[lmin:lmax]
4891	#print "latitudes_new", latitudes_new
4892	#print "longitudes_news",longitudes_news
4893
4894	self.failUnless(latitudes == latitudes_new and \
4895	longitudes == longitudes_news,
4896	'failed')
4897
4898	## 3rd test
4899	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(\
4900	latitudes,
4901	longitudes,
4902	1.1,1.9,12,17)
4903	#print "kmin",kmin;print "kmax",kmax
4904	#print "lmin",lmin;print "lmax",lmax
4905	latitudes_new = latitudes[kmin:kmax]
4906	longitudes_news = longitudes[lmin:lmax]
4907	#print "latitudes_new", latitudes_new
4908	#print "longitudes_news",longitudes_news
4909
4910	self.failUnless(latitudes_new == [2, 1] and \
4911	longitudes_news == [10, 20],
4912	'failed')
4913
4914
4915	## 4th test
4916	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4917	latitudes,longitudes,
4918	-0.1,1.9,-2,17)
4919	#print "kmin",kmin;print "kmax",kmax
4920	#print "lmin",lmin;print "lmax",lmax
4921	latitudes_new = latitudes[kmin:kmax]
4922	longitudes_news = longitudes[lmin:lmax]
4923	#print "latitudes_new", latitudes_new
4924	#print "longitudes_news",longitudes_news
4925
4926	self.failUnless(latitudes_new == [2, 1, 0] and \
4927	longitudes_news == [0, 10, 20],
4928	'failed')
4929	## 5th test
4930	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4931	latitudes,longitudes,
4932	0.1,1.9,2,17)
4933	#print "kmin",kmin;print "kmax",kmax
4934	#print "lmin",lmin;print "lmax",lmax
4935	latitudes_new = latitudes[kmin:kmax]
4936	longitudes_news = longitudes[lmin:lmax]
4937	#print "latitudes_new", latitudes_new
4938	#print "longitudes_news",longitudes_news
4939
4940	self.failUnless(latitudes_new == [2, 1, 0] and \
4941	longitudes_news == [0, 10, 20],
4942	'failed')
4943
4944	## 6th test
4945
4946	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4947	latitudes,longitudes,
4948	1.5,4,18,32)
4949	#print "kmin",kmin;print "kmax",kmax
4950	#print "lmin",lmin;print "lmax",lmax
4951	latitudes_new = latitudes[kmin:kmax]
4952	longitudes_news = longitudes[lmin:lmax]
4953	#print "latitudes_new", latitudes_new
4954	#print "longitudes_news",longitudes_news
4955
4956	self.failUnless(latitudes_new == [3, 2, 1] and \
4957	longitudes_news == [10, 20, 30],
4958	'failed')
4959
4960
4961	## 7th test
4962	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4963	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4964	latitudes,longitudes,
4965	1.5,1.5,15,15)
4966	#print "kmin",kmin;print "kmax",kmax
4967	#print "lmin",lmin;print "lmax",lmax
4968	latitudes_new = latitudes[kmin:kmax]
4969	longitudes_news = longitudes[lmin:lmax]
4970	m2d = m2d[kmin:kmax,lmin:lmax]
4971	#print "m2d", m2d
4972	#print "latitudes_new", latitudes_new
4973	#print "longitudes_news",longitudes_news
4974
4975	self.failUnless(latitudes_new == [2, 1] and \
4976	longitudes_news == [10, 20],
4977	'failed')
4978
4979	self.failUnless(m2d == [[5,6],[9,10]],
4980	'failed')
4981
4982	def test_get_min_max_indexes_lat_ascending(self):
4983	latitudes = [0,1,2,3]
4984	longitudes = [0,10,20,30]
4985
4986	# k - lat
4987	# l - lon
4988	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4989	latitudes,longitudes,
4990	-10,4,-10,31)
4991
4992	#print "kmin",kmin;print "kmax",kmax
4993	#print "lmin",lmin;print "lmax",lmax
4994	latitudes_new = latitudes[kmin:kmax]
4995	longitudes_news = longitudes[lmin:lmax]
4996	#print "latitudes_new", latitudes_new
4997	#print "longitudes_news",longitudes_news
4998	self.failUnless(latitudes == latitudes_new and \
4999	longitudes == longitudes_news,
5000	'failed')
5001
5002	## 3rd test
5003	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(\
5004	latitudes,
5005	longitudes,
5006	1.1,1.9,12,17)
5007	#print "kmin",kmin;print "kmax",kmax
5008	#print "lmin",lmin;print "lmax",lmax
5009	latitudes_new = latitudes[kmin:kmax]
5010	longitudes_news = longitudes[lmin:lmax]
5011	#print "latitudes_new", latitudes_new
5012	#print "longitudes_news",longitudes_news
5013
5014	self.failUnless(latitudes_new == [1, 2] and \
5015	longitudes_news == [10, 20],
5016	'failed')
5017
5018	def test_get_min_max_indexes2(self):
5019	latitudes = [-30,-35,-40,-45]
5020	longitudes = [148,149,150,151]
5021
5022	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
5023
5024	# k - lat
5025	# l - lon
5026	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5027	latitudes,longitudes,
5028	-37,-27,147,149.5)
5029
5030	#print "kmin",kmin;print "kmax",kmax
5031	#print "lmin",lmin;print "lmax",lmax
5032	#print "m2d", m2d
5033	#print "latitudes", latitudes
5034	#print "longitudes",longitudes
5035	#print "latitudes[kmax]", latitudes[kmax]
5036	latitudes_new = latitudes[kmin:kmax]
5037	longitudes_new = longitudes[lmin:lmax]
5038	m2d = m2d[kmin:kmax,lmin:lmax]
5039	#print "m2d", m2d
5040	#print "latitudes_new", latitudes_new
5041	#print "longitudes_new",longitudes_new
5042
5043	self.failUnless(latitudes_new == [-30, -35, -40] and \
5044	longitudes_new == [148, 149,150],
5045	'failed')
5046	self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
5047	'failed')
5048
5049	def test_get_min_max_indexes3(self):
5050	latitudes = [-30,-35,-40,-45,-50,-55,-60]
5051	longitudes = [148,149,150,151]
5052
5053	# k - lat
5054	# l - lon
5055	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5056	latitudes,longitudes,
5057	-43,-37,148.5,149.5)
5058
5059
5060	#print "kmin",kmin;print "kmax",kmax
5061	#print "lmin",lmin;print "lmax",lmax
5062	#print "latitudes", latitudes
5063	#print "longitudes",longitudes
5064	latitudes_new = latitudes[kmin:kmax]
5065	longitudes_news = longitudes[lmin:lmax]
5066	#print "latitudes_new", latitudes_new
5067	#print "longitudes_news",longitudes_news
5068
5069	self.failUnless(latitudes_new == [-35, -40, -45] and \
5070	longitudes_news == [148, 149,150],
5071	'failed')
5072
5073	def test_get_min_max_indexes4(self):
5074	latitudes = [-30,-35,-40,-45,-50,-55,-60]
5075	longitudes = [148,149,150,151]
5076
5077	# k - lat
5078	# l - lon
5079	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
5080	latitudes,longitudes)
5081
5082	#print "kmin",kmin;print "kmax",kmax
5083	#print "lmin",lmin;print "lmax",lmax
5084	#print "latitudes", latitudes
5085	#print "longitudes",longitudes
5086	latitudes_new = latitudes[kmin:kmax]
5087	longitudes_news = longitudes[lmin:lmax]
5088	#print "latitudes_new", latitudes_new
5089	#print "longitudes_news",longitudes_news
5090
5091	self.failUnless(latitudes_new == latitudes and \
5092	longitudes_news == longitudes,
5093	'failed')
5094
5095	def test_tsh2sww(self):
5096	import os
5097	import tempfile
5098
5099	tsh_file = tempfile.mktemp(".tsh")
5100	file = open(tsh_file,"w")
5101	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
5102	0 0.0 0.0 0.0 0.0 0.01 \n \
5103	1 1.0 0.0 10.0 10.0 0.02 \n \
5104	2 0.0 1.0 0.0 10.0 0.03 \n \
5105	3 0.5 0.25 8.0 12.0 0.04 \n \
5106	# Vert att title \n \
5107	elevation \n \
5108	stage \n \
5109	friction \n \
5110	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
5111	0 0 3 2 -1 -1 1 dsg\n\
5112	1 0 1 3 -1 0 -1 ole nielsen\n\
5113	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
5114	0 1 0 2 \n\
5115	1 0 2 3 \n\
5116	2 2 3 \n\
5117	3 3 1 1 \n\
5118	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
5119	0 216.0 -86.0 \n \
5120	1 160.0 -167.0 \n \
5121	2 114.0 -91.0 \n \
5122	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
5123	0 0 1 0 \n \
5124	1 1 2 0 \n \
5125	2 2 0 0 \n \
5126	0 # <x> <y> ...Mesh Holes... \n \
5127	0 # <x> <y> <attribute>...Mesh Regions... \n \
5128	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
5129	#Geo reference \n \
5130	56 \n \
5131	140 \n \
5132	120 \n")
5133	file.close()
5134
5135	#sww_file = tempfile.mktemp(".sww")
5136	#print "sww_file",sww_file
5137	#print "sww_file",tsh_file
5138	tsh2sww(tsh_file,
5139	verbose=self.verbose)
5140
5141	os.remove(tsh_file)
5142	os.remove(tsh_file[:-4] + '.sww')
5143
5144
5145
5146
5147	########## testing nbed class ##################
5148	def test_exposure_csv_loading(self):
5149	file_name = tempfile.mktemp(".csv")
5150	file = open(file_name,"w")
5151	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
5152	115.0, -21.0, splat, 0.0\n\
5153	114.0, -21.7, pow, 10.0\n\
5154	114.5, -21.4, bang, 40.0\n")
5155	file.close()
5156	exposure = Exposure_csv(file_name, title_check_list = ['speed','sound'])
5157	exposure.get_column("sound")
5158
5159	self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
5160	'FAILED!')
5161	self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
5162	'FAILED!')
5163
5164	os.remove(file_name)
5165
5166	def test_exposure_csv_loadingII(self):
5167
5168
5169	file_name = tempfile.mktemp(".txt")
5170	file = open(file_name,"w")
5171	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
5172	115.0, -21.0, splat, 0.0\n\
5173	114.0, -21.7, pow, 10.0\n\
5174	114.5, -21.4, bang, 40.0\n")
5175	file.close()
5176	exposure = Exposure_csv(file_name)
5177	exposure.get_column("sound")
5178
5179	self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
5180	'FAILED!')
5181	self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
5182	'FAILED!')
5183
5184	os.remove(file_name)
5185
5186	def test_exposure_csv_loading_title_check_list(self):
5187
5188	# I can't get cvs.reader to close the exposure file
5189	# The hacks below are to get around this.
5190	if sys.platform == 'win32':
5191	file_name = tempfile.gettempdir() + \
5192	"test_exposure_csv_loading_title_check_list.csv"
5193	else:
5194	file_name = tempfile.mktemp(".csv")
5195	file = open(file_name,"w")
5196	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
5197	115.0, -21.0, splat, 0.0\n\
5198	114.0, -21.7, pow, 10.0\n\
5199	114.5, -21.4, bang, 40.0\n")
5200	file.close()
5201	try:
5202	exposure = Exposure_csv(file_name, title_check_list = ['SOUND'])
5203	except IOError:
5204	pass
5205	else:
5206	self.failUnless(0 ==1, 'Assertion not thrown error!')
5207
5208	if not sys.platform == 'win32':
5209	os.remove(file_name)
5210
5211	def test_exposure_csv_cmp(self):
5212	file_name = tempfile.mktemp(".csv")
5213	file = open(file_name,"w")
5214	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
5215	115.0, -21.0, splat, 0.0\n\
5216	114.0, -21.7, pow, 10.0\n\
5217	114.5, -21.4, bang, 40.0\n")
5218	file.close()
5219
5220	e1 = Exposure_csv(file_name)
5221	e2 = Exposure_csv(file_name)
5222	os.remove(file_name)
5223
5224	self.failUnless(cmp(e1,e2)==0,
5225	'FAILED!')
5226
5227	self.failUnless(cmp(e1,"hey")==1,
5228	'FAILED!')
5229
5230	file_name = tempfile.mktemp(".csv")
5231	file = open(file_name,"w")
5232	# Note, this has less spaces in the title,
5233	# the instances will be the same.
5234	file.write("LATITUDE,LONGITUDE ,sound, speed \n\
5235	115.0, -21.0, splat, 0.0\n\
5236	114.0, -21.7, pow, 10.0\n\
5237	114.5, -21.4, bang, 40.0\n")
5238	file.close()
5239	e3 = Exposure_csv(file_name)
5240	os.remove(file_name)
5241
5242	self.failUnless(cmp(e3,e2)==0,
5243	'FAILED!')
5244
5245	file_name = tempfile.mktemp(".csv")
5246	file = open(file_name,"w")
5247	# Note, 40 changed to 44 .
5248	file.write("LATITUDE,LONGITUDE ,sound, speed \n\
5249	115.0, -21.0, splat, 0.0\n\
5250	114.0, -21.7, pow, 10.0\n\
5251	114.5, -21.4, bang, 44.0\n")
5252	file.close()
5253	e4 = Exposure_csv(file_name)
5254	os.remove(file_name)
5255	#print "e4",e4._attribute_dic
5256	#print "e2",e2._attribute_dic
5257	self.failUnless(cmp(e4,e2)<>0,
5258	'FAILED!')
5259
5260	file_name = tempfile.mktemp(".csv")
5261	file = open(file_name,"w")
5262	# Note, the first two columns are swapped.
5263	file.write("LONGITUDE,LATITUDE ,sound, speed \n\
5264	-21.0,115.0, splat, 0.0\n\
5265	-21.7,114.0, pow, 10.0\n\
5266	-21.4,114.5, bang, 40.0\n")
5267	file.close()
5268	e5 = Exposure_csv(file_name)
5269	os.remove(file_name)
5270
5271	self.failUnless(cmp(e3,e5)<>0,
5272	'FAILED!')
5273
5274	def test_exposure_csv_saving(self):
5275
5276
5277	file_name = tempfile.mktemp(".csv")
5278	file = open(file_name,"w")
5279	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
5280	115.0, -21.0, splat, 0.0\n\
5281	114.0, -21.7, pow, 10.0\n\
5282	114.5, -21.4, bang, 40.0\n")
5283	file.close()
5284	e1 = Exposure_csv(file_name)
5285
5286	file_name2 = tempfile.mktemp(".csv")
5287	e1.save(file_name = file_name2)
5288	e2 = Exposure_csv(file_name2)
5289
5290	self.failUnless(cmp(e1,e2)==0,
5291	'FAILED!')
5292	os.remove(file_name)
5293	os.remove(file_name2)
5294
5295	def test_exposure_csv_get_location(self):
5296	file_name = tempfile.mktemp(".csv")
5297	file = open(file_name,"w")
5298	file.write("LONGITUDE , LATITUDE, sound , speed \n\
5299	150.916666667, -34.5, splat, 0.0\n\
5300	150.0, -34.0, pow, 10.0\n")
5301	file.close()
5302	e1 = Exposure_csv(file_name)
5303
5304	gsd = e1.get_location()
5305
5306	points = gsd.get_data_points(absolute=True)
5307
5308	assert allclose(points[0][0], 308728.009)
5309	assert allclose(points[0][1], 6180432.601)
5310	assert allclose(points[1][0], 222908.705)
5311	assert allclose(points[1][1], 6233785.284)
5312	self.failUnless(gsd.get_geo_reference().get_zone() == 56,
5313	'Bad zone error!')
5314
5315	os.remove(file_name)
5316
5317	def test_exposure_csv_set_column_get_column(self):
5318	file_name = tempfile.mktemp(".csv")
5319	file = open(file_name,"w")
5320	file.write("LONGITUDE , LATITUDE, sound , speed \n\
5321	150.916666667, -34.5, splat, 0.0\n\
5322	150.0, -34.0, pow, 10.0\n")
5323	file.close()
5324	e1 = Exposure_csv(file_name)
5325	os.remove(file_name)
5326
5327	new_title = "feast"
5328	new_values = ["chicken","soup"]
5329	e1.set_column(new_title, new_values)
5330	returned_values = e1.get_column(new_title)
5331	self.failUnless(returned_values == new_values,
5332	' Error!')
5333
5334	file_name2 = tempfile.mktemp(".csv")
5335	e1.save(file_name = file_name2)
5336	e2 = Exposure_csv(file_name2)
5337	returned_values = e2.get_column(new_title)
5338	self.failUnless(returned_values == new_values,
5339	' Error!')
5340	os.remove(file_name2)
5341
5342	def test_exposure_csv_set_column_get_column_error_checking(self):
5343	file_name = tempfile.mktemp(".csv")
5344	file = open(file_name,"w")
5345	file.write("LONGITUDE , LATITUDE, sound , speed \n\
5346	150.916666667, -34.5, splat, 0.0\n\
5347	150.0, -34.0, pow, 10.0\n")
5348	file.close()
5349	e1 = Exposure_csv(file_name)
5350	os.remove(file_name)
5351
5352	new_title = "sound"
5353	new_values = [12.5,7.6]
5354	try:
5355	e1.set_column(new_title, new_values)
5356	except TitleValueError:
5357	pass
5358	else:
5359	self.failUnless(0 ==1, 'Error not thrown error!')
5360
5361	e1.set_column(new_title, new_values, overwrite=True)
5362	returned_values = e1.get_column(new_title)
5363	self.failUnless(returned_values == new_values,
5364	' Error!')
5365
5366	new2_title = "short list"
5367	new2_values = [12.5]
5368	try:
5369	e1.set_column(new2_title, new2_values)
5370	except DataMissingValuesError:
5371	pass
5372	else:
5373	self.failUnless(0 ==1, 'Error not thrown error!')
5374
5375	new2_title = "long list"
5376	new2_values = [12.5, 7,8]
5377	try:
5378	e1.set_column(new2_title, new2_values)
5379	except DataMissingValuesError:
5380	pass
5381	else:
5382	self.failUnless(0 ==1, 'Error not thrown error!')
5383	file_name2 = tempfile.mktemp(".csv")
5384	e1.save(file_name = file_name2)
5385	e2 = Exposure_csv(file_name2)
5386	returned_values = e2.get_column(new_title)
5387	for returned, new in map(None, returned_values, new_values):
5388	self.failUnless(returned == str(new), ' Error!')
5389	#self.failUnless(returned_values == new_values, ' Error!')
5390	os.remove(file_name2)
5391
5392	try:
5393	e1.get_column("toe jam")
5394	except TitleValueError:
5395	pass
5396	else:
5397	self.failUnless(0 ==1, 'Error not thrown error!')
5398
5399	def test_exposure_csv_loading_x_y(self):
5400
5401
5402	file_name = tempfile.mktemp(".csv")
5403	file = open(file_name,"w")
5404	file.write("x, y ,sound , speed \n\
5405	115.0, 7, splat, 0.0\n\
5406	114.0, 8.0, pow, 10.0\n\
5407	114.5, 9., bang, 40.0\n")
5408	file.close()
5409	e1 = Exposure_csv(file_name, is_x_y_locations=True)
5410	gsd = e1.get_location()
5411
5412	points = gsd.get_data_points(absolute=True)
5413
5414	assert allclose(points[0][0], 115)
5415	assert allclose(points[0][1], 7)
5416	assert allclose(points[1][0], 114)
5417	assert allclose(points[1][1], 8)
5418	assert allclose(points[2][0], 114.5)
5419	assert allclose(points[2][1], 9)
5420	self.failUnless(gsd.get_geo_reference().get_zone() == -1,
5421	'Bad zone error!')
5422
5423	os.remove(file_name)
5424
5425
5426	def test_exposure_csv_loading_x_y2(self):
5427
5428	csv_file = tempfile.mktemp(".csv")
5429	fd = open(csv_file,'wb')
5430	writer = csv.writer(fd)
5431	writer.writerow(['x','y','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
5432	writer.writerow([5.5,0.5,'199770','130000','Metal','Timber',20])
5433	writer.writerow([4.5,1.0,'150000','76000','Metal','Double Brick',20])
5434	writer.writerow([4.5,1.5,'150000','76000','Metal','Brick Veneer',20])
5435	fd.close()
5436
5437	e1 = Exposure_csv(csv_file)
5438	gsd = e1.get_location()
5439
5440	points = gsd.get_data_points(absolute=True)
5441	assert allclose(points[0][0], 5.5)
5442	assert allclose(points[0][1], 0.5)
5443	assert allclose(points[1][0], 4.5)
5444	assert allclose(points[1][1], 1.0)
5445	assert allclose(points[2][0], 4.5)
5446	assert allclose(points[2][1], 1.5)
5447	self.failUnless(gsd.get_geo_reference().get_zone() == -1,
5448	'Bad zone error!')
5449
5450	os.remove(csv_file)
5451
5452	#### TESTS FOR URS 2 SWW ###
5453
5454	def create_mux(self, points_num=None):
5455	# write all the mux stuff.
5456	time_step_count = 3
5457	time_step = 0.5
5458
5459	longitudes = [150.66667, 150.83334, 151., 151.16667]
5460	latitudes = [-34.5, -34.33333, -34.16667, -34]
5461
5462	if points_num == None:
5463	points_num = len(longitudes) * len(latitudes)
5464
5465	lonlatdeps = []
5466	quantities = ['HA','UA','VA']
5467	mux_names = [WAVEHEIGHT_MUX_LABEL,
5468	EAST_VELOCITY_LABEL,
5469	NORTH_VELOCITY_LABEL]
5470	quantities_init = [[],[],[]]
5471	# urs binary is latitude fastest
5472	for i,lon in enumerate(longitudes):
5473	for j,lat in enumerate(latitudes):
5474	_ , e, n = redfearn(lat, lon)
5475	lonlatdeps.append([lon, lat, n])
5476	quantities_init[0].append(e) # HA
5477	quantities_init[1].append(n ) # UA
5478	quantities_init[2].append(e) # VA
5479	#print "lonlatdeps",lonlatdeps
5480
5481	file_handle, base_name = tempfile.mkstemp("")
5482	os.close(file_handle)
5483	os.remove(base_name)
5484
5485	files = []
5486	for i,q in enumerate(quantities):
5487	quantities_init[i] = ensure_numeric(quantities_init[i])
5488	#print "HA_init", HA_init
5489	q_time = zeros((time_step_count, points_num), Float)
5490	for time in range(time_step_count):
5491	q_time[time,:] = quantities_init[i] #* time * 4
5492
5493	#Write C files
5494	columns = 3 # long, lat , depth
5495	file = base_name + mux_names[i]
5496	#print "base_name file",file
5497	f = open(file, 'wb')
5498	files.append(file)
5499	f.write(pack('i',points_num))
5500	f.write(pack('i',time_step_count))
5501	f.write(pack('f',time_step))
5502
5503	#write lat/long info
5504	for lonlatdep in lonlatdeps:
5505	for float in lonlatdep:
5506	f.write(pack('f',float))
5507
5508	# Write quantity info
5509	for time in range(time_step_count):
5510	for point_i in range(points_num):
5511	f.write(pack('f',q_time[time,point_i]))
5512	#print " mux_names[i]", mux_names[i]
5513	#print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
5514	f.close()
5515	return base_name, files
5516
5517	def write_mux(self,lat_long_points, time_step_count, time_step,
5518	depth=None, ha=None, ua=None, va=None
5519	):
5520	"""
5521	This will write 3 non-gridded mux files, for testing.
5522	If no quantities are passed in,
5523	na and va quantities will be the Easting values.
5524	Depth and ua will be the Northing value.
5525	"""
5526	#print "lat_long_points", lat_long_points
5527	#print "time_step_count",time_step_count
5528	#print "time_step",
5529
5530	points_num = len(lat_long_points)
5531	lonlatdeps = []
5532	quantities = ['HA','UA','VA']
5533
5534	mux_names = [WAVEHEIGHT_MUX_LABEL,
5535	EAST_VELOCITY_LABEL,
5536	NORTH_VELOCITY_LABEL]
5537	quantities_init = [[],[],[]]
5538	# urs binary is latitude fastest
5539	for point in lat_long_points:
5540	lat = point[0]
5541	lon = point[1]
5542	_ , e, n = redfearn(lat, lon)
5543	if depth is None:
5544	this_depth = n
5545	else:
5546	this_depth = depth
5547	if ha is None:
5548	this_ha = e
5549	else:
5550	this_ha = ha
5551	if ua is None:
5552	this_ua = n
5553	else:
5554	this_ua = ua
5555	if va is None:
5556	this_va = e
5557	else:
5558	this_va = va
5559	lonlatdeps.append([lon, lat, this_depth])
5560	quantities_init[0].append(this_ha) # HA
5561	quantities_init[1].append(this_ua) # UA
5562	quantities_init[2].append(this_va) # VA
5563
5564	file_handle, base_name = tempfile.mkstemp("")
5565	os.close(file_handle)
5566	os.remove(base_name)
5567
5568	files = []
5569	for i,q in enumerate(quantities):
5570	quantities_init[i] = ensure_numeric(quantities_init[i])
5571	#print "HA_init", HA_init
5572	q_time = zeros((time_step_count, points_num), Float)
5573	for time in range(time_step_count):
5574	q_time[time,:] = quantities_init[i] #* time * 4
5575
5576	#Write C files
5577	columns = 3 # long, lat , depth
5578	file = base_name + mux_names[i]
5579	#print "base_name file",file
5580	f = open(file, 'wb')
5581	files.append(file)
5582	f.write(pack('i',points_num))
5583	f.write(pack('i',time_step_count))
5584	f.write(pack('f',time_step))
5585
5586	#write lat/long info
5587	for lonlatdep in lonlatdeps:
5588	for float in lonlatdep:
5589	f.write(pack('f',float))
5590
5591	# Write quantity info
5592	for time in range(time_step_count):
5593	for point_i in range(points_num):
5594	f.write(pack('f',q_time[time,point_i]))
5595	#print " mux_names[i]", mux_names[i]
5596	#print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
5597	f.close()
5598	return base_name, files
5599
5600
5601	def delete_mux(self, files):
5602	for file in files:
5603	os.remove(file)
5604
5605	def test_urs2sww_test_fail(self):
5606	points_num = -100
5607	time_step_count = 45
5608	time_step = -7
5609	file_handle, base_name = tempfile.mkstemp("")
5610	os.close(file_handle)
5611	os.remove(base_name)
5612
5613	files = []
5614	quantities = ['HA','UA','VA']
5615
5616	mux_names = [WAVEHEIGHT_MUX_LABEL,
5617	EAST_VELOCITY_LABEL,
5618	NORTH_VELOCITY_LABEL]
5619	for i,q in enumerate(quantities):
5620	#Write C files
5621	columns = 3 # long, lat , depth
5622	file = base_name + mux_names[i]
5623	f = open(file, 'wb')
5624	files.append(file)
5625	f.write(pack('i',points_num))
5626	f.write(pack('i',time_step_count))
5627	f.write(pack('f',time_step))
5628
5629	f.close()
5630	tide = 1
5631	try:
5632	urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
5633	verbose=self.verbose)
5634	except ANUGAError:
5635	pass
5636	else:
5637	self.delete_mux(files)
5638	msg = 'Should have raised exception'
5639	raise msg
5640	sww_file = base_name + '.sww'
5641	self.delete_mux(files)
5642
5643	def test_urs2sww_test_fail2(self):
5644	base_name = 'Harry-high-pants'
5645	try:
5646	urs2sww(base_name)
5647	except IOError:
5648	pass
5649	else:
5650	self.delete_mux(files)
5651	msg = 'Should have raised exception'
5652	raise msg
5653
5654	def test_urs2sww(self):
5655	tide = 1
5656	base_name, files = self.create_mux()
5657	urs2sww(base_name
5658	#, origin=(0,0,0)
5659	, mean_stage=tide
5660	, remove_nc_files=True,
5661	verbose=self.verbose
5662	)
5663	sww_file = base_name + '.sww'
5664
5665	#Let's interigate the sww file
5666	# Note, the sww info is not gridded. It is point data.
5667	fid = NetCDFFile(sww_file)
5668
5669	x = fid.variables['x'][:]
5670	y = fid.variables['y'][:]
5671	geo_reference = Geo_reference(NetCDFObject=fid)
5672
5673
5674	#Check that first coordinate is correctly represented
5675	#Work out the UTM coordinates for first point
5676	zone, e, n = redfearn(-34.5, 150.66667)
5677
5678	assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
5679
5680	# Make x and y absolute
5681	points = geo_reference.get_absolute(map(None, x, y))
5682	points = ensure_numeric(points)
5683	x = points[:,0]
5684	y = points[:,1]
5685
5686	#Check first value
5687	stage = fid.variables['stage'][:]
5688	xmomentum = fid.variables['xmomentum'][:]
5689	ymomentum = fid.variables['ymomentum'][:]
5690	elevation = fid.variables['elevation'][:]
5691	assert allclose(stage[0,0], e +tide) #Meters
5692
5693	#Check the momentums - ua
5694	#momentum = velocity*(stage-elevation)
5695	# elevation = - depth
5696	#momentum = velocity_ua *(stage+depth)
5697	# = n*(e+tide+n) based on how I'm writing these files
5698	#
5699	answer_x = n*(e+tide+n)
5700	actual_x = xmomentum[0,0]
5701	#print "answer_x",answer_x
5702	#print "actual_x",actual_x
5703	assert allclose(answer_x, actual_x) #Meters
5704
5705	#Check the momentums - va
5706	#momentum = velocity*(stage-elevation)
5707	# -(-elevation) since elevation is inverted in mux files
5708	#momentum = velocity_va *(stage+elevation)
5709	# = e*(e+tide+n) based on how I'm writing these files
5710	answer_y = e(e+tide+n) -1 # talking into account mux file format
5711	actual_y = ymomentum[0,0]
5712	#print "answer_y",answer_y
5713	#print "actual_y",actual_y
5714	assert allclose(answer_y, actual_y) #Meters
5715
5716	assert allclose(answer_x, actual_x) #Meters
5717
5718	# check the stage values, first time step.
5719	# These arrays are equal since the Easting values were used as
5720	# the stage
5721	assert allclose(stage[0], x +tide) #Meters
5722
5723	# check the elevation values.
5724	# -ve since urs measures depth, sww meshers height,
5725	# these arrays are equal since the northing values were used as
5726	# the elevation
5727	assert allclose(-elevation, y) #Meters
5728
5729	fid.close()
5730	self.delete_mux(files)
5731	os.remove(sww_file)
5732
5733	def test_urs2sww_momentum(self):
5734	tide = 1
5735	time_step_count = 3
5736	time_step = 2
5737	#lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
5738	# This is gridded
5739	lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
5740	depth=20
5741	ha=2
5742	ua=5
5743	va=-10 #-ve added to take into account mux file format where south
5744	# is positive.
5745	base_name, files = self.write_mux(lat_long_points,
5746	time_step_count, time_step,
5747	depth=depth,
5748	ha=ha,
5749	ua=ua,
5750	va=va)
5751	# write_mux(self,lat_long_points, time_step_count, time_step,
5752	# depth=None, ha=None, ua=None, va=None
5753	urs2sww(base_name
5754	#, origin=(0,0,0)
5755	, mean_stage=tide
5756	, remove_nc_files=True,
5757	verbose=self.verbose
5758	)
5759	sww_file = base_name + '.sww'
5760
5761	#Let's interigate the sww file
5762	# Note, the sww info is not gridded. It is point data.
5763	fid = NetCDFFile(sww_file)
5764
5765	x = fid.variables['x'][:]
5766	y = fid.variables['y'][:]
5767	geo_reference = Geo_reference(NetCDFObject=fid)
5768
5769	#Check first value
5770	stage = fid.variables['stage'][:]
5771	xmomentum = fid.variables['xmomentum'][:]
5772	ymomentum = fid.variables['ymomentum'][:]
5773	elevation = fid.variables['elevation'][:]
5774	#assert allclose(stage[0,0], e + tide) #Meters
5775	#print "xmomentum", xmomentum
5776	#print "ymomentum", ymomentum
5777	#Check the momentums - ua
5778	#momentum = velocity*water height
5779	#water height = mux_depth + mux_height +tide
5780	#water height = mux_depth + mux_height +tide
5781	#momentum = velocity*(mux_depth + mux_height +tide)
5782	#
5783
5784	answer = 115
5785	actual = xmomentum[0,0]
5786	assert allclose(answer, actual) #Meters^2/ sec
5787	answer = 230
5788	actual = ymomentum[0,0]
5789	#print "answer",answer
5790	#print "actual",actual
5791	assert allclose(answer, actual) #Meters^2/ sec
5792
5793	# check the stage values, first time step.
5794	# These arrays are equal since the Easting values were used as
5795	# the stage
5796
5797	#assert allclose(stage[0], x +tide) #Meters
5798
5799	# check the elevation values.
5800	# -ve since urs measures depth, sww meshers height,
5801	# these arrays are equal since the northing values were used as
5802	# the elevation
5803	#assert allclose(-elevation, y) #Meters
5804
5805	fid.close()
5806	self.delete_mux(files)
5807	os.remove(sww_file)
5808
5809
5810	def test_urs2sww_origin(self):
5811	tide = 1
5812	base_name, files = self.create_mux()
5813	urs2sww(base_name
5814	, origin=(0,0,0)
5815	, mean_stage=tide
5816	, remove_nc_files=True