Context Navigation

source: anuga_core/source/anuga/pyvolution/test_data_manager.py @ 3546

Last change on this file since 3546 was 3532, checked in by duncan, 19 years ago
name change
File size: 154.3 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.pyvolution.data_manager import *
13	from anuga.pyvolution.shallow_water import *
14	from anuga.config import epsilon
15
16	# This is needed to run the tests of local functions
17	import anuga.pyvolution.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_sww_depth(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_sww_depth = 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))
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
2609	os.remove('test_va.nc')
2610	os.remove('test_ua.nc')
2611	os.remove('test_ha.nc')
2612	os.remove('test_e.nc')
2613
2614	#Read output file 'test.sww'
2615	fid = NetCDFFile('test.sww')
2616
2617
2618	#Check first value
2619	elevation = fid.variables['elevation'][:]
2620	stage = fid.variables['stage'][:]
2621	xmomentum = fid.variables['xmomentum'][:]
2622	ymomentum = fid.variables['ymomentum'][:]
2623
2624	#print ymomentum
2625	first_height = first_amp/100 - first_elevation
2626	third_height = third_amp/100 - third_elevation
2627	first_momentum=first_speed*first_height/100
2628	third_momentum=third_speed*third_height/100
2629
2630	assert allclose(ymomentum[0][0],first_momentum) #Meters
2631	assert allclose(ymomentum[0][2],third_momentum) #Meters
2632
2633	fid.close()
2634
2635	#Cleanup
2636	os.remove('test.sww')
2637
2638
2639
2640
2641	def test_ferret2sww_nz_origin(self):
2642	from Scientific.IO.NetCDF import NetCDFFile
2643	from anuga.coordinate_transforms.redfearn import redfearn
2644
2645	#Call conversion (with nonzero origin)
2646	ferret2sww(self.test_MOST_file, verbose=False,
2647	origin = (56, 100000, 200000))
2648
2649
2650	#Work out the UTM coordinates for first point
2651	zone, e, n = redfearn(-34.5, 150.66667)
2652
2653	#Read output file 'small.sww'
2654	#fid = NetCDFFile('small.sww', 'r')
2655	fid = NetCDFFile(self.test_MOST_file + '.sww')
2656
2657	x = fid.variables['x'][:]
2658	y = fid.variables['y'][:]
2659
2660	#Check that first coordinate is correctly represented
2661	assert allclose(x[0], e-100000)
2662	assert allclose(y[0], n-200000)
2663
2664	fid.close()
2665
2666	#Cleanup
2667	os.remove(self.test_MOST_file + '.sww')
2668
2669
2670
2671	def test_sww_extent(self):
2672	"""Not a test, rather a look at the sww format
2673	"""
2674
2675	import time, os
2676	from Numeric import array, zeros, allclose, Float, concatenate
2677	from Scientific.IO.NetCDF import NetCDFFile
2678
2679	self.domain.filename = 'datatest' + str(id(self))
2680	self.domain.format = 'sww'
2681	self.domain.smooth = True
2682	self.domain.reduction = mean
2683	self.domain.set_datadir('.')
2684
2685
2686	sww = get_dataobject(self.domain)
2687	sww.store_connectivity()
2688	sww.store_timestep('stage')
2689	self.domain.time = 2.
2690
2691	#Modify stage at second timestep
2692	stage = self.domain.quantities['stage'].vertex_values
2693	self.domain.set_quantity('stage', stage/2)
2694
2695	sww.store_timestep('stage')
2696
2697	file_and_extension_name = self.domain.filename + ".sww"
2698	#print "file_and_extension_name",file_and_extension_name
2699	[xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2700	extent_sww(file_and_extension_name )
2701
2702	assert allclose(xmin, 0.0)
2703	assert allclose(xmax, 1.0)
2704	assert allclose(ymin, 0.0)
2705	assert allclose(ymax, 1.0)
2706	assert allclose(stagemin, -0.85)
2707	assert allclose(stagemax, 0.15)
2708
2709
2710	#Cleanup
2711	os.remove(sww.filename)
2712
2713
2714
2715	def test_sww2domain1(self):
2716	################################################
2717	#Create a test domain, and evolve and save it.
2718	################################################
2719	from mesh_factory import rectangular
2720	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2721	Constant_height, Time_boundary, Transmissive_boundary
2722	from Numeric import array
2723
2724	#Create basic mesh
2725
2726	yiel=0.01
2727	points, vertices, boundary = rectangular(10,10)
2728
2729	#Create shallow water domain
2730	domain = Domain(points, vertices, boundary)
2731	domain.geo_reference = Geo_reference(56,11,11)
2732	domain.smooth = False
2733	domain.visualise = False
2734	domain.store = True
2735	domain.filename = 'bedslope'
2736	domain.default_order=2
2737	#Bed-slope and friction
2738	domain.set_quantity('elevation', lambda x,y: -x/3)
2739	domain.set_quantity('friction', 0.1)
2740	# Boundary conditions
2741	from math import sin, pi
2742	Br = Reflective_boundary(domain)
2743	Bt = Transmissive_boundary(domain)
2744	Bd = Dirichlet_boundary([0.2,0.,0.])
2745	Bw = Time_boundary(domain=domain,f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
2746
2747	#domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2748	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2749
2750	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2751	#Initial condition
2752	h = 0.05
2753	elevation = domain.quantities['elevation'].vertex_values
2754	domain.set_quantity('stage', elevation + h)
2755
2756	domain.check_integrity()
2757	#Evolution
2758	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2759	#domain.write_time()
2760	pass
2761
2762
2763	##########################################
2764	#Import the example's file as a new domain
2765	##########################################
2766	from anuga.pyvolution.data_manager import sww2domain
2767	from Numeric import allclose
2768	import os
2769
2770	filename = domain.datadir+os.sep+domain.filename+'.sww'
2771	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2772	#points, vertices, boundary = rectangular(15,15)
2773	#domain2.boundary = boundary
2774	###################
2775	##NOW TEST IT!!!
2776	###################
2777
2778	#os.remove(domain.filename + '.sww')
2779	os.remove(filename)
2780
2781	bits = ['vertex_coordinates']
2782	for quantity in ['elevation']+domain.quantities_to_be_stored:
2783	bits.append('get_quantity("%s").get_integral()' %quantity)
2784	bits.append('get_quantity("%s").get_values()' %quantity)
2785
2786	for bit in bits:
2787	#print 'testing that domain.'+bit+' has been restored'
2788	#print bit
2789	#print 'done'
2790	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2791
2792	######################################
2793	#Now evolve them both, just to be sure
2794	######################################x
2795	visualise = False
2796	#visualise = True
2797	domain.visualise = visualise
2798	domain.time = 0.
2799	from time import sleep
2800
2801	final = .1
2802	domain.set_quantity('friction', 0.1)
2803	domain.store = False
2804	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2805
2806
2807	for t in domain.evolve(yieldstep = yiel, finaltime = final):
2808	if visualise: sleep(1.)
2809	#domain.write_time()
2810	pass
2811
2812	final = final - (domain2.starttime-domain.starttime)
2813	#BUT since domain1 gets time hacked back to 0:
2814	final = final + (domain2.starttime-domain.starttime)
2815
2816	domain2.smooth = False
2817	domain2.visualise = visualise
2818	domain2.store = False
2819	domain2.default_order=2
2820	domain2.set_quantity('friction', 0.1)
2821	#Bed-slope and friction
2822	# Boundary conditions
2823	Bd2=Dirichlet_boundary([0.2,0.,0.])
2824	domain2.boundary = domain.boundary
2825	#print 'domain2.boundary'
2826	#print domain2.boundary
2827	domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2828	#domain2.set_boundary({'exterior': Bd})
2829
2830	domain2.check_integrity()
2831
2832	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2833	if visualise: sleep(1.)
2834	#domain2.write_time()
2835	pass
2836
2837	###################
2838	##NOW TEST IT!!!
2839	##################
2840
2841	bits = ['vertex_coordinates']
2842
2843	for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
2844	bits.append('get_quantity("%s").get_integral()' %quantity)
2845	bits.append('get_quantity("%s").get_values()' %quantity)
2846
2847	#print bits
2848	for bit in bits:
2849	#print bit
2850	#print eval('domain.'+bit)
2851	#print eval('domain2.'+bit)
2852
2853	#print eval('domain.'+bit+'-domain2.'+bit)
2854	msg = 'Values in the two domains are different for ' + bit
2855	assert allclose(eval('domain.'+bit),eval('domain2.'+bit),
2856	rtol=1.e-5, atol=3.e-8), msg
2857
2858
2859	def test_sww2domain2(self):
2860	##################################################################
2861	#Same as previous test, but this checks how NaNs are handled.
2862	##################################################################
2863
2864
2865	from mesh_factory import rectangular
2866	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2867	Constant_height, Time_boundary, Transmissive_boundary
2868	from Numeric import array
2869
2870	#Create basic mesh
2871	points, vertices, boundary = rectangular(2,2)
2872
2873	#Create shallow water domain
2874	domain = Domain(points, vertices, boundary)
2875	domain.smooth = False
2876	domain.visualise = False
2877	domain.store = True
2878	domain.set_name('test_file')
2879	domain.set_datadir('.')
2880	domain.default_order=2
2881	domain.quantities_to_be_stored=['stage']
2882
2883	domain.set_quantity('elevation', lambda x,y: -x/3)
2884	domain.set_quantity('friction', 0.1)
2885
2886	from math import sin, pi
2887	Br = Reflective_boundary(domain)
2888	Bt = Transmissive_boundary(domain)
2889	Bd = Dirichlet_boundary([0.2,0.,0.])
2890	Bw = Time_boundary(domain=domain,
2891	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
2892
2893	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2894
2895	h = 0.05
2896	elevation = domain.quantities['elevation'].vertex_values
2897	domain.set_quantity('stage', elevation + h)
2898
2899	domain.check_integrity()
2900
2901	for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2902	pass
2903	#domain.write_time()
2904
2905
2906
2907	##################################
2908	#Import the file as a new domain
2909	##################################
2910	from anuga.pyvolution.data_manager import sww2domain
2911	from Numeric import allclose
2912	import os
2913
2914	filename = domain.datadir+os.sep+domain.filename+'.sww'
2915
2916	#Fail because NaNs are present
2917	try:
2918	domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2919	except:
2920	#Now import it, filling NaNs to be 0
2921	filler = 0
2922	domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2923
2924	#Clean up
2925	os.remove(filename)
2926
2927
2928	bits = [ 'geo_reference.get_xllcorner()',
2929	'geo_reference.get_yllcorner()',
2930	'vertex_coordinates']
2931
2932	for quantity in ['elevation']+domain.quantities_to_be_stored:
2933	bits.append('get_quantity("%s").get_integral()' %quantity)
2934	bits.append('get_quantity("%s").get_values()' %quantity)
2935
2936	for bit in bits:
2937	# print 'testing that domain.'+bit+' has been restored'
2938	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2939
2940	assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
2941	assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
2942	assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
2943	assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
2944
2945
2946
2947	#def test_weed(self):
2948	from anuga.pyvolution.data_manager import weed
2949
2950	coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2951	volumes1 = [[0,1,2],[3,4,5]]
2952	boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2953	coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2954
2955	points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2956
2957	assert len(points2)==len(coordinates2)
2958	for i in range(len(coordinates2)):
2959	coordinate = tuple(coordinates2[i])
2960	assert points2.has_key(coordinate)
2961	points2[coordinate]=i
2962
2963	for triangle in volumes1:
2964	for coordinate in triangle:
2965	assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2966	assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2967
2968
2969	#FIXME This fails - smooth makes the comparism too hard for allclose
2970	def ztest_sww2domain3(self):
2971	################################################
2972	#DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2973	################################################
2974	from mesh_factory import rectangular
2975	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2976	Constant_height, Time_boundary, Transmissive_boundary
2977	from Numeric import array
2978	#Create basic mesh
2979
2980	yiel=0.01
2981	points, vertices, boundary = rectangular(10,10)
2982
2983	#Create shallow water domain
2984	domain = Domain(points, vertices, boundary)
2985	domain.geo_reference = Geo_reference(56,11,11)
2986	domain.smooth = True
2987	domain.visualise = False
2988	domain.store = True
2989	domain.filename = 'bedslope'
2990	domain.default_order=2
2991	#Bed-slope and friction
2992	domain.set_quantity('elevation', lambda x,y: -x/3)
2993	domain.set_quantity('friction', 0.1)
2994	# Boundary conditions
2995	from math import sin, pi
2996	Br = Reflective_boundary(domain)
2997	Bt = Transmissive_boundary(domain)
2998	Bd = Dirichlet_boundary([0.2,0.,0.])
2999	Bw = Time_boundary(domain=domain,
3000	f=lambda t: [(0.1sin(t2*pi)), 0.0, 0.0])
3001
3002	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
3003
3004	domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
3005	#Initial condition
3006	h = 0.05
3007	elevation = domain.quantities['elevation'].vertex_values
3008	domain.set_quantity('stage', elevation + h)
3009
3010
3011	domain.check_integrity()
3012	#Evolution
3013	for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
3014	# domain.write_time()
3015	pass
3016
3017
3018	##########################################
3019	#Import the example's file as a new domain
3020	##########################################
3021	from anuga.pyvolution.data_manager import sww2domain
3022	from Numeric import allclose
3023	import os
3024
3025	filename = domain.datadir+os.sep+domain.filename+'.sww'
3026	domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
3027	#points, vertices, boundary = rectangular(15,15)
3028	#domain2.boundary = boundary
3029	###################
3030	##NOW TEST IT!!!
3031	###################
3032
3033	os.remove(domain.filename + '.sww')
3034
3035	#FIXME smooth domain so that they can be compared
3036
3037
3038	bits = []#'vertex_coordinates']
3039	for quantity in ['elevation']+domain.quantities_to_be_stored:
3040	bits.append('quantities["%s"].get_integral()'%quantity)
3041
3042
3043	for bit in bits:
3044	#print 'testing that domain.'+bit+' has been restored'
3045	#print bit
3046	#print 'done'
3047	#print ('domain.'+bit), eval('domain.'+bit)
3048	#print ('domain2.'+bit), eval('domain2.'+bit)
3049	assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
3050	pass
3051
3052	######################################
3053	#Now evolve them both, just to be sure
3054	######################################x
3055	visualise = False
3056	visualise = True
3057	domain.visualise = visualise
3058	domain.time = 0.
3059	from time import sleep
3060
3061	final = .5
3062	domain.set_quantity('friction', 0.1)
3063	domain.store = False
3064	domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
3065
3066	for t in domain.evolve(yieldstep = yiel, finaltime = final):
3067	if visualise: sleep(.03)
3068	#domain.write_time()
3069	pass
3070
3071	domain2.smooth = True
3072	domain2.visualise = visualise
3073	domain2.store = False
3074	domain2.default_order=2
3075	domain2.set_quantity('friction', 0.1)
3076	#Bed-slope and friction
3077	# Boundary conditions
3078	Bd2=Dirichlet_boundary([0.2,0.,0.])
3079	Br2 = Reflective_boundary(domain2)
3080	domain2.boundary = domain.boundary
3081	#print 'domain2.boundary'
3082	#print domain2.boundary
3083	domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
3084	#domain2.boundary = domain.boundary
3085	#domain2.set_boundary({'exterior': Bd})
3086
3087	domain2.check_integrity()
3088
3089	for t in domain2.evolve(yieldstep = yiel, finaltime = final):
3090	if visualise: sleep(.03)
3091	#domain2.write_time()
3092	pass
3093
3094	###################
3095	##NOW TEST IT!!!
3096	##################
3097
3098	print '><><><><>>'
3099	bits = [ 'vertex_coordinates']
3100
3101	for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
3102	#bits.append('quantities["%s"].get_integral()'%quantity)
3103	bits.append('get_quantity("%s").get_values()' %quantity)
3104
3105	for bit in bits:
3106	print bit
3107	assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
3108
3109
3110	def test_decimate_dem(self):
3111	"""Test decimation of dem file
3112	"""
3113
3114	import os
3115	from Numeric import ones, allclose, Float, arange
3116	from Scientific.IO.NetCDF import NetCDFFile
3117
3118	#Write test dem file
3119	root = 'decdemtest'
3120
3121	filename = root + '.dem'
3122	fid = NetCDFFile(filename, 'w')
3123
3124	fid.institution = 'Geoscience Australia'
3125	fid.description = 'NetCDF DEM format for compact and portable ' +\
3126	'storage of spatial point data'
3127
3128	nrows = 15
3129	ncols = 18
3130
3131	fid.ncols = ncols
3132	fid.nrows = nrows
3133	fid.xllcorner = 2000.5
3134	fid.yllcorner = 3000.5
3135	fid.cellsize = 25
3136	fid.NODATA_value = -9999
3137
3138	fid.zone = 56
3139	fid.false_easting = 0.0
3140	fid.false_northing = 0.0
3141	fid.projection = 'UTM'
3142	fid.datum = 'WGS84'
3143	fid.units = 'METERS'
3144
3145	fid.createDimension('number_of_points', nrows*ncols)
3146
3147	fid.createVariable('elevation', Float, ('number_of_points',))
3148
3149	elevation = fid.variables['elevation']
3150
3151	elevation[:] = (arange(nrows*ncols))
3152
3153	fid.close()
3154
3155	#generate the elevation values expected in the decimated file
3156	ref_elevation = [( 0+ 1+ 2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3157	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3158	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3159	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3160	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3161	( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3162	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3163	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3164	(144+145+146+162+163+164+180+181+182) / 9.0,
3165	(148+149+150+166+167+168+184+185+186) / 9.0,
3166	(152+153+154+170+171+172+188+189+190) / 9.0,
3167	(156+157+158+174+175+176+192+193+194) / 9.0,
3168	(216+217+218+234+235+236+252+253+254) / 9.0,
3169	(220+221+222+238+239+240+256+257+258) / 9.0,
3170	(224+225+226+242+243+244+260+261+262) / 9.0,
3171	(228+229+230+246+247+248+264+265+266) / 9.0]
3172
3173	#generate a stencil for computing the decimated values
3174	stencil = ones((3,3), Float) / 9.0
3175
3176	decimate_dem(root, stencil=stencil, cellsize_new=100)
3177
3178	#Open decimated NetCDF file
3179	fid = NetCDFFile(root + '_100.dem', 'r')
3180
3181	# Get decimated elevation
3182	elevation = fid.variables['elevation']
3183
3184	#Check values
3185	assert allclose(elevation, ref_elevation)
3186
3187	#Cleanup
3188	fid.close()
3189
3190	os.remove(root + '.dem')
3191	os.remove(root + '_100.dem')
3192
3193	def test_decimate_dem_NODATA(self):
3194	"""Test decimation of dem file that includes NODATA values
3195	"""
3196
3197	import os
3198	from Numeric import ones, allclose, Float, arange, reshape
3199	from Scientific.IO.NetCDF import NetCDFFile
3200
3201	#Write test dem file
3202	root = 'decdemtest'
3203
3204	filename = root + '.dem'
3205	fid = NetCDFFile(filename, 'w')
3206
3207	fid.institution = 'Geoscience Australia'
3208	fid.description = 'NetCDF DEM format for compact and portable ' +\
3209	'storage of spatial point data'
3210
3211	nrows = 15
3212	ncols = 18
3213	NODATA_value = -9999
3214
3215	fid.ncols = ncols
3216	fid.nrows = nrows
3217	fid.xllcorner = 2000.5
3218	fid.yllcorner = 3000.5
3219	fid.cellsize = 25
3220	fid.NODATA_value = NODATA_value
3221
3222	fid.zone = 56
3223	fid.false_easting = 0.0
3224	fid.false_northing = 0.0
3225	fid.projection = 'UTM'
3226	fid.datum = 'WGS84'
3227	fid.units = 'METERS'
3228
3229	fid.createDimension('number_of_points', nrows*ncols)
3230
3231	fid.createVariable('elevation', Float, ('number_of_points',))
3232
3233	elevation = fid.variables['elevation']
3234
3235	#generate initial elevation values
3236	elevation_tmp = (arange(nrows*ncols))
3237	#add some NODATA values
3238	elevation_tmp[0] = NODATA_value
3239	elevation_tmp[95] = NODATA_value
3240	elevation_tmp[188] = NODATA_value
3241	elevation_tmp[189] = NODATA_value
3242	elevation_tmp[190] = NODATA_value
3243	elevation_tmp[209] = NODATA_value
3244	elevation_tmp[252] = NODATA_value
3245
3246	elevation[:] = elevation_tmp
3247
3248	fid.close()
3249
3250	#generate the elevation values expected in the decimated file
3251	ref_elevation = [NODATA_value,
3252	( 4+ 5+ 6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3253	( 8+ 9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3254	( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3255	( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3256	NODATA_value,
3257	( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3258	( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3259	(144+145+146+162+163+164+180+181+182) / 9.0,
3260	(148+149+150+166+167+168+184+185+186) / 9.0,
3261	NODATA_value,
3262	(156+157+158+174+175+176+192+193+194) / 9.0,
3263	NODATA_value,
3264	(220+221+222+238+239+240+256+257+258) / 9.0,
3265	(224+225+226+242+243+244+260+261+262) / 9.0,
3266	(228+229+230+246+247+248+264+265+266) / 9.0]
3267
3268	#generate a stencil for computing the decimated values
3269	stencil = ones((3,3), Float) / 9.0
3270
3271	decimate_dem(root, stencil=stencil, cellsize_new=100)
3272
3273	#Open decimated NetCDF file
3274	fid = NetCDFFile(root + '_100.dem', 'r')
3275
3276	# Get decimated elevation
3277	elevation = fid.variables['elevation']
3278
3279	#Check values
3280	assert allclose(elevation, ref_elevation)
3281
3282	#Cleanup
3283	fid.close()
3284
3285	os.remove(root + '.dem')
3286	os.remove(root + '_100.dem')
3287
3288	def xxxtestz_sww2ers_real(self):
3289	"""Test that sww information can be converted correctly to asc/prj
3290	format readable by e.g. ArcView
3291	"""
3292
3293	import time, os
3294	from Numeric import array, zeros, allclose, Float, concatenate
3295	from Scientific.IO.NetCDF import NetCDFFile
3296
3297	# the memory optimised least squares
3298	# cellsize = 20, # this one seems to hang
3299	# cellsize = 200000, # Ran 1 test in 269.703s
3300	#Ran 1 test in 267.344s
3301	# cellsize = 20000, # Ran 1 test in 460.922s
3302	# cellsize = 2000 #Ran 1 test in 5340.250s
3303	# cellsize = 200 #this one seems to hang, building matirx A
3304
3305	# not optimised
3306	# seems to hang
3307	# cellsize = 2000 # Ran 1 test in 5334.563s
3308	#Export to ascii/prj files
3309	sww2dem('karratha_100m',
3310	quantity = 'depth',
3311	cellsize = 200000,
3312	verbose = True)
3313
3314	def test_read_asc(self):
3315	"""Test conversion from dem in ascii format to native NetCDF xya format
3316	"""
3317
3318	import time, os
3319	from Numeric import array, zeros, allclose, Float, concatenate
3320	from Scientific.IO.NetCDF import NetCDFFile
3321
3322	from anuga.pyvolution.data_manager import _read_asc
3323	#Write test asc file
3324	filename = tempfile.mktemp(".000")
3325	fid = open(filename, 'w')
3326	fid.write("""ncols 7
3327	nrows 4
3328	xllcorner 2000.5
3329	yllcorner 3000.5
3330	cellsize 25
3331	NODATA_value -9999
3332	97.921 99.285 125.588 180.830 258.645 342.872 415.836
3333	473.157 514.391 553.893 607.120 678.125 777.283 883.038
3334	984.494 1040.349 1008.161 900.738 730.882 581.430 514.980
3335	502.645 516.230 504.739 450.604 388.500 338.097 514.980
3336	""")
3337	fid.close()
3338	bath_metadata, grid = _read_asc(filename, verbose=False)
3339	self.failUnless(bath_metadata['xllcorner'] == 2000.5, 'Failed')
3340	self.failUnless(bath_metadata['yllcorner'] == 3000.5, 'Failed')
3341	self.failUnless(bath_metadata['cellsize'] == 25, 'Failed')
3342	self.failUnless(bath_metadata['NODATA_value'] == -9999, 'Failed')
3343	self.failUnless(grid[0][0] == 97.921, 'Failed')
3344	self.failUnless(grid[3][6] == 514.980, 'Failed')
3345
3346	os.remove(filename)
3347
3348	def test_asc_csiro2sww(self):
3349	import tempfile
3350
3351	bath_dir = tempfile.mkdtemp()
3352	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3353	#bath_dir = 'bath_data_manager_test'
3354	#print "os.getcwd( )",os.getcwd( )
3355	elevation_dir = tempfile.mkdtemp()
3356	#elevation_dir = 'elev_expanded'
3357	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3358	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3359
3360	fid = open(bath_dir_filename, 'w')
3361	fid.write(""" ncols 3
3362	nrows 2
3363	xllcorner 148.00000
3364	yllcorner -38.00000
3365	cellsize 0.25
3366	nodata_value -9999.0
3367	9000.000 -1000.000 3000.0
3368	-1000.000 9000.000 -1000.000
3369	""")
3370	fid.close()
3371
3372	fid = open(elevation_dir_filename1, 'w')
3373	fid.write(""" ncols 3
3374	nrows 2
3375	xllcorner 148.00000
3376	yllcorner -38.00000
3377	cellsize 0.25
3378	nodata_value -9999.0
3379	9000.000 0.000 3000.0
3380	0.000 9000.000 0.000
3381	""")
3382	fid.close()
3383
3384	fid = open(elevation_dir_filename2, 'w')
3385	fid.write(""" ncols 3
3386	nrows 2
3387	xllcorner 148.00000
3388	yllcorner -38.00000
3389	cellsize 0.25
3390	nodata_value -9999.0
3391	9000.000 4000.000 4000.0
3392	4000.000 9000.000 4000.000
3393	""")
3394	fid.close()
3395
3396	ucur_dir = tempfile.mkdtemp()
3397	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3398	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3399
3400	fid = open(ucur_dir_filename1, 'w')
3401	fid.write(""" ncols 3
3402	nrows 2
3403	xllcorner 148.00000
3404	yllcorner -38.00000
3405	cellsize 0.25
3406	nodata_value -9999.0
3407	90.000 60.000 30.0
3408	10.000 10.000 10.000
3409	""")
3410	fid.close()
3411	fid = open(ucur_dir_filename2, 'w')
3412	fid.write(""" ncols 3
3413	nrows 2
3414	xllcorner 148.00000
3415	yllcorner -38.00000
3416	cellsize 0.25
3417	nodata_value -9999.0
3418	90.000 60.000 30.0
3419	10.000 10.000 10.000
3420	""")
3421	fid.close()
3422
3423	vcur_dir = tempfile.mkdtemp()
3424	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3425	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3426
3427	fid = open(vcur_dir_filename1, 'w')
3428	fid.write(""" ncols 3
3429	nrows 2
3430	xllcorner 148.00000
3431	yllcorner -38.00000
3432	cellsize 0.25
3433	nodata_value -9999.0
3434	90.000 60.000 30.0
3435	10.000 10.000 10.000
3436	""")
3437	fid.close()
3438	fid = open(vcur_dir_filename2, 'w')
3439	fid.write(""" ncols 3
3440	nrows 2
3441	xllcorner 148.00000
3442	yllcorner -38.00000
3443	cellsize 0.25
3444	nodata_value -9999.0
3445	90.000 60.000 30.0
3446	10.000 10.000 10.000
3447	""")
3448	fid.close()
3449
3450	sww_file = 'a_test.sww'
3451	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
3452
3453	# check the sww file
3454
3455	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
3456	x = fid.variables['x'][:]
3457	y = fid.variables['y'][:]
3458	z = fid.variables['z'][:]
3459	stage = fid.variables['stage'][:]
3460	xmomentum = fid.variables['xmomentum'][:]
3461	geo_ref = Geo_reference(NetCDFObject=fid)
3462	#print "geo_ref",geo_ref
3463	x_ref = geo_ref.get_xllcorner()
3464	y_ref = geo_ref.get_yllcorner()
3465	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
3466	assert allclose(x_ref, 587798.418) # (-38, 148)
3467	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3468
3469	#Zone: 55
3470	#Easting: 588095.674 Northing: 5821451.722
3471	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
3472	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3473
3474	#Zone: 55
3475	#Easting: 632145.632 Northing: 5820863.269
3476	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3477	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3478
3479	#Zone: 55
3480	#Easting: 609748.788 Northing: 5793447.860
3481	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
3482	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
3483
3484	assert allclose(z[0],9000.0 )
3485	assert allclose(stage[0][1],0.0 )
3486
3487	#(4000+1000)*60
3488	assert allclose(xmomentum[1][1],300000.0 )
3489
3490
3491	fid.close()
3492
3493	#tidy up
3494	os.remove(bath_dir_filename)
3495	os.rmdir(bath_dir)
3496
3497	os.remove(elevation_dir_filename1)
3498	os.remove(elevation_dir_filename2)
3499	os.rmdir(elevation_dir)
3500
3501	os.remove(ucur_dir_filename1)
3502	os.remove(ucur_dir_filename2)
3503	os.rmdir(ucur_dir)
3504
3505	os.remove(vcur_dir_filename1)
3506	os.remove(vcur_dir_filename2)
3507	os.rmdir(vcur_dir)
3508
3509
3510	# remove sww file
3511	os.remove(sww_file)
3512
3513	def test_asc_csiro2sww2(self):
3514	import tempfile
3515
3516	bath_dir = tempfile.mkdtemp()
3517	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3518	#bath_dir = 'bath_data_manager_test'
3519	#print "os.getcwd( )",os.getcwd( )
3520	elevation_dir = tempfile.mkdtemp()
3521	#elevation_dir = 'elev_expanded'
3522	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3523	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3524
3525	fid = open(bath_dir_filename, 'w')
3526	fid.write(""" ncols 3
3527	nrows 2
3528	xllcorner 148.00000
3529	yllcorner -38.00000
3530	cellsize 0.25
3531	nodata_value -9999.0
3532	9000.000 -1000.000 3000.0
3533	-1000.000 9000.000 -1000.000
3534	""")
3535	fid.close()
3536
3537	fid = open(elevation_dir_filename1, 'w')
3538	fid.write(""" ncols 3
3539	nrows 2
3540	xllcorner 148.00000
3541	yllcorner -38.00000
3542	cellsize 0.25
3543	nodata_value -9999.0
3544	9000.000 0.000 3000.0
3545	0.000 -9999.000 -9999.000
3546	""")
3547	fid.close()
3548
3549	fid = open(elevation_dir_filename2, 'w')
3550	fid.write(""" ncols 3
3551	nrows 2
3552	xllcorner 148.00000
3553	yllcorner -38.00000
3554	cellsize 0.25
3555	nodata_value -9999.0
3556	9000.000 4000.000 4000.0
3557	4000.000 9000.000 4000.000
3558	""")
3559	fid.close()
3560
3561	ucur_dir = tempfile.mkdtemp()
3562	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3563	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3564
3565	fid = open(ucur_dir_filename1, 'w')
3566	fid.write(""" ncols 3
3567	nrows 2
3568	xllcorner 148.00000
3569	yllcorner -38.00000
3570	cellsize 0.25
3571	nodata_value -9999.0
3572	90.000 60.000 30.0
3573	10.000 10.000 10.000
3574	""")
3575	fid.close()
3576	fid = open(ucur_dir_filename2, 'w')
3577	fid.write(""" ncols 3
3578	nrows 2
3579	xllcorner 148.00000
3580	yllcorner -38.00000
3581	cellsize 0.25
3582	nodata_value -9999.0
3583	90.000 60.000 30.0
3584	10.000 10.000 10.000
3585	""")
3586	fid.close()
3587
3588	vcur_dir = tempfile.mkdtemp()
3589	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3590	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3591
3592	fid = open(vcur_dir_filename1, 'w')
3593	fid.write(""" ncols 3
3594	nrows 2
3595	xllcorner 148.00000
3596	yllcorner -38.00000
3597	cellsize 0.25
3598	nodata_value -9999.0
3599	90.000 60.000 30.0
3600	10.000 10.000 10.000
3601	""")
3602	fid.close()
3603	fid = open(vcur_dir_filename2, 'w')
3604	fid.write(""" ncols 3
3605	nrows 2
3606	xllcorner 148.00000
3607	yllcorner -38.00000
3608	cellsize 0.25
3609	nodata_value -9999.0
3610	90.000 60.000 30.0
3611	10.000 10.000 10.000
3612	""")
3613	fid.close()
3614
3615	try:
3616	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3617	vcur_dir, sww_file)
3618	except:
3619	#tidy up
3620	os.remove(bath_dir_filename)
3621	os.rmdir(bath_dir)
3622
3623	os.remove(elevation_dir_filename1)
3624	os.remove(elevation_dir_filename2)
3625	os.rmdir(elevation_dir)
3626
3627	os.remove(ucur_dir_filename1)
3628	os.remove(ucur_dir_filename2)
3629	os.rmdir(ucur_dir)
3630
3631	os.remove(vcur_dir_filename1)
3632	os.remove(vcur_dir_filename2)
3633	os.rmdir(vcur_dir)
3634	else:
3635	#tidy up
3636	os.remove(bath_dir_filename)
3637	os.rmdir(bath_dir)
3638
3639	os.remove(elevation_dir_filename1)
3640	os.remove(elevation_dir_filename2)
3641	os.rmdir(elevation_dir)
3642	raise 'Should raise exception'
3643
3644	os.remove(ucur_dir_filename1)
3645	os.remove(ucur_dir_filename2)
3646	os.rmdir(ucur_dir)
3647
3648	os.remove(vcur_dir_filename1)
3649	os.remove(vcur_dir_filename2)
3650	os.rmdir(vcur_dir)
3651
3652
3653
3654	def test_asc_csiro2sww3(self):
3655	import tempfile
3656
3657	bath_dir = tempfile.mkdtemp()
3658	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3659	#bath_dir = 'bath_data_manager_test'
3660	#print "os.getcwd( )",os.getcwd( )
3661	elevation_dir = tempfile.mkdtemp()
3662	#elevation_dir = 'elev_expanded'
3663	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3664	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3665
3666	fid = open(bath_dir_filename, 'w')
3667	fid.write(""" ncols 3
3668	nrows 2
3669	xllcorner 148.00000
3670	yllcorner -38.00000
3671	cellsize 0.25
3672	nodata_value -9999.0
3673	9000.000 -1000.000 3000.0
3674	-1000.000 9000.000 -1000.000
3675	""")
3676	fid.close()
3677
3678	fid = open(elevation_dir_filename1, 'w')
3679	fid.write(""" ncols 3
3680	nrows 2
3681	xllcorner 148.00000
3682	yllcorner -38.00000
3683	cellsize 0.25
3684	nodata_value -9999.0
3685	9000.000 0.000 3000.0
3686	0.000 -9999.000 -9999.000
3687	""")
3688	fid.close()
3689
3690	fid = open(elevation_dir_filename2, 'w')
3691	fid.write(""" ncols 3
3692	nrows 2
3693	xllcorner 148.00000
3694	yllcorner -38.00000
3695	cellsize 0.25
3696	nodata_value -9999.0
3697	9000.000 4000.000 4000.0
3698	4000.000 9000.000 4000.000
3699	""")
3700	fid.close()
3701
3702	ucur_dir = tempfile.mkdtemp()
3703	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3704	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3705
3706	fid = open(ucur_dir_filename1, 'w')
3707	fid.write(""" ncols 3
3708	nrows 2
3709	xllcorner 148.00000
3710	yllcorner -38.00000
3711	cellsize 0.25
3712	nodata_value -9999.0
3713	90.000 60.000 30.0
3714	10.000 10.000 10.000
3715	""")
3716	fid.close()
3717	fid = open(ucur_dir_filename2, 'w')
3718	fid.write(""" ncols 3
3719	nrows 2
3720	xllcorner 148.00000
3721	yllcorner -38.00000
3722	cellsize 0.25
3723	nodata_value -9999.0
3724	90.000 60.000 30.0
3725	10.000 10.000 10.000
3726	""")
3727	fid.close()
3728
3729	vcur_dir = tempfile.mkdtemp()
3730	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3731	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3732
3733	fid = open(vcur_dir_filename1, 'w')
3734	fid.write(""" ncols 3
3735	nrows 2
3736	xllcorner 148.00000
3737	yllcorner -38.00000
3738	cellsize 0.25
3739	nodata_value -9999.0
3740	90.000 60.000 30.0
3741	10.000 10.000 10.000
3742	""")
3743	fid.close()
3744	fid = open(vcur_dir_filename2, 'w')
3745	fid.write(""" ncols 3
3746	nrows 2
3747	xllcorner 148.00000
3748	yllcorner -38.00000
3749	cellsize 0.25
3750	nodata_value -9999.0
3751	90.000 60.000 30.0
3752	10.000 10.000 10.000
3753	""")
3754	fid.close()
3755
3756	sww_file = 'a_test.sww'
3757	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3758	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3759	mean_stage = 100)
3760
3761	# check the sww file
3762
3763	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
3764	x = fid.variables['x'][:]
3765	y = fid.variables['y'][:]
3766	z = fid.variables['z'][:]
3767	stage = fid.variables['stage'][:]
3768	xmomentum = fid.variables['xmomentum'][:]
3769	geo_ref = Geo_reference(NetCDFObject=fid)
3770	#print "geo_ref",geo_ref
3771	x_ref = geo_ref.get_xllcorner()
3772	y_ref = geo_ref.get_yllcorner()
3773	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
3774	assert allclose(x_ref, 587798.418) # (-38, 148)
3775	assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3776
3777	#Zone: 55
3778	#Easting: 588095.674 Northing: 5821451.722
3779	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
3780	assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3781
3782	#Zone: 55
3783	#Easting: 632145.632 Northing: 5820863.269
3784	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3785	assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3786
3787	#Zone: 55
3788	#Easting: 609748.788 Northing: 5793447.860
3789	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
3790	assert allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
3791
3792	assert allclose(z[0],9000.0 )
3793	assert allclose(stage[0][4],100.0 )
3794	assert allclose(stage[0][5],100.0 )
3795
3796	#(100.0 - 9000)*10
3797	assert allclose(xmomentum[0][4], -89000.0 )
3798
3799	#(100.0 - -1000.000)*10
3800	assert allclose(xmomentum[0][5], 11000.0 )
3801
3802	fid.close()
3803
3804	#tidy up
3805	os.remove(bath_dir_filename)
3806	os.rmdir(bath_dir)
3807
3808	os.remove(elevation_dir_filename1)
3809	os.remove(elevation_dir_filename2)
3810	os.rmdir(elevation_dir)
3811
3812	os.remove(ucur_dir_filename1)
3813	os.remove(ucur_dir_filename2)
3814	os.rmdir(ucur_dir)
3815
3816	os.remove(vcur_dir_filename1)
3817	os.remove(vcur_dir_filename2)
3818	os.rmdir(vcur_dir)
3819
3820	# remove sww file
3821	os.remove(sww_file)
3822
3823
3824	def test_asc_csiro2sww4(self):
3825	"""
3826	Test specifying the extent
3827	"""
3828
3829	import tempfile
3830
3831	bath_dir = tempfile.mkdtemp()
3832	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3833	#bath_dir = 'bath_data_manager_test'
3834	#print "os.getcwd( )",os.getcwd( )
3835	elevation_dir = tempfile.mkdtemp()
3836	#elevation_dir = 'elev_expanded'
3837	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3838	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3839
3840	fid = open(bath_dir_filename, 'w')
3841	fid.write(""" ncols 4
3842	nrows 4
3843	xllcorner 148.00000
3844	yllcorner -38.00000
3845	cellsize 0.25
3846	nodata_value -9999.0
3847	-9000.000 -1000.000 -3000.0 -2000.000
3848	-1000.000 9000.000 -1000.000 -3000.000
3849	-4000.000 6000.000 2000.000 -5000.000
3850	-9000.000 -1000.000 -3000.0 -2000.000
3851	""")
3852	fid.close()
3853
3854	fid = open(elevation_dir_filename1, 'w')
3855	fid.write(""" ncols 4
3856	nrows 4
3857	xllcorner 148.00000
3858	yllcorner -38.00000
3859	cellsize 0.25
3860	nodata_value -9999.0
3861	-900.000 -100.000 -300.0 -200.000
3862	-100.000 900.000 -100.000 -300.000
3863	-400.000 600.000 200.000 -500.000
3864	-900.000 -100.000 -300.0 -200.000
3865	""")
3866	fid.close()
3867
3868	fid = open(elevation_dir_filename2, 'w')
3869	fid.write(""" ncols 4
3870	nrows 4
3871	xllcorner 148.00000
3872	yllcorner -38.00000
3873	cellsize 0.25
3874	nodata_value -9999.0
3875	-990.000 -110.000 -330.0 -220.000
3876	-110.000 990.000 -110.000 -330.000
3877	-440.000 660.000 220.000 -550.000
3878	-990.000 -110.000 -330.0 -220.000
3879	""")
3880	fid.close()
3881
3882	ucur_dir = tempfile.mkdtemp()
3883	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3884	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3885
3886	fid = open(ucur_dir_filename1, 'w')
3887	fid.write(""" ncols 4
3888	nrows 4
3889	xllcorner 148.00000
3890	yllcorner -38.00000
3891	cellsize 0.25
3892	nodata_value -9999.0
3893	-90.000 -10.000 -30.0 -20.000
3894	-10.000 90.000 -10.000 -30.000
3895	-40.000 60.000 20.000 -50.000
3896	-90.000 -10.000 -30.0 -20.000
3897	""")
3898	fid.close()
3899	fid = open(ucur_dir_filename2, 'w')
3900	fid.write(""" ncols 4
3901	nrows 4
3902	xllcorner 148.00000
3903	yllcorner -38.00000
3904	cellsize 0.25
3905	nodata_value -9999.0
3906	-90.000 -10.000 -30.0 -20.000
3907	-10.000 99.000 -11.000 -30.000
3908	-40.000 66.000 22.000 -50.000
3909	-90.000 -10.000 -30.0 -20.000
3910	""")
3911	fid.close()
3912
3913	vcur_dir = tempfile.mkdtemp()
3914	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3915	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3916
3917	fid = open(vcur_dir_filename1, 'w')
3918	fid.write(""" ncols 4
3919	nrows 4
3920	xllcorner 148.00000
3921	yllcorner -38.00000
3922	cellsize 0.25
3923	nodata_value -9999.0
3924	-90.000 -10.000 -30.0 -20.000
3925	-10.000 80.000 -20.000 -30.000
3926	-40.000 50.000 10.000 -50.000
3927	-90.000 -10.000 -30.0 -20.000
3928	""")
3929	fid.close()
3930	fid = open(vcur_dir_filename2, 'w')
3931	fid.write(""" ncols 4
3932	nrows 4
3933	xllcorner 148.00000
3934	yllcorner -38.00000
3935	cellsize 0.25
3936	nodata_value -9999.0
3937	-90.000 -10.000 -30.0 -20.000
3938	-10.000 88.000 -22.000 -30.000
3939	-40.000 55.000 11.000 -50.000
3940	-90.000 -10.000 -30.0 -20.000
3941	""")
3942	fid.close()
3943
3944	sww_file = tempfile.mktemp(".sww")
3945	#sww_file = 'a_test.sww'
3946	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3947	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3948	mean_stage = 100,
3949	minlat = -37.6, maxlat = -37.6,
3950	minlon = 148.3, maxlon = 148.3
3951	#,verbose = True
3952	)
3953
3954	# check the sww file
3955
3956	fid = NetCDFFile(sww_file, 'r') #Open existing file for read
3957	x = fid.variables['x'][:]
3958	y = fid.variables['y'][:]
3959	z = fid.variables['z'][:]
3960	stage = fid.variables['stage'][:]
3961	xmomentum = fid.variables['xmomentum'][:]
3962	ymomentum = fid.variables['ymomentum'][:]
3963	geo_ref = Geo_reference(NetCDFObject=fid)
3964	#print "geo_ref",geo_ref
3965	x_ref = geo_ref.get_xllcorner()
3966	y_ref = geo_ref.get_yllcorner()
3967	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
3968
3969	assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3970	assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3971	assert allclose(y_ref, 5820863.269 )# (-37.45, 148.5)
3972
3973	#Easting: 632145.632 Northing: 5820863.269
3974	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3975
3976	#print "x",x
3977	#print "y",y
3978	self.failUnless(len(x) == 4,'failed') # 2*2
3979	self.failUnless(len(x) == 4,'failed') # 2*2
3980
3981	#Zone: 55
3982	#Easting: 632145.632 Northing: 5820863.269
3983	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
3984	# magic number - y is close enough for me.
3985	assert allclose(x[3], 632145.63 - x_ref)
3986	assert allclose(y[3], 5820863.269 - y_ref + 5.22155314684e-005)
3987
3988	assert allclose(z[0],9000.0 ) #z is elevation info
3989	#print "z",z
3990	# 2 time steps, 4 points
3991	self.failUnless(xmomentum.shape == (2,4), 'failed')
3992	self.failUnless(ymomentum.shape == (2,4), 'failed')
3993
3994	#(100.0 - -1000.000)*10
3995	#assert allclose(xmomentum[0][5], 11000.0 )
3996
3997	fid.close()
3998
3999	# is the sww file readable?
4000	#Lets see if we can convert it to a dem!
4001	#print "sww_file",sww_file
4002	#dem_file = tempfile.mktemp(".dem")
4003	domain = sww2domain(sww_file) ###, dem_file)
4004	domain.check_integrity()
4005
4006	#tidy up
4007	os.remove(bath_dir_filename)
4008	os.rmdir(bath_dir)
4009
4010	os.remove(elevation_dir_filename1)
4011	os.remove(elevation_dir_filename2)
4012	os.rmdir(elevation_dir)
4013
4014	os.remove(ucur_dir_filename1)
4015	os.remove(ucur_dir_filename2)
4016	os.rmdir(ucur_dir)
4017
4018	os.remove(vcur_dir_filename1)
4019	os.remove(vcur_dir_filename2)
4020	os.rmdir(vcur_dir)
4021
4022
4023
4024
4025	# remove sww file
4026	os.remove(sww_file)
4027
4028	# remove dem file
4029	#os.remove(dem_file)
4030
4031	def test_get_min_max_indexes(self):
4032	latitudes = [3,2,1,0]
4033	longitudes = [0,10,20,30]
4034
4035	# k - lat
4036	# l - lon
4037	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4038	latitudes,longitudes,
4039	-10,4,-10,31)
4040
4041	#print "kmin",kmin;print "kmax",kmax
4042	#print "lmin",lmin;print "lmax",lmax
4043	latitudes_new = latitudes[kmin:kmax]
4044	longitudes_news = longitudes[lmin:lmax]
4045	#print "latitudes_new", latitudes_new
4046	#print "longitudes_news",longitudes_news
4047	self.failUnless(latitudes == latitudes_new and \
4048	longitudes == longitudes_news,
4049	'failed')
4050
4051	## 2nd test
4052	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4053	latitudes,longitudes,
4054	0.5,2.5,5,25)
4055	#print "kmin",kmin;print "kmax",kmax
4056	#print "lmin",lmin;print "lmax",lmax
4057	latitudes_new = latitudes[kmin:kmax]
4058	longitudes_news = longitudes[lmin:lmax]
4059	#print "latitudes_new", latitudes_new
4060	#print "longitudes_news",longitudes_news
4061
4062	self.failUnless(latitudes == latitudes_new and \
4063	longitudes == longitudes_news,
4064	'failed')
4065
4066	## 3rd test
4067	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(\
4068	latitudes,
4069	longitudes,
4070	1.1,1.9,12,17)
4071	#print "kmin",kmin;print "kmax",kmax
4072	#print "lmin",lmin;print "lmax",lmax
4073	latitudes_new = latitudes[kmin:kmax]
4074	longitudes_news = longitudes[lmin:lmax]
4075	#print "latitudes_new", latitudes_new
4076	#print "longitudes_news",longitudes_news
4077
4078	self.failUnless(latitudes_new == [2, 1] and \
4079	longitudes_news == [10, 20],
4080	'failed')
4081
4082
4083	## 4th test
4084	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4085	latitudes,longitudes,
4086	-0.1,1.9,-2,17)
4087	#print "kmin",kmin;print "kmax",kmax
4088	#print "lmin",lmin;print "lmax",lmax
4089	latitudes_new = latitudes[kmin:kmax]
4090	longitudes_news = longitudes[lmin:lmax]
4091	#print "latitudes_new", latitudes_new
4092	#print "longitudes_news",longitudes_news
4093
4094	self.failUnless(latitudes_new == [2, 1, 0] and \
4095	longitudes_news == [0, 10, 20],
4096	'failed')
4097	## 5th test
4098	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4099	latitudes,longitudes,
4100	0.1,1.9,2,17)
4101	#print "kmin",kmin;print "kmax",kmax
4102	#print "lmin",lmin;print "lmax",lmax
4103	latitudes_new = latitudes[kmin:kmax]
4104	longitudes_news = longitudes[lmin:lmax]
4105	#print "latitudes_new", latitudes_new
4106	#print "longitudes_news",longitudes_news
4107
4108	self.failUnless(latitudes_new == [2, 1, 0] and \
4109	longitudes_news == [0, 10, 20],
4110	'failed')
4111
4112	## 6th test
4113
4114	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4115	latitudes,longitudes,
4116	1.5,4,18,32)
4117	#print "kmin",kmin;print "kmax",kmax
4118	#print "lmin",lmin;print "lmax",lmax
4119	latitudes_new = latitudes[kmin:kmax]
4120	longitudes_news = longitudes[lmin:lmax]
4121	#print "latitudes_new", latitudes_new
4122	#print "longitudes_news",longitudes_news
4123
4124	self.failUnless(latitudes_new == [3, 2, 1] and \
4125	longitudes_news == [10, 20, 30],
4126	'failed')
4127
4128
4129	## 7th test
4130	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4131	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4132	latitudes,longitudes,
4133	1.5,1.5,15,15)
4134	#print "kmin",kmin;print "kmax",kmax
4135	#print "lmin",lmin;print "lmax",lmax
4136	latitudes_new = latitudes[kmin:kmax]
4137	longitudes_news = longitudes[lmin:lmax]
4138	m2d = m2d[kmin:kmax,lmin:lmax]
4139	#print "m2d", m2d
4140	#print "latitudes_new", latitudes_new
4141	#print "longitudes_news",longitudes_news
4142
4143	self.failUnless(latitudes_new == [2, 1] and \
4144	longitudes_news == [10, 20],
4145	'failed')
4146
4147	self.failUnless(m2d == [[5,6],[9,10]],
4148	'failed')
4149
4150	def test_get_min_max_indexes2(self):
4151	latitudes = [-30,-35,-40,-45]
4152	longitudes = [148,149,150,151]
4153
4154	m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4155
4156	# k - lat
4157	# l - lon
4158	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4159	latitudes,longitudes,
4160	-37,-27,147,149.5)
4161
4162	#print "kmin",kmin;print "kmax",kmax
4163	#print "lmin",lmin;print "lmax",lmax
4164	#print "m2d", m2d
4165	#print "latitudes", latitudes
4166	#print "longitudes",longitudes
4167	#print "latitudes[kmax]", latitudes[kmax]
4168	latitudes_new = latitudes[kmin:kmax]
4169	longitudes_new = longitudes[lmin:lmax]
4170	m2d = m2d[kmin:kmax,lmin:lmax]
4171	#print "m2d", m2d
4172	#print "latitudes_new", latitudes_new
4173	#print "longitudes_new",longitudes_new
4174
4175	self.failUnless(latitudes_new == [-30, -35, -40] and \
4176	longitudes_new == [148, 149,150],
4177	'failed')
4178	self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
4179	'failed')
4180
4181	def test_get_min_max_indexes3(self):
4182	latitudes = [-30,-35,-40,-45,-50,-55,-60]
4183	longitudes = [148,149,150,151]
4184
4185	# k - lat
4186	# l - lon
4187	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4188	latitudes,longitudes,
4189	-43,-37,148.5,149.5)
4190
4191
4192	#print "kmin",kmin;print "kmax",kmax
4193	#print "lmin",lmin;print "lmax",lmax
4194	#print "latitudes", latitudes
4195	#print "longitudes",longitudes
4196	latitudes_new = latitudes[kmin:kmax]
4197	longitudes_news = longitudes[lmin:lmax]
4198	#print "latitudes_new", latitudes_new
4199	#print "longitudes_news",longitudes_news
4200
4201	self.failUnless(latitudes_new == [-35, -40, -45] and \
4202	longitudes_news == [148, 149,150],
4203	'failed')
4204
4205	def test_get_min_max_indexes4(self):
4206	latitudes = [-30,-35,-40,-45,-50,-55,-60]
4207	longitudes = [148,149,150,151]
4208
4209	# k - lat
4210	# l - lon
4211	kmin, kmax, lmin, lmax = anuga.pyvolution.data_manager._get_min_max_indexes(
4212	latitudes,longitudes)
4213
4214	#print "kmin",kmin;print "kmax",kmax
4215	#print "lmin",lmin;print "lmax",lmax
4216	#print "latitudes", latitudes
4217	#print "longitudes",longitudes
4218	latitudes_new = latitudes[kmin:kmax]
4219	longitudes_news = longitudes[lmin:lmax]
4220	#print "latitudes_new", latitudes_new
4221	#print "longitudes_news",longitudes_news
4222
4223	self.failUnless(latitudes_new == latitudes and \
4224	longitudes_news == longitudes,
4225	'failed')
4226
4227	def test_tsh2sww(self):
4228	import os
4229	import tempfile
4230
4231	tsh_file = tempfile.mktemp(".tsh")
4232	file = open(tsh_file,"w")
4233	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
4234	0 0.0 0.0 0.0 0.0 0.01 \n \
4235	1 1.0 0.0 10.0 10.0 0.02 \n \
4236	2 0.0 1.0 0.0 10.0 0.03 \n \
4237	3 0.5 0.25 8.0 12.0 0.04 \n \
4238	# Vert att title \n \
4239	elevation \n \
4240	stage \n \
4241	friction \n \
4242	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
4243	0 0 3 2 -1 -1 1 dsg\n\
4244	1 0 1 3 -1 0 -1 ole nielsen\n\
4245	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
4246	0 1 0 2 \n\
4247	1 0 2 3 \n\
4248	2 2 3 \n\
4249	3 3 1 1 \n\
4250	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
4251	0 216.0 -86.0 \n \
4252	1 160.0 -167.0 \n \
4253	2 114.0 -91.0 \n \
4254	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
4255	0 0 1 0 \n \
4256	1 1 2 0 \n \
4257	2 2 0 0 \n \
4258	0 # <x> <y> ...Mesh Holes... \n \
4259	0 # <x> <y> <attribute>...Mesh Regions... \n \
4260	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
4261	#Geo reference \n \
4262	56 \n \
4263	140 \n \
4264	120 \n")
4265	file.close()
4266
4267	#sww_file = tempfile.mktemp(".sww")
4268	#print "sww_file",sww_file
4269	#print "sww_file",tsh_file
4270	tsh2sww(tsh_file)
4271
4272	os.remove(tsh_file)
4273	os.remove(tsh_file[:-4] + '.sww')
4274
4275
4276
4277	def NOT_test_dem2pts(self):
4278	"""Test conversion from dem in ascii format to native NetCDF xya format
4279	"""
4280
4281	import time, os
4282	from Numeric import array, zeros, allclose, Float, concatenate
4283	from Scientific.IO.NetCDF import NetCDFFile
4284
4285	#Write test asc file
4286	root = 'demtest'
4287
4288	filename = root+'.asc'
4289	fid = open(filename, 'w')
4290	fid.write("""ncols 5
4291	nrows 6
4292	xllcorner 2000.5
4293	yllcorner 3000.5
4294	cellsize 25
4295	NODATA_value -9999
4296	""")
4297	#Create linear function
4298
4299	ref_points = []
4300	ref_elevation = []
4301	for i in range(6):
4302	y = (6-i)*25.0
4303	for j in range(5):
4304	x = j*25.0
4305	z = x+2*y
4306
4307	ref_points.append( [x,y] )
4308	ref_elevation.append(z)
4309	fid.write('%f ' %z)
4310	fid.write('\n')
4311
4312	fid.close()
4313
4314	#Write prj file with metadata
4315	metafilename = root+'.prj'
4316	fid = open(metafilename, 'w')
4317
4318
4319	fid.write("""Projection UTM
4320	Zone 56
4321	Datum WGS84
4322	Zunits NO
4323	Units METERS
4324	Spheroid WGS84
4325	Xshift 0.0000000000
4326	Yshift 10000000.0000000000
4327	Parameters
4328	""")
4329	fid.close()
4330
4331	#Convert to NetCDF pts
4332	convert_dem_from_ascii2netcdf(root)
4333	dem2pts(root)
4334
4335	#Check contents
4336	#Get NetCDF
4337	fid = NetCDFFile(root+'.pts', 'r')
4338
4339	# Get the variables
4340	#print fid.variables.keys()
4341	points = fid.variables['points']
4342	elevation = fid.variables['elevation']
4343
4344	#Check values
4345
4346	#print points[:]
4347	#print ref_points
4348	assert allclose(points, ref_points)
4349
4350	#print attributes[:]
4351	#print ref_elevation
4352	assert allclose(elevation, ref_elevation)
4353
4354	#Cleanup
4355	fid.close()
4356
4357
4358	os.remove(root + '.pts')
4359	os.remove(root + '.dem')
4360	os.remove(root + '.asc')
4361	os.remove(root + '.prj')
4362
4363
4364
4365	def NOT_test_dem2pts_bounding_box(self):
4366	"""Test conversion from dem in ascii format to native NetCDF xya format
4367	"""
4368
4369	import time, os
4370	from Numeric import array, zeros, allclose, Float, concatenate
4371	from Scientific.IO.NetCDF import NetCDFFile
4372
4373	#Write test asc file
4374	root = 'demtest'
4375
4376	filename = root+'.asc'
4377	fid = open(filename, 'w')
4378	fid.write("""ncols 5
4379	nrows 6
4380	xllcorner 2000.5
4381	yllcorner 3000.5
4382	cellsize 25
4383	NODATA_value -9999
4384	""")
4385	#Create linear function
4386
4387	ref_points = []
4388	ref_elevation = []
4389	for i in range(6):
4390	y = (6-i)*25.0
4391	for j in range(5):
4392	x = j*25.0
4393	z = x+2*y
4394
4395	ref_points.append( [x,y] )
4396	ref_elevation.append(z)
4397	fid.write('%f ' %z)
4398	fid.write('\n')
4399
4400	fid.close()
4401
4402	#Write prj file with metadata
4403	metafilename = root+'.prj'
4404	fid = open(metafilename, 'w')
4405
4406
4407	fid.write("""Projection UTM
4408	Zone 56
4409	Datum WGS84
4410	Zunits NO
4411	Units METERS
4412	Spheroid WGS84
4413	Xshift 0.0000000000
4414	Yshift 10000000.0000000000
4415	Parameters
4416	""")
4417	fid.close()
4418
4419	#Convert to NetCDF pts
4420	convert_dem_from_ascii2netcdf(root)
4421	dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4422	northing_min=3035.0, northing_max=3125.5)
4423
4424	#Check contents
4425	#Get NetCDF
4426	fid = NetCDFFile(root+'.pts', 'r')
4427
4428	# Get the variables
4429	#print fid.variables.keys()
4430	points = fid.variables['points']
4431	elevation = fid.variables['elevation']
4432
4433	#Check values
4434	assert fid.xllcorner[0] == 2010.0
4435	assert fid.yllcorner[0] == 3035.0
4436
4437	#create new reference points
4438	ref_points = []
4439	ref_elevation = []
4440	for i in range(4):
4441	y = (4-i)*25.0 + 25.0
4442	y_new = y + 3000.5 - 3035.0
4443	for j in range(4):
4444	x = j*25.0 + 25.0
4445	x_new = x + 2000.5 - 2010.0
4446	z = x+2*y
4447
4448	ref_points.append( [x_new,y_new] )
4449	ref_elevation.append(z)
4450
4451	#print points[:]
4452	#print ref_points
4453	assert allclose(points, ref_points)
4454
4455	#print attributes[:]
4456	#print ref_elevation
4457	assert allclose(elevation, ref_elevation)
4458
4459	#Cleanup
4460	fid.close()
4461
4462
4463	os.remove(root + '.pts')
4464	os.remove(root + '.dem')
4465	os.remove(root + '.asc')
4466	os.remove(root + '.prj')
4467
4468
4469
4470	def NOT_test_dem2pts_remove_Nullvalues(self):
4471	"""Test conversion from dem in ascii format to native NetCDF xya format
4472	"""
4473
4474	import time, os
4475	from Numeric import array, zeros, allclose, Float, concatenate
4476	from Scientific.IO.NetCDF import NetCDFFile
4477
4478	#Write test asc file
4479	root = 'demtest'
4480
4481	filename = root+'.asc'
4482	fid = open(filename, 'w')
4483	fid.write("""ncols 5
4484	nrows 6
4485	xllcorner 2000.5
4486	yllcorner 3000.5
4487	cellsize 25
4488	NODATA_value -9999
4489	""")
4490	#Create linear function
4491	# ref_ will write all the values
4492	# new_ref_ will write the values except for NODATA_values
4493	ref_points = []
4494	ref_elevation = []
4495	new_ref_pts = []
4496	new_ref_elev = []
4497	NODATA_value = -9999
4498	for i in range(6):
4499	y = (6-i)*25.0
4500	for j in range(5):
4501	x = j*25.0
4502	z = x+2*y
4503	if j == 4: z = NODATA_value # column
4504	if i == 2 and j == 2: z = NODATA_value # random
4505	if i == 5 and j == 1: z = NODATA_value
4506	if i == 1: z = NODATA_value # row
4507	if i == 3 and j == 1: z = NODATA_value # two pts/row
4508	if i == 3 and j == 3: z = NODATA_value
4509
4510
4511	if z <> NODATA_value:
4512	new_ref_elev.append(z)
4513	new_ref_pts.append( [x,y] )
4514
4515	ref_points.append( [x,y] )
4516	ref_elevation.append(z)
4517
4518	fid.write('%f ' %z)
4519	fid.write('\n')
4520
4521	fid.close()
4522
4523
4524	#Write prj file with metadata
4525	metafilename = root+'.prj'
4526	fid = open(metafilename, 'w')
4527
4528
4529	fid.write("""Projection UTM
4530	Zone 56
4531	Datum WGS84
4532	Zunits NO
4533	Units METERS
4534	Spheroid WGS84
4535	Xshift 0.0000000000
4536	Yshift 10000000.0000000000
4537	Parameters
4538	""")
4539	fid.close()
4540
4541	#Convert to NetCDF pts
4542	convert_dem_from_ascii2netcdf(root)
4543	dem2pts(root)
4544
4545	#Check contents
4546	#Get NetCDF
4547	fid = NetCDFFile(root+'.pts', 'r')
4548
4549	# Get the variables
4550	#print fid.variables.keys()
4551	points = fid.variables['points']
4552	elevation = fid.variables['elevation']
4553
4554	#Check values
4555	#print 'points', points[:]
4556	assert len(points) == len(new_ref_pts), 'length of returned points not correct'
4557	assert allclose(points, new_ref_pts), 'points do not align'
4558
4559	#print 'elevation', elevation[:]
4560	assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
4561	assert allclose(elevation, new_ref_elev), 'elevations do not align'
4562
4563	#Cleanup
4564	fid.close()
4565
4566
4567	os.remove(root + '.pts')
4568	os.remove(root + '.dem')
4569	os.remove(root + '.asc')
4570	os.remove(root + '.prj')
4571
4572	def NOT_test_dem2pts_bounding_box_Nullvalues(self):
4573	"""Test conversion from dem in ascii format to native NetCDF xya format
4574	"""
4575
4576	import time, os
4577	from Numeric import array, zeros, allclose, Float, concatenate
4578	from Scientific.IO.NetCDF import NetCDFFile
4579
4580	#Write test asc file
4581	root = 'demtest'
4582
4583	filename = root+'.asc'
4584	fid = open(filename, 'w')
4585	fid.write("""ncols 5
4586	nrows 6
4587	xllcorner 2000.5
4588	yllcorner 3000.5
4589	cellsize 25
4590	NODATA_value -9999
4591	""")
4592	#Create linear function
4593
4594	ref_points = []
4595	ref_elevation = []
4596	new_ref_pts1 = []
4597	new_ref_elev1 = []
4598	NODATA_value = -9999
4599	for i in range(6):
4600	y = (6-i)*25.0
4601	for j in range(5):
4602	x = j*25.0
4603	z = x+2*y
4604	if j == 4: z = NODATA_value # column
4605	if i == 2 and j == 2: z = NODATA_value # random
4606	if i == 5 and j == 1: z = NODATA_value
4607	if i == 1: z = NODATA_value # row
4608	if i == 3 and j == 1: z = NODATA_value # two pts/row
4609	if i == 3 and j == 3: z = NODATA_value
4610
4611	if z <> NODATA_value:
4612	new_ref_elev1.append(z)
4613	new_ref_pts1.append( [x,y] )
4614
4615	ref_points.append( [x,y] )
4616	ref_elevation.append(z)
4617	fid.write('%f ' %z)
4618	fid.write('\n')
4619
4620	fid.close()
4621
4622	#Write prj file with metadata
4623	metafilename = root+'.prj'
4624	fid = open(metafilename, 'w')
4625
4626
4627	fid.write("""Projection UTM
4628	Zone 56
4629	Datum WGS84
4630	Zunits NO
4631	Units METERS
4632	Spheroid WGS84
4633	Xshift 0.0000000000
4634	Yshift 10000000.0000000000
4635	Parameters
4636	""")
4637	fid.close()
4638
4639	#Convert to NetCDF pts
4640	convert_dem_from_ascii2netcdf(root)
4641	dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4642	northing_min=3035.0, northing_max=3125.5)
4643
4644	#Check contents
4645	#Get NetCDF
4646	fid = NetCDFFile(root+'.pts', 'r')
4647
4648	# Get the variables
4649	#print fid.variables.keys()
4650	points = fid.variables['points']
4651	elevation = fid.variables['elevation']
4652
4653	#Check values
4654	assert fid.xllcorner[0] == 2010.0
4655	assert fid.yllcorner[0] == 3035.0
4656
4657	#create new reference points
4658	ref_points = []
4659	ref_elevation = []
4660	new_ref_pts2 = []
4661	new_ref_elev2 = []
4662	for i in range(4):
4663	y = (4-i)*25.0 + 25.0
4664	y_new = y + 3000.5 - 3035.0
4665	for j in range(4):
4666	x = j*25.0 + 25.0
4667	x_new = x + 2000.5 - 2010.0
4668	z = x+2*y
4669
4670	if j == 3: z = NODATA_value # column
4671	if i == 1 and j == 1: z = NODATA_value # random
4672	if i == 4 and j == 0: z = NODATA_value
4673	if i == 0: z = NODATA_value # row
4674	if i == 2 and j == 0: z = NODATA_value # two pts/row
4675	if i == 2 and j == 2: z = NODATA_value
4676
4677	if z <> NODATA_value:
4678	new_ref_elev2.append(z)
4679	new_ref_pts2.append( [x_new,y_new] )
4680
4681
4682	ref_points.append( [x_new,y_new] )
4683	ref_elevation.append(z)
4684
4685	#print points[:]
4686	#print ref_points
4687	#assert allclose(points, ref_points)
4688
4689	#print attributes[:]
4690	#print ref_elevation
4691	#assert allclose(elevation, ref_elevation)
4692
4693
4694	assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
4695	assert allclose(points, new_ref_pts2), 'points do not align'
4696
4697	#print 'elevation', elevation[:]
4698	assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
4699	assert allclose(elevation, new_ref_elev2), 'elevations do not align'
4700	#Cleanup
4701	fid.close()
4702
4703
4704	os.remove(root + '.pts')
4705	os.remove(root + '.dem')
4706	os.remove(root + '.asc')
4707	os.remove(root + '.prj')
4708
4709
4710	########## testing nbed class ##################
4711	def test_exposure_csv_loading(self):
4712
4713
4714	file_name = tempfile.mktemp(".xya")
4715	file = open(file_name,"w")
4716	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
4717	115.0, -21.0, splat, 0.0\n\
4718	114.0, -21.7, pow, 10.0\n\
4719	114.5, -21.4, bang, 40.0\n")
4720	file.close()
4721	exposure = Exposure_csv(file_name)
4722	exposure.get_column("sound")
4723
4724	self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
4725	'FAILED!')
4726	self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
4727	'FAILED!')
4728
4729	os.remove(file_name)
4730
4731	def test_exposure_csv_loading(self):
4732
4733
4734	file_name = tempfile.mktemp(".xya")
4735	file = open(file_name,"w")
4736	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
4737	115.0, -21.0, splat, 0.0\n\
4738	114.0, -21.7, pow, 10.0\n\
4739	114.5, -21.4, bang, 40.0\n")
4740	file.close()
4741	exposure = Exposure_csv(file_name)
4742	exposure.get_column("sound")
4743
4744	self.failUnless(exposure._attribute_dic['sound'][2]==' bang',
4745	'FAILED!')
4746	self.failUnless(exposure._attribute_dic['speed'][2]==' 40.0',
4747	'FAILED!')
4748
4749	os.remove(file_name)
4750
4751	def test_exposure_csv_cmp(self):
4752	file_name = tempfile.mktemp(".xya")
4753	file = open(file_name,"w")
4754	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
4755	115.0, -21.0, splat, 0.0\n\
4756	114.0, -21.7, pow, 10.0\n\
4757	114.5, -21.4, bang, 40.0\n")
4758	file.close()
4759
4760	e1 = Exposure_csv(file_name)
4761	e2 = Exposure_csv(file_name)
4762	os.remove(file_name)
4763
4764	self.failUnless(cmp(e1,e2)==0,
4765	'FAILED!')
4766
4767	self.failUnless(cmp(e1,"hey")==1,
4768	'FAILED!')
4769
4770	file_name = tempfile.mktemp(".xya")
4771	file = open(file_name,"w")
4772	# Note, this has less spaces in the title,
4773	# the instances will be the same.
4774	file.write("LATITUDE,LONGITUDE ,sound, speed \n\
4775	115.0, -21.0, splat, 0.0\n\
4776	114.0, -21.7, pow, 10.0\n\
4777	114.5, -21.4, bang, 40.0\n")
4778	file.close()
4779	e3 = Exposure_csv(file_name)
4780	os.remove(file_name)
4781
4782	self.failUnless(cmp(e3,e2)==0,
4783	'FAILED!')
4784
4785	file_name = tempfile.mktemp(".xya")
4786	file = open(file_name,"w")
4787	# Note, 40 changed to 44 .
4788	file.write("LATITUDE,LONGITUDE ,sound, speed \n\
4789	115.0, -21.0, splat, 0.0\n\
4790	114.0, -21.7, pow, 10.0\n\
4791	114.5, -21.4, bang, 44.0\n")
4792	file.close()
4793	e4 = Exposure_csv(file_name)
4794	os.remove(file_name)
4795	#print "e4",e4._attribute_dic
4796	#print "e2",e2._attribute_dic
4797	self.failUnless(cmp(e4,e2)<>0,
4798	'FAILED!')
4799
4800	file_name = tempfile.mktemp(".xya")
4801	file = open(file_name,"w")
4802	# Note, the first two columns are swapped.
4803	file.write("LONGITUDE,LATITUDE ,sound, speed \n\
4804	-21.0,115.0, splat, 0.0\n\
4805	-21.7,114.0, pow, 10.0\n\
4806	-21.4,114.5, bang, 40.0\n")
4807	file.close()
4808	e5 = Exposure_csv(file_name)
4809	os.remove(file_name)
4810
4811	self.failUnless(cmp(e3,e5)<>0,
4812	'FAILED!')
4813
4814	def test_exposure_csv_saving(self):
4815
4816
4817	file_name = tempfile.mktemp(".xya")
4818	file = open(file_name,"w")
4819	file.write("LATITUDE, LONGITUDE ,sound , speed \n\
4820	115.0, -21.0, splat, 0.0\n\
4821	114.0, -21.7, pow, 10.0\n\
4822	114.5, -21.4, bang, 40.0\n")
4823	file.close()
4824	e1 = Exposure_csv(file_name)
4825
4826	file_name2 = tempfile.mktemp(".xya")
4827	e1.save(file_name = file_name2)
4828	e2 = Exposure_csv(file_name2)
4829
4830	self.failUnless(cmp(e1,e2)==0,
4831	'FAILED!')
4832	os.remove(file_name)
4833	os.remove(file_name2)
4834
4835	def test_exposure_csv_get_location(self):
4836	file_name = tempfile.mktemp(".xya")
4837	file = open(file_name,"w")
4838	file.write("LONGITUDE , LATITUDE, sound , speed \n\
4839	150.916666667, -34.5, splat, 0.0\n\
4840	150.0, -34.0, pow, 10.0\n")
4841	file.close()
4842	e1 = Exposure_csv(file_name)
4843
4844	gsd = e1.get_location()
4845
4846	points = gsd.get_data_points(absolute=True)
4847
4848	assert allclose(points[0][0], 308728.009)
4849	assert allclose(points[0][1], 6180432.601)
4850	assert allclose(points[1][0], 222908.705)
4851	assert allclose(points[1][1], 6233785.284)
4852	self.failUnless(gsd.get_geo_reference().get_zone() == 56,
4853	'Bad zone error!')
4854
4855	os.remove(file_name)
4856
4857	def test_exposure_csv_set_column_get_column(self):
4858	file_name = tempfile.mktemp(".xya")
4859	file = open(file_name,"w")
4860	file.write("LONGITUDE , LATITUDE, sound , speed \n\
4861	150.916666667, -34.5, splat, 0.0\n\
4862	150.0, -34.0, pow, 10.0\n")
4863	file.close()
4864	e1 = Exposure_csv(file_name)
4865	os.remove(file_name)
4866
4867	new_title = "feast"
4868	new_values = ["chicken","soup"]
4869	e1.set_column(new_title, new_values)
4870	returned_values = e1.get_column(new_title)
4871	self.failUnless(returned_values == new_values,
4872	' Error!')
4873
4874	file_name2 = tempfile.mktemp(".xya")
4875	e1.save(file_name = file_name2)
4876	e2 = Exposure_csv(file_name2)
4877	returned_values = e2.get_column(new_title)
4878	self.failUnless(returned_values == new_values,
4879	' Error!')
4880	os.remove(file_name2)
4881
4882	def test_exposure_csv_set_column_get_column_error_checking(self):
4883	file_name = tempfile.mktemp(".xya")
4884	file = open(file_name,"w")
4885	file.write("LONGITUDE , LATITUDE, sound , speed \n\
4886	150.916666667, -34.5, splat, 0.0\n\
4887	150.0, -34.0, pow, 10.0\n")
4888	file.close()
4889	e1 = Exposure_csv(file_name)
4890	os.remove(file_name)
4891
4892	new_title = "sound"
4893	new_values = [12.5,7.6]
4894	try:
4895	e1.set_column(new_title, new_values)
4896	except TitleValueError:
4897	pass
4898	else:
4899	self.failUnless(0 ==1, 'Error not thrown error!')
4900
4901	e1.set_column(new_title, new_values, overwrite=True)
4902	returned_values = e1.get_column(new_title)
4903	self.failUnless(returned_values == new_values,
4904	' Error!')
4905
4906	new2_title = "short list"
4907	new2_values = [12.5]
4908	try:
4909	e1.set_column(new2_title, new2_values)
4910	except DataMissingValuesError:
4911	pass
4912	else:
4913	self.failUnless(0 ==1, 'Error not thrown error!')
4914
4915	new2_title = "long list"
4916	new2_values = [12.5, 7,8]
4917	try:
4918	e1.set_column(new2_title, new2_values)
4919	except DataMissingValuesError:
4920	pass
4921	else:
4922	self.failUnless(0 ==1, 'Error not thrown error!')
4923	file_name2 = tempfile.mktemp(".xya")
4924	e1.save(file_name = file_name2)
4925	e2 = Exposure_csv(file_name2)
4926	returned_values = e2.get_column(new_title)
4927	for returned, new in map(None, returned_values, new_values):
4928	self.failUnless(returned == str(new), ' Error!')
4929	#self.failUnless(returned_values == new_values, ' Error!')
4930	os.remove(file_name2)
4931
4932	try:
4933	e1.get_column("toe jam")
4934	except TitleValueError:
4935	pass
4936	else:
4937	self.failUnless(0 ==1, 'Error not thrown error!')
4938
4939	def test_exposure_csv_loading_x_y(self):
4940
4941
4942	file_name = tempfile.mktemp(".xya")
4943	file = open(file_name,"w")
4944	file.write("x, y ,sound , speed \n\
4945	115.0, 7, splat, 0.0\n\
4946	114.0, 8.0, pow, 10.0\n\
4947	114.5, 9., bang, 40.0\n")
4948	file.close()
4949	e1 = Exposure_csv(file_name, is_x_y_locations=True)
4950	gsd = e1.get_location()
4951
4952	points = gsd.get_data_points(absolute=True)
4953
4954	assert allclose(points[0][0], 115)
4955	assert allclose(points[0][1], 7)
4956	assert allclose(points[1][0], 114)
4957	assert allclose(points[1][1], 8)
4958	assert allclose(points[2][0], 114.5)
4959	assert allclose(points[2][1], 9)
4960	self.failUnless(gsd.get_geo_reference().get_zone() == -1,
4961	'Bad zone error!')
4962
4963	os.remove(file_name)
4964
4965
4966	def test_exposure_csv_loading_x_y2(self):
4967
4968	csv_file = tempfile.mktemp(".csv")
4969	fd = open(csv_file,'wb')
4970	writer = csv.writer(fd)
4971	writer.writerow(['x','y','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
4972	writer.writerow([5.5,0.5,'199770','130000','Metal','Timber',20])
4973	writer.writerow([4.5,1.0,'150000','76000','Metal','Double Brick',20])
4974	writer.writerow([4.5,1.5,'150000','76000','Metal','Brick Veneer',20])
4975	fd.close()
4976
4977	e1 = Exposure_csv(csv_file)
4978	gsd = e1.get_location()
4979
4980	points = gsd.get_data_points(absolute=True)
4981	assert allclose(points[0][0], 5.5)
4982	assert allclose(points[0][1], 0.5)
4983	assert allclose(points[1][0], 4.5)
4984	assert allclose(points[1][1], 1.0)
4985	assert allclose(points[2][0], 4.5)
4986	assert allclose(points[2][1], 1.5)
4987	self.failUnless(gsd.get_geo_reference().get_zone() == -1,
4988	'Bad zone error!')
4989
4990	os.remove(csv_file)
4991
4992	#-------------------------------------------------------------
4993	if __name__ == "__main__":
4994	#suite = unittest.makeSuite(Test_Data_Manager,'test_sww_minimum_sww_depth')
4995	suite = unittest.makeSuite(Test_Data_Manager,'test')
4996	runner = unittest.TextTestRunner()
4997	runner.run(suite)
4998
4999

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

Download in other formats: