Context Navigation

source: anuga_core/source/anuga/shallow_water/test_data_manager.py @ 4615

Last change on this file since 4615 was 4615, checked in by duncan, 17 years ago
extra unit test
File size: 231.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: