Context Navigation

source: trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py @ 7776

Last change on this file since 7776 was 7776, checked in by hudson, 13 years ago
Removed redundant data_manager class. Unit tests are running, but may fail.
File size: 76.9 KB

Rev	Line
[7742]	1
	2	from anuga.shallow_water.shallow_water_domain import Domain
	3	from ferret2sww import ferret2sww
	4	from anuga.utilities.numerical_tools import ensure_numeric, mean
	5	from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
	6	import Transmissive_boundary
	7	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
	8	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
[7776]	9	from anuga.file.sww import SWW_file
[7770]	10	from anuga.file.sww import extent_sww
[7742]	11	from anuga.config import netcdf_float, epsilon, g
	12	from Scientific.IO.NetCDF import NetCDFFile
[7758]	13	from anuga.file_conversion.file_conversion import tsh2sww, \
	14	pmesh_to_domain_instance
	15
[7776]	16
[7770]	17	from anuga.file.mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
	18	NORTH_VELOCITY_LABEL
	19
[7742]	20	import sys
	21	import unittest
	22	import numpy as num
	23	import copy
	24	import os
	25
	26
	27	class Test_File_Conversion(unittest.TestCase):
	28	""" A suite of tests to test file conversion functions.
	29	These tests are quite coarse-grained: converting a file
	30	and checking that its headers and some of its contents
	31	are correct.
	32	"""
	33	verbose = False
	34
	35	def set_verbose(self):
	36	Test_File_Conversion.verbose = True
	37
	38	def setUp(self):
	39	import time
	40
	41	self.verbose = Test_File_Conversion.verbose
	42	# Create basic mesh
	43	points, vertices, boundary = rectangular(2, 2)
	44
	45	# Create shallow water domain
	46	domain = Domain(points, vertices, boundary)
	47	domain.default_order = 2
	48
	49	# Set some field values
	50	domain.set_quantity('elevation', lambda x,y: -x)
	51	domain.set_quantity('friction', 0.03)
	52
	53
	54	######################
	55	# Boundary conditions
	56	B = Transmissive_boundary(domain)
	57	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	58
	59
	60	######################
	61	#Initial condition - with jumps
	62	bed = domain.quantities['elevation'].vertex_values
	63	stage = num.zeros(bed.shape, num.float)
	64
	65	h = 0.3
	66	for i in range(stage.shape[0]):
	67	if i % 2 == 0:
	68	stage[i,:] = bed[i,:] + h
	69	else:
	70	stage[i,:] = bed[i,:]
	71
	72	domain.set_quantity('stage', stage)
	73
	74
	75	domain.distribute_to_vertices_and_edges()
	76	self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
	77
	78
	79	self.domain = domain
	80
	81	C = domain.get_vertex_coordinates()
	82	self.X = C[:,0:6:2].copy()
	83	self.Y = C[:,1:6:2].copy()
	84
	85	self.F = bed
	86
	87	#Write A testfile (not realistic. Values aren't realistic)
	88	self.test_MOST_file = 'most_small'
	89
	90	longitudes = [150.66667, 150.83334, 151., 151.16667]
	91	latitudes = [-34.5, -34.33333, -34.16667, -34]
	92
	93	long_name = 'LON'
	94	lat_name = 'LAT'
	95
	96	nx = 4
	97	ny = 4
	98	six = 6
	99
	100
	101	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
	102	fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
	103
	104	fid.createDimension(long_name,nx)
	105	fid.createVariable(long_name,netcdf_float,(long_name,))
	106	fid.variables[long_name].point_spacing='uneven'
	107	fid.variables[long_name].units='degrees_east'
	108	fid.variables[long_name].assignValue(longitudes)
	109
	110	fid.createDimension(lat_name,ny)
	111	fid.createVariable(lat_name,netcdf_float,(lat_name,))
	112	fid.variables[lat_name].point_spacing='uneven'
	113	fid.variables[lat_name].units='degrees_north'
	114	fid.variables[lat_name].assignValue(latitudes)
	115
	116	fid.createDimension('TIME',six)
	117	fid.createVariable('TIME',netcdf_float,('TIME',))
	118	fid.variables['TIME'].point_spacing='uneven'
	119	fid.variables['TIME'].units='seconds'
	120	fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
	121
	122
	123	name = ext[1:3].upper()
	124	if name == 'E.': name = 'ELEVATION'
	125	fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
	126	fid.variables[name].units='CENTIMETERS'
	127	fid.variables[name].missing_value=-1.e+034
	128
	129	fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
	130	[-0.1214216, 0, 0, 0],
	131	[0, 0, 0, 0],
	132	[0, 0, 0, 0]],
	133	[[0.3400644, 2.291054e-005, -23.33335, -6.50198],
	134	[-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
	135	[2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
	136	[0, 2.291054e-005, 2.291054e-005, 0]],
	137	[[0.3400644, 0.0001374632, -23.31503, -6.50198],
	138	[-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
	139	[0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
	140	[0, 0.0001374632, 0.0001374632, 0]],
	141	[[0.3400644, 0.0002520159, -23.29672, -6.50198],
	142	[-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
	143	[0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
	144	[0, 0.0002520159, 0.0002520159, 0]],
	145	[[0.3400644, 0.0003665686, -23.27842, -6.50198],
	146	[-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
	147	[0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
	148	[0, 0.0003665686, 0.0003665686, 0]],
	149	[[0.3400644, 0.0004811212, -23.26012, -6.50198],
	150	[-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
	151	[0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
	152	[0, 0.0004811212, 0.0004811212, 0]]])
	153
	154
	155	fid.close()
	156
	157
	158
	159
	160	def tearDown(self):
	161	import os
	162	for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
	163	#print 'Trying to remove', self.test_MOST_file + ext
	164	os.remove(self.test_MOST_file + ext)
	165
	166	def test_ferret2sww1(self):
	167	"""Test that georeferencing etc works when converting from
	168	ferret format (lat/lon) to sww format (UTM)
	169	"""
	170	from Scientific.IO.NetCDF import NetCDFFile
	171	import os, sys
	172
	173	#The test file has
	174	# LON = 150.66667, 150.83334, 151, 151.16667
	175	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	176	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
	177	#
	178	# First value (index=0) in small_ha.nc is 0.3400644 cm,
	179	# Fourth value (index==3) is -6.50198 cm
	180
	181
	182
	183	#Read
	184	from anuga.coordinate_transforms.redfearn import redfearn
	185	#fid = NetCDFFile(self.test_MOST_file)
	186	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	187	first_value = fid.variables['HA'][:][0,0,0]
	188	fourth_value = fid.variables['HA'][:][0,0,3]
	189	fid.close()
	190
	191
	192	#Call conversion (with zero origin)
	193	#ferret2sww('small', verbose=False,
	194	# origin = (56, 0, 0))
	195	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	196	origin = (56, 0, 0))
	197
	198	#Work out the UTM coordinates for first point
	199	zone, e, n = redfearn(-34.5, 150.66667)
	200	#print zone, e, n
	201
	202	#Read output file 'small.sww'
	203	#fid = NetCDFFile('small.sww')
	204	fid = NetCDFFile(self.test_MOST_file + '.sww')
	205
	206	x = fid.variables['x'][:]
	207	y = fid.variables['y'][:]
	208
	209	#Check that first coordinate is correctly represented
	210	assert num.allclose(x[0], e)
	211	assert num.allclose(y[0], n)
	212
	213	#Check first value
	214	stage = fid.variables['stage'][:]
	215	xmomentum = fid.variables['xmomentum'][:]
	216	ymomentum = fid.variables['ymomentum'][:]
	217
	218	#print ymomentum
	219
	220	assert num.allclose(stage[0,0], first_value/100) #Meters
	221
	222	#Check fourth value
	223	assert num.allclose(stage[0,3], fourth_value/100) #Meters
	224
	225	fid.close()
	226
	227	#Cleanup
	228	import os
	229	os.remove(self.test_MOST_file + '.sww')
	230
	231
	232
	233	def test_ferret2sww_zscale(self):
	234	"""Test that zscale workse
	235	"""
	236	from Scientific.IO.NetCDF import NetCDFFile
	237	import os, sys
	238
	239	#The test file has
	240	# LON = 150.66667, 150.83334, 151, 151.16667
	241	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	242	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
	243	#
	244	# First value (index=0) in small_ha.nc is 0.3400644 cm,
	245	# Fourth value (index==3) is -6.50198 cm
	246
	247
	248	#Read
	249	from anuga.coordinate_transforms.redfearn import redfearn
	250	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	251	first_value = fid.variables['HA'][:][0,0,0]
	252	fourth_value = fid.variables['HA'][:][0,0,3]
	253	fid.close()
	254
	255	#Call conversion (with no scaling)
	256	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	257	origin = (56, 0, 0))
	258
	259	#Work out the UTM coordinates for first point
	260	fid = NetCDFFile(self.test_MOST_file + '.sww')
	261
	262	#Check values
	263	stage_1 = fid.variables['stage'][:]
	264	xmomentum_1 = fid.variables['xmomentum'][:]
	265	ymomentum_1 = fid.variables['ymomentum'][:]
	266
	267	assert num.allclose(stage_1[0,0], first_value/100) #Meters
	268	assert num.allclose(stage_1[0,3], fourth_value/100) #Meters
	269
	270	fid.close()
	271
	272	#Call conversion (with scaling)
	273	ferret2sww(self.test_MOST_file,
	274	zscale = 5,
	275	verbose=self.verbose,
	276	origin = (56, 0, 0))
	277
	278	#Work out the UTM coordinates for first point
	279	fid = NetCDFFile(self.test_MOST_file + '.sww')
	280
	281	#Check values
	282	stage_5 = fid.variables['stage'][:]
	283	xmomentum_5 = fid.variables['xmomentum'][:]
	284	ymomentum_5 = fid.variables['ymomentum'][:]
	285	elevation = fid.variables['elevation'][:]
	286
	287	assert num.allclose(stage_5[0,0], 5*first_value/100) #Meters
	288	assert num.allclose(stage_5[0,3], 5*fourth_value/100) #Meters
	289
	290	assert num.allclose(5*stage_1, stage_5)
	291
	292	# Momentum will also be changed due to new depth
	293
	294	depth_1 = stage_1-elevation
	295	depth_5 = stage_5-elevation
	296
	297
	298	for i in range(stage_1.shape[0]):
	299	for j in range(stage_1.shape[1]):
	300	if depth_1[i,j] > epsilon:
	301
	302	scale = depth_5[i,j]/depth_1[i,j]
	303	ref_xmomentum = xmomentum_1[i,j] * scale
	304	ref_ymomentum = ymomentum_1[i,j] * scale
	305
	306	#print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
	307
	308	assert num.allclose(xmomentum_5[i,j], ref_xmomentum)
	309	assert num.allclose(ymomentum_5[i,j], ref_ymomentum)
	310
	311
	312
	313	fid.close()
	314
	315
	316	#Cleanup
	317	import os
	318	os.remove(self.test_MOST_file + '.sww')
	319
	320
	321
	322	def test_ferret2sww_2(self):
	323	"""Test that georeferencing etc works when converting from
	324	ferret format (lat/lon) to sww format (UTM)
	325	"""
	326	from Scientific.IO.NetCDF import NetCDFFile
	327
	328	#The test file has
	329	# LON = 150.66667, 150.83334, 151, 151.16667
	330	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	331	# TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
	332	#
	333	# First value (index=0) in small_ha.nc is 0.3400644 cm,
	334	# Fourth value (index==3) is -6.50198 cm
	335
	336
	337	from anuga.coordinate_transforms.redfearn import redfearn
	338	#fid = NetCDFFile('small_ha.nc')
	339	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	340
	341	#Pick a coordinate and a value
	342
	343	time_index = 1
	344	lat_index = 0
	345	lon_index = 2
	346
	347	test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
	348	test_time = fid.variables['TIME'][:][time_index]
	349	test_lat = fid.variables['LAT'][:][lat_index]
	350	test_lon = fid.variables['LON'][:][lon_index]
	351
	352	linear_point_index = lat_index*4 + lon_index
	353	fid.close()
	354
	355	#Call conversion (with zero origin)
	356	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	357	origin = (56, 0, 0))
	358
	359
	360	#Work out the UTM coordinates for test point
	361	zone, e, n = redfearn(test_lat, test_lon)
	362
	363	#Read output file 'small.sww'
	364	fid = NetCDFFile(self.test_MOST_file + '.sww')
	365
	366	x = fid.variables['x'][:]
	367	y = fid.variables['y'][:]
	368
	369	#Check that test coordinate is correctly represented
	370	assert num.allclose(x[linear_point_index], e)
	371	assert num.allclose(y[linear_point_index], n)
	372
	373	#Check test value
	374	stage = fid.variables['stage'][:]
	375
	376	assert num.allclose(stage[time_index, linear_point_index], test_value/100)
	377
	378	fid.close()
	379
	380	#Cleanup
	381	import os
	382	os.remove(self.test_MOST_file + '.sww')
	383
	384
	385	def test_ferret2sww_lat_long(self):
	386	# Test that min lat long works
	387
	388	#The test file has
	389	# LON = 150.66667, 150.83334, 151, 151.16667
	390	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	391
	392	#Read
	393	from anuga.coordinate_transforms.redfearn import redfearn
	394	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	395	first_value = fid.variables['HA'][:][0,0,0]
	396	fourth_value = fid.variables['HA'][:][0,0,3]
	397	fid.close()
	398
	399
	400	#Call conversion (with zero origin)
	401	#ferret2sww('small', verbose=self.verbose,
	402	# origin = (56, 0, 0))
	403	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	404	origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
	405
	406	#Work out the UTM coordinates for first point
	407	zone, e, n = redfearn(-34.5, 150.66667)
	408	#print zone, e, n
	409
	410	#Read output file 'small.sww'
	411	#fid = NetCDFFile('small.sww')
	412	fid = NetCDFFile(self.test_MOST_file + '.sww')
	413
	414	x = fid.variables['x'][:]
	415	y = fid.variables['y'][:]
	416	#Check that first coordinate is correctly represented
	417	assert 16 == len(x)
	418
	419	fid.close()
	420
	421	#Cleanup
	422	import os
	423	os.remove(self.test_MOST_file + '.sww')
	424
	425
	426	def test_ferret2sww_lat_longII(self):
	427	# Test that min lat long works
	428
	429	#The test file has
	430	# LON = 150.66667, 150.83334, 151, 151.16667
	431	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	432
	433	#Read
	434	from anuga.coordinate_transforms.redfearn import redfearn
	435	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	436	first_value = fid.variables['HA'][:][0,0,0]
	437	fourth_value = fid.variables['HA'][:][0,0,3]
	438	fid.close()
	439
	440
	441	#Call conversion (with zero origin)
	442	#ferret2sww('small', verbose=False,
	443	# origin = (56, 0, 0))
	444	ferret2sww(self.test_MOST_file, verbose=False,
	445	origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
	446
	447	#Work out the UTM coordinates for first point
	448	zone, e, n = redfearn(-34.5, 150.66667)
	449	#print zone, e, n
	450
	451	#Read output file 'small.sww'
	452	#fid = NetCDFFile('small.sww')
	453	fid = NetCDFFile(self.test_MOST_file + '.sww')
	454
	455	x = fid.variables['x'][:]
	456	y = fid.variables['y'][:]
	457	#Check that first coordinate is correctly represented
	458	assert 12 == len(x)
	459
	460	fid.close()
	461
	462	#Cleanup
	463	import os
	464	os.remove(self.test_MOST_file + '.sww')
	465
	466	def test_ferret2sww3(self):
	467	"""Elevation included
	468	"""
	469	from Scientific.IO.NetCDF import NetCDFFile
	470
	471	#The test file has
	472	# LON = 150.66667, 150.83334, 151, 151.16667
	473	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	474	# ELEVATION = [-1 -2 -3 -4
	475	# -5 -6 -7 -8
	476	# ...
	477	# ... -16]
	478	# where the top left corner is -1m,
	479	# and the ll corner is -13.0m
	480	#
	481	# First value (index=0) in small_ha.nc is 0.3400644 cm,
	482	# Fourth value (index==3) is -6.50198 cm
	483
	484	from anuga.coordinate_transforms.redfearn import redfearn
	485	import os
	486	fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
	487	fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
	488	fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
	489	fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
	490
	491	h1_list = [150.66667,150.83334,151.]
	492	h2_list = [-34.5,-34.33333]
	493
	494	long_name = 'LON'
	495	lat_name = 'LAT'
	496	time_name = 'TIME'
	497
	498	nx = 3
	499	ny = 2
	500
	501	for fid in [fid1,fid2,fid3]:
	502	fid.createDimension(long_name,nx)
	503	fid.createVariable(long_name,netcdf_float,(long_name,))
	504	fid.variables[long_name].point_spacing='uneven'
	505	fid.variables[long_name].units='degrees_east'
	506	fid.variables[long_name].assignValue(h1_list)
	507
	508	fid.createDimension(lat_name,ny)
	509	fid.createVariable(lat_name,netcdf_float,(lat_name,))
	510	fid.variables[lat_name].point_spacing='uneven'
	511	fid.variables[lat_name].units='degrees_north'
	512	fid.variables[lat_name].assignValue(h2_list)
	513
	514	fid.createDimension(time_name,2)
	515	fid.createVariable(time_name,netcdf_float,(time_name,))
	516	fid.variables[time_name].point_spacing='uneven'
	517	fid.variables[time_name].units='seconds'
	518	fid.variables[time_name].assignValue([0.,1.])
	519	#if fid == fid3: break
	520
	521
	522	for fid in [fid4]:
	523	fid.createDimension(long_name,nx)
	524	fid.createVariable(long_name,netcdf_float,(long_name,))
	525	fid.variables[long_name].point_spacing='uneven'
	526	fid.variables[long_name].units='degrees_east'
	527	fid.variables[long_name].assignValue(h1_list)
	528
	529	fid.createDimension(lat_name,ny)
	530	fid.createVariable(lat_name,netcdf_float,(lat_name,))
	531	fid.variables[lat_name].point_spacing='uneven'
	532	fid.variables[lat_name].units='degrees_north'
	533	fid.variables[lat_name].assignValue(h2_list)
	534
	535	name = {}
	536	name[fid1]='HA'
	537	name[fid2]='UA'
	538	name[fid3]='VA'
	539	name[fid4]='ELEVATION'
	540
	541	units = {}
	542	units[fid1]='cm'
	543	units[fid2]='cm/s'
	544	units[fid3]='cm/s'
	545	units[fid4]='m'
	546
	547	values = {}
	548	values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
	549	values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
	550	values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
	551	values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
	552
	553	for fid in [fid1,fid2,fid3]:
	554	fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
	555	fid.variables[name[fid]].point_spacing='uneven'
	556	fid.variables[name[fid]].units=units[fid]
	557	fid.variables[name[fid]].assignValue(values[fid])
	558	fid.variables[name[fid]].missing_value = -99999999.
	559	#if fid == fid3: break
	560
	561	for fid in [fid4]:
	562	fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
	563	fid.variables[name[fid]].point_spacing='uneven'
	564	fid.variables[name[fid]].units=units[fid]
	565	fid.variables[name[fid]].assignValue(values[fid])
	566	fid.variables[name[fid]].missing_value = -99999999.
	567
	568
	569	fid1.sync(); fid1.close()
	570	fid2.sync(); fid2.close()
	571	fid3.sync(); fid3.close()
	572	fid4.sync(); fid4.close()
	573
	574	fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
	575	fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
	576	fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
	577
	578
	579	first_amp = fid1.variables['HA'][:][0,0,0]
	580	third_amp = fid1.variables['HA'][:][0,0,2]
	581	first_elevation = fid2.variables['ELEVATION'][0,0]
	582	third_elevation= fid2.variables['ELEVATION'][:][0,2]
	583	first_speed = fid3.variables['VA'][0,0,0]
	584	third_speed = fid3.variables['VA'][:][0,0,2]
	585
	586	fid1.close()
	587	fid2.close()
	588	fid3.close()
	589
	590	#Call conversion (with zero origin)
	591	ferret2sww('test', verbose=self.verbose,
	592	origin = (56, 0, 0), inverted_bathymetry=False)
	593
	594	os.remove('test_va.nc')
	595	os.remove('test_ua.nc')
	596	os.remove('test_ha.nc')
	597	os.remove('test_e.nc')
	598
	599	#Read output file 'test.sww'
	600	fid = NetCDFFile('test.sww')
	601
	602
	603	#Check first value
	604	elevation = fid.variables['elevation'][:]
	605	stage = fid.variables['stage'][:]
	606	xmomentum = fid.variables['xmomentum'][:]
	607	ymomentum = fid.variables['ymomentum'][:]
	608
	609	#print ymomentum
	610	first_height = first_amp/100 - first_elevation
	611	third_height = third_amp/100 - third_elevation
	612	first_momentum=first_speed*first_height/100
	613	third_momentum=third_speed*third_height/100
	614
	615	assert num.allclose(ymomentum[0][0],first_momentum) #Meters
	616	assert num.allclose(ymomentum[0][2],third_momentum) #Meters
	617
	618	fid.close()
	619
	620	#Cleanup
	621	os.remove('test.sww')
	622
	623
	624
	625	def test_ferret2sww4(self):
	626	"""Like previous but with augmented variable names as
	627	in files produced by ferret as opposed to MOST
	628	"""
	629	from Scientific.IO.NetCDF import NetCDFFile
	630
	631	#The test file has
	632	# LON = 150.66667, 150.83334, 151, 151.16667
	633	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	634	# ELEVATION = [-1 -2 -3 -4
	635	# -5 -6 -7 -8
	636	# ...
	637	# ... -16]
	638	# where the top left corner is -1m,
	639	# and the ll corner is -13.0m
	640	#
	641	# First value (index=0) in small_ha.nc is 0.3400644 cm,
	642	# Fourth value (index==3) is -6.50198 cm
	643
	644	from anuga.coordinate_transforms.redfearn import redfearn
	645	import os
	646	fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
	647	fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
	648	fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
	649	fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
	650
	651	h1_list = [150.66667,150.83334,151.]
	652	h2_list = [-34.5,-34.33333]
	653
	654	# long_name = 'LON961_1261'
	655	# lat_name = 'LAT481_841'
	656	# time_name = 'TIME1'
	657
	658	long_name = 'LON'
	659	lat_name = 'LAT'
	660	time_name = 'TIME'
	661
	662	nx = 3
	663	ny = 2
	664
	665	for fid in [fid1,fid2,fid3]:
	666	fid.createDimension(long_name,nx)
	667	fid.createVariable(long_name,netcdf_float,(long_name,))
	668	fid.variables[long_name].point_spacing='uneven'
	669	fid.variables[long_name].units='degrees_east'
	670	fid.variables[long_name].assignValue(h1_list)
	671
	672	fid.createDimension(lat_name,ny)
	673	fid.createVariable(lat_name,netcdf_float,(lat_name,))
	674	fid.variables[lat_name].point_spacing='uneven'
	675	fid.variables[lat_name].units='degrees_north'
	676	fid.variables[lat_name].assignValue(h2_list)
	677
	678	fid.createDimension(time_name,2)
	679	fid.createVariable(time_name,netcdf_float,(time_name,))
	680	fid.variables[time_name].point_spacing='uneven'
	681	fid.variables[time_name].units='seconds'
	682	fid.variables[time_name].assignValue([0.,1.])
	683	#if fid == fid3: break
	684
	685
	686	for fid in [fid4]:
	687	fid.createDimension(long_name,nx)
	688	fid.createVariable(long_name,netcdf_float,(long_name,))
	689	fid.variables[long_name].point_spacing='uneven'
	690	fid.variables[long_name].units='degrees_east'
	691	fid.variables[long_name].assignValue(h1_list)
	692
	693	fid.createDimension(lat_name,ny)
	694	fid.createVariable(lat_name,netcdf_float,(lat_name,))
	695	fid.variables[lat_name].point_spacing='uneven'
	696	fid.variables[lat_name].units='degrees_north'
	697	fid.variables[lat_name].assignValue(h2_list)
	698
	699	name = {}
	700	name[fid1]='HA'
	701	name[fid2]='UA'
	702	name[fid3]='VA'
	703	name[fid4]='ELEVATION'
	704
	705	units = {}
	706	units[fid1]='cm'
	707	units[fid2]='cm/s'
	708	units[fid3]='cm/s'
	709	units[fid4]='m'
	710
	711	values = {}
	712	values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
	713	values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
	714	values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
	715	values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
	716
	717	for fid in [fid1,fid2,fid3]:
	718	fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
	719	fid.variables[name[fid]].point_spacing='uneven'
	720	fid.variables[name[fid]].units=units[fid]
	721	fid.variables[name[fid]].assignValue(values[fid])
	722	fid.variables[name[fid]].missing_value = -99999999.
	723	#if fid == fid3: break
	724
	725	for fid in [fid4]:
	726	fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
	727	fid.variables[name[fid]].point_spacing='uneven'
	728	fid.variables[name[fid]].units=units[fid]
	729	fid.variables[name[fid]].assignValue(values[fid])
	730	fid.variables[name[fid]].missing_value = -99999999.
	731
	732
	733	fid1.sync(); fid1.close()
	734	fid2.sync(); fid2.close()
	735	fid3.sync(); fid3.close()
	736	fid4.sync(); fid4.close()
	737
	738	fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
	739	fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
	740	fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
	741
	742
	743	first_amp = fid1.variables['HA'][:][0,0,0]
	744	third_amp = fid1.variables['HA'][:][0,0,2]
	745	first_elevation = fid2.variables['ELEVATION'][0,0]
	746	third_elevation= fid2.variables['ELEVATION'][:][0,2]
	747	first_speed = fid3.variables['VA'][0,0,0]
	748	third_speed = fid3.variables['VA'][:][0,0,2]
	749
	750	fid1.close()
	751	fid2.close()
	752	fid3.close()
	753
	754	#Call conversion (with zero origin)
	755	ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
	756	, inverted_bathymetry=False)
	757
	758	os.remove('test_va.nc')
	759	os.remove('test_ua.nc')
	760	os.remove('test_ha.nc')
	761	os.remove('test_e.nc')
	762
	763	#Read output file 'test.sww'
	764	fid = NetCDFFile('test.sww')
	765
	766
	767	#Check first value
	768	elevation = fid.variables['elevation'][:]
	769	stage = fid.variables['stage'][:]
	770	xmomentum = fid.variables['xmomentum'][:]
	771	ymomentum = fid.variables['ymomentum'][:]
	772
	773	#print ymomentum
	774	first_height = first_amp/100 - first_elevation
	775	third_height = third_amp/100 - third_elevation
	776	first_momentum=first_speed*first_height/100
	777	third_momentum=third_speed*third_height/100
	778
	779	assert num.allclose(ymomentum[0][0],first_momentum) #Meters
	780	assert num.allclose(ymomentum[0][2],third_momentum) #Meters
	781
	782	fid.close()
	783
	784	#Cleanup
	785	os.remove('test.sww')
	786
	787
	788
	789
	790	def test_ferret2sww_nz_origin(self):
	791	from Scientific.IO.NetCDF import NetCDFFile
	792	from anuga.coordinate_transforms.redfearn import redfearn
	793
	794	#Call conversion (with nonzero origin)
	795	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	796	origin = (56, 100000, 200000))
	797
	798
	799	#Work out the UTM coordinates for first point
	800	zone, e, n = redfearn(-34.5, 150.66667)
	801
	802	#Read output file 'small.sww'
	803	#fid = NetCDFFile('small.sww', netcdf_mode_r)
	804	fid = NetCDFFile(self.test_MOST_file + '.sww')
	805
	806	x = fid.variables['x'][:]
	807	y = fid.variables['y'][:]
	808
	809	#Check that first coordinate is correctly represented
	810	assert num.allclose(x[0], e-100000)
	811	assert num.allclose(y[0], n-200000)
	812
	813	fid.close()
	814
	815	#Cleanup
	816	os.remove(self.test_MOST_file + '.sww')
	817
	818
	819	def test_ferret2sww_lat_longII(self):
	820	# Test that min lat long works
	821
	822	#The test file has
	823	# LON = 150.66667, 150.83334, 151, 151.16667
	824	# LAT = -34.5, -34.33333, -34.16667, -34 ;
	825
	826	#Read
	827	from anuga.coordinate_transforms.redfearn import redfearn
	828	fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
	829	first_value = fid.variables['HA'][:][0,0,0]
	830	fourth_value = fid.variables['HA'][:][0,0,3]
	831	fid.close()
	832
	833
	834	#Call conversion (with zero origin)
	835	#ferret2sww('small', verbose=self.verbose,
	836	# origin = (56, 0, 0))
	837	try:
	838	ferret2sww(self.test_MOST_file, verbose=self.verbose,
	839	origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
	840	except AssertionError:
	841	pass
	842	else:
	843	self.failUnless(0 ==1, 'Bad input did not throw exception error!')
	844
	845	def test_sww_extent(self):
	846	"""Not a test, rather a look at the sww format
	847	"""
	848
	849	import time, os
	850	from Scientific.IO.NetCDF import NetCDFFile
	851
	852	self.domain.set_name('datatest' + str(id(self)))
	853	self.domain.format = 'sww'
	854	self.domain.smooth = True
	855	self.domain.reduction = mean
	856	self.domain.set_datadir('.')
	857	#self.domain.tight_slope_limiters = 1
	858
	859
	860	sww = SWW_file(self.domain)
	861	sww.store_connectivity()
	862	sww.store_timestep()
	863	self.domain.time = 2.
	864
	865	#Modify stage at second timestep
	866	stage = self.domain.quantities['stage'].vertex_values
	867	self.domain.set_quantity('stage', stage/2)
	868
	869	sww.store_timestep()
	870
	871	file_and_extension_name = self.domain.get_name() + ".sww"
	872	#print "file_and_extension_name",file_and_extension_name
	873	[xmin, xmax, ymin, ymax, stagemin, stagemax] = \
	874	extent_sww(file_and_extension_name )
	875
	876	assert num.allclose(xmin, 0.0)
	877	assert num.allclose(xmax, 1.0)
	878	assert num.allclose(ymin, 0.0)
	879	assert num.allclose(ymax, 1.0)
	880
	881	# FIXME (Ole): Revisit these numbers
	882	#assert num.allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
	883	#assert num.allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
	884
	885
	886	#Cleanup
	887	os.remove(sww.filename)
	888
	889
	890
	891	#-------------------------------------------------------------
	892
	893	if __name__ == "__main__":
	894	suite = unittest.makeSuite(Test_File_Conversion, 'test_sww')
	895
	896	# FIXME(Ole): When Ross has implemented logging, we can
	897	# probably get rid of all this:
	898	if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
	899	Test_File_Conversion.verbose=True
	900	saveout = sys.stdout
	901	filename = ".temp_verbose"
	902	fid = open(filename, 'w')
	903	sys.stdout = fid
	904	else:
	905	pass
	906	runner = unittest.TextTestRunner() #verbosity=2)
	907	runner.run(suite)
	908
	909	# Cleaning up
	910	if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
	911	sys.stdout = saveout
	912	fid.close()
	913	os.remove(filename)
	914
	915
	916	def test_asc_csiro2sww(self):
	917	import tempfile
	918
	919	bath_dir = tempfile.mkdtemp()
	920	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
	921	#bath_dir = 'bath_data_manager_test'
	922	#print "os.getcwd( )",os.getcwd( )
	923	elevation_dir = tempfile.mkdtemp()
	924	#elevation_dir = 'elev_expanded'
	925	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
	926	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
	927
	928	fid = open(bath_dir_filename, 'w')
	929	fid.write(""" ncols 3
	930	nrows 2
	931	xllcorner 148.00000
	932	yllcorner -38.00000
	933	cellsize 0.25
	934	nodata_value -9999.0
	935	9000.000 -1000.000 3000.0
	936	-1000.000 9000.000 -1000.000
	937	""")
	938	fid.close()
	939
	940	fid = open(elevation_dir_filename1, 'w')
	941	fid.write(""" ncols 3
	942	nrows 2
	943	xllcorner 148.00000
	944	yllcorner -38.00000
	945	cellsize 0.25
	946	nodata_value -9999.0
	947	9000.000 0.000 3000.0
	948	0.000 9000.000 0.000
	949	""")
	950	fid.close()
	951
	952	fid = open(elevation_dir_filename2, 'w')
	953	fid.write(""" ncols 3
	954	nrows 2
	955	xllcorner 148.00000
	956	yllcorner -38.00000
	957	cellsize 0.25
	958	nodata_value -9999.0
	959	9000.000 4000.000 4000.0
	960	4000.000 9000.000 4000.000
	961	""")
	962	fid.close()
	963
	964	ucur_dir = tempfile.mkdtemp()
	965	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
	966	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
	967
	968	fid = open(ucur_dir_filename1, 'w')
	969	fid.write(""" ncols 3
	970	nrows 2
	971	xllcorner 148.00000
	972	yllcorner -38.00000
	973	cellsize 0.25
	974	nodata_value -9999.0
	975	90.000 60.000 30.0
	976	10.000 10.000 10.000
	977	""")
	978	fid.close()
	979	fid = open(ucur_dir_filename2, 'w')
	980	fid.write(""" ncols 3
	981	nrows 2
	982	xllcorner 148.00000
	983	yllcorner -38.00000
	984	cellsize 0.25
	985	nodata_value -9999.0
	986	90.000 60.000 30.0
	987	10.000 10.000 10.000
	988	""")
	989	fid.close()
	990
	991	vcur_dir = tempfile.mkdtemp()
	992	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
	993	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
	994
	995	fid = open(vcur_dir_filename1, 'w')
	996	fid.write(""" ncols 3
	997	nrows 2
	998	xllcorner 148.00000
	999	yllcorner -38.00000
	1000	cellsize 0.25
	1001	nodata_value -9999.0
	1002	90.000 60.000 30.0
	1003	10.000 10.000 10.000
	1004	""")
	1005	fid.close()
	1006	fid = open(vcur_dir_filename2, 'w')
	1007	fid.write(""" ncols 3
	1008	nrows 2
	1009	xllcorner 148.00000
	1010	yllcorner -38.00000
	1011	cellsize 0.25
	1012	nodata_value -9999.0
	1013	90.000 60.000 30.0
	1014	10.000 10.000 10.000
	1015	""")
	1016	fid.close()
	1017
	1018	sww_file = 'a_test.sww'
	1019	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
	1020	verbose=self.verbose)
	1021
	1022	# check the sww file
	1023
	1024	fid = NetCDFFile(sww_file, netcdf_mode_r) # Open existing file for read
	1025	x = fid.variables['x'][:]
	1026	y = fid.variables['y'][:]
	1027	z = fid.variables['elevation'][:]
	1028	stage = fid.variables['stage'][:]
	1029	xmomentum = fid.variables['xmomentum'][:]
	1030	geo_ref = Geo_reference(NetCDFObject=fid)
	1031	#print "geo_ref",geo_ref
	1032	x_ref = geo_ref.get_xllcorner()
	1033	y_ref = geo_ref.get_yllcorner()
	1034	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
	1035	assert num.allclose(x_ref, 587798.418) # (-38, 148)
	1036	assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
	1037
	1038	#Zone: 55
	1039	#Easting: 588095.674 Northing: 5821451.722
	1040	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
	1041	assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
	1042
	1043	#Zone: 55
	1044	#Easting: 632145.632 Northing: 5820863.269
	1045	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
	1046	assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
	1047
	1048	#Zone: 55
	1049	#Easting: 609748.788 Northing: 5793447.860
	1050	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
	1051	assert num.allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
	1052
	1053	assert num.allclose(z[0],9000.0 )
	1054	assert num.allclose(stage[0][1],0.0 )
	1055
	1056	#(4000+1000)*60
	1057	assert num.allclose(xmomentum[1][1],300000.0 )
	1058
	1059
	1060	fid.close()
	1061
	1062	#tidy up
	1063	os.remove(bath_dir_filename)
	1064	os.rmdir(bath_dir)
	1065
	1066	os.remove(elevation_dir_filename1)
	1067	os.remove(elevation_dir_filename2)
	1068	os.rmdir(elevation_dir)
	1069
	1070	os.remove(ucur_dir_filename1)
	1071	os.remove(ucur_dir_filename2)
	1072	os.rmdir(ucur_dir)
	1073
	1074	os.remove(vcur_dir_filename1)
	1075	os.remove(vcur_dir_filename2)
	1076	os.rmdir(vcur_dir)
	1077
	1078
	1079	# remove sww file
	1080	os.remove(sww_file)
	1081
	1082	def test_asc_csiro2sww2(self):
	1083	import tempfile
	1084
	1085	bath_dir = tempfile.mkdtemp()
	1086	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
	1087	#bath_dir = 'bath_data_manager_test'
	1088	#print "os.getcwd( )",os.getcwd( )
	1089	elevation_dir = tempfile.mkdtemp()
	1090	#elevation_dir = 'elev_expanded'
	1091	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
	1092	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
	1093
	1094	fid = open(bath_dir_filename, 'w')
	1095	fid.write(""" ncols 3
	1096	nrows 2
	1097	xllcorner 148.00000
	1098	yllcorner -38.00000
	1099	cellsize 0.25
	1100	nodata_value -9999.0
	1101	9000.000 -1000.000 3000.0
	1102	-1000.000 9000.000 -1000.000
	1103	""")
	1104	fid.close()
	1105
	1106	fid = open(elevation_dir_filename1, 'w')
	1107	fid.write(""" ncols 3
	1108	nrows 2
	1109	xllcorner 148.00000
	1110	yllcorner -38.00000
	1111	cellsize 0.25
	1112	nodata_value -9999.0
	1113	9000.000 0.000 3000.0
	1114	0.000 -9999.000 -9999.000
	1115	""")
	1116	fid.close()
	1117
	1118	fid = open(elevation_dir_filename2, 'w')
	1119	fid.write(""" ncols 3
	1120	nrows 2
	1121	xllcorner 148.00000
	1122	yllcorner -38.00000
	1123	cellsize 0.25
	1124	nodata_value -9999.0
	1125	9000.000 4000.000 4000.0
	1126	4000.000 9000.000 4000.000
	1127	""")
	1128	fid.close()
	1129
	1130	ucur_dir = tempfile.mkdtemp()
	1131	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
	1132	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
	1133
	1134	fid = open(ucur_dir_filename1, 'w')
	1135	fid.write(""" ncols 3
	1136	nrows 2
	1137	xllcorner 148.00000
	1138	yllcorner -38.00000
	1139	cellsize 0.25
	1140	nodata_value -9999.0
	1141	90.000 60.000 30.0
	1142	10.000 10.000 10.000
	1143	""")
	1144	fid.close()
	1145	fid = open(ucur_dir_filename2, 'w')
	1146	fid.write(""" ncols 3
	1147	nrows 2
	1148	xllcorner 148.00000
	1149	yllcorner -38.00000
	1150	cellsize 0.25
	1151	nodata_value -9999.0
	1152	90.000 60.000 30.0
	1153	10.000 10.000 10.000
	1154	""")
	1155	fid.close()
	1156
	1157	vcur_dir = tempfile.mkdtemp()
	1158	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
	1159	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
	1160
	1161	fid = open(vcur_dir_filename1, 'w')
	1162	fid.write(""" ncols 3
	1163	nrows 2
	1164	xllcorner 148.00000
	1165	yllcorner -38.00000
	1166	cellsize 0.25
	1167	nodata_value -9999.0
	1168	90.000 60.000 30.0
	1169	10.000 10.000 10.000
	1170	""")
	1171	fid.close()
	1172	fid = open(vcur_dir_filename2, 'w')
	1173	fid.write(""" ncols 3
	1174	nrows 2
	1175	xllcorner 148.00000
	1176	yllcorner -38.00000
	1177	cellsize 0.25
	1178	nodata_value -9999.0
	1179	90.000 60.000 30.0
	1180	10.000 10.000 10.000
	1181	""")
	1182	fid.close()
	1183
	1184	try:
	1185	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
	1186	vcur_dir, sww_file,
	1187	verbose=self.verbose)
	1188	except:
	1189	#tidy up
	1190	os.remove(bath_dir_filename)
	1191	os.rmdir(bath_dir)
	1192
	1193	os.remove(elevation_dir_filename1)
	1194	os.remove(elevation_dir_filename2)
	1195	os.rmdir(elevation_dir)
	1196
	1197	os.remove(ucur_dir_filename1)
	1198	os.remove(ucur_dir_filename2)
	1199	os.rmdir(ucur_dir)
	1200
	1201	os.remove(vcur_dir_filename1)
	1202	os.remove(vcur_dir_filename2)
	1203	os.rmdir(vcur_dir)
	1204	else:
	1205	#tidy up
	1206	os.remove(bath_dir_filename)
	1207	os.rmdir(bath_dir)
	1208
	1209	os.remove(elevation_dir_filename1)
	1210	os.remove(elevation_dir_filename2)
	1211	os.rmdir(elevation_dir)
	1212	raise 'Should raise exception'
	1213
	1214	os.remove(ucur_dir_filename1)
	1215	os.remove(ucur_dir_filename2)
	1216	os.rmdir(ucur_dir)
	1217
	1218	os.remove(vcur_dir_filename1)
	1219	os.remove(vcur_dir_filename2)
	1220	os.rmdir(vcur_dir)
	1221
	1222
	1223
	1224	def test_asc_csiro2sww3(self):
	1225	import tempfile
	1226
	1227	bath_dir = tempfile.mkdtemp()
	1228	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
	1229	#bath_dir = 'bath_data_manager_test'
	1230	#print "os.getcwd( )",os.getcwd( )
	1231	elevation_dir = tempfile.mkdtemp()
	1232	#elevation_dir = 'elev_expanded'
	1233	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
	1234	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
	1235
	1236	fid = open(bath_dir_filename, 'w')
	1237	fid.write(""" ncols 3
	1238	nrows 2
	1239	xllcorner 148.00000
	1240	yllcorner -38.00000
	1241	cellsize 0.25
	1242	nodata_value -9999.0
	1243	9000.000 -1000.000 3000.0
	1244	-1000.000 9000.000 -1000.000
	1245	""")
	1246	fid.close()
	1247
	1248	fid = open(elevation_dir_filename1, 'w')
	1249	fid.write(""" ncols 3
	1250	nrows 2
	1251	xllcorner 148.00000
	1252	yllcorner -38.00000
	1253	cellsize 0.25
	1254	nodata_value -9999.0
	1255	9000.000 0.000 3000.0
	1256	0.000 -9999.000 -9999.000
	1257	""")
	1258	fid.close()
	1259
	1260	fid = open(elevation_dir_filename2, 'w')
	1261	fid.write(""" ncols 3
	1262	nrows 2
	1263	xllcorner 148.00000
	1264	yllcorner -38.00000
	1265	cellsize 0.25
	1266	nodata_value -9999.0
	1267	9000.000 4000.000 4000.0
	1268	4000.000 9000.000 4000.000
	1269	""")
	1270	fid.close()
	1271
	1272	ucur_dir = tempfile.mkdtemp()
	1273	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
	1274	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
	1275
	1276	fid = open(ucur_dir_filename1, 'w')
	1277	fid.write(""" ncols 3
	1278	nrows 2
	1279	xllcorner 148.00000
	1280	yllcorner -38.00000
	1281	cellsize 0.25
	1282	nodata_value -9999.0
	1283	90.000 60.000 30.0
	1284	10.000 10.000 10.000
	1285	""")
	1286	fid.close()
	1287	fid = open(ucur_dir_filename2, 'w')
	1288	fid.write(""" ncols 3
	1289	nrows 2
	1290	xllcorner 148.00000
	1291	yllcorner -38.00000
	1292	cellsize 0.25
	1293	nodata_value -9999.0
	1294	90.000 60.000 30.0
	1295	10.000 10.000 10.000
	1296	""")
	1297	fid.close()
	1298
	1299	vcur_dir = tempfile.mkdtemp()
	1300	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
	1301	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
	1302
	1303	fid = open(vcur_dir_filename1, 'w')
	1304	fid.write(""" ncols 3
	1305	nrows 2
	1306	xllcorner 148.00000
	1307	yllcorner -38.00000
	1308	cellsize 0.25
	1309	nodata_value -9999.0
	1310	90.000 60.000 30.0
	1311	10.000 10.000 10.000
	1312	""")
	1313	fid.close()
	1314	fid = open(vcur_dir_filename2, 'w')
	1315	fid.write(""" ncols 3
	1316	nrows 2
	1317	xllcorner 148.00000
	1318	yllcorner -38.00000
	1319	cellsize 0.25
	1320	nodata_value -9999.0
	1321	90.000 60.000 30.0
	1322	10.000 10.000 10.000
	1323	""")
	1324	fid.close()
	1325
	1326	sww_file = 'a_test.sww'
	1327	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
	1328	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
	1329	mean_stage = 100,
	1330	verbose=self.verbose)
	1331
	1332	# check the sww file
	1333
	1334	fid = NetCDFFile(sww_file, netcdf_mode_r) # Open existing file for read
	1335	x = fid.variables['x'][:]
	1336	y = fid.variables['y'][:]
	1337	z = fid.variables['elevation'][:]
	1338	stage = fid.variables['stage'][:]
	1339	xmomentum = fid.variables['xmomentum'][:]
	1340	geo_ref = Geo_reference(NetCDFObject=fid)
	1341	#print "geo_ref",geo_ref
	1342	x_ref = geo_ref.get_xllcorner()
	1343	y_ref = geo_ref.get_yllcorner()
	1344	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
	1345	assert num.allclose(x_ref, 587798.418) # (-38, 148)
	1346	assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
	1347
	1348	#Zone: 55
	1349	#Easting: 588095.674 Northing: 5821451.722
	1350	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 0 ' 0.00000 ''
	1351	assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
	1352
	1353	#Zone: 55
	1354	#Easting: 632145.632 Northing: 5820863.269
	1355	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
	1356	assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
	1357
	1358	#Zone: 55
	1359	#Easting: 609748.788 Northing: 5793447.860
	1360	#Latitude: -38 0 ' 0.00000 '' Longitude: 148 15 ' 0.00000 ''
	1361	assert num.allclose((x[4],y[4]), (609748.788 - x_ref, 5793447.86 - y_ref))
	1362
	1363	assert num.allclose(z[0],9000.0 )
	1364	assert num.allclose(stage[0][4],100.0 )
	1365	assert num.allclose(stage[0][5],100.0 )
	1366
	1367	#(100.0 - 9000)*10
	1368	assert num.allclose(xmomentum[0][4], -89000.0 )
	1369
	1370	#(100.0 - -1000.000)*10
	1371	assert num.allclose(xmomentum[0][5], 11000.0 )
	1372
	1373	fid.close()
	1374
	1375	#tidy up
	1376	os.remove(bath_dir_filename)
	1377	os.rmdir(bath_dir)
	1378
	1379	os.remove(elevation_dir_filename1)
	1380	os.remove(elevation_dir_filename2)
	1381	os.rmdir(elevation_dir)
	1382
	1383	os.remove(ucur_dir_filename1)
	1384	os.remove(ucur_dir_filename2)
	1385	os.rmdir(ucur_dir)
	1386
	1387	os.remove(vcur_dir_filename1)
	1388	os.remove(vcur_dir_filename2)
	1389	os.rmdir(vcur_dir)
	1390
	1391	# remove sww file
	1392	os.remove(sww_file)
	1393
	1394
	1395	def test_asc_csiro2sww4(self):
	1396	"""
	1397	Test specifying the extent
	1398	"""
	1399
	1400	import tempfile
	1401
	1402	bath_dir = tempfile.mkdtemp()
	1403	bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
	1404	#bath_dir = 'bath_data_manager_test'
	1405	#print "os.getcwd( )",os.getcwd( )
	1406	elevation_dir = tempfile.mkdtemp()
	1407	#elevation_dir = 'elev_expanded'
	1408	elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
	1409	elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
	1410
	1411	fid = open(bath_dir_filename, 'w')
	1412	fid.write(""" ncols 4
	1413	nrows 4
	1414	xllcorner 148.00000
	1415	yllcorner -38.00000
	1416	cellsize 0.25
	1417	nodata_value -9999.0
	1418	-9000.000 -1000.000 -3000.0 -2000.000
	1419	-1000.000 9000.000 -1000.000 -3000.000
	1420	-4000.000 6000.000 2000.000 -5000.000
	1421	-9000.000 -1000.000 -3000.0 -2000.000
	1422	""")
	1423	fid.close()
	1424
	1425	fid = open(elevation_dir_filename1, 'w')
	1426	fid.write(""" ncols 4
	1427	nrows 4
	1428	xllcorner 148.00000
	1429	yllcorner -38.00000
	1430	cellsize 0.25
	1431	nodata_value -9999.0
	1432	-900.000 -100.000 -300.0 -200.000
	1433	-100.000 900.000 -100.000 -300.000
	1434	-400.000 600.000 200.000 -500.000
	1435	-900.000 -100.000 -300.0 -200.000
	1436	""")
	1437	fid.close()
	1438
	1439	fid = open(elevation_dir_filename2, 'w')
	1440	fid.write(""" ncols 4
	1441	nrows 4
	1442	xllcorner 148.00000
	1443	yllcorner -38.00000
	1444	cellsize 0.25
	1445	nodata_value -9999.0
	1446	-990.000 -110.000 -330.0 -220.000
	1447	-110.000 990.000 -110.000 -330.000
	1448	-440.000 660.000 220.000 -550.000
	1449	-990.000 -110.000 -330.0 -220.000
	1450	""")
	1451	fid.close()
	1452
	1453	ucur_dir = tempfile.mkdtemp()
	1454	ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
	1455	ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
	1456
	1457	fid = open(ucur_dir_filename1, 'w')
	1458	fid.write(""" ncols 4
	1459	nrows 4
	1460	xllcorner 148.00000
	1461	yllcorner -38.00000
	1462	cellsize 0.25
	1463	nodata_value -9999.0
	1464	-90.000 -10.000 -30.0 -20.000
	1465	-10.000 90.000 -10.000 -30.000
	1466	-40.000 60.000 20.000 -50.000
	1467	-90.000 -10.000 -30.0 -20.000
	1468	""")
	1469	fid.close()
	1470	fid = open(ucur_dir_filename2, 'w')
	1471	fid.write(""" ncols 4
	1472	nrows 4
	1473	xllcorner 148.00000
	1474	yllcorner -38.00000
	1475	cellsize 0.25
	1476	nodata_value -9999.0
	1477	-90.000 -10.000 -30.0 -20.000
	1478	-10.000 99.000 -11.000 -30.000
	1479	-40.000 66.000 22.000 -50.000
	1480	-90.000 -10.000 -30.0 -20.000
	1481	""")
	1482	fid.close()
	1483
	1484	vcur_dir = tempfile.mkdtemp()
	1485	vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
	1486	vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
	1487
	1488	fid = open(vcur_dir_filename1, 'w')
	1489	fid.write(""" ncols 4
	1490	nrows 4
	1491	xllcorner 148.00000
	1492	yllcorner -38.00000
	1493	cellsize 0.25
	1494	nodata_value -9999.0
	1495	-90.000 -10.000 -30.0 -20.000
	1496	-10.000 80.000 -20.000 -30.000
	1497	-40.000 50.000 10.000 -50.000
	1498	-90.000 -10.000 -30.0 -20.000
	1499	""")
	1500	fid.close()
	1501	fid = open(vcur_dir_filename2, 'w')
	1502	fid.write(""" ncols 4
	1503	nrows 4
	1504	xllcorner 148.00000
	1505	yllcorner -38.00000
	1506	cellsize 0.25
	1507	nodata_value -9999.0
	1508	-90.000 -10.000 -30.0 -20.000
	1509	-10.000 88.000 -22.000 -30.000
	1510	-40.000 55.000 11.000 -50.000
	1511	-90.000 -10.000 -30.0 -20.000
	1512	""")
	1513	fid.close()
	1514
	1515	sww_file = tempfile.mktemp(".sww")
	1516	#sww_file = 'a_test.sww'
	1517	asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
	1518	sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
	1519	mean_stage = 100,
	1520	minlat = -37.6, maxlat = -37.6,
	1521	minlon = 148.3, maxlon = 148.3,
	1522	verbose=self.verbose
	1523	#,verbose = True
	1524	)
	1525
	1526	# check the sww file
	1527
	1528	fid = NetCDFFile(sww_file, netcdf_mode_r) # Open existing file for read
	1529	x = fid.variables['x'][:]
	1530	y = fid.variables['y'][:]
	1531	z = fid.variables['elevation'][:]
	1532	stage = fid.variables['stage'][:]
	1533	xmomentum = fid.variables['xmomentum'][:]
	1534	ymomentum = fid.variables['ymomentum'][:]
	1535	geo_ref = Geo_reference(NetCDFObject=fid)
	1536	#print "geo_ref",geo_ref
	1537	x_ref = geo_ref.get_xllcorner()
	1538	y_ref = geo_ref.get_yllcorner()
	1539	self.failUnless(geo_ref.get_zone() == 55, 'Failed')
	1540
	1541	assert num.allclose(fid.starttime, 0.0) # (-37.45, 148.25)
	1542	assert num.allclose(x_ref, 610120.388) # (-37.45, 148.25)
	1543	assert num.allclose(y_ref, 5820863.269 )# (-37.45, 148.5)
	1544
	1545	#Easting: 632145.632 Northing: 5820863.269
	1546	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
	1547
	1548	#print "x",x
	1549	#print "y",y
	1550	self.failUnless(len(x) == 4,'failed') # 2*2
	1551	self.failUnless(len(x) == 4,'failed') # 2*2
	1552
	1553	#Zone: 55
	1554	#Easting: 632145.632 Northing: 5820863.269
	1555	#Latitude: -37 45 ' 0.00000 '' Longitude: 148 30 ' 0.00000 ''
	1556	# magic number - y is close enough for me.
	1557	assert num.allclose(x[3], 632145.63 - x_ref)
	1558	assert num.allclose(y[3], 5820863.269 - y_ref + 5.22155314684e-005)
	1559
	1560	assert num.allclose(z[0],9000.0 ) #z is elevation info
	1561	#print "z",z
	1562	# 2 time steps, 4 points
	1563	self.failUnless(xmomentum.shape == (2,4), 'failed')
	1564	self.failUnless(ymomentum.shape == (2,4), 'failed')
	1565
	1566	#(100.0 - -1000.000)*10
	1567	#assert num.allclose(xmomentum[0][5], 11000.0 )
	1568
	1569	fid.close()
	1570
	1571	# is the sww file readable?
	1572	#Lets see if we can convert it to a dem!
	1573	# if you uncomment, remember to delete the file
	1574	#print "sww_file",sww_file
	1575	#dem_file = tempfile.mktemp(".dem")
[7770]	1576	domain = load_sww_as_domain(sww_file) ###, dem_file)
[7742]	1577	domain.check_integrity()
	1578
	1579	#tidy up
	1580	os.remove(bath_dir_filename)
	1581	os.rmdir(bath_dir)
	1582
	1583	os.remove(elevation_dir_filename1)
	1584	os.remove(elevation_dir_filename2)
	1585	os.rmdir(elevation_dir)
	1586
	1587	os.remove(ucur_dir_filename1)
	1588	os.remove(ucur_dir_filename2)
	1589	os.rmdir(ucur_dir)
	1590
	1591	os.remove(vcur_dir_filename1)
	1592	os.remove(vcur_dir_filename2)
	1593	os.rmdir(vcur_dir)
	1594
	1595	# remove sww file
	1596	os.remove(sww_file)
	1597
	1598
	1599
	1600	def test_get_min_max_indexes(self):
	1601	latitudes = [3,2,1,0]
	1602	longitudes = [0,10,20,30]
	1603
	1604	# k - lat
	1605	# l - lon
[7743]	1606	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1607	latitudes,longitudes,
	1608	-10,4,-10,31)
	1609
	1610	#print "kmin",kmin;print "kmax",kmax
	1611	#print "lmin",lmin;print "lmax",lmax
	1612	latitudes_new = latitudes[kmin:kmax]
	1613	longitudes_news = longitudes[lmin:lmax]
	1614	#print "latitudes_new", latitudes_new
	1615	#print "longitudes_news",longitudes_news
	1616	self.failUnless(latitudes == latitudes_new and \
	1617	longitudes == longitudes_news,
	1618	'failed')
	1619
	1620	## 2nd test
[7743]	1621	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1622	latitudes,longitudes,
	1623	0.5,2.5,5,25)
	1624	#print "kmin",kmin;print "kmax",kmax
	1625	#print "lmin",lmin;print "lmax",lmax
	1626	latitudes_new = latitudes[kmin:kmax]
	1627	longitudes_news = longitudes[lmin:lmax]
	1628	#print "latitudes_new", latitudes_new
	1629	#print "longitudes_news",longitudes_news
	1630
	1631	self.failUnless(latitudes == latitudes_new and \
	1632	longitudes == longitudes_news,
	1633	'failed')
	1634
	1635	## 3rd test
[7743]	1636	kmin, kmax, lmin, lmax = get_min_max_indices(\
[7742]	1637	latitudes,
	1638	longitudes,
	1639	1.1,1.9,12,17)
	1640	#print "kmin",kmin;print "kmax",kmax
	1641	#print "lmin",lmin;print "lmax",lmax
	1642	latitudes_new = latitudes[kmin:kmax]
	1643	longitudes_news = longitudes[lmin:lmax]
	1644	#print "latitudes_new", latitudes_new
	1645	#print "longitudes_news",longitudes_news
	1646
	1647	self.failUnless(latitudes_new == [2, 1] and \
	1648	longitudes_news == [10, 20],
	1649	'failed')
	1650
	1651
	1652	## 4th test
[7743]	1653	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1654	latitudes,longitudes,
	1655	-0.1,1.9,-2,17)
	1656	#print "kmin",kmin;print "kmax",kmax
	1657	#print "lmin",lmin;print "lmax",lmax
	1658	latitudes_new = latitudes[kmin:kmax]
	1659	longitudes_news = longitudes[lmin:lmax]
	1660	#print "latitudes_new", latitudes_new
	1661	#print "longitudes_news",longitudes_news
	1662
	1663	self.failUnless(latitudes_new == [2, 1, 0] and \
	1664	longitudes_news == [0, 10, 20],
	1665	'failed')
	1666	## 5th test
[7743]	1667	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1668	latitudes,longitudes,
	1669	0.1,1.9,2,17)
	1670	#print "kmin",kmin;print "kmax",kmax
	1671	#print "lmin",lmin;print "lmax",lmax
	1672	latitudes_new = latitudes[kmin:kmax]
	1673	longitudes_news = longitudes[lmin:lmax]
	1674	#print "latitudes_new", latitudes_new
	1675	#print "longitudes_news",longitudes_news
	1676
	1677	self.failUnless(latitudes_new == [2, 1, 0] and \
	1678	longitudes_news == [0, 10, 20],
	1679	'failed')
	1680
	1681	## 6th test
	1682
[7743]	1683	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1684	latitudes,longitudes,
	1685	1.5,4,18,32)
	1686	#print "kmin",kmin;print "kmax",kmax
	1687	#print "lmin",lmin;print "lmax",lmax
	1688	latitudes_new = latitudes[kmin:kmax]
	1689	longitudes_news = longitudes[lmin:lmax]
	1690	#print "latitudes_new", latitudes_new
	1691	#print "longitudes_news",longitudes_news
	1692
	1693	self.failUnless(latitudes_new == [3, 2, 1] and \
	1694	longitudes_news == [10, 20, 30],
	1695	'failed')
	1696
	1697
	1698	## 7th test
	1699	m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int) #array default#
[7743]	1700	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1701	latitudes,longitudes,
	1702	1.5,1.5,15,15)
	1703	#print "kmin",kmin;print "kmax",kmax
	1704	#print "lmin",lmin;print "lmax",lmax
	1705	latitudes_new = latitudes[kmin:kmax]
	1706	longitudes_news = longitudes[lmin:lmax]
	1707	m2d = m2d[kmin:kmax,lmin:lmax]
	1708	#print "m2d", m2d
	1709	#print "latitudes_new", latitudes_new
	1710	#print "longitudes_news",longitudes_news
	1711
	1712	self.failUnless(num.alltrue(latitudes_new == [2, 1]) and
	1713	num.alltrue(longitudes_news == [10, 20]),
	1714	'failed')
	1715
	1716	self.failUnless(num.alltrue(m2d == [[5,6],[9,10]]), 'failed')
	1717
	1718	def test_get_min_max_indexes_lat_ascending(self):
	1719	latitudes = [0,1,2,3]
	1720	longitudes = [0,10,20,30]
	1721
	1722	# k - lat
	1723	# l - lon
[7743]	1724	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1725	latitudes,longitudes,
	1726	-10,4,-10,31)
	1727
	1728	#print "kmin",kmin;print "kmax",kmax
	1729	#print "lmin",lmin;print "lmax",lmax
	1730	latitudes_new = latitudes[kmin:kmax]
	1731	longitudes_news = longitudes[lmin:lmax]
	1732	#print "latitudes_new", latitudes_new
	1733	#print "longitudes_news",longitudes_news
	1734	self.failUnless(latitudes == latitudes_new and \
	1735	longitudes == longitudes_news,
	1736	'failed')
	1737
	1738	## 3rd test
[7743]	1739	kmin, kmax, lmin, lmax = get_min_max_indices(\
[7742]	1740	latitudes,
	1741	longitudes,
	1742	1.1,1.9,12,17)
	1743	#print "kmin",kmin;print "kmax",kmax
	1744	#print "lmin",lmin;print "lmax",lmax
	1745	latitudes_new = latitudes[kmin:kmax]
	1746	longitudes_news = longitudes[lmin:lmax]
	1747	#print "latitudes_new", latitudes_new
	1748	#print "longitudes_news",longitudes_news
	1749
	1750	self.failUnless(latitudes_new == [1, 2] and \
	1751	longitudes_news == [10, 20],
	1752	'failed')
	1753
	1754	def test_get_min_max_indexes2(self):
	1755	latitudes = [-30,-35,-40,-45]
	1756	longitudes = [148,149,150,151]
	1757
	1758	m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int) #array default#
	1759
	1760	# k - lat
	1761	# l - lon
[7743]	1762	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1763	latitudes,longitudes,
	1764	-37,-27,147,149.5)
	1765
	1766	#print "kmin",kmin;print "kmax",kmax
	1767	#print "lmin",lmin;print "lmax",lmax
	1768	#print "m2d", m2d
	1769	#print "latitudes", latitudes
	1770	#print "longitudes",longitudes
	1771	#print "latitudes[kmax]", latitudes[kmax]
	1772	latitudes_new = latitudes[kmin:kmax]
	1773	longitudes_new = longitudes[lmin:lmax]
	1774	m2d = m2d[kmin:kmax,lmin:lmax]
	1775	#print "m2d", m2d
	1776	#print "latitudes_new", latitudes_new
	1777	#print "longitudes_new",longitudes_new
	1778
	1779	self.failUnless(latitudes_new == [-30, -35, -40] and
	1780	longitudes_new == [148, 149,150],
	1781	'failed')
	1782	self.failUnless(num.alltrue(m2d == [[0,1,2],[4,5,6],[8,9,10]]),
	1783	'failed')
	1784
	1785	def test_get_min_max_indexes3(self):
	1786	latitudes = [-30,-35,-40,-45,-50,-55,-60]
	1787	longitudes = [148,149,150,151]
	1788
	1789	# k - lat
	1790	# l - lon
[7743]	1791	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1792	latitudes,longitudes,
	1793	-43,-37,148.5,149.5)
	1794
	1795
	1796	#print "kmin",kmin;print "kmax",kmax
	1797	#print "lmin",lmin;print "lmax",lmax
	1798	#print "latitudes", latitudes
	1799	#print "longitudes",longitudes
	1800	latitudes_new = latitudes[kmin:kmax]
	1801	longitudes_news = longitudes[lmin:lmax]
	1802	#print "latitudes_new", latitudes_new
	1803	#print "longitudes_news",longitudes_news
	1804
	1805	self.failUnless(latitudes_new == [-35, -40, -45] and
	1806	longitudes_news == [148, 149,150],
	1807	'failed')
	1808
	1809	def test_get_min_max_indexes4(self):
	1810	latitudes = [-30,-35,-40,-45,-50,-55,-60]
	1811	longitudes = [148,149,150,151]
	1812
	1813	# k - lat
	1814	# l - lon
[7743]	1815	kmin, kmax, lmin, lmax = get_min_max_indices(
[7742]	1816	latitudes,longitudes)
	1817
	1818	#print "kmin",kmin;print "kmax",kmax
	1819	#print "lmin",lmin;print "lmax",lmax
	1820	#print "latitudes", latitudes
	1821	#print "longitudes",longitudes
	1822	latitudes_new = latitudes[kmin:kmax]
	1823	longitudes_news = longitudes[lmin:lmax]
	1824	#print "latitudes_new", latitudes_new
	1825	#print "longitudes_news",longitudes_news
	1826
	1827	self.failUnless(latitudes_new == latitudes and \
	1828	longitudes_news == longitudes,
	1829	'failed')
	1830
	1831	def test_tsh2sww(self):
	1832	import os
	1833	import tempfile
	1834
	1835	tsh_file = tempfile.mktemp(".tsh")
	1836	file = open(tsh_file,"w")
	1837	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
	1838	0 0.0 0.0 0.0 0.0 0.01 \n \
	1839	1 1.0 0.0 10.0 10.0 0.02 \n \
	1840	2 0.0 1.0 0.0 10.0 0.03 \n \
	1841	3 0.5 0.25 8.0 12.0 0.04 \n \
	1842	# Vert att title \n \
	1843	elevation \n \
	1844	stage \n \
	1845	friction \n \
	1846	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
	1847	0 0 3 2 -1 -1 1 dsg\n\
	1848	1 0 1 3 -1 0 -1 ole nielsen\n\
	1849	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
	1850	0 1 0 2 \n\
	1851	1 0 2 3 \n\
	1852	2 2 3 \n\
	1853	3 3 1 1 \n\
	1854	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
	1855	0 216.0 -86.0 \n \
	1856	1 160.0 -167.0 \n \
	1857	2 114.0 -91.0 \n \
	1858	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
	1859	0 0 1 0 \n \
	1860	1 1 2 0 \n \
	1861	2 2 0 0 \n \
	1862	0 # <x> <y> ...Mesh Holes... \n \
	1863	0 # <x> <y> <attribute>...Mesh Regions... \n \
	1864	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
	1865	#Geo reference \n \
	1866	56 \n \
	1867	140 \n \
	1868	120 \n")
	1869	file.close()
	1870
	1871	#sww_file = tempfile.mktemp(".sww")
	1872	#print "sww_file",sww_file
	1873	#print "sww_file",tsh_file
	1874	tsh2sww(tsh_file,
	1875	verbose=self.verbose)
	1876
	1877	os.remove(tsh_file)
	1878	os.remove(tsh_file[:-4] + '.sww')
	1879
	1880
	1881	def test_urs2sww_test_fail(self):
	1882	points_num = -100
	1883	time_step_count = 45
	1884	time_step = -7
	1885	file_handle, base_name = tempfile.mkstemp("")
	1886	os.close(file_handle)
	1887	os.remove(base_name)
	1888
	1889	files = []
	1890	quantities = ['HA','UA','VA']
	1891
	1892	mux_names = [WAVEHEIGHT_MUX_LABEL,
	1893	EAST_VELOCITY_LABEL,
	1894	NORTH_VELOCITY_LABEL]
	1895	for i,q in enumerate(quantities):
	1896	#Write C files
	1897	columns = 3 # long, lat , depth
	1898	file = base_name + mux_names[i]
	1899	f = open(file, 'wb')
	1900	files.append(file)
	1901	f.write(pack('i',points_num))
	1902	f.write(pack('i',time_step_count))
	1903	f.write(pack('f',time_step))
	1904
	1905	f.close()
	1906	tide = 1
	1907	try:
	1908	urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
	1909	verbose=self.verbose)
	1910	except ANUGAError:
	1911	pass
	1912	else:
	1913	self.delete_mux(files)
	1914	msg = 'Should have raised exception'
	1915	raise msg
	1916	sww_file = base_name + '.sww'
	1917	self.delete_mux(files)
	1918
	1919	def test_urs2sww_test_fail2(self):
	1920	base_name = 'Harry-high-pants'
	1921	try:
	1922	urs2sww(base_name)
	1923	except IOError:
	1924	pass
	1925	else:
	1926	self.delete_mux(files)
	1927	msg = 'Should have raised exception'
	1928	raise msg
	1929
	1930	def test_urs2sww(self):
	1931	tide = 1
	1932	base_name, files = self.create_mux()
	1933	urs2sww(base_name
	1934	#, origin=(0,0,0)
	1935	, mean_stage=tide
	1936	, remove_nc_files=True,
	1937	verbose=self.verbose
	1938	)
	1939	sww_file = base_name + '.sww'
	1940
	1941	#Let's interigate the sww file
	1942	# Note, the sww info is not gridded. It is point data.
	1943	fid = NetCDFFile(sww_file)
	1944
	1945	x = fid.variables['x'][:]
	1946	y = fid.variables['y'][:]
	1947	geo_reference = Geo_reference(NetCDFObject=fid)
	1948
	1949
	1950	#Check that first coordinate is correctly represented
	1951	#Work out the UTM coordinates for first point
	1952	zone, e, n = redfearn(-34.5, 150.66667)
	1953
	1954	assert num.allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
	1955
	1956	# Make x and y absolute
	1957	points = geo_reference.get_absolute(map(None, x, y))
	1958	points = ensure_numeric(points)
	1959	x = points[:,0]
	1960	y = points[:,1]
	1961
	1962	#Check first value
	1963	stage = fid.variables['stage'][:]
	1964	xmomentum = fid.variables['xmomentum'][:]
	1965	ymomentum = fid.variables['ymomentum'][:]
	1966	elevation = fid.variables['elevation'][:]
	1967	assert num.allclose(stage[0,0], e +tide) #Meters
	1968
	1969	#Check the momentums - ua
	1970	#momentum = velocity*(stage-elevation)
	1971	# elevation = - depth
	1972	#momentum = velocity_ua *(stage+depth)
	1973	# = n*(e+tide+n) based on how I'm writing these files
	1974	#
	1975	answer_x = n*(e+tide+n)
	1976	actual_x = xmomentum[0,0]
	1977	#print "answer_x",answer_x
	1978	#print "actual_x",actual_x
	1979	assert num.allclose(answer_x, actual_x) #Meters
	1980
	1981	#Check the momentums - va
	1982	#momentum = velocity*(stage-elevation)
	1983	# -(-elevation) since elevation is inverted in mux files
	1984	#momentum = velocity_va *(stage+elevation)
	1985	# = e*(e+tide+n) based on how I'm writing these files
	1986	answer_y = e(e+tide+n) -1 # talking into account mux file format
	1987	actual_y = ymomentum[0,0]
	1988	#print "answer_y",answer_y
	1989	#print "actual_y",actual_y
	1990	assert num.allclose(answer_y, actual_y) #Meters
	1991
	1992	assert num.allclose(answer_x, actual_x) #Meters
	1993
	1994	# check the stage values, first time step.
	1995	# These arrays are equal since the Easting values were used as
	1996	# the stage
	1997	assert num.allclose(stage[0], x +tide) #Meters
	1998
	1999	# check the elevation values.
	2000	# -ve since urs measures depth, sww meshers height,
	2001	# these arrays are equal since the northing values were used as
	2002	# the elevation
	2003	assert num.allclose(-elevation, y) #Meters
	2004
	2005	fid.close()
	2006	self.delete_mux(files)
	2007	os.remove(sww_file)
	2008
	2009	def test_urs2sww_momentum(self):
	2010	tide = 1
	2011	time_step_count = 3
	2012	time_step = 2
	2013	#lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
	2014	# This is gridded
	2015	lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
	2016	depth=20
	2017	ha=2
	2018	ua=5
	2019	va=-10 #-ve added to take into account mux file format where south
	2020	# is positive.
	2021	base_name, files = self.write_mux(lat_long_points,
	2022	time_step_count, time_step,
	2023	depth=depth,
	2024	ha=ha,
	2025	ua=ua,
	2026	va=va)
	2027	# write_mux(self,lat_long_points, time_step_count, time_step,
	2028	# depth=None, ha=None, ua=None, va=None
	2029	urs2sww(base_name
	2030	#, origin=(0,0,0)
	2031	, mean_stage=tide
	2032	, remove_nc_files=True,
	2033	verbose=self.verbose
	2034	)
	2035	sww_file = base_name + '.sww'
	2036
	2037	#Let's interigate the sww file
	2038	# Note, the sww info is not gridded. It is point data.
	2039	fid = NetCDFFile(sww_file)
	2040
	2041	x = fid.variables['x'][:]
	2042	y = fid.variables['y'][:]
	2043	geo_reference = Geo_reference(NetCDFObject=fid)
	2044
	2045	#Check first value
	2046	stage = fid.variables['stage'][:]
	2047	xmomentum = fid.variables['xmomentum'][:]
	2048	ymomentum = fid.variables['ymomentum'][:]
	2049	elevation = fid.variables['elevation'][:]
	2050	#assert allclose(stage[0,0], e + tide) #Meters
	2051	#print "xmomentum", xmomentum
	2052	#print "ymomentum", ymomentum
	2053	#Check the momentums - ua
	2054	#momentum = velocity*water height
	2055	#water height = mux_depth + mux_height +tide
	2056	#water height = mux_depth + mux_height +tide
	2057	#momentum = velocity*(mux_depth + mux_height +tide)
	2058	#
	2059
	2060	answer = 115
	2061	actual = xmomentum[0,0]
	2062	assert num.allclose(answer, actual) #Meters^2/ sec
	2063	answer = 230
	2064	actual = ymomentum[0,0]
	2065	#print "answer",answer
	2066	#print "actual",actual
	2067	assert num.allclose(answer, actual) #Meters^2/ sec
	2068
	2069	# check the stage values, first time step.
	2070	# These arrays are equal since the Easting values were used as
	2071	# the stage
	2072
	2073	#assert allclose(stage[0], x +tide) #Meters
	2074
	2075	# check the elevation values.
	2076	# -ve since urs measures depth, sww meshers height,
	2077	# these arrays are equal since the northing values were used as
	2078	# the elevation
	2079	#assert allclose(-elevation, y) #Meters
	2080
	2081	fid.close()
	2082	self.delete_mux(files)
	2083	os.remove(sww_file)
	2084
	2085
	2086	def test_urs2sww_origin(self):
	2087	tide = 1
	2088	base_name, files = self.create_mux()
	2089	urs2sww(base_name
	2090	, origin=(0,0,0)
	2091	, mean_stage=tide
	2092	, remove_nc_files=True,
	2093	verbose=self.verbose
	2094	)
	2095	sww_file = base_name + '.sww'
	2096
	2097	#Let's interigate the sww file
	2098	# Note, the sww info is not gridded. It is point data.
	2099	fid = NetCDFFile(sww_file)
	2100
	2101	# x and y are absolute
	2102	x = fid.variables['x'][:]
	2103	y = fid.variables['y'][:]
	2104	geo_reference = Geo_reference(NetCDFObject=fid)
	2105
	2106
	2107	time = fid.variables['time'][:]
	2108	#print "time", time
	2109	assert num.allclose([0.,0.5,1.], time)
	2110	assert fid.starttime == 0.0
	2111	#Check that first coordinate is correctly represented
	2112	#Work out the UTM coordinates for first point
	2113	zone, e, n = redfearn(-34.5, 150.66667)
	2114
	2115	assert num.allclose([x[0],y[0]], [e,n])
	2116
	2117
	2118	#Check first value
	2119	stage = fid.variables['stage'][:]
	2120	xmomentum = fid.variables['xmomentum'][:]
	2121	ymomentum = fid.variables['ymomentum'][:]
	2122	elevation = fid.variables['elevation'][:]
	2123	assert num.allclose(stage[0,0], e +tide) #Meters
	2124
	2125	#Check the momentums - ua
	2126	#momentum = velocity*(stage-elevation)
	2127	#momentum = velocity*(stage+elevation)
	2128	# -(-elevation) since elevation is inverted in mux files
	2129	# = n*(e+tide+n) based on how I'm writing these files
	2130	answer = n*(e+tide+n)
	2131	actual = xmomentum[0,0]
	2132	assert num.allclose(answer, actual) #Meters
	2133
	2134	# check the stage values, first time step.
	2135	# These arrays are equal since the Easting values were used as
	2136	# the stage
	2137	assert num.allclose(stage[0], x +tide) #Meters
	2138
	2139	# check the elevation values.
	2140	# -ve since urs measures depth, sww meshers height,
	2141	# these arrays are equal since the northing values were used as
	2142	# the elevation
	2143	assert num.allclose(-elevation, y) #Meters
	2144
	2145	fid.close()
	2146	self.delete_mux(files)
	2147	os.remove(sww_file)
	2148
	2149	def test_urs2sww_minmaxlatlong(self):
	2150
	2151	#longitudes = [150.66667, 150.83334, 151., 151.16667]
	2152	#latitudes = [-34.5, -34.33333, -34.16667, -34]
	2153
	2154	tide = 1
	2155	base_name, files = self.create_mux()
	2156	urs2sww(base_name,
	2157	minlat=-34.5,
	2158	maxlat=-34,
	2159	minlon= 150.66667,
	2160	maxlon= 151.16667,
	2161	mean_stage=tide,
	2162	remove_nc_files=True,
	2163	verbose=self.verbose
	2164	)
	2165	sww_file = base_name + '.sww'
	2166
	2167	#Let's interigate the sww file
	2168	# Note, the sww info is not gridded. It is point data.
	2169	fid = NetCDFFile(sww_file)
	2170
	2171
	2172	# Make x and y absolute
	2173	x = fid.variables['x'][:]
	2174	y = fid.variables['y'][:]
	2175	geo_reference = Geo_reference(NetCDFObject=fid)
	2176	points = geo_reference.get_absolute(map(None, x, y))
	2177	points = ensure_numeric(points)
	2178	x = points[:,0]
	2179	y = points[:,1]
	2180
	2181	#Check that first coordinate is correctly represented
	2182	#Work out the UTM coordinates for first point
	2183	zone, e, n = redfearn(-34.5, 150.66667)
	2184	assert num.allclose([x[0],y[0]], [e,n])
	2185
	2186
	2187	#Check first value
	2188	stage = fid.variables['stage'][:]
	2189	xmomentum = fid.variables['xmomentum'][:]
	2190	ymomentum = fid.variables['ymomentum'][:]
	2191	elevation = fid.variables['elevation'][:]
	2192	assert num.allclose(stage[0,0], e +tide) #Meters
	2193
	2194	#Check the momentums - ua
	2195	#momentum = velocity*(stage-elevation)
	2196	#momentum = velocity*(stage+elevation)
	2197	# -(-elevation) since elevation is inverted in mux files
	2198	# = n*(e+tide+n) based on how I'm writing these files
	2199	answer = n*(e+tide+n)
	2200	actual = xmomentum[0,0]
	2201	assert num.allclose(answer, actual) #Meters
	2202
	2203	# check the stage values, first time step.
	2204	# These arrays are equal since the Easting values were used as
	2205	# the stage
	2206	assert num.allclose(stage[0], x +tide) #Meters
	2207
	2208	# check the elevation values.
	2209	# -ve since urs measures depth, sww meshers height,
	2210	# these arrays are equal since the northing values were used as
	2211	# the elevation
	2212	assert num.allclose(-elevation, y) #Meters
	2213
	2214	fid.close()
	2215	self.delete_mux(files)
	2216	os.remove(sww_file)
	2217
	2218	def test_urs2sww_minmaxmintmaxt(self):
	2219
	2220	#longitudes = [150.66667, 150.83334, 151., 151.16667]
	2221	#latitudes = [-34.5, -34.33333, -34.16667, -34]
	2222
	2223	tide = 1
	2224	base_name, files = self.create_mux()
	2225
	2226	urs2sww(base_name,
	2227	mint=0.25,
	2228	maxt=0.75,
	2229	mean_stage=tide,
	2230	remove_nc_files=True,
	2231	verbose=self.verbose)
	2232	sww_file = base_name + '.sww'
	2233
	2234	#Let's interigate the sww file
	2235	# Note, the sww info is not gridded. It is point data.
	2236	fid = NetCDFFile(sww_file)
	2237
	2238
	2239	time = fid.variables['time'][:]
	2240	assert num.allclose(time, [0.0]) # the time is relative
	2241	assert fid.starttime == 0.5
	2242
	2243	fid.close()
	2244	self.delete_mux(files)
	2245	#print "sww_file", sww_file
	2246	os.remove(sww_file)
	2247
	2248
	2249
[7744]	2250	def test_read_asc(self):
	2251	"""Test conversion from dem in ascii format to native NetCDF format
	2252	"""
	2253
	2254	import time, os
	2255
	2256	from file_conversion import _read_asc
	2257	#Write test asc file
	2258	filename = tempfile.mktemp(".000")
	2259	fid = open(filename, 'w')
	2260	fid.write("""ncols 7
	2261	nrows 4
	2262	xllcorner 2000.5
	2263	yllcorner 3000.5
	2264	cellsize 25
	2265	NODATA_value -9999
	2266	97.921 99.285 125.588 180.830 258.645 342.872 415.836
	2267	473.157 514.391 553.893 607.120 678.125 777.283 883.038
	2268	984.494 1040.349 1008.161 900.738 730.882 581.430 514.980
	2269	502.645 516.230 504.739 450.604 388.500 338.097 514.980
	2270	""")
	2271	fid.close()
	2272	bath_metadata, grid = _read_asc(filename, verbose=self.verbose)
	2273	self.failUnless(bath_metadata['xllcorner'] == 2000.5, 'Failed')
	2274	self.failUnless(bath_metadata['yllcorner'] == 3000.5, 'Failed')
	2275	self.failUnless(bath_metadata['cellsize'] == 25, 'Failed')
	2276	self.failUnless(bath_metadata['NODATA_value'] == -9999, 'Failed')
	2277	self.failUnless(grid[0][0] == 97.921, 'Failed')
	2278	self.failUnless(grid[3][6] == 514.980, 'Failed')
	2279
	2280	os.remove(filename)
	2281
	2282
	2283	#### TESTS FOR URS 2 SWW ###
	2284
	2285	def create_mux(self, points_num=None):
	2286	# write all the mux stuff.
	2287	time_step_count = 3
	2288	time_step = 0.5
	2289
	2290	longitudes = [150.66667, 150.83334, 151., 151.16667]
	2291	latitudes = [-34.5, -34.33333, -34.16667, -34]
	2292
	2293	if points_num == None:
	2294	points_num = len(longitudes) * len(latitudes)
	2295
	2296	lonlatdeps = []
	2297	quantities = ['HA','UA','VA']
	2298	mux_names = [WAVEHEIGHT_MUX_LABEL,
	2299	EAST_VELOCITY_LABEL,
	2300	NORTH_VELOCITY_LABEL]
	2301	quantities_init = [[],[],[]]
	2302	# urs binary is latitude fastest
	2303	for i,lon in enumerate(longitudes):
	2304	for j,lat in enumerate(latitudes):
	2305	_ , e, n = redfearn(lat, lon)
	2306	lonlatdeps.append([lon, lat, n])
	2307	quantities_init[0].append(e) # HA
	2308	quantities_init[1].append(n ) # UA
	2309	quantities_init[2].append(e) # VA
	2310	#print "lonlatdeps",lonlatdeps
	2311
	2312	file_handle, base_name = tempfile.mkstemp("")
	2313	os.close(file_handle)
	2314	os.remove(base_name)
	2315
	2316	files = []
	2317	for i,q in enumerate(quantities):
	2318	quantities_init[i] = ensure_numeric(quantities_init[i])
	2319	#print "HA_init", HA_init
	2320	q_time = num.zeros((time_step_count, points_num), num.float)
	2321	for time in range(time_step_count):
	2322	q_time[time,:] = quantities_init[i] #* time * 4
	2323
	2324	#Write C files
	2325	columns = 3 # long, lat , depth
	2326	file = base_name + mux_names[i]
	2327	#print "base_name file",file
	2328	f = open(file, 'wb')
	2329	files.append(file)
	2330	f.write(pack('i',points_num))
	2331	f.write(pack('i',time_step_count))
	2332	f.write(pack('f',time_step))
	2333
	2334	#write lat/long info
	2335	for lonlatdep in lonlatdeps:
	2336	for float in lonlatdep:
	2337	f.write(pack('f',float))
	2338
	2339	# Write quantity info
	2340	for time in range(time_step_count):
	2341	for point_i in range(points_num):
	2342	f.write(pack('f',q_time[time,point_i]))
	2343	#print " mux_names[i]", mux_names[i]
	2344	#print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
	2345	f.close()
	2346	return base_name, files
	2347
[7742]	2348	#-------------------------------------------------------------
	2349
	2350	if __name__ == "__main__":
	2351	suite = unittest.makeSuite(Test_File_Conversion,'test')
	2352	runner = unittest.TextTestRunner()
	2353	runner.run(suite)
	2354

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

Download in other formats: