Context Navigation

source: inundation/pyvolution/test_data_manager.py @ 2553

Last change on this file since 2553 was 2553, checked in by sexton, 18 years ago
maybe this is it for dem2pts ...
File size: 128.1 KB

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

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

Download in other formats: