Context Navigation

test_data_manager.py @ 3975

Last change on this file since 3975 was 3975, checked in by duncan, 17 years ago
extend tests of urs2sww.
File size: 152.4 KB

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