Context Navigation

test_data_manager.py @ 4376

Last change on this file since 4376 was 4376, checked in by ole, 16 years ago
Prepared tests for making new slope limiters default. Still a few unprepared tests remaining and one example in test_data_manager where small timesteps are generated.
File size: 188.3 KB

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