Context Navigation

test_data_manager.py @ 4705

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