Context Navigation

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

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

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

Download in other formats: