Context Navigation

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

Last change on this file since 3720 was 3720, checked in by duncan, 18 years ago
added urs2sww to convert urs tsunami info to an sww boundary
File size: 161.2 KB

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

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

Download in other formats: