Context Navigation

test_data_manager.py @ 4704

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