Context Navigation

test_data_manager.py @ 4699

Last change on this file since 4699 was 4699, checked in by ole, 16 years ago
Cosmetics while pondering over storage of extrema
File size: 231.7 KB

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