Context Navigation

source: inundation/pyvolution/test_data_manager.py @ 2334

Last change on this file since 2334 was 2305, checked in by ole, 19 years ago
Fixed problem with leaking files from tests as per ticket:42
File size: 118.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 util import mean
8	import tempfile
9	import os
10	from Scientific.IO.NetCDF import NetCDFFile
11
12	from data_manager import *
13	from shallow_water import *
14	from config import epsilon
15	import data_manager
16
17	from coordinate_transforms.geo_reference import Geo_reference
18
19	class Test_Data_Manager(unittest.TestCase):
20	def setUp(self):
21	import time
22	from mesh_factory import rectangular
23
24	#Create basic mesh
25	points, vertices, boundary = rectangular(2, 2)
26
27	#Create shallow water domain
28	domain = Domain(points, vertices, boundary)
29	domain.default_order=2
30
31
32	#Set some field values
33	domain.set_quantity('elevation', lambda x,y: -x)
34	domain.set_quantity('friction', 0.03)
35
36
37	######################
38	# Boundary conditions
39	B = Transmissive_boundary(domain)
40	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
41
42
43	######################
44	#Initial condition - with jumps
45
46
47	bed = domain.quantities['elevation'].vertex_values
48	stage = zeros(bed.shape, Float)
49
50	h = 0.3
51	for i in range(stage.shape[0]):
52	if i % 2 == 0:
53	stage[i,:] = bed[i,:] + h
54	else:
55	stage[i,:] = bed[i,:]
56
57	domain.set_quantity('stage', stage)
58	self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
59
60	domain.distribute_to_vertices_and_edges()
61
62
63	self.domain = domain
64
65	C = domain.get_vertex_coordinates()
66	self.X = C[:,0:6:2].copy()
67	self.Y = C[:,1:6:2].copy()
68
69	self.F = bed
70
71
72
73
74	#Write A testfile (not realistic. Values aren't realistic)
75
76	self.test_MOST_file = 'most_small'
77
78	longitudes = [150.66667, 150.83334, 151., 151.16667]
79	latitudes = [-34.5, -34.33333, -34.16667, -34]
80
81	long_name = 'LON'
82	lat_name = 'LAT'
83
84	nx = 4
85	ny = 4
86	six = 6
87
88
89	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
90	fid = NetCDFFile(self.test_MOST_file + ext, 'w')
91
92	fid.createDimension(long_name,nx)
93	fid.createVariable(long_name,'d',(long_name,))
94	fid.variables[long_name].point_spacing='uneven'
95	fid.variables[long_name].units='degrees_east'
96	fid.variables[long_name].assignValue(longitudes)
97
98	fid.createDimension(lat_name,ny)
99	fid.createVariable(lat_name,'d',(lat_name,))
100	fid.variables[lat_name].point_spacing='uneven'
101	fid.variables[lat_name].units='degrees_north'
102	fid.variables[lat_name].assignValue(latitudes)
103
104	fid.createDimension('TIME',six)
105	fid.createVariable('TIME','d',('TIME',))
106	fid.variables['TIME'].point_spacing='uneven'
107	fid.variables['TIME'].units='seconds'
108	fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
109
110
111	name = ext[1:3].upper()
112	if name == 'E.': name = 'ELEVATION'
113	fid.createVariable(name,'d',('TIME', lat_name, long_name))
114	fid.variables[name].units='CENTIMETERS'
115	fid.variables[name].missing_value=-1.e+034
116
117	fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
118	[-0.1214216, 0, 0, 0],
119	[0, 0, 0, 0],
120	[0, 0, 0, 0]],
121	[[0.3400644, 2.291054e-005, -23.33335, -6.50198],
122	[-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
123	[2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
124	[0, 2.291054e-005, 2.291054e-005, 0]],
125	[[0.3400644, 0.0001374632, -23.31503, -6.50198],
126	[-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
127	[0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
128	[0, 0.0001374632, 0.0001374632, 0]],
129	[[0.3400644, 0.0002520159, -23.29672, -6.50198],
130	[-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
131	[0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
132	[0, 0.0002520159, 0.0002520159, 0]],
133	[[0.3400644, 0.0003665686, -23.27842, -6.50198],
134	[-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
135	[0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
136	[0, 0.0003665686, 0.0003665686, 0]],
137	[[0.3400644, 0.0004811212, -23.26012, -6.50198],
138	[-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
139	[0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
140	[0, 0.0004811212, 0.0004811212, 0]]])
141
142
143	fid.close()
144
145
146
147
148	def tearDown(self):
149	import os
150	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
151	#print 'Trying to remove', self.test_MOST_file + ext
152	os.remove(self.test_MOST_file + ext)
153
154	def test_sww_constant(self):
155	"""Test that constant sww information can be written correctly
156	(non smooth)
157	"""
158
159	import time, os
160	from Numeric import array, zeros, allclose, Float, concatenate
161	from Scientific.IO.NetCDF import NetCDFFile
162
163	self.domain.filename = 'datatest' + str(id(self))
164	self.domain.format = 'sww'
165	self.domain.smooth = False
166
167	sww = get_dataobject(self.domain)
168	sww.store_connectivity()
169
170	#Check contents
171	#Get NetCDF
172	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
173
174	# Get the variables
175	x = fid.variables['x']
176	y = fid.variables['y']
177	z = fid.variables['elevation']
178
179	volumes = fid.variables['volumes']
180
181
182	assert allclose (x[:], self.X.flat)
183	assert allclose (y[:], self.Y.flat)
184	assert allclose (z[:], self.F.flat)
185
186	V = volumes
187
188	P = len(self.domain)
189	for k in range(P):
190	assert V[k, 0] == 3*k
191	assert V[k, 1] == 3*k+1
192	assert V[k, 2] == 3*k+2
193
194
195	fid.close()
196
197	#Cleanup
198	os.remove(sww.filename)
199
200
201	def test_sww_constant_smooth(self):
202	"""Test that constant sww information can be written correctly
203	(non smooth)
204	"""
205
206	import time, os
207	from Numeric import array, zeros, allclose, Float, concatenate
208	from Scientific.IO.NetCDF import NetCDFFile
209
210	self.domain.filename = 'datatest' + str(id(self))
211	self.domain.format = 'sww'
212	self.domain.smooth = True
213
214	sww = get_dataobject(self.domain)
215	sww.store_connectivity()
216
217	#Check contents
218	#Get NetCDF
219	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
220
221	# Get the variables
222	x = fid.variables['x']
223	y = fid.variables['y']
224	z = fid.variables['elevation']
225
226	volumes = fid.variables['volumes']
227
228	X = x[:]
229	Y = y[:]
230
231	assert allclose([X[0], Y[0]], array([0.0, 0.0]))
232	assert allclose([X[1], Y[1]], array([0.0, 0.5]))
233	assert allclose([X[2], Y[2]], array([0.0, 1.0]))
234
235	assert allclose([X[4], Y[4]], array([0.5, 0.5]))
236
237	assert allclose([X[7], Y[7]], array([1.0, 0.5]))
238
239	Z = z[:]
240	assert Z[4] == -0.5
241
242	V = volumes
243	assert V[2,0] == 4
244	assert V[2,1] == 5
245	assert V[2,2] == 1
246
247	assert V[4,0] == 6
248	assert V[4,1] == 7
249	assert V[4,2] == 3
250
251
252	fid.close()
253
254	#Cleanup
255	os.remove(sww.filename)
256
257
258
259	def test_sww_variable(self):
260	"""Test that sww information can be written correctly
261	"""
262
263	import time, os
264	from Numeric import array, zeros, allclose, Float, concatenate
265	from Scientific.IO.NetCDF import NetCDFFile
266
267	self.domain.filename = 'datatest' + str(id(self))
268	self.domain.format = 'sww'
269	self.domain.smooth = True
270	self.domain.reduction = mean
271
272	sww = get_dataobject(self.domain)
273	sww.store_connectivity()
274	sww.store_timestep('stage')
275
276	#Check contents
277	#Get NetCDF
278	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
279
280
281	# Get the variables
282	x = fid.variables['x']
283	y = fid.variables['y']
284	z = fid.variables['elevation']
285	time = fid.variables['time']
286	stage = fid.variables['stage']
287
288
289	Q = self.domain.quantities['stage']
290	Q0 = Q.vertex_values[:,0]
291	Q1 = Q.vertex_values[:,1]
292	Q2 = Q.vertex_values[:,2]
293
294	A = stage[0,:]
295	#print A[0], (Q2[0,0] + Q1[1,0])/2
296	assert allclose(A[0], (Q2[0] + Q1[1])/2)
297	assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
298	assert allclose(A[2], Q0[3])
299	assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
300
301	#Center point
302	assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
303	Q0[5] + Q2[6] + Q1[7])/6)
304
305
306
307	fid.close()
308
309	#Cleanup
310	os.remove(sww.filename)
311
312
313	def test_sww_variable2(self):
314	"""Test that sww information can be written correctly
315	multiple timesteps. Use average as reduction operator
316	"""
317
318	import time, os
319	from Numeric import array, zeros, allclose, Float, concatenate
320	from Scientific.IO.NetCDF import NetCDFFile
321
322	self.domain.filename = 'datatest' + str(id(self))
323	self.domain.format = 'sww'
324	self.domain.smooth = True
325
326	self.domain.reduction = mean
327
328	sww = get_dataobject(self.domain)
329	sww.store_connectivity()
330	sww.store_timestep('stage')
331	self.domain.evolve_to_end(finaltime = 0.01)
332	sww.store_timestep('stage')
333
334
335	#Check contents
336	#Get NetCDF
337	fid = NetCDFFile(sww.filename, 'r') #Open existing file for append
338
339	# Get the variables
340	x = fid.variables['x']
341	y = fid.variables['y']
342	z = fid.variables['elevation']
343	time = fid.variables['time']
344	stage = fid.variables['stage']
345
346	#Check values
347	Q = self.domain.quantities['stage']
348	Q0 = Q.vertex_values[:,0]
349	Q1 = Q.vertex_values[:,1]
350	Q2 = Q.vertex_values[:,2]
351
352	A = stage[1,:]
353	assert allclose(A[0], (Q2[0] + Q1[1])/2)
354	assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
355	assert allclose(A[2], Q0[3])
356	assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
357
358	#Center point
359	assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
360	Q0[5] + Q2[6] + Q1[7])/6)
361
362
363	fid.close()
364
365	#Cleanup
366	os.remove(sww.filename)
367
368	def test_sww_variable3(self):
369	"""Test that sww information can be written correctly
370	multiple timesteps using a different reduction operator (min)
371	"""
372
373	import time, os
374	from Numeric import array, zeros, allclose, Float, concatenate
375	from Scientific.IO.NetCDF import NetCDFFile
376
377	self.domain.filename = 'datatest' + str(id(self))
378	self.domain.format = 'sww'
379	self.domain.smooth = True
380	self.domain.reduction = min
381
382	sww = get_dataobject(self.domain)
383	sww.store_connectivity()
384	sww.store_timestep('stage')
385
386	self.domain.evolve_to_end(finaltime = 0.01)
387	sww.store_timestep('stage')
388
389
390	#Check contents
391	#Get NetCDF
392	fid = NetCDFFile(sww.filename, 'r')
393
394
395	# Get the variables
396	x = fid.variables['x']
397	y = fid.variables['y']
398	z = fid.variables['elevation']
399	time = fid.variables['time']
400	stage = fid.variables['stage']
401
402	#Check values
403	Q = self.domain.quantities['stage']
404	Q0 = Q.vertex_values[:,0]
405	Q1 = Q.vertex_values[:,1]
406	Q2 = Q.vertex_values[:,2]
407
408	A = stage[1,:]
409	assert allclose(A[0], min(Q2[0], Q1[1]))
410	assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
411	assert allclose(A[2], Q0[3])
412	assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
413
414	#Center point
415	assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
416	Q0[5], Q2[6], Q1[7]))
417
418
419	fid.close()
420
421	#Cleanup
422	os.remove(sww.filename)
423
424
425	def test_sync(self):
426	"""Test info stored at each timestep is as expected (incl initial condition)
427	"""
428
429	import time, os, config
430	from Numeric import array, zeros, allclose, Float, concatenate
431	from Scientific.IO.NetCDF import NetCDFFile
432
433	self.domain.filename = 'synctest'
434	self.domain.format = 'sww'
435	self.domain.smooth = False
436	self.domain.store = True
437	self.domain.beta_h = 0
438
439	#Evolution
440	for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
441	stage = self.domain.quantities['stage'].vertex_values
442
443	#Get NetCDF
444	fid = NetCDFFile(self.domain.writer.filename, 'r')
445	stage_file = fid.variables['stage']
446
447	if t == 0.0:
448	assert allclose(stage, self.initial_stage)
449	assert allclose(stage_file[:], stage.flat)
450	else:
451	assert not allclose(stage, self.initial_stage)
452	assert not allclose(stage_file[:], stage.flat)
453
454	fid.close()
455
456	os.remove(self.domain.writer.filename)
457
458
459
460	def test_sww_DSG(self):
461	"""Not a test, rather a look at the sww format
462	"""
463
464	import time, os
465	from Numeric import array, zeros, allclose, Float, concatenate
466	from Scientific.IO.NetCDF import NetCDFFile
467
468	self.domain.filename = 'datatest' + str(id(self))
469	self.domain.format = 'sww'
470	self.domain.smooth = True
471	self.domain.reduction = mean
472
473	sww = get_dataobject(self.domain)
474	sww.store_connectivity()
475	sww.store_timestep('stage')
476
477	#Check contents
478	#Get NetCDF
479	fid = NetCDFFile(sww.filename, 'r')
480
481	# Get the variables
482	x = fid.variables['x']
483	y = fid.variables['y']
484	z = fid.variables['elevation']
485
486	volumes = fid.variables['volumes']
487	time = fid.variables['time']
488
489	# 2D
490	stage = fid.variables['stage']
491
492	X = x[:]
493	Y = y[:]
494	Z = z[:]
495	V = volumes[:]
496	T = time[:]
497	S = stage[:,:]
498
499	# print "****************************"
500	# print "X ",X
501	# print "****************************"
502	# print "Y ",Y
503	# print "****************************"
504	# print "Z ",Z
505	# print "****************************"
506	# print "V ",V
507	# print "****************************"
508	# print "Time ",T
509	# print "****************************"
510	# print "Stage ",S
511	# print "****************************"
512
513
514	fid.close()
515
516	#Cleanup
517	os.remove(sww.filename)
518
519
520	#def test_write_pts(self):
521	# #Obsolete
522	#
523	# #Get (enough) datapoints
524	#
525	# from Numeric import array
526	# points = array([[ 0.66666667, 0.66666667],
527	# [ 1.33333333, 1.33333333],
528	# [ 2.66666667, 0.66666667],
529	# [ 0.66666667, 2.66666667],
530	# [ 0.0, 1.0],
531	# [ 0.0, 3.0],
532	# [ 1.0, 0.0],
533	# [ 1.0, 1.0],
534	# [ 1.0, 2.0],
535	# [ 1.0, 3.0],
536	# [ 2.0, 1.0],
537	# [ 3.0, 0.0],
538	# [ 3.0, 1.0]])
539	#
540	# z = points[:,0] + 2*points[:,1]
541	#
542	# ptsfile = 'testptsfile.pts'
543	# write_ptsfile(ptsfile, points, z,
544	# attribute_name = 'linear_combination')
545	#
546	# #Check contents
547	# #Get NetCDF
548	# from Scientific.IO.NetCDF import NetCDFFile
549	# fid = NetCDFFile(ptsfile, 'r')
550	#
551	# # Get the variables
552	# #print fid.variables.keys()
553	# points1 = fid.variables['points']
554	# z1 = fid.variables['linear_combination']
555	#
556	# #Check values#
557	#
558	# #print points[:]
559	# #print ref_points
560	# assert allclose(points, points1)
561	#
562	# #print attributes[:]
563	# #print ref_elevation
564	# assert allclose(z, z1)
565	#
566	# #Cleanup
567	# fid.close()
568	#
569	# import os
570	# os.remove(ptsfile)
571
572
573
574
575	def test_dem2pts(self):
576	"""Test conversion from dem in ascii format to native NetCDF xya format
577	"""
578
579	import time, os
580	from Numeric import array, zeros, allclose, Float, concatenate
581	from Scientific.IO.NetCDF import NetCDFFile
582
583	#Write test asc file
584	root = 'demtest'
585
586	filename = root+'.asc'
587	fid = open(filename, 'w')
588	fid.write("""ncols 5
589	nrows 6
590	xllcorner 2000.5
591	yllcorner 3000.5
592	cellsize 25
593	NODATA_value -9999
594	""")
595	#Create linear function
596
597	ref_points = []
598	ref_elevation = []
599	for i in range(6):
600	y = (6-i)*25.0
601	for j in range(5):
602	x = j*25.0
603	z = x+2*y
604
605	ref_points.append( [x,y] )
606	ref_elevation.append(z)
607	fid.write('%f ' %z)
608	fid.write('\n')
609
610	fid.close()
611
612	#Write prj file with metadata
613	metafilename = root+'.prj'
614	fid = open(metafilename, 'w')
615
616
617	fid.write("""Projection UTM
618	Zone 56
619	Datum WGS84
620	Zunits NO
621	Units METERS
622	Spheroid WGS84
623	Xshift 0.0000000000
624	Yshift 10000000.0000000000
625	Parameters
626	""")
627	fid.close()
628
629	#Convert to NetCDF pts
630	convert_dem_from_ascii2netcdf(root)
631	dem2pts(root)
632
633	#Check contents
634	#Get NetCDF
635	fid = NetCDFFile(root+'.pts', 'r')
636
637	# Get the variables
638	#print fid.variables.keys()
639	points = fid.variables['points']
640	elevation = fid.variables['elevation']
641
642	#Check values
643
644	#print points[:]
645	#print ref_points
646	assert allclose(points, ref_points)
647
648	#print attributes[:]
649	#print ref_elevation
650	assert allclose(elevation, ref_elevation)
651
652	#Cleanup
653	fid.close()
654
655
656	os.remove(root + '.pts')
657	os.remove(root + '.dem')
658	os.remove(root + '.asc')
659	os.remove(root + '.prj')
660
661
662
663	def test_dem2pts_bounding_box(self):
664	"""Test conversion from dem in ascii format to native NetCDF xya format
665	"""
666
667	import time, os
668	from Numeric import array, zeros, allclose, Float, concatenate
669	from Scientific.IO.NetCDF import NetCDFFile
670
671	#Write test asc file
672	root = 'demtest'
673
674	filename = root+'.asc'
675	fid = open(filename, 'w')
676	fid.write("""ncols 5
677	nrows 6
678	xllcorner 2000.5
679	yllcorner 3000.5
680	cellsize 25
681	NODATA_value -9999
682	""")
683	#Create linear function
684
685	ref_points = []
686	ref_elevation = []
687	for i in range(6):
688	y = (6-i)*25.0
689	for j in range(5):
690	x = j*25.0
691	z = x+2*y
692
693	ref_points.append( [x,y] )
694	ref_elevation.append(z)
695	fid.write('%f ' %z)
696	fid.write('\n')
697
698	fid.close()
699
700	#Write prj file with metadata
701	metafilename = root+'.prj'
702	fid = open(metafilename, 'w')
703
704
705	fid.write("""Projection UTM
706	Zone 56
707	Datum WGS84
708	Zunits NO
709	Units METERS
710	Spheroid WGS84
711	Xshift 0.0000000000
712	Yshift 10000000.0000000000
713	Parameters
714	""")
715	fid.close()
716
717	#Convert to NetCDF pts
718	convert_dem_from_ascii2netcdf(root)
719	dem2pts(root, easting_min=2010.0, easting_max=2110.0,
720	northing_min=3035.0, northing_max=3125.5)
721
722	#Check contents
723	#Get NetCDF
724	fid = NetCDFFile(root+'.pts', 'r')
725
726	# Get the variables
727	#print fid.variables.keys()
728	points = fid.variables['points']
729	elevation = fid.variables['elevation']
730
731	#Check values
732	assert fid.xllcorner[0] == 2010.0
733	assert fid.yllcorner[0] == 3035.0
734
735	#create new reference points
736	ref_points = []
737	ref_elevation = []
738	for i in range(4):
739	y = (4-i)*25.0 + 25.0
740	y_new = y + 3000.5 - 3035.0
741	for j in range(4):
742	x = j*25.0 + 25.0
743	x_new = x + 2000.5 - 2010.0
744	z = x+2*y
745
746	ref_points.append( [x_new,y_new] )
747	ref_elevation.append(z)
748
749	#print points[:]
750	#print ref_points
751	assert allclose(points, ref_points)
752
753	#print attributes[:]
754	#print ref_elevation
755	assert allclose(elevation, ref_elevation)
756
757	#Cleanup
758	fid.close()
759
760
761	os.remove(root + '.pts')
762	os.remove(root + '.dem')
763	os.remove(root + '.asc')
764	os.remove(root + '.prj')
765
766
767
768	def test_hecras_cross_sections2pts(self):
769	"""Test conversion from HECRAS cross sections in ascii format
770	to native NetCDF pts format
771	"""
772
773	import time, os
774	from Numeric import array, zeros, allclose, Float, concatenate
775	from Scientific.IO.NetCDF import NetCDFFile
776
777	#Write test asc file
778	root = 'hecrastest'
779
780	filename = root+'.sdf'
781	fid = open(filename, 'w')
782	fid.write("""
783	# RAS export file created on Mon 15Aug2005 11:42
784	# by HEC-RAS Version 3.1.1
785
786	BEGIN HEADER:
787	UNITS: METRIC
788	DTM TYPE: TIN
789	DTM: v:\1\cit\perth_topo\river_tin
790	STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
791	CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
792	MAP PROJECTION: UTM
793	PROJECTION ZONE: 50
794	DATUM: AGD66
795	VERTICAL DATUM:
796	NUMBER OF REACHES: 19
797	NUMBER OF CROSS-SECTIONS: 2
798	END HEADER:
799
800
801	BEGIN CROSS-SECTIONS:
802
803	CROSS-SECTION:
804	STREAM ID:Southern-Wungong
805	REACH ID:Southern-Wungong
806	STATION:21410
807	CUT LINE:
808	407546.08 , 6437277.542
809	407329.32 , 6437489.482
810	407283.11 , 6437541.232
811	SURFACE LINE:
812	407546.08, 6437277.54, 52.14
813	407538.88, 6437284.58, 51.07
814	407531.68, 6437291.62, 50.56
815	407524.48, 6437298.66, 49.58
816	407517.28, 6437305.70, 49.09
817	407510.08, 6437312.74, 48.76
818	END:
819
820	CROSS-SECTION:
821	STREAM ID:Swan River
822	REACH ID:Swan Mouth
823	STATION:840.*
824	CUT LINE:
825	381178.0855 , 6452559.0685
826	380485.4755 , 6453169.272
827	SURFACE LINE:
828	381178.09, 6452559.07, 4.17
829	381169.49, 6452566.64, 4.26
830	381157.78, 6452576.96, 4.34
831	381155.97, 6452578.56, 4.35
832	381143.72, 6452589.35, 4.43
833	381136.69, 6452595.54, 4.58
834	381114.74, 6452614.88, 4.41
835	381075.53, 6452649.43, 4.17
836	381071.47, 6452653.00, 3.99
837	381063.46, 6452660.06, 3.67
838	381054.41, 6452668.03, 3.67
839	END:
840	END CROSS-SECTIONS:
841	""")
842
843	fid.close()
844
845
846	#Convert to NetCDF pts
847	hecras_cross_sections2pts(root)
848
849	#Check contents
850	#Get NetCDF
851	fid = NetCDFFile(root+'.pts', 'r')
852
853	# Get the variables
854	#print fid.variables.keys()
855	points = fid.variables['points']
856	elevation = fid.variables['elevation']
857
858	#Check values
859	ref_points = [[407546.08, 6437277.54],
860	[407538.88, 6437284.58],
861	[407531.68, 6437291.62],
862	[407524.48, 6437298.66],
863	[407517.28, 6437305.70],
864	[407510.08, 6437312.74]]
865
866	ref_points += [[381178.09, 6452559.07],
867	[381169.49, 6452566.64],
868	[381157.78, 6452576.96],
869	[381155.97, 6452578.56],
870	[381143.72, 6452589.35],
871	[381136.69, 6452595.54],
872	[381114.74, 6452614.88],
873	[381075.53, 6452649.43],
874	[381071.47, 6452653.00],
875	[381063.46, 6452660.06],
876	[381054.41, 6452668.03]]
877
878
879	ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
880	ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
881
882	#print points[:]
883	#print ref_points
884	assert allclose(points, ref_points)
885
886	#print attributes[:]
887	#print ref_elevation
888	assert allclose(elevation, ref_elevation)
889
890	#Cleanup
891	fid.close()
892
893
894	os.remove(root + '.sdf')
895	os.remove(root + '.pts')
896
897
898
899	def test_sww2dem_asc_elevation(self):
900	"""Test that sww information can be converted correctly to asc/prj
901	format readable by e.g. ArcView
902	"""
903
904	import time, os
905	from Numeric import array, zeros, allclose, Float, concatenate
906	from Scientific.IO.NetCDF import NetCDFFile
907
908	#Setup
909	self.domain.filename = 'datatest'
910
911	prjfile = self.domain.filename + '_elevation.prj'
912	ascfile = self.domain.filename + '_elevation.asc'
913	swwfile = self.domain.filename + '.sww'
914
915	self.domain.set_datadir('.')
916	self.domain.format = 'sww'
917	self.domain.smooth = True
918	self.domain.set_quantity('elevation', lambda x,y: -x-y)
919
920	self.domain.geo_reference = Geo_reference(56,308500,6189000)
921
922	sww = get_dataobject(self.domain)
923	sww.store_connectivity()
924	sww.store_timestep('stage')
925
926	self.domain.evolve_to_end(finaltime = 0.01)
927	sww.store_timestep('stage')
928
929	cellsize = 0.25
930	#Check contents
931	#Get NetCDF
932
933	fid = NetCDFFile(sww.filename, 'r')
934
935	# Get the variables
936	x = fid.variables['x'][:]
937	y = fid.variables['y'][:]
938	z = fid.variables['elevation'][:]
939	time = fid.variables['time'][:]
940	stage = fid.variables['stage'][:]
941
942
943	#Export to ascii/prj files
944	sww2dem(self.domain.filename,
945	quantity = 'elevation',
946	cellsize = cellsize,
947	verbose = False,
948	format = 'asc')
949
950	#Check prj (meta data)
951	prjid = open(prjfile)
952	lines = prjid.readlines()
953	prjid.close()
954
955	L = lines[0].strip().split()
956	assert L[0].strip().lower() == 'projection'
957	assert L[1].strip().lower() == 'utm'
958
959	L = lines[1].strip().split()
960	assert L[0].strip().lower() == 'zone'
961	assert L[1].strip().lower() == '56'
962
963	L = lines[2].strip().split()
964	assert L[0].strip().lower() == 'datum'
965	assert L[1].strip().lower() == 'wgs84'
966
967	L = lines[3].strip().split()
968	assert L[0].strip().lower() == 'zunits'
969	assert L[1].strip().lower() == 'no'
970
971	L = lines[4].strip().split()
972	assert L[0].strip().lower() == 'units'
973	assert L[1].strip().lower() == 'meters'
974
975	L = lines[5].strip().split()
976	assert L[0].strip().lower() == 'spheroid'
977	assert L[1].strip().lower() == 'wgs84'
978
979	L = lines[6].strip().split()
980	assert L[0].strip().lower() == 'xshift'
981	assert L[1].strip().lower() == '500000'
982
983	L = lines[7].strip().split()
984	assert L[0].strip().lower() == 'yshift'
985	assert L[1].strip().lower() == '10000000'
986
987	L = lines[8].strip().split()
988	assert L[0].strip().lower() == 'parameters'
989
990
991	#Check asc file
992	ascid = open(ascfile)
993	lines = ascid.readlines()
994	ascid.close()
995
996	L = lines[0].strip().split()
997	assert L[0].strip().lower() == 'ncols'
998	assert L[1].strip().lower() == '5'
999
1000	L = lines[1].strip().split()
1001	assert L[0].strip().lower() == 'nrows'
1002	assert L[1].strip().lower() == '5'
1003
1004	L = lines[2].strip().split()
1005	assert L[0].strip().lower() == 'xllcorner'
1006	assert allclose(float(L[1].strip().lower()), 308500)
1007
1008	L = lines[3].strip().split()
1009	assert L[0].strip().lower() == 'yllcorner'
1010	assert allclose(float(L[1].strip().lower()), 6189000)
1011
1012	L = lines[4].strip().split()
1013	assert L[0].strip().lower() == 'cellsize'
1014	assert allclose(float(L[1].strip().lower()), cellsize)
1015
1016	L = lines[5].strip().split()
1017	assert L[0].strip() == 'NODATA_value'
1018	assert L[1].strip().lower() == '-9999'
1019
1020	#Check grid values
1021	for j in range(5):
1022	L = lines[6+j].strip().split()
1023	y = (4-j) * cellsize
1024	for i in range(5):
1025	assert allclose(float(L[i]), -i*cellsize - y)
1026
1027
1028	fid.close()
1029
1030	#Cleanup
1031	os.remove(prjfile)
1032	os.remove(ascfile)
1033	os.remove(swwfile)
1034
1035
1036
1037	def test_sww2dem_larger(self):
1038	"""Test that sww information can be converted correctly to asc/prj
1039	format readable by e.g. ArcView. Here:
1040
1041	ncols 11
1042	nrows 11
1043	xllcorner 308500
1044	yllcorner 6189000
1045	cellsize 10.000000
1046	NODATA_value -9999
1047	-100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1048	-90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1049	-80 -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1050	-70 -80 -90 -100 -110 -120 -130 -140 -150 -160 -170
1051	-60 -70 -80 -90 -100 -110 -120 -130 -140 -150 -160
1052	-50 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150
1053	-40 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140
1054	-30 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130
1055	-20 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120
1056	-10 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110
1057	0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100
1058
1059	"""
1060
1061	import time, os
1062	from Numeric import array, zeros, allclose, Float, concatenate
1063	from Scientific.IO.NetCDF import NetCDFFile
1064
1065	#Setup
1066
1067	from mesh_factory import rectangular
1068
1069	#Create basic mesh (100m x 100m)
1070	points, vertices, boundary = rectangular(2, 2, 100, 100)
1071
1072	#Create shallow water domain
1073	domain = Domain(points, vertices, boundary)
1074	domain.default_order = 2
1075
1076	domain.filename = 'datatest'
1077
1078	prjfile = domain.filename + '_elevation.prj'
1079	ascfile = domain.filename + '_elevation.asc'
1080	swwfile = domain.filename + '.sww'
1081
1082	domain.set_datadir('.')
1083	domain.format = 'sww'
1084	domain.smooth = True
1085	domain.geo_reference = Geo_reference(56, 308500, 6189000)
1086
1087	#
1088	domain.set_quantity('elevation', lambda x,y: -x-y)
1089	domain.set_quantity('stage', 0)
1090
1091	B = Transmissive_boundary(domain)
1092	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1093
1094
1095	#
1096	sww = get_dataobject(domain)
1097	sww.store_connectivity()
1098	sww.store_timestep('stage')
1099
1100	domain.evolve_to_end(finaltime = 0.01)
1101	sww.store_timestep('stage')
1102
1103	cellsize = 10 #10m grid
1104
1105
1106	#Check contents
1107	#Get NetCDF
1108
1109	fid = NetCDFFile(sww.filename, 'r')
1110
1111	# Get the variables
1112	x = fid.variables['x'][:]
1113	y = fid.variables['y'][:]
1114	z = fid.variables['elevation'][:]
1115	time = fid.variables['time'][:]
1116	stage = fid.variables['stage'][:]
1117
1118
1119	#Export to ascii/prj files
1120	sww2dem(domain.filename,
1121	quantity = 'elevation',
1122	cellsize = cellsize,
1123	verbose = False,
1124	format = 'asc')
1125
1126
1127	#Check prj (meta data)
1128	prjid = open(prjfile)
1129	lines = prjid.readlines()
1130	prjid.close()
1131
1132	L = lines[0].strip().split()
1133	assert L[0].strip().lower() == 'projection'
1134	assert L[1].strip().lower() == 'utm'
1135
1136	L = lines[1].strip().split()
1137	assert L[0].strip().lower() == 'zone'
1138	assert L[1].strip().lower() == '56'
1139
1140	L = lines[2].strip().split()
1141	assert L[0].strip().lower() == 'datum'
1142	assert L[1].strip().lower() == 'wgs84'
1143
1144	L = lines[3].strip().split()
1145	assert L[0].strip().lower() == 'zunits'
1146	assert L[1].strip().lower() == 'no'
1147
1148	L = lines[4].strip().split()
1149	assert L[0].strip().lower() == 'units'
1150	assert L[1].strip().lower() == 'meters'
1151
1152	L = lines[5].strip().split()
1153	assert L[0].strip().lower() == 'spheroid'
1154	assert L[1].strip().lower() == 'wgs84'
1155
1156	L = lines[6].strip().split()
1157	assert L[0].strip().lower() == 'xshift'
1158	assert L[1].strip().lower() == '500000'
1159
1160	L = lines[7].strip().split()
1161	assert L[0].strip().lower() == 'yshift'
1162	assert L[1].strip().lower() == '10000000'
1163
1164	L = lines[8].strip().split()
1165	assert L[0].strip().lower() == 'parameters'
1166
1167
1168	#Check asc file
1169	ascid = open(ascfile)
1170	lines = ascid.readlines()
1171	ascid.close()
1172
1173	L = lines[0].strip().split()
1174	assert L[0].strip().lower() == 'ncols'
1175	assert L[1].strip().lower() == '11'
1176
1177	L = lines[1].strip().split()
1178	assert L[0].strip().lower() == 'nrows'
1179	assert L[1].strip().lower() == '11'
1180
1181	L = lines[2].strip().split()
1182	assert L[0].strip().lower() == 'xllcorner'
1183	assert allclose(float(L[1].strip().lower()), 308500)
1184
1185	L = lines[3].strip().split()
1186	assert L[0].strip().lower() == 'yllcorner'
1187	assert allclose(float(L[1].strip().lower()), 6189000)
1188
1189	L = lines[4].strip().split()
1190	assert L[0].strip().lower() == 'cellsize'
1191	assert allclose(float(L[1].strip().lower()), cellsize)
1192
1193	L = lines[5].strip().split()
1194	assert L[0].strip() == 'NODATA_value'
1195	assert L[1].strip().lower() == '-9999'
1196
1197	#Check grid values (FIXME: Use same strategy for other sww2dem tests)
1198	for i, line in enumerate(lines[6:]):
1199	for j, value in enumerate( line.split() ):
1200	#assert allclose(float(value), -(10-i+j)*cellsize)
1201	assert float(value) == -(10-i+j)*cellsize
1202
1203
1204	fid.close()
1205
1206	#Cleanup
1207	os.remove(prjfile)
1208	os.remove(ascfile)
1209	os.remove(swwfile)
1210
1211
1212
1213	def test_sww2dem_boundingbox(self):
1214	"""Test that sww information can be converted correctly to asc/prj
1215	format readable by e.g. ArcView.
1216	This will test that mesh can be restricted by bounding box
1217
1218	Original extent is 100m x 100m:
1219
1220	Eastings: 308500 - 308600
1221	Northings: 6189000 - 6189100
1222
1223	Bounding box changes this to the 50m x 50m square defined by
1224
1225	Eastings: 308530 - 308570
1226	Northings: 6189050 - 6189100
1227
1228	The cropped values should be
1229
1230	-130 -140 -150 -160 -170
1231	-120 -130 -140 -150 -160
1232	-110 -120 -130 -140 -150
1233	-100 -110 -120 -130 -140
1234	-90 -100 -110 -120 -130
1235	-80 -90 -100 -110 -120
1236
1237	and the new lower reference point should be
1238	Eastings: 308530
1239	Northings: 6189050
1240
1241	Original dataset is the same as in test_sww2dem_larger()
1242
1243	"""
1244
1245	import time, os
1246	from Numeric import array, zeros, allclose, Float, concatenate
1247	from Scientific.IO.NetCDF import NetCDFFile
1248
1249	#Setup
1250
1251	from mesh_factory import rectangular
1252
1253	#Create basic mesh (100m x 100m)
1254	points, vertices, boundary = rectangular(2, 2, 100, 100)
1255
1256	#Create shallow water domain
1257	domain = Domain(points, vertices, boundary)
1258	domain.default_order = 2
1259
1260	domain.filename = 'datatest'
1261
1262	prjfile = domain.filename + '_elevation.prj'
1263	ascfile = domain.filename + '_elevation.asc'
1264	swwfile = domain.filename + '.sww'
1265
1266	domain.set_datadir('.')
1267	domain.format = 'sww'
1268	domain.smooth = True
1269	domain.geo_reference = Geo_reference(56, 308500, 6189000)
1270
1271	#
1272	domain.set_quantity('elevation', lambda x,y: -x-y)
1273	domain.set_quantity('stage', 0)
1274
1275	B = Transmissive_boundary(domain)
1276	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1277
1278
1279	#
1280	sww = get_dataobject(domain)
1281	sww.store_connectivity()
1282	sww.store_timestep('stage')
1283
1284	domain.evolve_to_end(finaltime = 0.01)
1285	sww.store_timestep('stage')
1286
1287	cellsize = 10 #10m grid
1288
1289
1290	#Check contents
1291	#Get NetCDF
1292
1293	fid = NetCDFFile(sww.filename, 'r')
1294
1295	# Get the variables
1296	x = fid.variables['x'][:]
1297	y = fid.variables['y'][:]
1298	z = fid.variables['elevation'][:]
1299	time = fid.variables['time'][:]
1300	stage = fid.variables['stage'][:]
1301
1302
1303	#Export to ascii/prj files
1304	sww2dem(domain.filename,
1305	quantity = 'elevation',
1306	cellsize = cellsize,
1307	easting_min = 308530,
1308	easting_max = 308570,
1309	northing_min = 6189050,
1310	northing_max = 6189100,
1311	verbose = False,
1312	format = 'asc')
1313
1314	fid.close()
1315
1316
1317	#Check prj (meta data)
1318	prjid = open(prjfile)
1319	lines = prjid.readlines()
1320	prjid.close()
1321
1322	L = lines[0].strip().split()
1323	assert L[0].strip().lower() == 'projection'
1324	assert L[1].strip().lower() == 'utm'
1325
1326	L = lines[1].strip().split()
1327	assert L[0].strip().lower() == 'zone'
1328	assert L[1].strip().lower() == '56'
1329
1330	L = lines[2].strip().split()
1331	assert L[0].strip().lower() == 'datum'
1332	assert L[1].strip().lower() == 'wgs84'
1333
1334	L = lines[3].strip().split()
1335	assert L[0].strip().lower() == 'zunits'
1336	assert L[1].strip().lower() == 'no'
1337
1338	L = lines[4].strip().split()
1339	assert L[0].strip().lower() == 'units'
1340	assert L[1].strip().lower() == 'meters'
1341
1342	L = lines[5].strip().split()
1343	assert L[0].strip().lower() == 'spheroid'
1344	assert L[1].strip().lower() == 'wgs84'
1345
1346	L = lines[6].strip().split()
1347	assert L[0].strip().lower() == 'xshift'
1348	assert L[1].strip().lower() == '500000'
1349
1350	L = lines[7].strip().split()
1351	assert L[0].strip().lower() == 'yshift'
1352	assert L[1].strip().lower() == '10000000'
1353
1354	L = lines[8].strip().split()
1355	assert L[0].strip().lower() == 'parameters'
1356
1357
1358	#Check asc file
1359	ascid = open(ascfile)
1360	lines = ascid.readlines()
1361	ascid.close()
1362
1363	L = lines[0].strip().split()
1364	assert L[0].strip().lower() == 'ncols'
1365	assert L[1].strip().lower() == '5'
1366
1367	L = lines[1].strip().split()
1368	assert L[0].strip().lower() == 'nrows'
1369	assert L[1].strip().lower() == '6'
1370
1371	L = lines[2].strip().split()
1372	assert L[0].strip().lower() == 'xllcorner'
1373	assert allclose(float(L[1].strip().lower()), 308530)
1374
1375	L = lines[3].strip().split()
1376	assert L[0].strip().lower() == 'yllcorner'
1377	assert allclose(float(L[1].strip().lower()), 6189050)
1378
1379	L = lines[4].strip().split()
1380	assert L[0].strip().lower() == 'cellsize'
1381	assert allclose(float(L[1].strip().lower()), cellsize)
1382
1383	L = lines[5].strip().split()
1384	assert L[0].strip() == 'NODATA_value'
1385	assert L[1].strip().lower() == '-9999'
1386
1387	#Check grid values
1388	for i, line in enumerate(lines[6:]):
1389	for j, value in enumerate( line.split() ):
1390	#assert float(value) == -(10-i+j)*cellsize
1391	assert float(value) == -(10-i+j+3)*cellsize
1392
1393
1394
1395	#Cleanup
1396	os.remove(prjfile)
1397	os.remove(ascfile)
1398	os.remove(swwfile)
1399
1400
1401
1402	def test_sww2dem_asc_stage_reduction(self):
1403	"""Test that sww information can be converted correctly to asc/prj
1404	format readable by e.g. ArcView
1405
1406	This tests the reduction of quantity stage using min
1407	"""
1408
1409	import time, os
1410	from Numeric import array, zeros, allclose, Float, concatenate
1411	from Scientific.IO.NetCDF import NetCDFFile
1412
1413	#Setup
1414	self.domain.filename = 'datatest'
1415
1416	prjfile = self.domain.filename + '_stage.prj'
1417	ascfile = self.domain.filename + '_stage.asc'
1418	swwfile = self.domain.filename + '.sww'
1419
1420	self.domain.set_datadir('.')
1421	self.domain.format = 'sww'
1422	self.domain.smooth = True
1423	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1424
1425	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1426
1427
1428	sww = get_dataobject(self.domain)
1429	sww.store_connectivity()
1430	sww.store_timestep('stage')
1431
1432	self.domain.evolve_to_end(finaltime = 0.01)
1433	sww.store_timestep('stage')
1434
1435	cellsize = 0.25
1436	#Check contents
1437	#Get NetCDF
1438
1439	fid = NetCDFFile(sww.filename, 'r')
1440
1441	# Get the variables
1442	x = fid.variables['x'][:]
1443	y = fid.variables['y'][:]
1444	z = fid.variables['elevation'][:]
1445	time = fid.variables['time'][:]
1446	stage = fid.variables['stage'][:]
1447
1448
1449	#Export to ascii/prj files
1450	sww2dem(self.domain.filename,
1451	quantity = 'stage',
1452	cellsize = cellsize,
1453	reduction = min,
1454	format = 'asc')
1455
1456
1457	#Check asc file
1458	ascid = open(ascfile)
1459	lines = ascid.readlines()
1460	ascid.close()
1461
1462	L = lines[0].strip().split()
1463	assert L[0].strip().lower() == 'ncols'
1464	assert L[1].strip().lower() == '5'
1465
1466	L = lines[1].strip().split()
1467	assert L[0].strip().lower() == 'nrows'
1468	assert L[1].strip().lower() == '5'
1469
1470	L = lines[2].strip().split()
1471	assert L[0].strip().lower() == 'xllcorner'
1472	assert allclose(float(L[1].strip().lower()), 308500)
1473
1474	L = lines[3].strip().split()
1475	assert L[0].strip().lower() == 'yllcorner'
1476	assert allclose(float(L[1].strip().lower()), 6189000)
1477
1478	L = lines[4].strip().split()
1479	assert L[0].strip().lower() == 'cellsize'
1480	assert allclose(float(L[1].strip().lower()), cellsize)
1481
1482	L = lines[5].strip().split()
1483	assert L[0].strip() == 'NODATA_value'
1484	assert L[1].strip().lower() == '-9999'
1485
1486
1487	#Check grid values (where applicable)
1488	for j in range(5):
1489	if j%2 == 0:
1490	L = lines[6+j].strip().split()
1491	jj = 4-j
1492	for i in range(5):
1493	if i%2 == 0:
1494	index = jj/2 + i/2*3
1495	val0 = stage[0,index]
1496	val1 = stage[1,index]
1497
1498	#print i, j, index, ':', L[i], val0, val1
1499	assert allclose(float(L[i]), min(val0, val1))
1500
1501
1502	fid.close()
1503
1504	#Cleanup
1505	os.remove(prjfile)
1506	os.remove(ascfile)
1507	#os.remove(swwfile)
1508
1509
1510
1511	def test_sww2dem_asc_derived_quantity(self):
1512	"""Test that sww information can be converted correctly to asc/prj
1513	format readable by e.g. ArcView
1514
1515	This tests the use of derived quantities
1516	"""
1517
1518	import time, os
1519	from Numeric import array, zeros, allclose, Float, concatenate
1520	from Scientific.IO.NetCDF import NetCDFFile
1521
1522	#Setup
1523	self.domain.filename = 'datatest'
1524
1525	prjfile = self.domain.filename + '_depth.prj'
1526	ascfile = self.domain.filename + '_depth.asc'
1527	swwfile = self.domain.filename + '.sww'
1528
1529	self.domain.set_datadir('.')
1530	self.domain.format = 'sww'
1531	self.domain.smooth = True
1532	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1533	self.domain.set_quantity('stage', 0.0)
1534
1535	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1536
1537
1538	sww = get_dataobject(self.domain)
1539	sww.store_connectivity()
1540	sww.store_timestep('stage')
1541
1542	self.domain.evolve_to_end(finaltime = 0.01)
1543	sww.store_timestep('stage')
1544
1545	cellsize = 0.25
1546	#Check contents
1547	#Get NetCDF
1548
1549	fid = NetCDFFile(sww.filename, 'r')
1550
1551	# Get the variables
1552	x = fid.variables['x'][:]
1553	y = fid.variables['y'][:]
1554	z = fid.variables['elevation'][:]
1555	time = fid.variables['time'][:]
1556	stage = fid.variables['stage'][:]
1557
1558
1559	#Export to ascii/prj files
1560	sww2dem(self.domain.filename,
1561	basename_out = 'datatest_depth',
1562	quantity = 'stage - elevation',
1563	cellsize = cellsize,
1564	reduction = min,
1565	format = 'asc',
1566	verbose = False)
1567
1568
1569	#Check asc file
1570	ascid = open(ascfile)
1571	lines = ascid.readlines()
1572	ascid.close()
1573
1574	L = lines[0].strip().split()
1575	assert L[0].strip().lower() == 'ncols'
1576	assert L[1].strip().lower() == '5'
1577
1578	L = lines[1].strip().split()
1579	assert L[0].strip().lower() == 'nrows'
1580	assert L[1].strip().lower() == '5'
1581
1582	L = lines[2].strip().split()
1583	assert L[0].strip().lower() == 'xllcorner'
1584	assert allclose(float(L[1].strip().lower()), 308500)
1585
1586	L = lines[3].strip().split()
1587	assert L[0].strip().lower() == 'yllcorner'
1588	assert allclose(float(L[1].strip().lower()), 6189000)
1589
1590	L = lines[4].strip().split()
1591	assert L[0].strip().lower() == 'cellsize'
1592	assert allclose(float(L[1].strip().lower()), cellsize)
1593
1594	L = lines[5].strip().split()
1595	assert L[0].strip() == 'NODATA_value'
1596	assert L[1].strip().lower() == '-9999'
1597
1598
1599	#Check grid values (where applicable)
1600	for j in range(5):
1601	if j%2 == 0:
1602	L = lines[6+j].strip().split()
1603	jj = 4-j
1604	for i in range(5):
1605	if i%2 == 0:
1606	index = jj/2 + i/2*3
1607	val0 = stage[0,index] - z[index]
1608	val1 = stage[1,index] - z[index]
1609
1610	#print i, j, index, ':', L[i], val0, val1
1611	assert allclose(float(L[i]), min(val0, val1))
1612
1613
1614	fid.close()
1615
1616	#Cleanup
1617	os.remove(prjfile)
1618	os.remove(ascfile)
1619	#os.remove(swwfile)
1620
1621
1622
1623
1624
1625	def test_sww2dem_asc_missing_points(self):
1626	"""Test that sww information can be converted correctly to asc/prj
1627	format readable by e.g. ArcView
1628
1629	This test includes the writing of missing values
1630	"""
1631
1632	import time, os
1633	from Numeric import array, zeros, allclose, Float, concatenate
1634	from Scientific.IO.NetCDF import NetCDFFile
1635
1636	#Setup mesh not coinciding with rectangle.
1637	#This will cause missing values to occur in gridded data
1638
1639
1640	points = [ [1.0, 1.0],
1641	[0.5, 0.5], [1.0, 0.5],
1642	[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1643
1644	vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1645
1646	#Create shallow water domain
1647	domain = Domain(points, vertices)
1648	domain.default_order=2
1649
1650
1651	#Set some field values
1652	domain.set_quantity('elevation', lambda x,y: -x-y)
1653	domain.set_quantity('friction', 0.03)
1654
1655
1656	######################
1657	# Boundary conditions
1658	B = Transmissive_boundary(domain)
1659	domain.set_boundary( {'exterior': B} )
1660
1661
1662	######################
1663	#Initial condition - with jumps
1664
1665	bed = domain.quantities['elevation'].vertex_values
1666	stage = zeros(bed.shape, Float)
1667
1668	h = 0.3
1669	for i in range(stage.shape[0]):
1670	if i % 2 == 0:
1671	stage[i,:] = bed[i,:] + h
1672	else:
1673	stage[i,:] = bed[i,:]
1674
1675	domain.set_quantity('stage', stage)
1676	domain.distribute_to_vertices_and_edges()
1677
1678	domain.filename = 'datatest'
1679
1680	prjfile = domain.filename + '_elevation.prj'
1681	ascfile = domain.filename + '_elevation.asc'
1682	swwfile = domain.filename + '.sww'
1683
1684	domain.set_datadir('.')
1685	domain.format = 'sww'
1686	domain.smooth = True
1687
1688	domain.geo_reference = Geo_reference(56,308500,6189000)
1689
1690	sww = get_dataobject(domain)
1691	sww.store_connectivity()
1692	sww.store_timestep('stage')
1693
1694	cellsize = 0.25
1695	#Check contents
1696	#Get NetCDF
1697
1698	fid = NetCDFFile(swwfile, 'r')
1699
1700	# Get the variables
1701	x = fid.variables['x'][:]
1702	y = fid.variables['y'][:]
1703	z = fid.variables['elevation'][:]
1704	time = fid.variables['time'][:]
1705
1706	try:
1707	geo_reference = Geo_reference(NetCDFObject=fid)
1708	except AttributeError, e:
1709	geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1710
1711	#Export to ascii/prj files
1712	sww2dem(domain.filename,
1713	quantity = 'elevation',
1714	cellsize = cellsize,
1715	verbose = False,
1716	format = 'asc')
1717
1718
1719	#Check asc file
1720	ascid = open(ascfile)
1721	lines = ascid.readlines()
1722	ascid.close()
1723
1724	L = lines[0].strip().split()
1725	assert L[0].strip().lower() == 'ncols'
1726	assert L[1].strip().lower() == '5'
1727
1728	L = lines[1].strip().split()
1729	assert L[0].strip().lower() == 'nrows'
1730	assert L[1].strip().lower() == '5'
1731
1732	L = lines[2].strip().split()
1733	assert L[0].strip().lower() == 'xllcorner'
1734	assert allclose(float(L[1].strip().lower()), 308500)
1735
1736	L = lines[3].strip().split()
1737	assert L[0].strip().lower() == 'yllcorner'
1738	assert allclose(float(L[1].strip().lower()), 6189000)
1739
1740	L = lines[4].strip().split()
1741	assert L[0].strip().lower() == 'cellsize'
1742	assert allclose(float(L[1].strip().lower()), cellsize)
1743
1744	L = lines[5].strip().split()
1745	assert L[0].strip() == 'NODATA_value'
1746	assert L[1].strip().lower() == '-9999'
1747
1748	#Check grid values
1749	for j in range(5):
1750	L = lines[6+j].strip().split()
1751	assert len(L) == 5
1752	y = (4-j) * cellsize
1753
1754	for i in range(5):
1755	#print i
1756	if i+j >= 4:
1757	assert allclose(float(L[i]), -i*cellsize - y)
1758	else:
1759	#Missing values
1760	assert allclose(float(L[i]), -9999)
1761
1762
1763
1764	fid.close()
1765
1766	#Cleanup
1767	os.remove(prjfile)
1768	os.remove(ascfile)
1769	os.remove(swwfile)
1770
1771	def test_sww2ers_simple(self):
1772	"""Test that sww information can be converted correctly to asc/prj
1773	format readable by e.g. ArcView
1774	"""
1775
1776	import time, os
1777	from Numeric import array, zeros, allclose, Float, concatenate
1778	from Scientific.IO.NetCDF import NetCDFFile
1779
1780
1781	NODATA_value = 1758323
1782
1783	#Setup
1784	self.domain.filename = 'datatest'
1785
1786	headerfile = self.domain.filename + '.ers'
1787	swwfile = self.domain.filename + '.sww'
1788
1789	self.domain.set_datadir('.')
1790	self.domain.format = 'sww'
1791	self.domain.smooth = True
1792	self.domain.set_quantity('elevation', lambda x,y: -x-y)
1793
1794	self.domain.geo_reference = Geo_reference(56,308500,6189000)
1795
1796	sww = get_dataobject(self.domain)
1797	sww.store_connectivity()
1798	sww.store_timestep('stage')
1799
1800	self.domain.evolve_to_end(finaltime = 0.01)
1801	sww.store_timestep('stage')
1802
1803	cellsize = 0.25
1804	#Check contents
1805	#Get NetCDF
1806
1807	fid = NetCDFFile(sww.filename, 'r')
1808
1809	# Get the variables
1810	x = fid.variables['x'][:]
1811	y = fid.variables['y'][:]
1812	z = fid.variables['elevation'][:]
1813	time = fid.variables['time'][:]
1814	stage = fid.variables['stage'][:]
1815
1816
1817	#Export to ers files
1818	#sww2ers(self.domain.filename,
1819	# quantity = 'elevation',
1820	# cellsize = cellsize,
1821	# verbose = False)
1822
1823	sww2dem(self.domain.filename,
1824	quantity = 'elevation',
1825	cellsize = cellsize,
1826	NODATA_value = NODATA_value,
1827	verbose = False,
1828	format = 'ers')
1829
1830	#Check header data
1831	from ermapper_grids import read_ermapper_header, read_ermapper_data
1832
1833	header = read_ermapper_header(self.domain.filename + '_elevation.ers')
1834	#print header
1835	assert header['projection'].lower() == '"utm-56"'
1836	assert header['datum'].lower() == '"wgs84"'
1837	assert header['units'].lower() == '"meters"'
1838	assert header['value'].lower() == '"elevation"'
1839	assert header['xdimension'] == '0.25'
1840	assert header['ydimension'] == '0.25'
1841	assert float(header['eastings']) == 308500.0 #xllcorner
1842	assert float(header['northings']) == 6189000.0 #yllcorner
1843	assert int(header['nroflines']) == 5
1844	assert int(header['nrofcellsperline']) == 5
1845	assert int(header['nullcellvalue']) == NODATA_value
1846	#FIXME - there is more in the header
1847
1848
1849	#Check grid data
1850	grid = read_ermapper_data(self.domain.filename + '_elevation')
1851
1852	#FIXME (Ole): Why is this the desired reference grid for -x-y?
1853	ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
1854	-1, -1.25, -1.5, -1.75, -2.0,
1855	-0.75, -1.0, -1.25, -1.5, -1.75,
1856	-0.5, -0.75, -1.0, -1.25, -1.5,
1857	-0.25, -0.5, -0.75, -1.0, -1.25]
1858
1859
1860	#print grid
1861	assert allclose(grid, ref_grid)
1862
1863	fid.close()
1864
1865	#Cleanup
1866	#FIXME the file clean-up doesn't work (eg Permission Denied Error)
1867	#Done (Ole) - it was because sww2ers didn't close it's sww file
1868	os.remove(sww.filename)
1869	os.remove(self.domain.filename + '_elevation')
1870	os.remove(self.domain.filename + '_elevation.ers')
1871
1872
1873
1874	def test_ferret2sww1(self):
1875	"""Test that georeferencing etc works when converting from
1876	ferret format (lat/lon) to sww format (UTM)
1877	"""
1878	from Scientific.IO.NetCDF import NetCDFFile
1879	import os, sys
1880
1881	#The test file has
1882	# LON = 150.66667, 150.83334, 151, 151.16667
1883	# LAT = -34.5, -34.33333, -34.16667, -34 ;
1884	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1885	#
1886	# First value (index=0) in small_ha.nc is 0.3400644 cm,
1887	# Fourth value (index==3) is -6.50198 cm
1888
1889
1890
1891	#Read
1892	from coordinate_transforms.redfearn import redfearn
1893	#fid = NetCDFFile(self.test_MOST_file)
1894	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
1895	first_value = fid.variables['HA'][:][0,0,0]
1896	fourth_value = fid.variables['HA'][:][0,0,3]
1897	fid.close()
1898
1899
1900	#Call conversion (with zero origin)
1901	#ferret2sww('small', verbose=False,
1902	# origin = (56, 0, 0))
1903	ferret2sww(self.test_MOST_file, verbose=False,
1904	origin = (56, 0, 0))
1905
1906	#Work out the UTM coordinates for first point
1907	zone, e, n = redfearn(-34.5, 150.66667)
1908	#print zone, e, n
1909
1910	#Read output file 'small.sww'
1911	#fid = NetCDFFile('small.sww')
1912	fid = NetCDFFile(self.test_MOST_file + '.sww')
1913
1914	x = fid.variables['x'][:]
1915	y = fid.variables['y'][:]
1916
1917	#Check that first coordinate is correctly represented
1918	assert allclose(x[0], e)
1919	assert allclose(y[0], n)
1920
1921	#Check first value
1922	stage = fid.variables['stage'][:]
1923	xmomentum = fid.variables['xmomentum'][:]
1924	ymomentum = fid.variables['ymomentum'][:]
1925
1926	#print ymomentum
1927
1928	assert allclose(stage[0,0], first_value/100) #Meters
1929
1930	#Check fourth value
1931	assert allclose(stage[0,3], fourth_value/100) #Meters
1932
1933	fid.close()
1934
1935	#Cleanup
1936	import os
1937	os.remove(self.test_MOST_file + '.sww')
1938
1939
1940	def test_ferret2sww_2(self):
1941	"""Test that georeferencing etc works when converting from
1942	ferret format (lat/lon) to sww format (UTM)
1943	"""
1944	from Scientific.IO.NetCDF import NetCDFFile
1945
1946	#The test file has
1947	# LON = 150.66667, 150.83334, 151, 151.16667
1948	# LAT = -34.5, -34.33333, -34.16667, -34 ;
1949	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
1950	#
1951	# First value (index=0) in small_ha.nc is 0.3400644 cm,
1952	# Fourth value (index==3) is -6.50198 cm
1953
1954
1955	from coordinate_transforms.redfearn import redfearn
1956
1957	#fid = NetCDFFile('small_ha.nc')
1958	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
1959
1960	#Pick a coordinate and a value
1961
1962	time_index = 1
1963	lat_index = 0
1964	lon_index = 2
1965
1966	test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
1967	test_time = fid.variables['TIME'][:][time_index]
1968	test_lat = fid.variables['LAT'][:][lat_index]
1969	test_lon = fid.variables['LON'][:][lon_index]
1970
1971	linear_point_index = lat_index*4 + lon_index
1972	fid.close()
1973
1974	#Call conversion (with zero origin)
1975	ferret2sww(self.test_MOST_file, verbose=False,
1976	origin = (56, 0, 0))
1977
1978
1979	#Work out the UTM coordinates for test point
1980	zone, e, n = redfearn(test_lat, test_lon)
1981
1982	#Read output file 'small.sww'
1983	fid = NetCDFFile(self.test_MOST_file + '.sww')
1984
1985	x = fid.variables['x'][:]
1986	y = fid.variables['y'][:]
1987
1988	#Check that test coordinate is correctly represented
1989	assert allclose(x[linear_point_index], e)
1990	assert allclose(y[linear_point_index], n)
1991
1992	#Check test value
1993	stage = fid.variables['stage'][:]
1994
1995	assert allclose(stage[time_index, linear_point_index], test_value/100)
1996
1997	fid.close()
1998
1999	#Cleanup
2000	import os
2001	os.remove(self.test_MOST_file + '.sww')
2002
2003
2004
2005	def test_ferret2sww3(self):
2006	"""Elevation included
2007	"""
2008	from Scientific.IO.NetCDF import NetCDFFile
2009
2010	#The test file has
2011	# LON = 150.66667, 150.83334, 151, 151.16667
2012	# LAT = -34.5, -34.33333, -34.16667, -34 ;
2013	# ELEVATION = [-1 -2 -3 -4
2014	# -5 -6 -7 -8
2015	# ...
2016	# ... -16]
2017	# where the top left corner is -1m,
2018	# and the ll corner is -13.0m
2019	#
2020	# First value (index=0) in small_ha.nc is 0.3400644 cm,
2021	# Fourth value (index==3) is -6.50198 cm
2022
2023	from coordinate_transforms.redfearn import redfearn
2024	import os
2025	fid1 = NetCDFFile('test_ha.nc','w')
2026	fid2 = NetCDFFile('test_ua.nc','w')
2027	fid3 = NetCDFFile('test_va.nc','w')
2028	fid4 = NetCDFFile('test_e.nc','w')
2029
2030	h1_list = [150.66667,150.83334,151.]
2031	h2_list = [-34.5,-34.33333]
2032
2033	long_name = 'LON'
2034	lat_name = 'LAT'
2035
2036	nx = 3
2037	ny = 2
2038
2039	for fid in [fid1,fid2,fid3]:
2040	fid.createDimension(long_name,nx)
2041	fid.createVariable(long_name,'d',(long_name,))
2042	fid.variables[long_name].point_spacing='uneven'
2043	fid.variables[long_name].units='degrees_east'
2044	fid.variables[long_name].assignValue(h1_list)
2045
2046	fid.createDimension(lat_name,ny)
2047	fid.createVariable(lat_name,'d',(lat_name,))
2048	fid.variables[lat_name].point_spacing='uneven'
2049	fid.variables[lat_name].units='degrees_north'
2050	fid.variables[lat_name].assignValue(h2_list)
2051
2052	fid.createDimension('TIME',2)
2053	fid.createVariable('TIME','d',('TIME',))
2054	fid.variables['TIME'].point_spacing='uneven'
2055	fid.variables['TIME'].units='seconds'
2056	fid.variables['TIME'].assignValue([0.,1.])
2057	if fid == fid3: break
2058
2059
2060	for fid in [fid4]:
2061	fid.createDimension(long_name,nx)
2062	fid.createVariable(long_name,'d',(long_name,))
2063	fid.variables[long_name].point_spacing='uneven'
2064	fid.variables[long_name].units='degrees_east'
2065	fid.variables[long_name].assignValue(h1_list)
2066
2067	fid.createDimension(lat_name,ny)
2068	fid.createVariable(lat_name,'d',(lat_name,))
2069	fid.variables[lat_name].point_spacing='uneven'
2070	fid.variables[lat_name].units='degrees_north'
2071	fid.variables[lat_name].assignValue(h2_list)
2072
2073	name = {}
2074	name[fid1]='HA'
2075	name[fid2]='UA'
2076	name[fid3]='VA'
2077	name[fid4]='ELEVATION'
2078
2079	units = {}
2080	units[fid1]='cm'
2081	units[fid2]='cm/s'
2082	units[fid3]='cm/s'
2083	units[fid4]='m'
2084
2085	values = {}
2086	values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2087	values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2088	values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2089	values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2090
2091	for fid in [fid1,fid2,fid3]:
2092	fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
2093	fid.variables[name[fid]].point_spacing='uneven'
2094	fid.variables[name[fid]].units=units[fid]
2095	fid.variables[name[fid]].assignValue(values[fid])
2096	fid.variables[name[fid]].missing_value = -99999999.
2097	if fid == fid3: break
2098
2099	for fid in [fid4]:
2100	fid.createVariable(name[fid],'d',(lat_name,long_name))
2101	fid.variables[name[fid]].point_spacing='uneven'
2102	fid.variables[name[fid]].units=units[fid]
2103	fid.variables[name[fid]].assignValue(values[fid])
2104	fid.variables[name[fid]].missing_value = -99999999.
2105
2106
2107	fid1.sync(); fid1.close()
2108	fid2.sync(); fid2.close()
2109	fid3.sync(); fid3.close()
2110	fid4.sync(); fid4.close()
2111
2112	fid1 = NetCDFFile('test_ha.nc','r')
2113	fid2 = NetCDFFile('test_e.nc','r')
2114	fid3 = NetCDFFile('test_va.nc','r')
2115
2116
2117	first_amp = fid1.variables['HA'][:][0,0,0]
2118	third_amp = fid1.variables['HA'][:][0,0,2]
2119	first_elevation = fid2.variables['ELEVATION'][0,0]
2120	third_elevation= fid2.variables['ELEVATION'][:][0,2]
2121	first_speed = fid3.variables['VA'][0,0,0]
2122	third_speed = fid3.variables['VA'][:][0,0,2]
2123
2124	fid1.close()
2125	fid2.close()
2126	fid3.close()
2127
2128	#Call conversion (with zero origin)
2129	ferret2sww('test', verbose=False,
2130	origin = (56, 0, 0))
2131
2132	os.remove('test_va.nc')
2133	os.remove('test_ua.nc')
2134	os.remove('test_ha.nc')
2135	os.remove('test_e.nc')
2136
2137	#Read output file 'test.sww'
2138	fid = NetCDFFile('test.sww')
2139
2140
2141	#Check first value
2142	elevation = fid.variables['elevation'][:]
2143	stage = fid.variables['stage'][:]
2144	xmomentum = fid.variables['xmomentum'][:]
2145	ymomentum = fid.variables['ymomentum'][:]
2146
2147	#print ymomentum
2148	first_height = first_amp/100 - first_elevation
2149	third_height = third_amp/100 - third_elevation
2150	first_momentum=first_speed*first_height/100
2151	third_momentum=third_speed*third_height/100
2152
2153	assert allclose(ymomentum[0][0],first_momentum) #Meters
2154	assert allclose(ymomentum[0][2],third_momentum) #Meters
2155
2156	fid.close()
2157
2158	#Cleanup
2159	os.remove('test.sww')
2160
2161
2162
2163
2164	def test_ferret2sww_nz_origin(self):
2165	from Scientific.IO.NetCDF import NetCDFFile
2166	from coordinate_transforms.redfearn import redfearn
2167
2168	#Call conversion (with nonzero origin)
2169	ferret2sww(self.test_MOST_file, verbose=False,
2170	origin = (56, 100000, 200000))
2171
2172
2173	#Work out the UTM coordinates for first point
2174	zone, e, n = redfearn(-34.5, 150.66667)
2175
2176	#Read output file 'small.sww'
2177	#fid = NetCDFFile('small.sww', 'r')
2178	fid = NetCDFFile(self.test_MOST_file + '.sww')
2179
2180	x = fid.variables['x'][:]
2181	y = fid.variables['y'][:]
2182
2183	#Check that first coordinate is correctly represented
2184	assert allclose(x[0], e-100000)
2185	assert allclose(y[0], n-200000)
2186
2187	fid.close()
2188
2189	#Cleanup
2190	os.remove(self.test_MOST_file + '.sww')
2191
2192
2193
2194	def test_sww_extent(self):
2195	"""Not a test, rather a look at the sww format
2196	"""
2197
2198	import time, os
2199	from Numeric import array, zeros, allclose, Float, concatenate
2200	from Scientific.IO.NetCDF import NetCDFFile
2201
2202	self.domain.filename = 'datatest' + str(id(self))
2203	self.domain.format = 'sww'
2204	self.domain.smooth = True
2205	self.domain.reduction = mean
2206	self.domain.set_datadir('.')
2207
2208
2209	sww = get_dataobject(self.domain)
2210	sww.store_connectivity()
2211	sww.store_timestep('stage')
2212	self.domain.time = 2.
2213
2214	#Modify stage at second timestep
2215	stage = self.domain.quantities['stage'].vertex_values
2216	self.domain.set_quantity('stage', stage/2)
2217
2218	sww.store_timestep('stage')
2219
2220	file_and_extension_name = self.domain.filename + ".sww"
2221	#print "file_and_extension_name",file_and_extension_name
2222	[xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2223	extent_sww(file_and_extension_name )
2224
2225	assert allclose(xmin, 0.0)
2226	assert allclose(xmax, 1.0)
2227	assert allclose(ymin, 0.0)
2228	assert allclose(ymax, 1.0)
2229	assert allclose(stagemin, -0.85)
2230	assert allclose(stagemax, 0.15)
2231
2232
2233	#Cleanup
2234	os.remove(sww.filename)
2235
2236
2237
2238	def test_sww2domain1(self):
2239	################################################
2240	#Create a test domain, and evolve and save it.
2241	################################################
2242	from mesh_factory import rectangular
2243	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2244	Constant_height, Time_boundary, Transmissive_boundary
2245	from Numeric import array
2246
2247	#Create basic mesh
2248
2249	yiel=0.01
2250	points, vertices, boundary = rectangular(10,10)
2251
2252	#Create shallow water domain
2253	domain = Domain(points, vertices, boundary)
2254	domain.geo_reference = Geo_reference(56,11,11)
2255	domain.smooth = False
2256	domain.visualise = False
2257	domain.store = True
2258	domain.filename = 'bedslope'
2259	domain.default_order=2
2260	#Bed-slope and friction
2261	domain.set_quantity('elevation', lambda x,y: -x/3)
2262	domain.set_quantity('friction', 0.1)
2263	# Boundary conditions
2264	from math import sin, pi
2265	Br = Reflective_boundary(domain)
2266	Bt = Transmissive_boundary(domain)
2267	Bd = Dirichlet_boundary([0.2,0.,0.])
2268	Bw = Time_boundary(domain=domain,
2269	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
2270
2271	#domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2272	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2273
2274	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2275	#Initial condition
2276	h = 0.05
2277	elevation = domain.quantities['elevation'].vertex_values
2278	domain.set_quantity('stage', elevation + h)
2279	#elevation = domain.get_quantity('elevation')
2280	#domain.set_quantity('stage', elevation + h)
2281
2282	domain.check_integrity()
2283	#Evolution
2284	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2285	# domain.write_time()
2286	pass
2287
2288
2289	##########################################
2290	#Import the example's file as a new domain
2291	##########################################
2292	from data_manager import sww2domain
2293	from Numeric import allclose
2294	import os
2295
2296	filename = domain.datadir+os.sep+domain.filename+'.sww'
2297	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2298	#points, vertices, boundary = rectangular(15,15)
2299	#domain2.boundary = boundary
2300	###################
2301	##NOW TEST IT!!!
2302	###################
2303
2304	#os.remove(domain.filename + '.sww')
2305	os.remove(filename)
2306
2307	bits = ['vertex_coordinates']
2308	for quantity in ['elevation']+domain.quantities_to_be_stored:
2309	bits.append('quantities["%s"].get_integral()'%quantity)
2310	bits.append('get_quantity("%s")'%quantity)
2311
2312	for bit in bits:
2313	#print 'testing that domain.'+bit+' has been restored'
2314	#print bit
2315	#print 'done'
2316	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2317
2318	######################################
2319	#Now evolve them both, just to be sure
2320	######################################x
2321	visualise = False
2322	#visualise = True
2323	domain.visualise = visualise
2324	domain.time = 0.
2325	from time import sleep
2326
2327	final = .1
2328	domain.set_quantity('friction', 0.1)
2329	domain.store = False
2330	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2331
2332
2333	for t in domain.evolve(yieldstep = yiel, finaltime = final):
2334	if visualise: sleep(1.)
2335	#domain.write_time()
2336	pass
2337
2338	final = final - (domain2.starttime-domain.starttime)
2339	#BUT since domain1 gets time hacked back to 0:
2340	final = final + (domain2.starttime-domain.starttime)
2341
2342	domain2.smooth = False
2343	domain2.visualise = visualise
2344	domain2.store = False
2345	domain2.default_order=2
2346	domain2.set_quantity('friction', 0.1)
2347	#Bed-slope and friction
2348	# Boundary conditions
2349	Bd2=Dirichlet_boundary([0.2,0.,0.])
2350	domain2.boundary = domain.boundary
2351	#print 'domain2.boundary'
2352	#print domain2.boundary
2353	domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2354	#domain2.set_boundary({'exterior': Bd})
2355
2356	domain2.check_integrity()
2357
2358	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2359	if visualise: sleep(1.)
2360	#domain2.write_time()
2361	pass
2362
2363	###################
2364	##NOW TEST IT!!!
2365	##################
2366
2367	bits = [ 'vertex_coordinates']
2368
2369	for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2370	bits.append('quantities["%s"].get_integral()'%quantity)
2371	bits.append('get_quantity("%s")'%quantity)
2372
2373	for bit in bits:
2374	#print bit
2375	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2376
2377
2378	def test_sww2domain2(self):
2379	##################################################################
2380	#Same as previous test, but this checks how NaNs are handled.
2381	##################################################################
2382
2383
2384	from mesh_factory import rectangular
2385	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2386	Constant_height, Time_boundary, Transmissive_boundary
2387	from Numeric import array
2388
2389	#Create basic mesh
2390	points, vertices, boundary = rectangular(2,2)
2391
2392	#Create shallow water domain
2393	domain = Domain(points, vertices, boundary)
2394	domain.smooth = False
2395	domain.visualise = False
2396	domain.store = True
2397	domain.set_name('test_file')
2398	domain.set_datadir('.')
2399	domain.default_order=2
2400	domain.quantities_to_be_stored=['stage']
2401
2402	domain.set_quantity('elevation', lambda x,y: -x/3)
2403	domain.set_quantity('friction', 0.1)
2404
2405	from math import sin, pi
2406	Br = Reflective_boundary(domain)
2407	Bt = Transmissive_boundary(domain)
2408	Bd = Dirichlet_boundary([0.2,0.,0.])
2409	Bw = Time_boundary(domain=domain,
2410	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
2411
2412	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2413
2414	h = 0.05
2415	elevation = domain.quantities['elevation'].vertex_values
2416	domain.set_quantity('stage', elevation + h)
2417
2418	domain.check_integrity()
2419
2420	for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2421	pass
2422	#domain.write_time()
2423
2424
2425
2426	##################################
2427	#Import the file as a new domain
2428	##################################
2429	from data_manager import sww2domain
2430	from Numeric import allclose
2431	import os
2432
2433	filename = domain.datadir+os.sep+domain.filename+'.sww'
2434
2435	#Fail because NaNs are present
2436	try:
2437	domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2438	except:
2439	#Now import it, filling NaNs to be 0
2440	filler = 0
2441	domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2442
2443	#Clean up
2444	os.remove(filename)
2445
2446
2447	bits = [ 'geo_reference.get_xllcorner()',
2448	'geo_reference.get_yllcorner()',
2449	'vertex_coordinates']
2450
2451	for quantity in ['elevation']+domain.quantities_to_be_stored:
2452	bits.append('quantities["%s"].get_integral()'%quantity)
2453	bits.append('get_quantity("%s")'%quantity)
2454
2455	for bit in bits:
2456	# print 'testing that domain.'+bit+' has been restored'
2457	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2458
2459	assert max(max(domain2.get_quantity('xmomentum')))==filler
2460	assert min(min(domain2.get_quantity('xmomentum')))==filler
2461	assert max(max(domain2.get_quantity('ymomentum')))==filler
2462	assert min(min(domain2.get_quantity('ymomentum')))==filler
2463
2464
2465
2466	#def test_weed(self):
2467	from data_manager import weed
2468
2469	coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2470	volumes1 = [[0,1,2],[3,4,5]]
2471	boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2472	coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2473
2474	points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2475
2476	assert len(points2)==len(coordinates2)
2477	for i in range(len(coordinates2)):
2478	coordinate = tuple(coordinates2[i])
2479	assert points2.has_key(coordinate)
2480	points2[coordinate]=i
2481
2482	for triangle in volumes1:
2483	for coordinate in triangle:
2484	assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2485	assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2486
2487
2488	#FIXME This fails - smooth makes the comparism too hard for allclose
2489	def ztest_sww2domain3(self):
2490	################################################
2491	#DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2492	################################################
2493	from mesh_factory import rectangular
2494	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2495	Constant_height, Time_boundary, Transmissive_boundary
2496	from Numeric import array
2497	#Create basic mesh
2498
2499	yiel=0.01
2500	points, vertices, boundary = rectangular(10,10)
2501
2502	#Create shallow water domain
2503	domain = Domain(points, vertices, boundary)
2504	domain.geo_reference = Geo_reference(56,11,11)
2505	domain.smooth = True
2506	domain.visualise = False
2507	domain.store = True
2508	domain.filename = 'bedslope'
2509	domain.default_order=2
2510	#Bed-slope and friction
2511	domain.set_quantity('elevation', lambda x,y: -x/3)
2512	domain.set_quantity('friction', 0.1)
2513	# Boundary conditions
2514	from math import sin, pi
2515	Br = Reflective_boundary(domain)
2516	Bt = Transmissive_boundary(domain)
2517	Bd = Dirichlet_boundary([0.2,0.,0.])
2518	Bw = Time_boundary(domain=domain,
2519	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
2520
2521	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2522
2523	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2524	#Initial condition
2525	h = 0.05
2526	elevation = domain.quantities['elevation'].vertex_values
2527	domain.set_quantity('stage', elevation + h)
2528	#elevation = domain.get_quantity('elevation')
2529	#domain.set_quantity('stage', elevation + h)
2530
2531	domain.check_integrity()
2532	#Evolution
2533	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2534	# domain.write_time()
2535	pass
2536
2537
2538	##########################################
2539	#Import the example's file as a new domain
2540	##########################################
2541	from data_manager import sww2domain
2542	from Numeric import allclose
2543	import os
2544
2545	filename = domain.datadir+os.sep+domain.filename+'.sww'
2546	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2547	#points, vertices, boundary = rectangular(15,15)
2548	#domain2.boundary = boundary
2549	###################
2550	##NOW TEST IT!!!
2551	###################
2552
2553	os.remove(domain.filename + '.sww')
2554
2555	#FIXME smooth domain so that they can be compared
2556
2557
2558	bits = []#'vertex_coordinates']
2559	for quantity in ['elevation']+domain.quantities_to_be_stored:
2560	bits.append('quantities["%s"].get_integral()'%quantity)
2561	#bits.append('get_quantity("%s")'%quantity)
2562
2563	for bit in bits:
2564	#print 'testing that domain.'+bit+' has been restored'
2565	#print bit
2566	#print 'done'
2567	#print ('domain.'+bit), eval('domain.'+bit)
2568	#print ('domain2.'+bit), eval('domain2.'+bit)
2569	assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2570	pass
2571
2572	######################################
2573	#Now evolve them both, just to be sure
2574	######################################x
2575	visualise = False
2576	visualise = True
2577	domain.visualise = visualise
2578	domain.time = 0.
2579	from time import sleep
2580
2581	final = .5
2582	domain.set_quantity('friction', 0.1)
2583	domain.store = False
2584	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2585
2586	for t in domain.evolve(yieldstep = yiel, finaltime = final):
2587	if visualise: sleep(.03)
2588	#domain.write_time()
2589	pass
2590
2591	domain2.smooth = True
2592	domain2.visualise = visualise
2593	domain2.store = False
2594	domain2.default_order=2
2595	domain2.set_quantity('friction', 0.1)
2596	#Bed-slope and friction
2597	# Boundary conditions
2598	Bd2=Dirichlet_boundary([0.2,0.,0.])
2599	Br2 = Reflective_boundary(domain2)
2600	domain2.boundary = domain.boundary
2601	#print 'domain2.boundary'
2602	#print domain2.boundary
2603	domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2604	#domain2.boundary = domain.boundary
2605	#domain2.set_boundary({'exterior': Bd})
2606
2607	domain2.check_integrity()
2608
2609	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2610	if visualise: sleep(.03)
2611	#domain2.write_time()
2612	pass
2613
2614	###################
2615	##NOW TEST IT!!!
2616	##################
2617
2618	bits = [ 'vertex_coordinates']
2619
2620	for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2621	#bits.append('quantities["%s"].get_integral()'%quantity)
2622	bits.append('get_quantity("%s")'%quantity)
2623
2624	for bit in bits:
2625	print bit
2626	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2627
2628
2629	def test_decimate_dem(self):
2630	"""Test decimation of dem file
2631	"""
2632
2633	import os
2634	from Numeric import ones, allclose, Float, arange
2635	from Scientific.IO.NetCDF import NetCDFFile
2636
2637	#Write test dem file
2638	root = 'decdemtest'
2639
2640	filename = root + '.dem'
2641	fid = NetCDFFile(filename, 'w')
2642
2643	fid.institution = 'Geoscience Australia'
2644	fid.description = 'NetCDF DEM format for compact and portable ' +\
2645	'storage of spatial point data'
2646
2647	nrows = 15
2648	ncols = 18
2649
2650	fid.ncols = ncols
2651	fid.nrows = nrows
2652	fid.xllcorner = 2000.5
2653	fid.yllcorner = 3000.5
2654	fid.cellsize = 25
2655	fid.NODATA_value = -9999
2656
2657	fid.zone = 56
2658	fid.false_easting = 0.0
2659	fid.false_northing = 0.0
2660	fid.projection = 'UTM'
2661	fid.datum = 'WGS84'
2662	fid.units = 'METERS'
2663
2664	fid.createDimension('number_of_points', nrows*ncols)
2665
2666	fid.createVariable('elevation', Float, ('number_of_points',))
2667
2668	elevation = fid.variables['elevation']
2669
2670	elevation[:] = (arange(nrows*ncols))
2671
2672	fid.close()
2673
2674	#generate the elevation values expected in the decimated file
2675	ref_elevation = [( 0+ 1+ 2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2676	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2677	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2678	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2679	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2680	( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2681	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2682	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2683	(144+145+146+162+163+164+180+181+182) / 9.0,
2684	(148+149+150+166+167+168+184+185+186) / 9.0,
2685	(152+153+154+170+171+172+188+189+190) / 9.0,
2686	(156+157+158+174+175+176+192+193+194) / 9.0,
2687	(216+217+218+234+235+236+252+253+254) / 9.0,
2688	(220+221+222+238+239+240+256+257+258) / 9.0,
2689	(224+225+226+242+243+244+260+261+262) / 9.0,
2690	(228+229+230+246+247+248+264+265+266) / 9.0]
2691
2692	#generate a stencil for computing the decimated values
2693	stencil = ones((3,3), Float) / 9.0
2694
2695	decimate_dem(root, stencil=stencil, cellsize_new=100)
2696
2697	#Open decimated NetCDF file
2698	fid = NetCDFFile(root + '_100.dem', 'r')
2699
2700	# Get decimated elevation
2701	elevation = fid.variables['elevation']
2702
2703	#Check values
2704	assert allclose(elevation, ref_elevation)
2705
2706	#Cleanup
2707	fid.close()
2708
2709	os.remove(root + '.dem')
2710	os.remove(root + '_100.dem')
2711
2712	def test_decimate_dem_NODATA(self):
2713	"""Test decimation of dem file that includes NODATA values
2714	"""
2715
2716	import os
2717	from Numeric import ones, allclose, Float, arange, reshape
2718	from Scientific.IO.NetCDF import NetCDFFile
2719
2720	#Write test dem file
2721	root = 'decdemtest'
2722
2723	filename = root + '.dem'
2724	fid = NetCDFFile(filename, 'w')
2725
2726	fid.institution = 'Geoscience Australia'
2727	fid.description = 'NetCDF DEM format for compact and portable ' +\
2728	'storage of spatial point data'
2729
2730	nrows = 15
2731	ncols = 18
2732	NODATA_value = -9999
2733
2734	fid.ncols = ncols
2735	fid.nrows = nrows
2736	fid.xllcorner = 2000.5
2737	fid.yllcorner = 3000.5
2738	fid.cellsize = 25
2739	fid.NODATA_value = NODATA_value
2740
2741	fid.zone = 56
2742	fid.false_easting = 0.0
2743	fid.false_northing = 0.0
2744	fid.projection = 'UTM'
2745	fid.datum = 'WGS84'
2746	fid.units = 'METERS'
2747
2748	fid.createDimension('number_of_points', nrows*ncols)
2749
2750	fid.createVariable('elevation', Float, ('number_of_points',))
2751
2752	elevation = fid.variables['elevation']
2753
2754	#generate initial elevation values
2755	elevation_tmp = (arange(nrows*ncols))
2756	#add some NODATA values
2757	elevation_tmp[0] = NODATA_value
2758	elevation_tmp[95] = NODATA_value
2759	elevation_tmp[188] = NODATA_value
2760	elevation_tmp[189] = NODATA_value
2761	elevation_tmp[190] = NODATA_value
2762	elevation_tmp[209] = NODATA_value
2763	elevation_tmp[252] = NODATA_value
2764
2765	elevation[:] = elevation_tmp
2766
2767	fid.close()
2768
2769	#generate the elevation values expected in the decimated file
2770	ref_elevation = [NODATA_value,
2771	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2772	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2773	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2774	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2775	NODATA_value,
2776	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2777	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2778	(144+145+146+162+163+164+180+181+182) / 9.0,
2779	(148+149+150+166+167+168+184+185+186) / 9.0,
2780	NODATA_value,
2781	(156+157+158+174+175+176+192+193+194) / 9.0,
2782	NODATA_value,
2783	(220+221+222+238+239+240+256+257+258) / 9.0,
2784	(224+225+226+242+243+244+260+261+262) / 9.0,
2785	(228+229+230+246+247+248+264+265+266) / 9.0]
2786
2787	#generate a stencil for computing the decimated values
2788	stencil = ones((3,3), Float) / 9.0
2789
2790	decimate_dem(root, stencil=stencil, cellsize_new=100)
2791
2792	#Open decimated NetCDF file
2793	fid = NetCDFFile(root + '_100.dem', 'r')
2794
2795	# Get decimated elevation
2796	elevation = fid.variables['elevation']
2797
2798	#Check values
2799	assert allclose(elevation, ref_elevation)
2800
2801	#Cleanup
2802	fid.close()
2803
2804	os.remove(root + '.dem')
2805	os.remove(root + '_100.dem')
2806
2807	def xxxtestz_sww2ers_real(self):
2808	"""Test that sww information can be converted correctly to asc/prj
2809	format readable by e.g. ArcView
2810	"""
2811
2812	import time, os
2813	from Numeric import array, zeros, allclose, Float, concatenate
2814	from Scientific.IO.NetCDF import NetCDFFile
2815
2816	# the memory optimised least squares
2817	# cellsize = 20, # this one seems to hang
2818	# cellsize = 200000, # Ran 1 test in 269.703s
2819	#Ran 1 test in 267.344s
2820	# cellsize = 20000, # Ran 1 test in 460.922s
2821	# cellsize = 2000 #Ran 1 test in 5340.250s
2822	# cellsize = 200 #this one seems to hang, building matirx A
2823
2824	# not optimised
2825	# seems to hang
2826	# cellsize = 2000 # Ran 1 test in 5334.563s
2827	#Export to ascii/prj files
2828	sww2dem('karratha_100m',
2829	quantity = 'depth',
2830	cellsize = 200000,
2831	verbose = True)
2832
2833	def test_read_asc(self):
2834	"""Test conversion from dem in ascii format to native NetCDF xya format
2835	"""
2836
2837	import time, os
2838	from Numeric import array, zeros, allclose, Float, concatenate
2839	from Scientific.IO.NetCDF import NetCDFFile
2840
2841	import data_manager
2842	#Write test asc file
2843	filename = tempfile.mktemp(".000")
2844	fid = open(filename, 'w')
2845	fid.write("""ncols 7
2846	nrows 4
2847	xllcorner 2000.5
2848	yllcorner 3000.5
2849	cellsize 25
2850	NODATA_value -9999
2851	97.921 99.285 125.588 180.830 258.645 342.872 415.836
2852	473.157 514.391 553.893 607.120 678.125 777.283 883.038
2853	984.494 1040.349 1008.161 900.738 730.882 581.430 514.980
2854	502.645 516.230 504.739 450.604 388.500 338.097 514.980
2855	""")
2856	fid.close()
2857	bath_metadata, grid = \
2858	data_manager._read_asc(filename, verbose=False)
2859	self.failUnless(bath_metadata['xllcorner'] == 2000.5, 'Failed')
2860	self.failUnless(bath_metadata['yllcorner'] == 3000.5, 'Failed')
2861	self.failUnless(bath_metadata['cellsize'] == 25, 'Failed')
2862	self.failUnless(bath_metadata['NODATA_value'] == -9999, 'Failed')
2863	self.failUnless(grid[0][0] == 97.921, 'Failed')
2864	self.failUnless(grid[3][6] == 514.980, 'Failed')
2865
2866	os.remove(filename)
2867
2868	def test_asc_csiro2sww(self):
2869	import tempfile
2870
2871	bath_dir = tempfile.mkdtemp()
2872	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
2873	#bath_dir = 'bath_data_manager_test'
2874	#print "os.getcwd( )",os.getcwd( )
2875	elevation_dir = tempfile.mkdtemp()
2876	#elevation_dir = 'elev_expanded'
2877	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
2878	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
2879
2880	fid = open(bath_dir_filename, 'w')
2881	fid.write(""" ncols 3
2882	nrows 2
2883	xllcorner 148.00000
2884	yllcorner -38.00000
2885	cellsize 0.25
2886	nodata_value -9999.0
2887	9000.000 -1000.000 3000.0
2888	-1000.000 9000.000 -1000.000
2889	""")
2890	fid.close()
2891
2892	fid = open(elevation_dir_filename1, 'w')
2893	fid.write(""" ncols 3
2894	nrows 2
2895	xllcorner 148.00000
2896	yllcorner -38.00000
2897	cellsize 0.25
2898	nodata_value -9999.0
2899	9000.000 0.000 3000.0
2900	0.000 9000.000 0.000
2901	""")
2902	fid.close()
2903
2904	fid = open(elevation_dir_filename2, 'w')
2905	fid.write(""" ncols 3
2906	nrows 2
2907	xllcorner 148.00000
2908	yllcorner -38.00000
2909	cellsize 0.25
2910	nodata_value -9999.0
2911	9000.000 4000.000 4000.0
2912	4000.000 9000.000 4000.000
2913	""")
2914	fid.close()
2915
2916	ucur_dir = tempfile.mkdtemp()
2917	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
2918	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
2919
2920	fid = open(ucur_dir_filename1, 'w')
2921	fid.write(""" ncols 3
2922	nrows 2
2923	xllcorner 148.00000
2924	yllcorner -38.00000
2925	cellsize 0.25
2926	nodata_value -9999.0
2927	90.000 60.000 30.0
2928	10.000 10.000 10.000
2929	""")
2930	fid.close()
2931	fid = open(ucur_dir_filename2, 'w')
2932	fid.write(""" ncols 3
2933	nrows 2
2934	xllcorner 148.00000
2935	yllcorner -38.00000
2936	cellsize 0.25
2937	nodata_value -9999.0
2938	90.000 60.000 30.0
2939	10.000 10.000 10.000
2940	""")
2941	fid.close()
2942
2943	vcur_dir = tempfile.mkdtemp()
2944	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
2945	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
2946
2947	fid = open(vcur_dir_filename1, 'w')
2948	fid.write(""" ncols 3
2949	nrows 2
2950	xllcorner 148.00000
2951	yllcorner -38.00000
2952	cellsize 0.25
2953	nodata_value -9999.0
2954	90.000 60.000 30.0
2955	10.000 10.000 10.000
2956	""")
2957	fid.close()
2958	fid = open(vcur_dir_filename2, 'w')
2959	fid.write(""" ncols 3
2960	nrows 2
2961	xllcorner 148.00000
2962	yllcorner -38.00000
2963	cellsize 0.25
2964	nodata_value -9999.0
2965	90.000 60.000 30.0
2966	10.000 10.000 10.000
2967	""")
2968	fid.close()
2969
2970	sww_file = 'a_test.sww'
2971	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
2972
2973	# check the sww file
2974
2975	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
2976	x = fid.variables['x'][:]
2977	y = fid.variables['y'][:]
2978	z = fid.variables['z'][:]
2979	stage = fid.variables['stage'][:]
2980	xmomentum = fid.variables['xmomentum'][:]
2981	geo_ref = Geo_reference(NetCDFObject=fid)
2982	#print "geo_ref",geo_ref
2983	x_ref = geo_ref.get_xllcorner()
2984	y_ref = geo_ref.get_yllcorner()
2985	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
2986	assert allclose(x_ref, 587798.418) # (-38, 148)
2987	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
2988
2989	#Zone: 55
2990	#Easting: 588095.674 Northing: 5821451.722
2991	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
2992	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
2993
2994	#Zone: 55
2995	#Easting: 632145.632 Northing: 5820863.269
2996	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
2997	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
2998
2999	#Zone: 55
3000	#Easting: 609748.788 Northing: 5793447.860
3001	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
3002	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
3003
3004	assert allclose(z[0],9000.0 )
3005	assert allclose(stage[0][1],0.0 )
3006
3007	#(4000+1000)*60
3008	assert allclose(xmomentum[1][1],300000.0 )
3009
3010
3011	fid.close()
3012
3013	#tidy up
3014	os.remove(bath_dir_filename)
3015	os.rmdir(bath_dir)
3016
3017	os.remove(elevation_dir_filename1)
3018	os.remove(elevation_dir_filename2)
3019	os.rmdir(elevation_dir)
3020
3021	os.remove(ucur_dir_filename1)
3022	os.remove(ucur_dir_filename2)
3023	os.rmdir(ucur_dir)
3024
3025	os.remove(vcur_dir_filename1)
3026	os.remove(vcur_dir_filename2)
3027	os.rmdir(vcur_dir)
3028
3029
3030	# remove sww file
3031	os.remove(sww_file)
3032
3033	def test_asc_csiro2sww2(self):
3034	import tempfile
3035
3036	bath_dir = tempfile.mkdtemp()
3037	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3038	#bath_dir = 'bath_data_manager_test'
3039	#print "os.getcwd( )",os.getcwd( )
3040	elevation_dir = tempfile.mkdtemp()
3041	#elevation_dir = 'elev_expanded'
3042	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3043	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3044
3045	fid = open(bath_dir_filename, 'w')
3046	fid.write(""" ncols 3
3047	nrows 2
3048	xllcorner 148.00000
3049	yllcorner -38.00000
3050	cellsize 0.25
3051	nodata_value -9999.0
3052	9000.000 -1000.000 3000.0
3053	-1000.000 9000.000 -1000.000
3054	""")
3055	fid.close()
3056
3057	fid = open(elevation_dir_filename1, 'w')
3058	fid.write(""" ncols 3
3059	nrows 2
3060	xllcorner 148.00000
3061	yllcorner -38.00000
3062	cellsize 0.25
3063	nodata_value -9999.0
3064	9000.000 0.000 3000.0
3065	0.000 -9999.000 -9999.000
3066	""")
3067	fid.close()
3068
3069	fid = open(elevation_dir_filename2, 'w')
3070	fid.write(""" ncols 3
3071	nrows 2
3072	xllcorner 148.00000
3073	yllcorner -38.00000
3074	cellsize 0.25
3075	nodata_value -9999.0
3076	9000.000 4000.000 4000.0
3077	4000.000 9000.000 4000.000
3078	""")
3079	fid.close()
3080
3081	ucur_dir = tempfile.mkdtemp()
3082	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3083	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3084
3085	fid = open(ucur_dir_filename1, 'w')
3086	fid.write(""" ncols 3
3087	nrows 2
3088	xllcorner 148.00000
3089	yllcorner -38.00000
3090	cellsize 0.25
3091	nodata_value -9999.0
3092	90.000 60.000 30.0
3093	10.000 10.000 10.000
3094	""")
3095	fid.close()
3096	fid = open(ucur_dir_filename2, 'w')
3097	fid.write(""" ncols 3
3098	nrows 2
3099	xllcorner 148.00000
3100	yllcorner -38.00000
3101	cellsize 0.25
3102	nodata_value -9999.0
3103	90.000 60.000 30.0
3104	10.000 10.000 10.000
3105	""")
3106	fid.close()
3107
3108	vcur_dir = tempfile.mkdtemp()
3109	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3110	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3111
3112	fid = open(vcur_dir_filename1, 'w')
3113	fid.write(""" ncols 3
3114	nrows 2
3115	xllcorner 148.00000
3116	yllcorner -38.00000
3117	cellsize 0.25
3118	nodata_value -9999.0
3119	90.000 60.000 30.0
3120	10.000 10.000 10.000
3121	""")
3122	fid.close()
3123	fid = open(vcur_dir_filename2, 'w')
3124	fid.write(""" ncols 3
3125	nrows 2
3126	xllcorner 148.00000
3127	yllcorner -38.00000
3128	cellsize 0.25
3129	nodata_value -9999.0
3130	90.000 60.000 30.0
3131	10.000 10.000 10.000
3132	""")
3133	fid.close()
3134
3135	try:
3136	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3137	vcur_dir, sww_file)
3138	except:
3139	#tidy up
3140	os.remove(bath_dir_filename)
3141	os.rmdir(bath_dir)
3142
3143	os.remove(elevation_dir_filename1)
3144	os.remove(elevation_dir_filename2)
3145	os.rmdir(elevation_dir)
3146
3147	os.remove(ucur_dir_filename1)
3148	os.remove(ucur_dir_filename2)
3149	os.rmdir(ucur_dir)
3150
3151	os.remove(vcur_dir_filename1)
3152	os.remove(vcur_dir_filename2)
3153	os.rmdir(vcur_dir)
3154	else:
3155	#tidy up
3156	os.remove(bath_dir_filename)
3157	os.rmdir(bath_dir)
3158
3159	os.remove(elevation_dir_filename1)
3160	os.remove(elevation_dir_filename2)
3161	os.rmdir(elevation_dir)
3162	raise 'Should raise exception'
3163
3164	os.remove(ucur_dir_filename1)
3165	os.remove(ucur_dir_filename2)
3166	os.rmdir(ucur_dir)
3167
3168	os.remove(vcur_dir_filename1)
3169	os.remove(vcur_dir_filename2)
3170	os.rmdir(vcur_dir)
3171
3172
3173
3174	def test_asc_csiro2sww3(self):
3175	import tempfile
3176
3177	bath_dir = tempfile.mkdtemp()
3178	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3179	#bath_dir = 'bath_data_manager_test'
3180	#print "os.getcwd( )",os.getcwd( )
3181	elevation_dir = tempfile.mkdtemp()
3182	#elevation_dir = 'elev_expanded'
3183	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3184	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3185
3186	fid = open(bath_dir_filename, 'w')
3187	fid.write(""" ncols 3
3188	nrows 2
3189	xllcorner 148.00000
3190	yllcorner -38.00000
3191	cellsize 0.25
3192	nodata_value -9999.0
3193	9000.000 -1000.000 3000.0
3194	-1000.000 9000.000 -1000.000
3195	""")
3196	fid.close()
3197
3198	fid = open(elevation_dir_filename1, 'w')
3199	fid.write(""" ncols 3
3200	nrows 2
3201	xllcorner 148.00000
3202	yllcorner -38.00000
3203	cellsize 0.25
3204	nodata_value -9999.0
3205	9000.000 0.000 3000.0
3206	0.000 -9999.000 -9999.000
3207	""")
3208	fid.close()
3209
3210	fid = open(elevation_dir_filename2, 'w')
3211	fid.write(""" ncols 3
3212	nrows 2
3213	xllcorner 148.00000
3214	yllcorner -38.00000
3215	cellsize 0.25
3216	nodata_value -9999.0
3217	9000.000 4000.000 4000.0
3218	4000.000 9000.000 4000.000
3219	""")
3220	fid.close()
3221
3222	ucur_dir = tempfile.mkdtemp()
3223	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3224	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3225
3226	fid = open(ucur_dir_filename1, 'w')
3227	fid.write(""" ncols 3
3228	nrows 2
3229	xllcorner 148.00000
3230	yllcorner -38.00000
3231	cellsize 0.25
3232	nodata_value -9999.0
3233	90.000 60.000 30.0
3234	10.000 10.000 10.000
3235	""")
3236	fid.close()
3237	fid = open(ucur_dir_filename2, 'w')
3238	fid.write(""" ncols 3
3239	nrows 2
3240	xllcorner 148.00000
3241	yllcorner -38.00000
3242	cellsize 0.25
3243	nodata_value -9999.0
3244	90.000 60.000 30.0
3245	10.000 10.000 10.000
3246	""")
3247	fid.close()
3248
3249	vcur_dir = tempfile.mkdtemp()
3250	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3251	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3252
3253	fid = open(vcur_dir_filename1, 'w')
3254	fid.write(""" ncols 3
3255	nrows 2
3256	xllcorner 148.00000
3257	yllcorner -38.00000
3258	cellsize 0.25
3259	nodata_value -9999.0
3260	90.000 60.000 30.0
3261	10.000 10.000 10.000
3262	""")
3263	fid.close()
3264	fid = open(vcur_dir_filename2, 'w')
3265	fid.write(""" ncols 3
3266	nrows 2
3267	xllcorner 148.00000
3268	yllcorner -38.00000
3269	cellsize 0.25
3270	nodata_value -9999.0
3271	90.000 60.000 30.0
3272	10.000 10.000 10.000
3273	""")
3274	fid.close()
3275
3276	sww_file = 'a_test.sww'
3277	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3278	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3279	mean_stage = 100)
3280
3281	# check the sww file
3282
3283	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
3284	x = fid.variables['x'][:]
3285	y = fid.variables['y'][:]
3286	z = fid.variables['z'][:]
3287	stage = fid.variables['stage'][:]
3288	xmomentum = fid.variables['xmomentum'][:]
3289	geo_ref = Geo_reference(NetCDFObject=fid)
3290	#print "geo_ref",geo_ref
3291	x_ref = geo_ref.get_xllcorner()
3292	y_ref = geo_ref.get_yllcorner()
3293	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
3294	assert allclose(x_ref, 587798.418) # (-38, 148)
3295	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3296
3297	#Zone: 55
3298	#Easting: 588095.674 Northing: 5821451.722
3299	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
3300	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3301
3302	#Zone: 55
3303	#Easting: 632145.632 Northing: 5820863.269
3304	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3305	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3306
3307	#Zone: 55
3308	#Easting: 609748.788 Northing: 5793447.860
3309	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
3310	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
3311
3312	assert allclose(z[0],9000.0 )
3313	assert allclose(stage[0][4],100.0 )
3314	assert allclose(stage[0][5],100.0 )
3315
3316	#(100.0 - 9000)*10
3317	assert allclose(xmomentum[0][4], -89000.0 )
3318
3319	#(100.0 - -1000.000)*10
3320	assert allclose(xmomentum[0][5], 11000.0 )
3321
3322	fid.close()
3323
3324	#tidy up
3325	os.remove(bath_dir_filename)
3326	os.rmdir(bath_dir)
3327
3328	os.remove(elevation_dir_filename1)
3329	os.remove(elevation_dir_filename2)
3330	os.rmdir(elevation_dir)
3331
3332	os.remove(ucur_dir_filename1)
3333	os.remove(ucur_dir_filename2)
3334	os.rmdir(ucur_dir)
3335
3336	os.remove(vcur_dir_filename1)
3337	os.remove(vcur_dir_filename2)
3338	os.rmdir(vcur_dir)
3339
3340	# remove sww file
3341	os.remove(sww_file)
3342
3343
3344	def test_asc_csiro2sww4(self):
3345	"""
3346	Test specifying the extent
3347	"""
3348
3349	import tempfile
3350
3351	bath_dir = tempfile.mkdtemp()
3352	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3353	#bath_dir = 'bath_data_manager_test'
3354	#print "os.getcwd( )",os.getcwd( )
3355	elevation_dir = tempfile.mkdtemp()
3356	#elevation_dir = 'elev_expanded'
3357	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3358	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3359
3360	fid = open(bath_dir_filename, 'w')
3361	fid.write(""" ncols 4
3362	nrows 4
3363	xllcorner 148.00000
3364	yllcorner -38.00000
3365	cellsize 0.25
3366	nodata_value -9999.0
3367	-9000.000 -1000.000 -3000.0 -2000.000
3368	-1000.000 9000.000 -1000.000 -3000.000
3369	-4000.000 6000.000 2000.000 -5000.000
3370	-9000.000 -1000.000 -3000.0 -2000.000
3371	""")
3372	fid.close()
3373
3374	fid = open(elevation_dir_filename1, 'w')
3375	fid.write(""" ncols 4
3376	nrows 4
3377	xllcorner 148.00000
3378	yllcorner -38.00000
3379	cellsize 0.25
3380	nodata_value -9999.0
3381	-900.000 -100.000 -300.0 -200.000
3382	-100.000 900.000 -100.000 -300.000
3383	-400.000 600.000 200.000 -500.000
3384	-900.000 -100.000 -300.0 -200.000
3385	""")
3386	fid.close()
3387
3388	fid = open(elevation_dir_filename2, 'w')
3389	fid.write(""" ncols 4
3390	nrows 4
3391	xllcorner 148.00000
3392	yllcorner -38.00000
3393	cellsize 0.25
3394	nodata_value -9999.0
3395	-990.000 -110.000 -330.0 -220.000
3396	-110.000 990.000 -110.000 -330.000
3397	-440.000 660.000 220.000 -550.000
3398	-990.000 -110.000 -330.0 -220.000
3399	""")
3400	fid.close()
3401
3402	ucur_dir = tempfile.mkdtemp()
3403	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3404	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3405
3406	fid = open(ucur_dir_filename1, 'w')
3407	fid.write(""" ncols 4
3408	nrows 4
3409	xllcorner 148.00000
3410	yllcorner -38.00000
3411	cellsize 0.25
3412	nodata_value -9999.0
3413	-90.000 -10.000 -30.0 -20.000
3414	-10.000 90.000 -10.000 -30.000
3415	-40.000 60.000 20.000 -50.000
3416	-90.000 -10.000 -30.0 -20.000
3417	""")
3418	fid.close()
3419	fid = open(ucur_dir_filename2, 'w')
3420	fid.write(""" ncols 4
3421	nrows 4
3422	xllcorner 148.00000
3423	yllcorner -38.00000
3424	cellsize 0.25
3425	nodata_value -9999.0
3426	-90.000 -10.000 -30.0 -20.000
3427	-10.000 99.000 -11.000 -30.000
3428	-40.000 66.000 22.000 -50.000
3429	-90.000 -10.000 -30.0 -20.000
3430	""")
3431	fid.close()
3432
3433	vcur_dir = tempfile.mkdtemp()
3434	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3435	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3436
3437	fid = open(vcur_dir_filename1, 'w')
3438	fid.write(""" ncols 4
3439	nrows 4
3440	xllcorner 148.00000
3441	yllcorner -38.00000
3442	cellsize 0.25
3443	nodata_value -9999.0
3444	-90.000 -10.000 -30.0 -20.000
3445	-10.000 80.000 -20.000 -30.000
3446	-40.000 50.000 10.000 -50.000
3447	-90.000 -10.000 -30.0 -20.000
3448	""")
3449	fid.close()
3450	fid = open(vcur_dir_filename2, 'w')
3451	fid.write(""" ncols 4
3452	nrows 4
3453	xllcorner 148.00000
3454	yllcorner -38.00000
3455	cellsize 0.25
3456	nodata_value -9999.0
3457	-90.000 -10.000 -30.0 -20.000
3458	-10.000 88.000 -22.000 -30.000
3459	-40.000 55.000 11.000 -50.000
3460	-90.000 -10.000 -30.0 -20.000
3461	""")
3462	fid.close()
3463
3464	sww_file = tempfile.mktemp(".sww")
3465	#sww_file = 'a_test.sww'
3466	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3467	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3468	mean_stage = 100,
3469	minlat = -37.6, maxlat = -37.6,
3470	minlon = 148.3, maxlon = 148.3
3471	#,verbose = True
3472	)
3473
3474	# check the sww file
3475
3476	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
3477	x = fid.variables['x'][:]
3478	y = fid.variables['y'][:]
3479	z = fid.variables['z'][:]
3480	stage = fid.variables['stage'][:]
3481	xmomentum = fid.variables['xmomentum'][:]
3482	ymomentum = fid.variables['ymomentum'][:]
3483	geo_ref = Geo_reference(NetCDFObject=fid)
3484	#print "geo_ref",geo_ref
3485	x_ref = geo_ref.get_xllcorner()
3486	y_ref = geo_ref.get_yllcorner()
3487	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
3488
3489	assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3490	assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3491	assert allclose(y_ref, 5820863.269 )# (-37.45, 148.5)
3492
3493	#Easting: 632145.632 Northing: 5820863.269
3494	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3495
3496	#print "x",x
3497	#print "y",y
3498	self.failUnless(len(x) == 4,'failed') # 2*2
3499	self.failUnless(len(x) == 4,'failed') # 2*2
3500
3501	#Zone: 55
3502	#Easting: 632145.632 Northing: 5820863.269
3503	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3504	# magic number - y is close enough for me.
3505	assert allclose(x[3], 632145.63 - x_ref)
3506	assert allclose(y[3], 5820863.269 - y_ref + 5.22155314684e-005)
3507
3508	assert allclose(z[0],9000.0 ) #z is elevation info
3509	#print "z",z
3510	# 2 time steps, 4 points
3511	self.failUnless(xmomentum.shape == (2,4), 'failed')
3512	self.failUnless(ymomentum.shape == (2,4), 'failed')
3513
3514	#(100.0 - -1000.000)*10
3515	#assert allclose(xmomentum[0][5], 11000.0 )
3516
3517	fid.close()
3518
3519	# is the sww file readable?
3520	#Lets see if we can convert it to a dem!
3521	#print "sww_file",sww_file
3522	#dem_file = tempfile.mktemp(".dem")
3523	domain = sww2domain(sww_file) ###, dem_file)
3524	domain.check_integrity()
3525
3526	#tidy up
3527	os.remove(bath_dir_filename)
3528	os.rmdir(bath_dir)
3529
3530	os.remove(elevation_dir_filename1)
3531	os.remove(elevation_dir_filename2)
3532	os.rmdir(elevation_dir)
3533
3534	os.remove(ucur_dir_filename1)
3535	os.remove(ucur_dir_filename2)
3536	os.rmdir(ucur_dir)
3537
3538	os.remove(vcur_dir_filename1)
3539	os.remove(vcur_dir_filename2)
3540	os.rmdir(vcur_dir)
3541
3542
3543
3544
3545	# remove sww file
3546	os.remove(sww_file)
3547
3548	# remove dem file
3549	#os.remove(dem_file)
3550
3551	def test_get_min_max_indexes(self):
3552	latitudes = [3,2,1,0]
3553	longitudes = [0,10,20,30]
3554
3555	# k - lat
3556	# l - lon
3557	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3558	latitudes,longitudes,
3559	-10,4,-10,31)
3560
3561	#print "kmin",kmin;print "kmax",kmax
3562	#print "lmin",lmin;print "lmax",lmax
3563	latitudes_new = latitudes[kmin:kmax]
3564	longitudes_news = longitudes[lmin:lmax]
3565	#print "latitudes_new", latitudes_new
3566	#print "longitudes_news",longitudes_news
3567	self.failUnless(latitudes == latitudes_new and \
3568	longitudes == longitudes_news,
3569	'failed')
3570
3571	## 2nd test
3572	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3573	latitudes,longitudes,
3574	0.5,2.5,5,25)
3575	#print "kmin",kmin;print "kmax",kmax
3576	#print "lmin",lmin;print "lmax",lmax
3577	latitudes_new = latitudes[kmin:kmax]
3578	longitudes_news = longitudes[lmin:lmax]
3579	#print "latitudes_new", latitudes_new
3580	#print "longitudes_news",longitudes_news
3581
3582	self.failUnless(latitudes == latitudes_new and \
3583	longitudes == longitudes_news,
3584	'failed')
3585
3586	## 3rd test
3587	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3588	longitudes,
3589	1.1,1.9,12,17)
3590	#print "kmin",kmin;print "kmax",kmax
3591	#print "lmin",lmin;print "lmax",lmax
3592	latitudes_new = latitudes[kmin:kmax]
3593	longitudes_news = longitudes[lmin:lmax]
3594	#print "latitudes_new", latitudes_new
3595	#print "longitudes_news",longitudes_news
3596
3597	self.failUnless(latitudes_new == [2, 1] and \
3598	longitudes_news == [10, 20],
3599	'failed')
3600
3601
3602	## 4th test
3603	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3604	latitudes,longitudes,
3605	-0.1,1.9,-2,17)
3606	#print "kmin",kmin;print "kmax",kmax
3607	#print "lmin",lmin;print "lmax",lmax
3608	latitudes_new = latitudes[kmin:kmax]
3609	longitudes_news = longitudes[lmin:lmax]
3610	#print "latitudes_new", latitudes_new
3611	#print "longitudes_news",longitudes_news
3612
3613	self.failUnless(latitudes_new == [2, 1, 0] and \
3614	longitudes_news == [0, 10, 20],
3615	'failed')
3616	## 5th test
3617	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3618	latitudes,longitudes,
3619	0.1,1.9,2,17)
3620	#print "kmin",kmin;print "kmax",kmax
3621	#print "lmin",lmin;print "lmax",lmax
3622	latitudes_new = latitudes[kmin:kmax]
3623	longitudes_news = longitudes[lmin:lmax]
3624	#print "latitudes_new", latitudes_new
3625	#print "longitudes_news",longitudes_news
3626
3627	self.failUnless(latitudes_new == [2, 1, 0] and \
3628	longitudes_news == [0, 10, 20],
3629	'failed')
3630
3631	## 6th test
3632
3633	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3634	latitudes,longitudes,
3635	1.5,4,18,32)
3636	#print "kmin",kmin;print "kmax",kmax
3637	#print "lmin",lmin;print "lmax",lmax
3638	latitudes_new = latitudes[kmin:kmax]
3639	longitudes_news = longitudes[lmin:lmax]
3640	#print "latitudes_new", latitudes_new
3641	#print "longitudes_news",longitudes_news
3642
3643	self.failUnless(latitudes_new == [3, 2, 1] and \
3644	longitudes_news == [10, 20, 30],
3645	'failed')
3646
3647
3648	## 7th test
3649	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3650	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3651	latitudes,longitudes,
3652	1.5,1.5,15,15)
3653	#print "kmin",kmin;print "kmax",kmax
3654	#print "lmin",lmin;print "lmax",lmax
3655	latitudes_new = latitudes[kmin:kmax]
3656	longitudes_news = longitudes[lmin:lmax]
3657	m2d = m2d[kmin:kmax,lmin:lmax]
3658	#print "m2d", m2d
3659	#print "latitudes_new", latitudes_new
3660	#print "longitudes_news",longitudes_news
3661
3662	self.failUnless(latitudes_new == [2, 1] and \
3663	longitudes_news == [10, 20],
3664	'failed')
3665
3666	self.failUnless(m2d == [[5,6],[9,10]],
3667	'failed')
3668
3669	def test_get_min_max_indexes2(self):
3670	latitudes = [-30,-35,-40,-45]
3671	longitudes = [148,149,150,151]
3672
3673	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3674
3675	# k - lat
3676	# l - lon
3677	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3678	latitudes,longitudes,
3679	-37,-27,147,149.5)
3680
3681	#print "kmin",kmin;print "kmax",kmax
3682	#print "lmin",lmin;print "lmax",lmax
3683	#print "m2d", m2d
3684	#print "latitudes", latitudes
3685	#print "longitudes",longitudes
3686	#print "latitudes[kmax]", latitudes[kmax]
3687	latitudes_new = latitudes[kmin:kmax]
3688	longitudes_new = longitudes[lmin:lmax]
3689	m2d = m2d[kmin:kmax,lmin:lmax]
3690	#print "m2d", m2d
3691	#print "latitudes_new", latitudes_new
3692	#print "longitudes_new",longitudes_new
3693
3694	self.failUnless(latitudes_new == [-30, -35, -40] and \
3695	longitudes_new == [148, 149,150],
3696	'failed')
3697	self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
3698	'failed')
3699
3700	def test_get_min_max_indexes3(self):
3701	latitudes = [-30,-35,-40,-45,-50,-55,-60]
3702	longitudes = [148,149,150,151]
3703
3704	# k - lat
3705	# l - lon
3706	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3707	latitudes,longitudes,
3708	-43,-37,148.5,149.5)
3709
3710
3711	#print "kmin",kmin;print "kmax",kmax
3712	#print "lmin",lmin;print "lmax",lmax
3713	#print "latitudes", latitudes
3714	#print "longitudes",longitudes
3715	latitudes_new = latitudes[kmin:kmax]
3716	longitudes_news = longitudes[lmin:lmax]
3717	#print "latitudes_new", latitudes_new
3718	#print "longitudes_news",longitudes_news
3719
3720	self.failUnless(latitudes_new == [-35, -40, -45] and \
3721	longitudes_news == [148, 149,150],
3722	'failed')
3723
3724	def test_get_min_max_indexes4(self):
3725	latitudes = [-30,-35,-40,-45,-50,-55,-60]
3726	longitudes = [148,149,150,151]
3727
3728	# k - lat
3729	# l - lon
3730	kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3731	latitudes,longitudes)
3732
3733	#print "kmin",kmin;print "kmax",kmax
3734	#print "lmin",lmin;print "lmax",lmax
3735	#print "latitudes", latitudes
3736	#print "longitudes",longitudes
3737	latitudes_new = latitudes[kmin:kmax]
3738	longitudes_news = longitudes[lmin:lmax]
3739	#print "latitudes_new", latitudes_new
3740	#print "longitudes_news",longitudes_news
3741
3742	self.failUnless(latitudes_new == latitudes and \
3743	longitudes_news == longitudes,
3744	'failed')
3745
3746	def test_tsh2sww(self):
3747	import os
3748	import tempfile
3749
3750	tsh_file = tempfile.mktemp(".tsh")
3751	file = open(tsh_file,"w")
3752	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
3753	0 0.0 0.0 0.0 0.0 0.01 \n \
3754	1 1.0 0.0 10.0 10.0 0.02 \n \
3755	2 0.0 1.0 0.0 10.0 0.03 \n \
3756	3 0.5 0.25 8.0 12.0 0.04 \n \
3757	# Vert att title \n \
3758	elevation \n \
3759	stage \n \
3760	friction \n \
3761	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
3762	0 0 3 2 -1 -1 1 dsg\n\
3763	1 0 1 3 -1 0 -1 ole nielsen\n\
3764	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
3765	0 1 0 2 \n\
3766	1 0 2 3 \n\
3767	2 2 3 \n\
3768	3 3 1 1 \n\
3769	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
3770	0 216.0 -86.0 \n \
3771	1 160.0 -167.0 \n \
3772	2 114.0 -91.0 \n \
3773	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
3774	0 0 1 0 \n \
3775	1 1 2 0 \n \
3776	2 2 0 0 \n \
3777	0 # <x> <y> ...Mesh Holes... \n \
3778	0 # <x> <y> <attribute>...Mesh Regions... \n \
3779	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
3780	#Geo reference \n \
3781	56 \n \
3782	140 \n \
3783	120 \n")
3784	file.close()
3785
3786	#sww_file = tempfile.mktemp(".sww")
3787	#print "sww_file",sww_file
3788	#print "sww_file",tsh_file
3789	tsh2sww(tsh_file)
3790
3791	os.remove(tsh_file)
3792	os.remove(tsh_file[:-4] + '.sww')
3793	#-------------------------------------------------------------
3794	if __name__ == "__main__":
3795	#suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
3796	suite = unittest.makeSuite(Test_Data_Manager,'test')
3797	#suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
3798	#suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')
3799	#suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
3800	#suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
3801	#suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
3802	runner = unittest.TextTestRunner()
3803	runner.run(suite)
3804

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

Download in other formats: