Context Navigation

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

Last change on this file since 3694 was 3694, checked in by duncan, 17 years ago
In ferret2sww set default of inverted_bathymetry to True
File size: 153.7 KB

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

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

Download in other formats: