Context Navigation

source: inundation/pyvolution/test_data_manager.py @ 3354

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

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

Download in other formats: