Context Navigation

source: inundation/pyvolution/test_data_manager.py @ 2267

Last change on this file since 2267 was 2263, checked in by ole, 19 years ago
Work towards a test_all.py at the inundation directory level. Most tests pass except a few in least_squares, data_manager and cg_solve...
File size: 117.2 KB

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

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

Download in other formats: