Context Navigation

test_data_manager.py @ 4348

Last change on this file since 4348 was 4348, checked in by duncan, 16 years ago
Changing the definition of a mux file. South is now the positive axis, not North.
File size: 194.0 KB

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