Context Navigation

source: anuga_core/source/anuga/shallow_water/data_manager.py @ 4731

Last change on this file since 4731 was 4706, checked in by ole, 18 years ago
Replaced ':' character with '.' as Windows version of netcdf doesn't allow it.
File size: 188.3 KB

Line
1	"""datamanager.py - input output for AnuGA
2
3
4	This module takes care of reading and writing datafiles such as topograhies,
5	model output, etc
6
7
8	Formats used within AnuGA:
9
10	.sww: Netcdf format for storing model output f(t,x,y)
11	.tms: Netcdf format for storing time series f(t)
12
13	.csv: ASCII format for storing arbitrary points and associated attributes
14	.pts: NetCDF format for storing arbitrary points and associated attributes
15
16	.asc: ASCII format of regular DEMs as output from ArcView
17	.prj: Associated ArcView file giving more meta data for asc format
18	.ers: ERMapper header format of regular DEMs for ArcView
19
20	.dem: NetCDF representation of regular DEM data
21
22	.tsh: ASCII format for storing meshes and associated boundary and region info
23	.msh: NetCDF format for storing meshes and associated boundary and region info
24
25	.nc: Native ferret NetCDF format
26	.geo: Houdinis ascii geometry format (?)
27
28
29	A typical dataflow can be described as follows
30
31	Manually created files:
32	ASC, PRJ: Digital elevation models (gridded)
33	TSH: Triangular meshes (e.g. created from anuga.pmesh)
34	NC Model outputs for use as boundary conditions (e.g from MOST)
35
36
37	AUTOMATICALLY CREATED FILES:
38
39	ASC, PRJ -> DEM -> PTS: Conversion of DEM's to native pts file
40
41	NC -> SWW: Conversion of MOST bundary files to boundary sww
42
43	PTS + TSH -> TSH with elevation: Least squares fit
44
45	TSH -> SWW: Conversion of TSH to sww viewable using Swollen
46
47	TSH + Boundary SWW -> SWW: Simluation using abstract_2d_finite_volumes
48
49	"""
50
51	import exceptions
52	class TitleValueError(exceptions.Exception): pass
53	class DataMissingValuesError(exceptions.Exception): pass
54	class DataFileNotOpenError(exceptions.Exception): pass
55	class DataTimeError(exceptions.Exception): pass
56	class DataDomainError(exceptions.Exception): pass
57	class NewQuantity(exceptions.Exception): pass
58
59
60
61	import csv
62	import os, sys
63	import shutil
64	from struct import unpack
65	import array as p_array
66	#import time, os
67	from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
68
69
70	from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
71	searchsorted, zeros, allclose, around, reshape, transpose, sort, \
72	NewAxis, ArrayType, compress, take, arange, argmax, alltrue,shape,Float
73
74	import string
75
76	from Scientific.IO.NetCDF import NetCDFFile
77	#from shutil import copy
78	from os.path import exists, basename, join
79	from os import getcwd
80
81
82	from anuga.coordinate_transforms.redfearn import redfearn, \
83	convert_from_latlon_to_utm
84	from anuga.coordinate_transforms.geo_reference import Geo_reference, \
85	write_NetCDF_georeference, ensure_geo_reference
86	from anuga.geospatial_data.geospatial_data import Geospatial_data,\
87	ensure_absolute
88	from anuga.config import minimum_storable_height as default_minimum_storable_height
89	from anuga.config import max_float
90	from anuga.utilities.numerical_tools import ensure_numeric, mean
91	from anuga.caching.caching import myhash
92	from anuga.utilities.anuga_exceptions import ANUGAError
93	from anuga.shallow_water import Domain
94	from anuga.abstract_2d_finite_volumes.pmesh2domain import \
95	pmesh_to_domain_instance
96	from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
97	remove_lone_verts, sww2timeseries, get_centroid_values
98	from anuga.load_mesh.loadASCII import export_mesh_file
99	# formula mappings
100
101	quantity_formula = {'momentum':'(xmomentum2 + ymomentum2)**0.5',
102	'depth':'stage-elevation',
103	'speed': \
104	'(xmomentum2 + ymomentum2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
105
106
107
108	def make_filename(s):
109	"""Transform argument string into a Sexsuitable filename
110	"""
111
112	s = s.strip()
113	s = s.replace(' ', '_')
114	s = s.replace('(', '')
115	s = s.replace(')', '')
116	s = s.replace('__', '_')
117
118	return s
119
120
121	def check_dir(path, verbose=None):
122	"""Check that specified path exists.
123	If path does not exist it will be created if possible
124
125	USAGE:
126	checkdir(path, verbose):
127
128	ARGUMENTS:
129	path -- Directory
130	verbose -- Flag verbose output (default: None)
131
132	RETURN VALUE:
133	Verified path including trailing separator
134
135	"""
136
137	import os.path
138
139	if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
140	unix = 0
141	else:
142	unix = 1
143
144
145	if path[-1] != os.sep:
146	path = path + os.sep # Add separator for directories
147
148	path = os.path.expanduser(path) # Expand ~ or ~user in pathname
149	if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
150	try:
151	exitcode=os.mkdir(path)
152
153	# Change access rights if possible
154	#
155	if unix:
156	exitcode=os.system('chmod 775 '+path)
157	else:
158	pass # FIXME: What about acces rights under Windows?
159
160	if verbose: print 'MESSAGE: Directory', path, 'created.'
161
162	except:
163	print 'WARNING: Directory', path, 'could not be created.'
164	if unix:
165	path = '/tmp/'
166	else:
167	path = 'C:'
168
169	print 'Using directory %s instead' %path
170
171	return(path)
172
173
174
175	def del_dir(path):
176	"""Recursively delete directory path and all its contents
177	"""
178
179	import os
180
181	if os.path.isdir(path):
182	for file in os.listdir(path):
183	X = os.path.join(path, file)
184
185
186	if os.path.isdir(X) and not os.path.islink(X):
187	del_dir(X)
188	else:
189	try:
190	os.remove(X)
191	except:
192	print "Could not remove file %s" %X
193
194	os.rmdir(path)
195
196
197	# ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007
198	def rmgeneric(path, __func__,verbose=False):
199	ERROR_STR= """Error removing %(path)s, %(error)s """
200
201	try:
202	__func__(path)
203	if verbose: print 'Removed ', path
204	except OSError, (errno, strerror):
205	print ERROR_STR % {'path' : path, 'error': strerror }
206
207	def removeall(path,verbose=False):
208
209	if not os.path.isdir(path):
210	return
211
212	files=os.listdir(path)
213
214	for x in files:
215	fullpath=os.path.join(path, x)
216	if os.path.isfile(fullpath):
217	f=os.remove
218	rmgeneric(fullpath, f)
219	elif os.path.isdir(fullpath):
220	removeall(fullpath)
221	f=os.rmdir
222	rmgeneric(fullpath, f,verbose)
223
224
225
226	def create_filename(datadir, filename, format, size=None, time=None):
227
228	import os
229	#from anuga.config import data_dir
230
231	FN = check_dir(datadir) + filename
232
233	if size is not None:
234	FN += '_size%d' %size
235
236	if time is not None:
237	FN += '_time%.2f' %time
238
239	FN += '.' + format
240	return FN
241
242
243	def get_files(datadir, filename, format, size):
244	"""Get all file (names) with given name, size and format
245	"""
246
247	import glob
248
249	import os
250	#from anuga.config import data_dir
251
252	dir = check_dir(datadir)
253
254	pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
255	return glob.glob(pattern)
256
257
258
259	#Generic class for storing output to e.g. visualisation or checkpointing
260	class Data_format:
261	"""Generic interface to data formats
262	"""
263
264
265	def __init__(self, domain, extension, mode = 'w'):
266	assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
267	''' 'w' (write)'''+\
268	''' 'r' (read)''' +\
269	''' 'a' (append)'''
270
271	#Create filename
272	self.filename = create_filename(domain.get_datadir(),
273	domain.get_name(), extension)
274
275	#print 'F', self.filename
276	self.timestep = 0
277	self.domain = domain
278
279
280
281	# Exclude ghosts in case this is a parallel domain
282	self.number_of_nodes = domain.number_of_full_nodes
283	self.number_of_volumes = domain.number_of_full_triangles
284	#self.number_of_volumes = len(domain)
285
286
287
288
289	#FIXME: Should we have a general set_precision function?
290
291
292
293	#Class for storing output to e.g. visualisation
294	class Data_format_sww(Data_format):
295	"""Interface to native NetCDF format (.sww) for storing model output
296
297	There are two kinds of data
298
299	1: Constant data: Vertex coordinates and field values. Stored once
300	2: Variable data: Conserved quantities. Stored once per timestep.
301
302	All data is assumed to reside at vertex locations.
303	"""
304
305
306	def __init__(self, domain, mode = 'w',\
307	max_size = 2000000000,
308	recursion = False):
309	from Scientific.IO.NetCDF import NetCDFFile
310	from Numeric import Int, Float, Float32
311
312	self.precision = Float32 #Use single precision for quantities
313	if hasattr(domain, 'max_size'):
314	self.max_size = domain.max_size #file size max is 2Gig
315	else:
316	self.max_size = max_size
317	self.recursion = recursion
318	self.mode = mode
319
320	Data_format.__init__(self, domain, 'sww', mode)
321
322	if hasattr(domain, 'minimum_storable_height'):
323	self.minimum_storable_height = domain.minimum_storable_height
324	else:
325	self.minimum_storable_height = default_minimum_storable_height
326
327	# NetCDF file definition
328	fid = NetCDFFile(self.filename, mode)
329
330	if mode == 'w':
331	description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
332	self.writer = Write_sww()
333	self.writer.store_header(fid,
334	domain.starttime,
335	self.number_of_volumes,
336	self.domain.number_of_full_nodes,
337	description=description,
338	smoothing=domain.smooth,
339	order=domain.default_order)
340
341	# Extra optional information
342	if hasattr(domain, 'texture'):
343	fid.texture = domain.texture
344
345	if domain.quantities_to_be_monitored is not None:
346	fid.createDimension('singleton', 1)
347	fid.createDimension('two', 2)
348
349	poly = domain.monitor_polygon
350	if poly is not None:
351	N = len(poly)
352	fid.createDimension('polygon_length', N)
353	fid.createVariable('extrema.polygon',
354	self.precision,
355	('polygon_length',
356	'two'))
357	fid.variables['extrema.polygon'][:] = poly
358
359
360	interval = domain.monitor_time_interval
361	if interval is not None:
362	fid.createVariable('extrema.time_interval',
363	self.precision,
364	('two',))
365	fid.variables['extrema.time_interval'][:] = interval
366
367
368	for q in domain.quantities_to_be_monitored:
369	#print 'doing', q
370	fid.createVariable(q+'.extrema', self.precision,
371	('numbers_in_range',))
372	fid.createVariable(q+'.min_location', self.precision,
373	('numbers_in_range',))
374	fid.createVariable(q+'.max_location', self.precision,
375	('numbers_in_range',))
376	fid.createVariable(q+'.min_time', self.precision,
377	('singleton',))
378	fid.createVariable(q+'.max_time', self.precision,
379	('singleton',))
380
381
382	fid.close()
383
384
385	def store_connectivity(self):
386	"""Specialisation of store_connectivity for net CDF format
387
388	Writes x,y,z coordinates of triangles constituting
389	the bed elevation.
390	"""
391
392	from Scientific.IO.NetCDF import NetCDFFile
393
394	from Numeric import concatenate, Int
395
396	domain = self.domain
397
398	#Get NetCDF
399	fid = NetCDFFile(self.filename, 'a') #Open existing file for append
400
401	# Get the variables
402	x = fid.variables['x']
403	y = fid.variables['y']
404	z = fid.variables['elevation']
405
406	volumes = fid.variables['volumes']
407
408	# Get X, Y and bed elevation Z
409	Q = domain.quantities['elevation']
410	X,Y,Z,V = Q.get_vertex_values(xy=True,
411	precision=self.precision)
412
413	#
414	points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
415	self.writer.store_triangulation(fid,
416	points,
417	V.astype(volumes.typecode()),
418	Z,
419	points_georeference= \
420	domain.geo_reference)
421
422	# Close
423	fid.close()
424
425
426	def store_timestep(self, names):
427	"""Store time and named quantities to file
428	"""
429	from Scientific.IO.NetCDF import NetCDFFile
430	import types
431	from time import sleep
432	from os import stat
433
434	from Numeric import choose
435
436	# Get NetCDF
437	retries = 0
438	file_open = False
439	while not file_open and retries < 10:
440	try:
441	fid = NetCDFFile(self.filename, 'a') # Open existing file
442	except IOError:
443	# This could happen if someone was reading the file.
444	# In that case, wait a while and try again
445	msg = 'Warning (store_timestep): File %s could not be opened'\
446	%self.filename
447	msg += ' - trying step %s again' %self.domain.time
448	print msg
449	retries += 1
450	sleep(1)
451	else:
452	file_open = True
453
454	if not file_open:
455	msg = 'File %s could not be opened for append' %self.filename
456	raise DataFileNotOpenError, msg
457
458
459
460	# Check to see if the file is already too big:
461	time = fid.variables['time']
462	i = len(time)+1
463	file_size = stat(self.filename)[6]
464	file_size_increase = file_size/i
465	if file_size + file_size_increase > self.max_size(2*self.recursion):
466	# In order to get the file name and start time correct,
467	# I change the domain.filename and domain.starttime.
468	# This is the only way to do this without changing
469	# other modules (I think).
470
471	# Write a filename addon that won't break swollens reader
472	# (10.sww is bad)
473	filename_ext = '_time_%s'%self.domain.time
474	filename_ext = filename_ext.replace('.', '_')
475
476	# Remember the old filename, then give domain a
477	# name with the extension
478	old_domain_filename = self.domain.get_name()
479	if not self.recursion:
480	self.domain.set_name(old_domain_filename+filename_ext)
481
482
483	# Change the domain starttime to the current time
484	old_domain_starttime = self.domain.starttime
485	self.domain.starttime = self.domain.time
486
487	# Build a new data_structure.
488	next_data_structure=\
489	Data_format_sww(self.domain, mode=self.mode,\
490	max_size = self.max_size,\
491	recursion = self.recursion+1)
492	if not self.recursion:
493	print ' file_size = %s'%file_size
494	print ' saving file to %s'%next_data_structure.filename
495	#set up the new data_structure
496	self.domain.writer = next_data_structure
497
498	#FIXME - could be cleaner to use domain.store_timestep etc.
499	next_data_structure.store_connectivity()
500	next_data_structure.store_timestep(names)
501	fid.sync()
502	fid.close()
503
504	#restore the old starttime and filename
505	self.domain.starttime = old_domain_starttime
506	self.domain.set_name(old_domain_filename)
507	else:
508	self.recursion = False
509	domain = self.domain
510
511	# Get the variables
512	time = fid.variables['time']
513	stage = fid.variables['stage']
514	xmomentum = fid.variables['xmomentum']
515	ymomentum = fid.variables['ymomentum']
516	i = len(time)
517	if type(names) not in [types.ListType, types.TupleType]:
518	names = [names]
519
520	if 'stage' in names and 'xmomentum' in names and \
521	'ymomentum' in names:
522
523	# Get stage
524	Q = domain.quantities['stage']
525	A,_ = Q.get_vertex_values(xy = False,
526	precision = self.precision)
527	z = fid.variables['elevation']
528	stage = choose(A-z[:] >= self.minimum_storable_height,
529	(z[:], A))
530
531	# Get xmomentum
532	Q = domain.quantities['xmomentum']
533	xmomentum, _ = Q.get_vertex_values(xy = False,
534	precision = self.precision)
535
536	# Get ymomentum
537	Q = domain.quantities['ymomentum']
538	ymomentum, _ = Q.get_vertex_values(xy = False,
539	precision = self.precision)
540
541	# Write quantities to NetCDF
542	self.writer.store_quantities(fid,
543	time=self.domain.time,
544	precision=self.precision,
545	stage=stage,
546	xmomentum=xmomentum,
547	ymomentum=ymomentum)
548	else:
549	# This is producing a sww that is not standard.
550	# Store time
551	time[i] = self.domain.time
552
553	for name in names:
554	# Get quantity
555	Q = domain.quantities[name]
556	A,V = Q.get_vertex_values(xy = False,
557	precision = self.precision)
558
559	# FIXME: Make this general (see below)
560	if name == 'stage':
561	z = fid.variables['elevation']
562	A = choose(A-z[:] >= self.minimum_storable_height,
563	(z[:], A))
564	stage[i,:] = A.astype(self.precision)
565	elif name == 'xmomentum':
566	xmomentum[i,:] = A.astype(self.precision)
567	elif name == 'ymomentum':
568	ymomentum[i,:] = A.astype(self.precision)
569
570	#As in....
571	#eval( name + '[i,:] = A.astype(self.precision)' )
572	#FIXME (Ole): But we need a UNIT test for that before
573	# refactoring
574
575
576
577	# Update extrema if requested
578	domain = self.domain
579	if domain.quantities_to_be_monitored is not None:
580	for q, info in domain.quantities_to_be_monitored.items():
581
582	if info['min'] is not None:
583	fid.variables[q + '.extrema'][0] = info['min']
584	fid.variables[q + '.min_location'][:] =\
585	info['min_location']
586	fid.variables[q + '.min_time'][0] = info['min_time']
587
588	if info['max'] is not None:
589	fid.variables[q + '.extrema'][1] = info['max']
590	fid.variables[q + '.max_location'][:] =\
591	info['max_location']
592	fid.variables[q + '.max_time'][0] = info['max_time']
593
594
595
596	#Flush and close
597	fid.sync()
598	fid.close()
599
600
601
602	# Class for handling checkpoints data
603	class Data_format_cpt(Data_format):
604	"""Interface to native NetCDF format (.cpt)
605	"""
606
607
608	def __init__(self, domain, mode = 'w'):
609	from Scientific.IO.NetCDF import NetCDFFile
610	from Numeric import Int, Float, Float
611
612	self.precision = Float #Use full precision
613
614	Data_format.__init__(self, domain, 'sww', mode)
615
616
617	# NetCDF file definition
618	fid = NetCDFFile(self.filename, mode)
619
620	if mode == 'w':
621	#Create new file
622	fid.institution = 'Geoscience Australia'
623	fid.description = 'Checkpoint data'
624	#fid.smooth = domain.smooth
625	fid.order = domain.default_order
626
627	# dimension definitions
628	fid.createDimension('number_of_volumes', self.number_of_volumes)
629	fid.createDimension('number_of_vertices', 3)
630
631	#Store info at all vertices (no smoothing)
632	fid.createDimension('number_of_points', 3*self.number_of_volumes)
633	fid.createDimension('number_of_timesteps', None) #extensible
634
635	# variable definitions
636
637	#Mesh
638	fid.createVariable('x', self.precision, ('number_of_points',))
639	fid.createVariable('y', self.precision, ('number_of_points',))
640
641
642	fid.createVariable('volumes', Int, ('number_of_volumes',
643	'number_of_vertices'))
644
645	fid.createVariable('time', self.precision,
646	('number_of_timesteps',))
647
648	#Allocate space for all quantities
649	for name in domain.quantities.keys():
650	fid.createVariable(name, self.precision,
651	('number_of_timesteps',
652	'number_of_points'))
653
654	#Close
655	fid.close()
656
657
658	def store_checkpoint(self):
659	"""
660	Write x,y coordinates of triangles.
661	Write connectivity (
662	constituting
663	the bed elevation.
664	"""
665
666	from Scientific.IO.NetCDF import NetCDFFile
667
668	from Numeric import concatenate
669
670	domain = self.domain
671
672	#Get NetCDF
673	fid = NetCDFFile(self.filename, 'a') #Open existing file for append
674
675	# Get the variables
676	x = fid.variables['x']
677	y = fid.variables['y']
678
679	volumes = fid.variables['volumes']
680
681	# Get X, Y and bed elevation Z
682	Q = domain.quantities['elevation']
683	X,Y,Z,V = Q.get_vertex_values(xy=True,
684	precision = self.precision)
685
686
687
688	x[:] = X.astype(self.precision)
689	y[:] = Y.astype(self.precision)
690	z[:] = Z.astype(self.precision)
691
692	volumes[:] = V
693
694	#Close
695	fid.close()
696
697
698	def store_timestep(self, name):
699	"""Store time and named quantity to file
700	"""
701	from Scientific.IO.NetCDF import NetCDFFile
702	from time import sleep
703
704	#Get NetCDF
705	retries = 0
706	file_open = False
707	while not file_open and retries < 10:
708	try:
709	fid = NetCDFFile(self.filename, 'a') #Open existing file
710	except IOError:
711	#This could happen if someone was reading the file.
712	#In that case, wait a while and try again
713	msg = 'Warning (store_timestep): File %s could not be opened'\
714	%self.filename
715	msg += ' - trying again'
716	print msg
717	retries += 1
718	sleep(1)
719	else:
720	file_open = True
721
722	if not file_open:
723	msg = 'File %s could not be opened for append' %self.filename
724	raise DataFileNotOPenError, msg
725
726
727	domain = self.domain
728
729	# Get the variables
730	time = fid.variables['time']
731	stage = fid.variables['stage']
732	i = len(time)
733
734	#Store stage
735	time[i] = self.domain.time
736
737	# Get quantity
738	Q = domain.quantities[name]
739	A,V = Q.get_vertex_values(xy=False,
740	precision = self.precision)
741
742	stage[i,:] = A.astype(self.precision)
743
744	#Flush and close
745	fid.sync()
746	fid.close()
747
748
749	#### NED is national exposure database (name changed to NEXIS)
750
751	LAT_TITLE = 'LATITUDE'
752	LONG_TITLE = 'LONGITUDE'
753	X_TITLE = 'x'
754	Y_TITLE = 'y'
755	class Exposure_csv:
756	def __init__(self,file_name, latitude_title=LAT_TITLE,
757	longitude_title=LONG_TITLE, is_x_y_locations=None,
758	x_title=X_TITLE, y_title=Y_TITLE,
759	refine_polygon=None, title_check_list=None):
760	"""
761	This class is for handling the exposure csv file.
762	It reads the file in and converts the lats and longs to a geospatial
763	data object.
764	Use the methods to read and write columns.
765
766	The format of the csv files it reads is;
767	The first row is a title row.
768	comma's are the delimiters
769	each column is a 'set' of data
770
771	Feel free to use/expand it to read other csv files.
772
773
774	It is not for adding and deleting rows
775
776	Can geospatial handle string attributes? It's not made for them.
777	Currently it can't load and save string att's.
778
779	So just use geospatial to hold the x, y and georef? Bad, since
780	different att's are in diferent structures. Not so bad, the info
781	to write if the .csv file is saved is in attribute_dic
782
783	The location info is in the geospatial attribute.
784
785
786	"""
787	self._file_name = file_name
788	self._geospatial = None #
789
790	# self._attribute_dic is a dictionary.
791	#The keys are the column titles.
792	#The values are lists of column data
793
794	# self._title_index_dic is a dictionary.
795	#The keys are the column titles.
796	#The values are the index positions of file columns.
797	self._attribute_dic, self._title_index_dic = \
798	csv2dict(self._file_name, title_check_list=title_check_list)
799	try:
800	#Have code here that handles caps or lower
801	lats = self._attribute_dic[latitude_title]
802	longs = self._attribute_dic[longitude_title]
803
804	except KeyError:
805	# maybe a warning..
806	#Let's see if this works..
807	if False != is_x_y_locations:
808	is_x_y_locations = True
809	pass
810	else:
811	self._geospatial = Geospatial_data(latitudes = lats,
812	longitudes = longs)
813
814	if is_x_y_locations is True:
815	if self._geospatial is not None:
816	pass #fixme throw an error
817	try:
818	xs = self._attribute_dic[x_title]
819	ys = self._attribute_dic[y_title]
820	points = [[float(i),float(j)] for i,j in map(None,xs,ys)]
821	except KeyError:
822	# maybe a warning..
823	msg = "Could not find location information."
824	raise TitleValueError, msg
825	else:
826	self._geospatial = Geospatial_data(data_points=points)
827
828	# create a list of points that are in the refining_polygon
829	# described by a list of indexes representing the points
830
831	def __cmp__(self, other):
832	#print "self._attribute_dic",self._attribute_dic
833	#print "other._attribute_dic",other._attribute_dic
834	#print "self._title_index_dic", self._title_index_dic
835	#print "other._title_index_dic", other._title_index_dic
836
837	#check that a is an instance of this class
838	if isinstance(self, type(other)):
839	result = cmp(self._attribute_dic, other._attribute_dic)
840	if result <>0:
841	return result
842	# The order of the columns is important. Therefore..
843	result = cmp(self._title_index_dic, other._title_index_dic)
844	if result <>0:
845	return result
846	for self_ls, other_ls in map(None,self._attribute_dic, \
847	other._attribute_dic):
848	result = cmp(self._attribute_dic[self_ls],
849	other._attribute_dic[other_ls])
850	if result <>0:
851	return result
852	return 0
853	else:
854	return 1
855
856
857	def get_column(self, column_name, use_refind_polygon=False):
858	"""
859	Given a column name return a list of the column values
860
861	Note, the type of the values will be String!
862	do this to change a list of strings to a list of floats
863	time = [float(x) for x in time]
864
865	Not implemented:
866	if use_refind_polygon is True, only return values in the
867	refined polygon
868	"""
869	if not self._attribute_dic.has_key(column_name):
870	msg = 'Therer is no column called %s!' %column_name
871	raise TitleValueError, msg
872	return self._attribute_dic[column_name]
873
874
875	def get_value(self, value_column_name,
876	known_column_name,
877	known_values,
878	use_refind_polygon=False):
879	"""
880	Do linear interpolation on the known_colum, using the known_value,
881	to return a value of the column_value_name.
882	"""
883	pass
884
885
886	def get_location(self, use_refind_polygon=False):
887	"""
888	Return a geospatial object which describes the
889	locations of the location file.
890
891	Note, if there is not location info, this returns None.
892
893	Not implemented:
894	if use_refind_polygon is True, only return values in the
895	refined polygon
896	"""
897	return self._geospatial
898
899	def set_column(self, column_name, column_values, overwrite=False):
900	"""
901	Add a column to the 'end' (with the right most column being the end)
902	of the csv file.
903
904	Set overwrite to True if you want to overwrite a column.
905
906	Note, in column_name white space is removed and case is not checked.
907	Precondition
908	The column_name and column_values cannot have comma's in it.
909	"""
910
911	value_row_count = \
912	len(self._attribute_dic[self._title_index_dic.keys()[0]])
913	if len(column_values) <> value_row_count:
914	msg = 'The number of column values must equal the number of rows.'
915	raise DataMissingValuesError, msg
916
917	if self._attribute_dic.has_key(column_name):
918	if not overwrite:
919	msg = 'Column name %s already in use!' %column_name
920	raise TitleValueError, msg
921	else:
922	# New title. Add it to the title index.
923	self._title_index_dic[column_name] = len(self._title_index_dic)
924	self._attribute_dic[column_name] = column_values
925	#print "self._title_index_dic[column_name]",self._title_index_dic
926
927	def save(self, file_name=None):
928	"""
929	Save the exposure csv file
930	"""
931	if file_name is None:
932	file_name = self._file_name
933
934	fd = open(file_name,'wb')
935	writer = csv.writer(fd)
936
937	#Write the title to a cvs file
938	line = [None]* len(self._title_index_dic)
939	for title in self._title_index_dic.iterkeys():
940	line[self._title_index_dic[title]]= title
941	writer.writerow(line)
942
943	# Write the values to a cvs file
944	value_row_count = \
945	len(self._attribute_dic[self._title_index_dic.keys()[0]])
946	for row_i in range(value_row_count):
947	line = [None]* len(self._title_index_dic)
948	for title in self._title_index_dic.iterkeys():
949	line[self._title_index_dic[title]]= \
950	self._attribute_dic[title][row_i]
951	writer.writerow(line)
952
953
954	def csv2dict(file_name, title_check_list=None):
955	"""
956	Load in the csv as a dic, title as key and column info as value, .
957	Also, create a dic, title as key and column index as value,
958	to keep track of the column order.
959
960	WARNING: Vaules are returned as strings.
961	do this to change a list of strings to a list of floats
962	time = [float(x) for x in time]
963	"""
964
965	#
966	attribute_dic = {}
967	title_index_dic = {}
968	titles_stripped = [] # list of titles
969	reader = csv.reader(file(file_name))
970
971	# Read in and manipulate the title info
972	titles = reader.next()
973	for i,title in enumerate(titles):
974	titles_stripped.append(title.strip())
975	title_index_dic[title.strip()] = i
976	title_count = len(titles_stripped)
977	#print "title_index_dic",title_index_dic
978	if title_check_list is not None:
979	for title_check in title_check_list:
980	#msg = "Reading error. This row is not present ", title_check
981	#assert title_index_dic.has_key(title_check), msg
982	if not title_index_dic.has_key(title_check):
983	#reader.close()
984	msg = "Reading error. This row is not present ", \
985	title_check
986	raise IOError, msg
987
988
989
990	#create a dic of colum values, indexed by column title
991	for line in reader:
992	if len(line) <> title_count:
993	raise IOError #FIXME make this nicer
994	for i, value in enumerate(line):
995	attribute_dic.setdefault(titles_stripped[i],[]).append(value)
996
997	return attribute_dic, title_index_dic
998
999
1000	#Auxiliary
1001	def write_obj(filename,x,y,z):
1002	"""Store x,y,z vectors into filename (obj format)
1003	Vectors are assumed to have dimension (M,3) where
1004	M corresponds to the number elements.
1005	triangles are assumed to be disconnected
1006
1007	The three numbers in each vector correspond to three vertices,
1008
1009	e.g. the x coordinate of vertex 1 of element i is in x[i,1]
1010
1011	"""
1012	#print 'Writing obj to %s' % filename
1013
1014	import os.path
1015
1016	root, ext = os.path.splitext(filename)
1017	if ext == '.obj':
1018	FN = filename
1019	else:
1020	FN = filename + '.obj'
1021
1022
1023	outfile = open(FN, 'wb')
1024	outfile.write("# Triangulation as an obj file\n")
1025
1026	M, N = x.shape
1027	assert N==3 #Assuming three vertices per element
1028
1029	for i in range(M):
1030	for j in range(N):
1031	outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
1032
1033	for i in range(M):
1034	base = i*N
1035	outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
1036
1037	outfile.close()
1038
1039
1040	#########################################################
1041	#Conversion routines
1042	########################################################
1043
1044	def sww2obj(basefilename, size):
1045	"""Convert netcdf based data output to obj
1046	"""
1047	from Scientific.IO.NetCDF import NetCDFFile
1048
1049	from Numeric import Float, zeros
1050
1051	#Get NetCDF
1052	FN = create_filename('.', basefilename, 'sww', size)
1053	print 'Reading from ', FN
1054	fid = NetCDFFile(FN, 'r') #Open existing file for read
1055
1056
1057	# Get the variables
1058	x = fid.variables['x']
1059	y = fid.variables['y']
1060	z = fid.variables['elevation']
1061	time = fid.variables['time']
1062	stage = fid.variables['stage']
1063
1064	M = size #Number of lines
1065	xx = zeros((M,3), Float)
1066	yy = zeros((M,3), Float)
1067	zz = zeros((M,3), Float)
1068
1069	for i in range(M):
1070	for j in range(3):
1071	xx[i,j] = x[i+j*M]
1072	yy[i,j] = y[i+j*M]
1073	zz[i,j] = z[i+j*M]
1074
1075	#Write obj for bathymetry
1076	FN = create_filename('.', basefilename, 'obj', size)
1077	write_obj(FN,xx,yy,zz)
1078
1079
1080	#Now read all the data with variable information, combine with
1081	#x,y info and store as obj
1082
1083	for k in range(len(time)):
1084	t = time[k]
1085	print 'Processing timestep %f' %t
1086
1087	for i in range(M):
1088	for j in range(3):
1089	zz[i,j] = stage[k,i+j*M]
1090
1091
1092	#Write obj for variable data
1093	#FN = create_filename(basefilename, 'obj', size, time=t)
1094	FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
1095	write_obj(FN,xx,yy,zz)
1096
1097
1098	def dat2obj(basefilename):
1099	"""Convert line based data output to obj
1100	FIXME: Obsolete?
1101	"""
1102
1103	import glob, os
1104	from anuga.config import data_dir
1105
1106
1107	#Get bathymetry and x,y's
1108	lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
1109
1110	from Numeric import zeros, Float
1111
1112	M = len(lines) #Number of lines
1113	x = zeros((M,3), Float)
1114	y = zeros((M,3), Float)
1115	z = zeros((M,3), Float)
1116
1117	##i = 0
1118	for i, line in enumerate(lines):
1119	tokens = line.split()
1120	values = map(float,tokens)
1121
1122	for j in range(3):
1123	x[i,j] = values[j*3]
1124	y[i,j] = values[j*3+1]
1125	z[i,j] = values[j*3+2]
1126
1127	##i += 1
1128
1129
1130	#Write obj for bathymetry
1131	write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
1132
1133
1134	#Now read all the data files with variable information, combine with
1135	#x,y info
1136	#and store as obj
1137
1138	files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
1139
1140	for filename in files:
1141	print 'Processing %s' % filename
1142
1143	lines = open(data_dir+os.sep+filename,'r').readlines()
1144	assert len(lines) == M
1145	root, ext = os.path.splitext(filename)
1146
1147	#Get time from filename
1148	i0 = filename.find('_time=')
1149	if i0 == -1:
1150	#Skip bathymetry file
1151	continue
1152
1153	i0 += 6 #Position where time starts
1154	i1 = filename.find('.dat')
1155
1156	if i1 > i0:
1157	t = float(filename[i0:i1])
1158	else:
1159	raise DataTimeError, 'Hmmmm'
1160
1161
1162
1163	##i = 0
1164	for i, line in enumerate(lines):
1165	tokens = line.split()
1166	values = map(float,tokens)
1167
1168	for j in range(3):
1169	z[i,j] = values[j]
1170
1171	##i += 1
1172
1173	#Write obj for variable data
1174	write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
1175
1176
1177	def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
1178	"""Read netcdf filename1, pick timesteps first:step:last and save to
1179	nettcdf file filename2
1180	"""
1181	from Scientific.IO.NetCDF import NetCDFFile
1182
1183	#Get NetCDF
1184	infile = NetCDFFile(filename1, 'r') #Open existing file for read
1185	outfile = NetCDFFile(filename2, 'w') #Open new file
1186
1187
1188	#Copy dimensions
1189	for d in infile.dimensions:
1190	outfile.createDimension(d, infile.dimensions[d])
1191
1192	for name in infile.variables:
1193	var = infile.variables[name]
1194	outfile.createVariable(name, var.typecode(), var.dimensions)
1195
1196
1197	#Copy the static variables
1198	for name in infile.variables:
1199	if name == 'time' or name == 'stage':
1200	pass
1201	else:
1202	#Copy
1203	outfile.variables[name][:] = infile.variables[name][:]
1204
1205	#Copy selected timesteps
1206	time = infile.variables['time']
1207	stage = infile.variables['stage']
1208
1209	newtime = outfile.variables['time']
1210	newstage = outfile.variables['stage']
1211
1212	if last is None:
1213	last = len(time)
1214
1215	selection = range(first, last, step)
1216	for i, j in enumerate(selection):
1217	print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
1218	newtime[i] = time[j]
1219	newstage[i,:] = stage[j,:]
1220
1221	#Close
1222	infile.close()
1223	outfile.close()
1224
1225
1226	#Get data objects
1227	def get_dataobject(domain, mode='w'):
1228	"""Return instance of class of given format using filename
1229	"""
1230
1231	cls = eval('Data_format_%s' %domain.format)
1232	return cls(domain, mode)
1233
1234
1235
1236
1237	def dem2pts(basename_in, basename_out=None,
1238	easting_min=None, easting_max=None,
1239	northing_min=None, northing_max=None,
1240	use_cache=False, verbose=False,):
1241	"""Read Digitial Elevation model from the following NetCDF format (.dem)
1242
1243	Example:
1244
1245	ncols 3121
1246	nrows 1800
1247	xllcorner 722000
1248	yllcorner 5893000
1249	cellsize 25
1250	NODATA_value -9999
1251	138.3698 137.4194 136.5062 135.5558 ..........
1252
1253	Convert to NetCDF pts format which is
1254
1255	points: (Nx2) Float array
1256	elevation: N Float array
1257	"""
1258
1259
1260
1261	kwargs = {'basename_out': basename_out,
1262	'easting_min': easting_min,
1263	'easting_max': easting_max,
1264	'northing_min': northing_min,
1265	'northing_max': northing_max,
1266	'verbose': verbose}
1267
1268	if use_cache is True:
1269	from caching import cache
1270	result = cache(_dem2pts, basename_in, kwargs,
1271	dependencies = [basename_in + '.dem'],
1272	verbose = verbose)
1273
1274	else:
1275	result = apply(_dem2pts, [basename_in], kwargs)
1276
1277	return result
1278
1279
1280	def _dem2pts(basename_in, basename_out=None, verbose=False,
1281	easting_min=None, easting_max=None,
1282	northing_min=None, northing_max=None):
1283	"""Read Digitial Elevation model from the following NetCDF format (.dem)
1284
1285	Internal function. See public function dem2pts for details.
1286	"""
1287
1288	#FIXME: Can this be written feasibly using write_pts?
1289
1290	import os
1291	from Scientific.IO.NetCDF import NetCDFFile
1292	from Numeric import Float, zeros, reshape, sum
1293
1294	root = basename_in
1295
1296	#Get NetCDF
1297	infile = NetCDFFile(root + '.dem', 'r') #Open existing netcdf file for read
1298
1299	if verbose: print 'Reading DEM from %s' %(root + '.dem')
1300
1301	ncols = infile.ncols[0]
1302	nrows = infile.nrows[0]
1303	xllcorner = infile.xllcorner[0] #Easting of lower left corner
1304	yllcorner = infile.yllcorner[0] #Northing of lower left corner
1305	cellsize = infile.cellsize[0]
1306	NODATA_value = infile.NODATA_value[0]
1307	dem_elevation = infile.variables['elevation']
1308
1309	zone = infile.zone[0]
1310	false_easting = infile.false_easting[0]
1311	false_northing = infile.false_northing[0]
1312
1313	#Text strings
1314	projection = infile.projection
1315	datum = infile.datum
1316	units = infile.units
1317
1318
1319	#Get output file
1320	if basename_out == None:
1321	ptsname = root + '.pts'
1322	else:
1323	ptsname = basename_out + '.pts'
1324
1325	if verbose: print 'Store to NetCDF file %s' %ptsname
1326	# NetCDF file definition
1327	outfile = NetCDFFile(ptsname, 'w')
1328
1329	#Create new file
1330	outfile.institution = 'Geoscience Australia'
1331	outfile.description = 'NetCDF pts format for compact and portable storage ' +\
1332	'of spatial point data'
1333	#assign default values
1334	if easting_min is None: easting_min = xllcorner
1335	if easting_max is None: easting_max = xllcorner + ncols*cellsize
1336	if northing_min is None: northing_min = yllcorner
1337	if northing_max is None: northing_max = yllcorner + nrows*cellsize
1338
1339	#compute offsets to update georeferencing
1340	easting_offset = xllcorner - easting_min
1341	northing_offset = yllcorner - northing_min
1342
1343	#Georeferencing
1344	outfile.zone = zone
1345	outfile.xllcorner = easting_min #Easting of lower left corner
1346	outfile.yllcorner = northing_min #Northing of lower left corner
1347	outfile.false_easting = false_easting
1348	outfile.false_northing = false_northing
1349
1350	outfile.projection = projection
1351	outfile.datum = datum
1352	outfile.units = units
1353
1354
1355	#Grid info (FIXME: probably not going to be used, but heck)
1356	outfile.ncols = ncols
1357	outfile.nrows = nrows
1358
1359	dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
1360	totalnopoints = nrows*ncols
1361
1362	# calculating number of NODATA_values for each row in clipped region
1363	#FIXME: use array operations to do faster
1364	nn = 0
1365	k = 0
1366	i1_0 = 0
1367	j1_0 = 0
1368	thisj = 0
1369	thisi = 0
1370	for i in range(nrows):
1371	y = (nrows-i-1)*cellsize + yllcorner
1372	for j in range(ncols):
1373	x = j*cellsize + xllcorner
1374	if easting_min <= x <= easting_max and \
1375	northing_min <= y <= northing_max:
1376	thisj = j
1377	thisi = i
1378	if dem_elevation_r[i,j] == NODATA_value: nn += 1
1379
1380	if k == 0:
1381	i1_0 = i
1382	j1_0 = j
1383	k += 1
1384
1385	index1 = j1_0
1386	index2 = thisj
1387
1388	# dimension definitions
1389	nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
1390	ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
1391
1392	clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
1393	nopoints = clippednopoints-nn
1394
1395	clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
1396
1397	if verbose:
1398	print 'There are %d values in the elevation' %totalnopoints
1399	print 'There are %d values in the clipped elevation' %clippednopoints
1400	print 'There are %d NODATA_values in the clipped elevation' %nn
1401
1402	outfile.createDimension('number_of_points', nopoints)
1403	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1404
1405	# variable definitions
1406	outfile.createVariable('points', Float, ('number_of_points',
1407	'number_of_dimensions'))
1408	outfile.createVariable('elevation', Float, ('number_of_points',))
1409
1410	# Get handles to the variables
1411	points = outfile.variables['points']
1412	elevation = outfile.variables['elevation']
1413
1414	lenv = index2-index1+1
1415	#Store data
1416	global_index = 0
1417	#for i in range(nrows):
1418	for i in range(i1_0,thisi+1,1):
1419	if verbose and i%((nrows+10)/10)==0:
1420	print 'Processing row %d of %d' %(i, nrows)
1421
1422	lower_index = global_index
1423
1424	v = dem_elevation_r[i,index1:index2+1]
1425	no_NODATA = sum(v == NODATA_value)
1426	if no_NODATA > 0:
1427	newcols = lenv - no_NODATA #ncols_in_bounding_box - no_NODATA
1428	else:
1429	newcols = lenv #ncols_in_bounding_box
1430
1431	telev = zeros(newcols, Float)
1432	tpoints = zeros((newcols, 2), Float)
1433
1434	local_index = 0
1435
1436	y = (nrows-i-1)*cellsize + yllcorner
1437	#for j in range(ncols):
1438	for j in range(j1_0,index2+1,1):
1439
1440	x = j*cellsize + xllcorner
1441	if easting_min <= x <= easting_max and \
1442	northing_min <= y <= northing_max and \
1443	dem_elevation_r[i,j] <> NODATA_value:
1444	tpoints[local_index, :] = [x-easting_min,y-northing_min]
1445	telev[local_index] = dem_elevation_r[i, j]
1446	global_index += 1
1447	local_index += 1
1448
1449	upper_index = global_index
1450
1451	if upper_index == lower_index + newcols:
1452	points[lower_index:upper_index, :] = tpoints
1453	elevation[lower_index:upper_index] = telev
1454
1455	assert global_index == nopoints, 'index not equal to number of points'
1456
1457	infile.close()
1458	outfile.close()
1459
1460
1461
1462	def _read_hecras_cross_sections(lines):
1463	"""Return block of surface lines for each cross section
1464	Starts with SURFACE LINE,
1465	Ends with END CROSS-SECTION
1466	"""
1467
1468	points = []
1469
1470	reading_surface = False
1471	for i, line in enumerate(lines):
1472
1473	if len(line.strip()) == 0: #Ignore blanks
1474	continue
1475
1476	if lines[i].strip().startswith('SURFACE LINE'):
1477	reading_surface = True
1478	continue
1479
1480	if lines[i].strip().startswith('END') and reading_surface:
1481	yield points
1482	reading_surface = False
1483	points = []
1484
1485	if reading_surface:
1486	fields = line.strip().split(',')
1487	easting = float(fields[0])
1488	northing = float(fields[1])
1489	elevation = float(fields[2])
1490	points.append([easting, northing, elevation])
1491
1492
1493
1494
1495	def hecras_cross_sections2pts(basename_in,
1496	basename_out=None,
1497	verbose=False):
1498	"""Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
1499
1500	Example:
1501
1502
1503	# RAS export file created on Mon 15Aug2005 11:42
1504	# by HEC-RAS Version 3.1.1
1505
1506	BEGIN HEADER:
1507	UNITS: METRIC
1508	DTM TYPE: TIN
1509	DTM: v:\1\cit\perth_topo\river_tin
1510	STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
1511	CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
1512	MAP PROJECTION: UTM
1513	PROJECTION ZONE: 50
1514	DATUM: AGD66
1515	VERTICAL DATUM:
1516	NUMBER OF REACHES: 19
1517	NUMBER OF CROSS-SECTIONS: 14206
1518	END HEADER:
1519
1520
1521	Only the SURFACE LINE data of the following form will be utilised
1522
1523	CROSS-SECTION:
1524	STREAM ID:Southern-Wungong
1525	REACH ID:Southern-Wungong
1526	STATION:19040.*
1527	CUT LINE:
1528	405548.671603161 , 6438142.7594925
1529	405734.536092045 , 6438326.10404912
1530	405745.130459356 , 6438331.48627354
1531	405813.89633823 , 6438368.6272789
1532	SURFACE LINE:
1533	405548.67, 6438142.76, 35.37
1534	405552.24, 6438146.28, 35.41
1535	405554.78, 6438148.78, 35.44
1536	405555.80, 6438149.79, 35.44
1537	405559.37, 6438153.31, 35.45
1538	405560.88, 6438154.81, 35.44
1539	405562.93, 6438156.83, 35.42
1540	405566.50, 6438160.35, 35.38
1541	405566.99, 6438160.83, 35.37
1542	...
1543	END CROSS-SECTION
1544
1545	Convert to NetCDF pts format which is
1546
1547	points: (Nx2) Float array
1548	elevation: N Float array
1549	"""
1550
1551	import os
1552	from Scientific.IO.NetCDF import NetCDFFile
1553	from Numeric import Float, zeros, reshape
1554
1555	root = basename_in
1556
1557	#Get ASCII file
1558	infile = open(root + '.sdf', 'r') #Open SDF file for read
1559
1560	if verbose: print 'Reading DEM from %s' %(root + '.sdf')
1561
1562	lines = infile.readlines()
1563	infile.close()
1564
1565	if verbose: print 'Converting to pts format'
1566
1567	i = 0
1568	while lines[i].strip() == '' or lines[i].strip().startswith('#'):
1569	i += 1
1570
1571	assert lines[i].strip().upper() == 'BEGIN HEADER:'
1572	i += 1
1573
1574	assert lines[i].strip().upper().startswith('UNITS:')
1575	units = lines[i].strip().split()[1]
1576	i += 1
1577
1578	assert lines[i].strip().upper().startswith('DTM TYPE:')
1579	i += 1
1580
1581	assert lines[i].strip().upper().startswith('DTM:')
1582	i += 1
1583
1584	assert lines[i].strip().upper().startswith('STREAM')
1585	i += 1
1586
1587	assert lines[i].strip().upper().startswith('CROSS')
1588	i += 1
1589
1590	assert lines[i].strip().upper().startswith('MAP PROJECTION:')
1591	projection = lines[i].strip().split(':')[1]
1592	i += 1
1593
1594	assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
1595	zone = int(lines[i].strip().split(':')[1])
1596	i += 1
1597
1598	assert lines[i].strip().upper().startswith('DATUM:')
1599	datum = lines[i].strip().split(':')[1]
1600	i += 1
1601
1602	assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
1603	i += 1
1604
1605	assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
1606	i += 1
1607
1608	assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
1609	number_of_cross_sections = int(lines[i].strip().split(':')[1])
1610	i += 1
1611
1612
1613	#Now read all points
1614	points = []
1615	elevation = []
1616	for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
1617	for k, entry in enumerate(entries):
1618	points.append(entry[:2])
1619	elevation.append(entry[2])
1620
1621
1622	msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
1623	%(j+1, number_of_cross_sections)
1624	assert j+1 == number_of_cross_sections, msg
1625
1626	#Get output file
1627	if basename_out == None:
1628	ptsname = root + '.pts'
1629	else:
1630	ptsname = basename_out + '.pts'
1631
1632	geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
1633	geo = Geospatial_data(points, {"elevation":elevation},
1634	verbose=verbose, geo_reference=geo_ref)
1635	geo.export_points_file(ptsname)
1636
1637	def export_grid(basename_in, extra_name_out = None,
1638	quantities = None, # defaults to elevation
1639	timestep = None,
1640	reduction = None,
1641	cellsize = 10,
1642	NODATA_value = -9999,
1643	easting_min = None,
1644	easting_max = None,
1645	northing_min = None,
1646	northing_max = None,
1647	verbose = False,
1648	origin = None,
1649	datum = 'WGS84',
1650	format = 'ers'):
1651	"""
1652
1653	Wrapper for sww2dem. - see sww2dem to find out what most of the
1654	parameters do.
1655
1656	Quantities is a list of quantities. Each quantity will be
1657	calculated for each sww file.
1658
1659	This returns the basenames of the files returned, which is made up
1660	of the dir and all of the file name, except the extension.
1661
1662	This function returns the names of the files produced.
1663
1664	It will also produce as many output files as there are input sww files.
1665	"""
1666
1667	if quantities is None:
1668	quantities = ['elevation']
1669
1670	if type(quantities) is str:
1671	quantities = [quantities]
1672
1673	# How many sww files are there?
1674	dir, base = os.path.split(basename_in)
1675	#print "basename_in",basename_in
1676	#print "base",base
1677
1678	iterate_over = get_all_swwfiles(dir,base,verbose)
1679
1680	if dir == "":
1681	dir = "." # Unix compatibility
1682
1683	files_out = []
1684	#print 'sww_file',iterate_over
1685	for sww_file in iterate_over:
1686	for quantity in quantities:
1687	if extra_name_out is None:
1688	basename_out = sww_file + '_' + quantity
1689	else:
1690	basename_out = sww_file + '_' + quantity + '_' \
1691	+ extra_name_out
1692	# print "basename_out", basename_out
1693
1694	file_out = sww2dem(dir+sep+sww_file, dir+sep+basename_out,
1695	quantity,
1696	timestep,
1697	reduction,
1698	cellsize,
1699	NODATA_value,
1700	easting_min,
1701	easting_max,
1702	northing_min,
1703	northing_max,
1704	verbose,
1705	origin,
1706	datum,
1707	format)
1708	files_out.append(file_out)
1709	#print "basenames_out after",basenames_out
1710	return files_out
1711
1712
1713	def get_timeseries(production_dirs, output_dir, scenario_name, gauges_dir_name,
1714	plot_quantity, generate_fig = False,
1715	reportname = None, surface = False, time_min = None,
1716	time_max = None, title_on = False, verbose = True,
1717	nodes=None):
1718	"""
1719	nodes - number of processes used.
1720
1721	warning - this function has no tests
1722	"""
1723	if reportname == None:
1724	report = False
1725	else:
1726	report = True
1727
1728	if nodes is None:
1729	is_parallel = False
1730	else:
1731	is_parallel = True
1732
1733	# Generate figures
1734	swwfiles = {}
1735
1736	if is_parallel is True:
1737	for i in range(nodes):
1738	print 'Sending node %d of %d' %(i,nodes)
1739	swwfiles = {}
1740	if not reportname == None:
1741	reportname = report_name + '_%s' %(i)
1742	for label_id in production_dirs.keys():
1743	if label_id == 'boundaries':
1744	swwfile = best_boundary_sww
1745	else:
1746	file_loc = output_dir + label_id + sep
1747	sww_extra = '_P%s_%s' %(i,nodes)
1748	swwfile = file_loc + scenario_name + sww_extra + '.sww'
1749	print 'swwfile',swwfile
1750	swwfiles[swwfile] = label_id
1751
1752	texname, elev_output = sww2timeseries(swwfiles,
1753	gauges_dir_name,
1754	production_dirs,
1755	report = report,
1756	reportname = reportname,
1757	plot_quantity = plot_quantity,
1758	generate_fig = generate_fig,
1759	surface = surface,
1760	time_min = time_min,
1761	time_max = time_max,
1762	title_on = title_on,
1763	verbose = verbose)
1764	else:
1765	for label_id in production_dirs.keys():
1766	if label_id == 'boundaries':
1767	print 'boundaries'
1768	file_loc = project.boundaries_in_dir
1769	swwfile = project.boundaries_dir_name3 + '.sww'
1770	# swwfile = boundary_dir_filename
1771	else:
1772	file_loc = output_dir + label_id + sep
1773	swwfile = file_loc + scenario_name + '.sww'
1774	swwfiles[swwfile] = label_id
1775
1776	texname, elev_output = sww2timeseries(swwfiles,
1777	gauges_dir_name,
1778	production_dirs,
1779	report = report,
1780	reportname = reportname,
1781	plot_quantity = plot_quantity,
1782	generate_fig = generate_fig,
1783	surface = surface,
1784	time_min = time_min,
1785	time_max = time_max,
1786	title_on = title_on,
1787	verbose = verbose)
1788
1789
1790
1791	def sww2dem(basename_in, basename_out = None,
1792	quantity = None, # defaults to elevation
1793	timestep = None,
1794	reduction = None,
1795	cellsize = 10,
1796	NODATA_value = -9999,
1797	easting_min = None,
1798	easting_max = None,
1799	northing_min = None,
1800	northing_max = None,
1801	verbose = False,
1802	origin = None,
1803	datum = 'WGS84',
1804	format = 'ers'):
1805
1806	"""Read SWW file and convert to Digitial Elevation model format
1807	(.asc or .ers)
1808
1809	Example (ASC):
1810
1811	ncols 3121
1812	nrows 1800
1813	xllcorner 722000
1814	yllcorner 5893000
1815	cellsize 25
1816	NODATA_value -9999
1817	138.3698 137.4194 136.5062 135.5558 ..........
1818
1819	Also write accompanying file with same basename_in but extension .prj
1820	used to fix the UTM zone, datum, false northings and eastings.
1821
1822	The prj format is assumed to be as
1823
1824	Projection UTM
1825	Zone 56
1826	Datum WGS84
1827	Zunits NO
1828	Units METERS
1829	Spheroid WGS84
1830	Xshift 0.0000000000
1831	Yshift 10000000.0000000000
1832	Parameters
1833
1834
1835	The parameter quantity must be the name of an existing quantity or
1836	an expression involving existing quantities. The default is
1837	'elevation'. Quantity is not a list of quantities.
1838
1839	if timestep (an index) is given, output quantity at that timestep
1840
1841	if reduction is given use that to reduce quantity over all timesteps.
1842
1843	datum
1844
1845	format can be either 'asc' or 'ers'
1846	"""
1847
1848	import sys
1849	from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, \
1850	sometrue
1851	from Numeric import array2string
1852
1853	from anuga.utilities.polygon import inside_polygon, outside_polygon, \
1854	separate_points_by_polygon
1855	from anuga.abstract_2d_finite_volumes.util import \
1856	apply_expression_to_dictionary
1857
1858	msg = 'Format must be either asc or ers'
1859	assert format.lower() in ['asc', 'ers'], msg
1860
1861
1862	false_easting = 500000
1863	false_northing = 10000000
1864
1865	if quantity is None:
1866	quantity = 'elevation'
1867
1868	if reduction is None:
1869	reduction = max
1870
1871	if basename_out is None:
1872	basename_out = basename_in + '_%s' %quantity
1873
1874	if quantity_formula.has_key(quantity):
1875	quantity = quantity_formula[quantity]
1876
1877	swwfile = basename_in + '.sww'
1878	demfile = basename_out + '.' + format
1879	# Note the use of a .ers extension is optional (write_ermapper_grid will
1880	# deal with either option
1881
1882	#if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors
1883
1884	# Read sww file
1885	if verbose:
1886	print 'Reading from %s' %swwfile
1887	print 'Output directory is %s' %basename_out
1888
1889	from Scientific.IO.NetCDF import NetCDFFile
1890	fid = NetCDFFile(swwfile)
1891
1892	#Get extent and reference
1893	x = fid.variables['x'][:]
1894	y = fid.variables['y'][:]
1895	volumes = fid.variables['volumes'][:]
1896	if timestep is not None:
1897	times = fid.variables['time'][timestep]
1898	else:
1899	times = fid.variables['time'][:]
1900
1901	number_of_timesteps = fid.dimensions['number_of_timesteps']
1902	number_of_points = fid.dimensions['number_of_points']
1903
1904	if origin is None:
1905
1906	# Get geo_reference
1907	# sww files don't have to have a geo_ref
1908	try:
1909	geo_reference = Geo_reference(NetCDFObject=fid)
1910	except AttributeError, e:
1911	geo_reference = Geo_reference() # Default georef object
1912
1913	xllcorner = geo_reference.get_xllcorner()
1914	yllcorner = geo_reference.get_yllcorner()
1915	zone = geo_reference.get_zone()
1916	else:
1917	zone = origin[0]
1918	xllcorner = origin[1]
1919	yllcorner = origin[2]
1920
1921
1922
1923	# FIXME: Refactor using code from Interpolation_function.statistics
1924	# (in interpolate.py)
1925	# Something like print swwstats(swwname)
1926	if verbose:
1927	print '------------------------------------------------'
1928	print 'Statistics of SWW file:'
1929	print ' Name: %s' %swwfile
1930	print ' Reference:'
1931	print ' Lower left corner: [%f, %f]'\
1932	%(xllcorner, yllcorner)
1933	if timestep is not None:
1934	print ' Time: %f' %(times)
1935	else:
1936	print ' Start time: %f' %fid.starttime[0]
1937	print ' Extent:'
1938	print ' x [m] in [%f, %f], len(x) == %d'\
1939	%(min(x.flat), max(x.flat), len(x.flat))
1940	print ' y [m] in [%f, %f], len(y) == %d'\
1941	%(min(y.flat), max(y.flat), len(y.flat))
1942	if timestep is not None:
1943	print ' t [s] = %f, len(t) == %d' %(times, 1)
1944	else:
1945	print ' t [s] in [%f, %f], len(t) == %d'\
1946	%(min(times), max(times), len(times))
1947	print ' Quantities [SI units]:'
1948	for name in ['stage', 'xmomentum', 'ymomentum']:
1949	q = fid.variables[name][:].flat
1950	if timestep is not None:
1951	q = q[timesteplen(x):(timestep+1)len(x)]
1952	if verbose: print ' %s in [%f, %f]' %(name, min(q), max(q))
1953	for name in ['elevation']:
1954	q = fid.variables[name][:].flat
1955	if verbose: print ' %s in [%f, %f]' %(name, min(q), max(q))
1956
1957	# Get quantity and reduce if applicable
1958	if verbose: print 'Processing quantity %s' %quantity
1959
1960	# Turn NetCDF objects into Numeric arrays
1961	quantity_dict = {}
1962	for name in fid.variables.keys():
1963	quantity_dict[name] = fid.variables[name][:]
1964
1965
1966	# Convert quantity expression to quantities found in sww file
1967	q = apply_expression_to_dictionary(quantity, quantity_dict)
1968
1969	if len(q.shape) == 2:
1970	#q has a time component and needs to be reduced along
1971	#the temporal dimension
1972	if verbose: print 'Reducing quantity %s' %quantity
1973	q_reduced = zeros( number_of_points, Float )
1974
1975	if timestep is not None:
1976	for k in range(number_of_points):
1977	q_reduced[k] = q[timestep,k]
1978	else:
1979	for k in range(number_of_points):
1980	q_reduced[k] = reduction( q[:,k] )
1981
1982	q = q_reduced
1983
1984	#Post condition: Now q has dimension: number_of_points
1985	assert len(q.shape) == 1
1986	assert q.shape[0] == number_of_points
1987
1988
1989	if verbose:
1990	print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
1991
1992
1993	#Create grid and update xll/yll corner and x,y
1994
1995	#Relative extent
1996	if easting_min is None:
1997	xmin = min(x)
1998	else:
1999	xmin = easting_min - xllcorner
2000
2001	if easting_max is None:
2002	xmax = max(x)
2003	else:
2004	xmax = easting_max - xllcorner
2005
2006	if northing_min is None:
2007	ymin = min(y)
2008	else:
2009	ymin = northing_min - yllcorner
2010
2011	if northing_max is None:
2012	ymax = max(y)
2013	else:
2014	ymax = northing_max - yllcorner
2015
2016
2017
2018	if verbose: print 'Creating grid'
2019	ncols = int((xmax-xmin)/cellsize)+1
2020	nrows = int((ymax-ymin)/cellsize)+1
2021
2022
2023	#New absolute reference and coordinates
2024	newxllcorner = xmin+xllcorner
2025	newyllcorner = ymin+yllcorner
2026
2027	x = x+xllcorner-newxllcorner
2028	y = y+yllcorner-newyllcorner
2029
2030	vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2031	assert len(vertex_points.shape) == 2
2032
2033	grid_points = zeros ( (ncols*nrows, 2), Float )
2034
2035
2036	for i in xrange(nrows):
2037	if format.lower() == 'asc':
2038	yg = i*cellsize
2039	else:
2040	#this will flip the order of the y values for ers
2041	yg = (nrows-i)*cellsize
2042
2043	for j in xrange(ncols):
2044	xg = j*cellsize
2045	k = i*ncols + j
2046
2047	grid_points[k,0] = xg
2048	grid_points[k,1] = yg
2049
2050	#Interpolate
2051	from anuga.fit_interpolate.interpolate import Interpolate
2052
2053	# Remove loners from vertex_points, volumes here
2054	vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
2055	#export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
2056	#import sys; sys.exit()
2057	interp = Interpolate(vertex_points, volumes, verbose = verbose)
2058
2059	#Interpolate using quantity values
2060	if verbose: print 'Interpolating'
2061	grid_values = interp.interpolate(q, grid_points).flat
2062
2063
2064	if verbose:
2065	print 'Interpolated values are in [%f, %f]' %(min(grid_values),
2066	max(grid_values))
2067
2068	#Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2069	P = interp.mesh.get_boundary_polygon()
2070	outside_indices = outside_polygon(grid_points, P, closed=True)
2071
2072	for i in outside_indices:
2073	grid_values[i] = NODATA_value
2074
2075
2076
2077
2078	if format.lower() == 'ers':
2079	# setup ERS header information
2080	grid_values = reshape(grid_values,(nrows, ncols))
2081	header = {}
2082	header['datum'] = '"' + datum + '"'
2083	# FIXME The use of hardwired UTM and zone number needs to be made optional
2084	# FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
2085	header['projection'] = '"UTM-' + str(zone) + '"'
2086	header['coordinatetype'] = 'EN'
2087	if header['coordinatetype'] == 'LL':
2088	header['longitude'] = str(newxllcorner)
2089	header['latitude'] = str(newyllcorner)
2090	elif header['coordinatetype'] == 'EN':
2091	header['eastings'] = str(newxllcorner)
2092	header['northings'] = str(newyllcorner)
2093	header['nullcellvalue'] = str(NODATA_value)
2094	header['xdimension'] = str(cellsize)
2095	header['ydimension'] = str(cellsize)
2096	header['value'] = '"' + quantity + '"'
2097	#header['celltype'] = 'IEEE8ByteReal' #FIXME: Breaks unit test
2098
2099
2100	#Write
2101	if verbose: print 'Writing %s' %demfile
2102	import ermapper_grids
2103	ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
2104
2105	fid.close()
2106	else:
2107	#Write to Ascii format
2108
2109	#Write prj file
2110	prjfile = basename_out + '.prj'
2111
2112	if verbose: print 'Writing %s' %prjfile
2113	prjid = open(prjfile, 'w')
2114	prjid.write('Projection %s\n' %'UTM')
2115	prjid.write('Zone %d\n' %zone)
2116	prjid.write('Datum %s\n' %datum)
2117	prjid.write('Zunits NO\n')
2118	prjid.write('Units METERS\n')
2119	prjid.write('Spheroid %s\n' %datum)
2120	prjid.write('Xshift %d\n' %false_easting)
2121	prjid.write('Yshift %d\n' %false_northing)
2122	prjid.write('Parameters\n')
2123	prjid.close()
2124
2125
2126
2127	if verbose: print 'Writing %s' %demfile
2128
2129	ascid = open(demfile, 'w')
2130
2131	ascid.write('ncols %d\n' %ncols)
2132	ascid.write('nrows %d\n' %nrows)
2133	ascid.write('xllcorner %d\n' %newxllcorner)
2134	ascid.write('yllcorner %d\n' %newyllcorner)
2135	ascid.write('cellsize %f\n' %cellsize)
2136	ascid.write('NODATA_value %d\n' %NODATA_value)
2137
2138
2139	#Get bounding polygon from mesh
2140	#P = interp.mesh.get_boundary_polygon()
2141	#inside_indices = inside_polygon(grid_points, P)
2142
2143	for i in range(nrows):
2144	if verbose and i%((nrows+10)/10)==0:
2145	print 'Doing row %d of %d' %(i, nrows)
2146
2147	base_index = (nrows-i-1)*ncols
2148
2149	slice = grid_values[base_index:base_index+ncols]
2150	s = array2string(slice, max_line_width=sys.maxint)
2151	ascid.write(s[1:-1] + '\n')
2152
2153
2154	#print
2155	#for j in range(ncols):
2156	# index = base_index+j#
2157	# print grid_values[index],
2158	# ascid.write('%f ' %grid_values[index])
2159	#ascid.write('\n')
2160
2161	#Close
2162	ascid.close()
2163	fid.close()
2164	return basename_out
2165
2166	#Backwards compatibility
2167	def sww2asc(basename_in, basename_out = None,
2168	quantity = None,
2169	timestep = None,
2170	reduction = None,
2171	cellsize = 10,
2172	verbose = False,
2173	origin = None):
2174	print 'sww2asc will soon be obsoleted - please use sww2dem'
2175	sww2dem(basename_in,
2176	basename_out = basename_out,
2177	quantity = quantity,
2178	timestep = timestep,
2179	reduction = reduction,
2180	cellsize = cellsize,
2181	verbose = verbose,
2182	origin = origin,
2183	datum = 'WGS84',
2184	format = 'asc')
2185
2186	def sww2ers(basename_in, basename_out = None,
2187	quantity = None,
2188	timestep = None,
2189	reduction = None,
2190	cellsize = 10,
2191	verbose = False,
2192	origin = None,
2193	datum = 'WGS84'):
2194	print 'sww2ers will soon be obsoleted - please use sww2dem'
2195	sww2dem(basename_in,
2196	basename_out = basename_out,
2197	quantity = quantity,
2198	timestep = timestep,
2199	reduction = reduction,
2200	cellsize = cellsize,
2201	verbose = verbose,
2202	origin = origin,
2203	datum = datum,
2204	format = 'ers')
2205	################################# END COMPATIBILITY ##############
2206
2207
2208
2209	def sww2pts(basename_in, basename_out=None,
2210	data_points=None,
2211	quantity=None,
2212	timestep=None,
2213	reduction=None,
2214	NODATA_value=-9999,
2215	verbose=False,
2216	origin=None):
2217	#datum = 'WGS84')
2218
2219
2220	"""Read SWW file and convert to interpolated values at selected points
2221
2222	The parameter quantity' must be the name of an existing quantity or
2223	an expression involving existing quantities. The default is
2224	'elevation'.
2225
2226	if timestep (an index) is given, output quantity at that timestep
2227
2228	if reduction is given use that to reduce quantity over all timesteps.
2229
2230	data_points (Nx2 array) give locations of points where quantity is to be computed.
2231
2232	"""
2233
2234	import sys
2235	from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
2236	from Numeric import array2string
2237
2238	from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
2239	from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
2240
2241	from anuga.geospatial_data.geospatial_data import Geospatial_data
2242
2243	if quantity is None:
2244	quantity = 'elevation'
2245
2246	if reduction is None:
2247	reduction = max
2248
2249	if basename_out is None:
2250	basename_out = basename_in + '_%s' %quantity
2251
2252	swwfile = basename_in + '.sww'
2253	ptsfile = basename_out + '.pts'
2254
2255	# Read sww file
2256	if verbose: print 'Reading from %s' %swwfile
2257	from Scientific.IO.NetCDF import NetCDFFile
2258	fid = NetCDFFile(swwfile)
2259
2260	# Get extent and reference
2261	x = fid.variables['x'][:]
2262	y = fid.variables['y'][:]
2263	volumes = fid.variables['volumes'][:]
2264
2265	number_of_timesteps = fid.dimensions['number_of_timesteps']
2266	number_of_points = fid.dimensions['number_of_points']
2267	if origin is None:
2268
2269	# Get geo_reference
2270	# sww files don't have to have a geo_ref
2271	try:
2272	geo_reference = Geo_reference(NetCDFObject=fid)
2273	except AttributeError, e:
2274	geo_reference = Geo_reference() #Default georef object
2275
2276	xllcorner = geo_reference.get_xllcorner()
2277	yllcorner = geo_reference.get_yllcorner()
2278	zone = geo_reference.get_zone()
2279	else:
2280	zone = origin[0]
2281	xllcorner = origin[1]
2282	yllcorner = origin[2]
2283
2284
2285
2286	# FIXME: Refactor using code from file_function.statistics
2287	# Something like print swwstats(swwname)
2288	if verbose:
2289	x = fid.variables['x'][:]
2290	y = fid.variables['y'][:]
2291	times = fid.variables['time'][:]
2292	print '------------------------------------------------'
2293	print 'Statistics of SWW file:'
2294	print ' Name: %s' %swwfile
2295	print ' Reference:'
2296	print ' Lower left corner: [%f, %f]'\
2297	%(xllcorner, yllcorner)
2298	print ' Start time: %f' %fid.starttime[0]
2299	print ' Extent:'
2300	print ' x [m] in [%f, %f], len(x) == %d'\
2301	%(min(x.flat), max(x.flat), len(x.flat))
2302	print ' y [m] in [%f, %f], len(y) == %d'\
2303	%(min(y.flat), max(y.flat), len(y.flat))
2304	print ' t [s] in [%f, %f], len(t) == %d'\
2305	%(min(times), max(times), len(times))
2306	print ' Quantities [SI units]:'
2307	for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2308	q = fid.variables[name][:].flat
2309	print ' %s in [%f, %f]' %(name, min(q), max(q))
2310
2311
2312
2313	# Get quantity and reduce if applicable
2314	if verbose: print 'Processing quantity %s' %quantity
2315
2316	# Turn NetCDF objects into Numeric arrays
2317	quantity_dict = {}
2318	for name in fid.variables.keys():
2319	quantity_dict[name] = fid.variables[name][:]
2320
2321
2322
2323	# Convert quantity expression to quantities found in sww file
2324	q = apply_expression_to_dictionary(quantity, quantity_dict)
2325
2326
2327
2328	if len(q.shape) == 2:
2329	# q has a time component and needs to be reduced along
2330	# the temporal dimension
2331	if verbose: print 'Reducing quantity %s' %quantity
2332	q_reduced = zeros( number_of_points, Float )
2333
2334	for k in range(number_of_points):
2335	q_reduced[k] = reduction( q[:,k] )
2336
2337	q = q_reduced
2338
2339	# Post condition: Now q has dimension: number_of_points
2340	assert len(q.shape) == 1
2341	assert q.shape[0] == number_of_points
2342
2343
2344	if verbose:
2345	print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
2346
2347
2348	# Create grid and update xll/yll corner and x,y
2349	vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2350	assert len(vertex_points.shape) == 2
2351
2352	# Interpolate
2353	from anuga.fit_interpolate.interpolate import Interpolate
2354	interp = Interpolate(vertex_points, volumes, verbose = verbose)
2355
2356	# Interpolate using quantity values
2357	if verbose: print 'Interpolating'
2358	interpolated_values = interp.interpolate(q, data_points).flat
2359
2360
2361	if verbose:
2362	print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
2363	max(interpolated_values))
2364
2365	# Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2366	P = interp.mesh.get_boundary_polygon()
2367	outside_indices = outside_polygon(data_points, P, closed=True)
2368
2369	for i in outside_indices:
2370	interpolated_values[i] = NODATA_value
2371
2372
2373	# Store results
2374	G = Geospatial_data(data_points=data_points,
2375	attributes=interpolated_values)
2376
2377	G.export_points_file(ptsfile, absolute = True)
2378
2379	fid.close()
2380
2381
2382	def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2383	use_cache = False,
2384	verbose = False):
2385	"""Read Digitial Elevation model from the following ASCII format (.asc)
2386
2387	Example:
2388
2389	ncols 3121
2390	nrows 1800
2391	xllcorner 722000
2392	yllcorner 5893000
2393	cellsize 25
2394	NODATA_value -9999
2395	138.3698 137.4194 136.5062 135.5558 ..........
2396
2397	Convert basename_in + '.asc' to NetCDF format (.dem)
2398	mimicking the ASCII format closely.
2399
2400
2401	An accompanying file with same basename_in but extension .prj must exist
2402	and is used to fix the UTM zone, datum, false northings and eastings.
2403
2404	The prj format is assumed to be as
2405
2406	Projection UTM
2407	Zone 56
2408	Datum WGS84
2409	Zunits NO
2410	Units METERS
2411	Spheroid WGS84
2412	Xshift 0.0000000000
2413	Yshift 10000000.0000000000
2414	Parameters
2415	"""
2416
2417
2418
2419	kwargs = {'basename_out': basename_out, 'verbose': verbose}
2420
2421	if use_cache is True:
2422	from caching import cache
2423	result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
2424	dependencies = [basename_in + '.asc',
2425	basename_in + '.prj'],
2426	verbose = verbose)
2427
2428	else:
2429	result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
2430
2431	return result
2432
2433
2434
2435
2436
2437
2438	def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2439	verbose = False):
2440	"""Read Digitial Elevation model from the following ASCII format (.asc)
2441
2442	Internal function. See public function convert_dem_from_ascii2netcdf for details.
2443	"""
2444
2445	import os
2446	from Scientific.IO.NetCDF import NetCDFFile
2447	from Numeric import Float, array
2448
2449	#root, ext = os.path.splitext(basename_in)
2450	root = basename_in
2451
2452	###########################################
2453	# Read Meta data
2454	if verbose: print 'Reading METADATA from %s' %root + '.prj'
2455	metadatafile = open(root + '.prj')
2456	metalines = metadatafile.readlines()
2457	metadatafile.close()
2458
2459	L = metalines[0].strip().split()
2460	assert L[0].strip().lower() == 'projection'
2461	projection = L[1].strip() #TEXT
2462
2463	L = metalines[1].strip().split()
2464	assert L[0].strip().lower() == 'zone'
2465	zone = int(L[1].strip())
2466
2467	L = metalines[2].strip().split()
2468	assert L[0].strip().lower() == 'datum'
2469	datum = L[1].strip() #TEXT
2470
2471	L = metalines[3].strip().split()
2472	assert L[0].strip().lower() == 'zunits' #IGNORE
2473	zunits = L[1].strip() #TEXT
2474
2475	L = metalines[4].strip().split()
2476	assert L[0].strip().lower() == 'units'
2477	units = L[1].strip() #TEXT
2478
2479	L = metalines[5].strip().split()
2480	assert L[0].strip().lower() == 'spheroid' #IGNORE
2481	spheroid = L[1].strip() #TEXT
2482
2483	L = metalines[6].strip().split()
2484	assert L[0].strip().lower() == 'xshift'
2485	false_easting = float(L[1].strip())
2486
2487	L = metalines[7].strip().split()
2488	assert L[0].strip().lower() == 'yshift'
2489	false_northing = float(L[1].strip())
2490
2491	#print false_easting, false_northing, zone, datum
2492
2493
2494	###########################################
2495	#Read DEM data
2496
2497	datafile = open(basename_in + '.asc')
2498
2499	if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
2500	lines = datafile.readlines()
2501	datafile.close()
2502
2503	if verbose: print 'Got', len(lines), ' lines'
2504
2505	ncols = int(lines[0].split()[1].strip())
2506	nrows = int(lines[1].split()[1].strip())
2507	xllcorner = float(lines[2].split()[1].strip())
2508	yllcorner = float(lines[3].split()[1].strip())
2509	cellsize = float(lines[4].split()[1].strip())
2510	NODATA_value = int(lines[5].split()[1].strip())
2511
2512	assert len(lines) == nrows + 6
2513
2514
2515	##########################################
2516
2517
2518	if basename_out == None:
2519	netcdfname = root + '.dem'
2520	else:
2521	netcdfname = basename_out + '.dem'
2522
2523	if verbose: print 'Store to NetCDF file %s' %netcdfname
2524	# NetCDF file definition
2525	fid = NetCDFFile(netcdfname, 'w')
2526
2527	#Create new file
2528	fid.institution = 'Geoscience Australia'
2529	fid.description = 'NetCDF DEM format for compact and portable storage ' +\
2530	'of spatial point data'
2531
2532	fid.ncols = ncols
2533	fid.nrows = nrows
2534	fid.xllcorner = xllcorner
2535	fid.yllcorner = yllcorner
2536	fid.cellsize = cellsize
2537	fid.NODATA_value = NODATA_value
2538
2539	fid.zone = zone
2540	fid.false_easting = false_easting
2541	fid.false_northing = false_northing
2542	fid.projection = projection
2543	fid.datum = datum
2544	fid.units = units
2545
2546
2547	# dimension definitions
2548	fid.createDimension('number_of_rows', nrows)
2549	fid.createDimension('number_of_columns', ncols)
2550
2551	# variable definitions
2552	fid.createVariable('elevation', Float, ('number_of_rows',
2553	'number_of_columns'))
2554
2555	# Get handles to the variables
2556	elevation = fid.variables['elevation']
2557
2558	#Store data
2559	n = len(lines[6:])
2560	for i, line in enumerate(lines[6:]):
2561	fields = line.split()
2562	if verbose and i%((n+10)/10)==0:
2563	print 'Processing row %d of %d' %(i, nrows)
2564
2565	elevation[i, :] = array([float(x) for x in fields])
2566
2567	fid.close()
2568
2569
2570
2571
2572
2573	def ferret2sww(basename_in, basename_out = None,
2574	verbose = False,
2575	minlat = None, maxlat = None,
2576	minlon = None, maxlon = None,
2577	mint = None, maxt = None, mean_stage = 0,
2578	origin = None, zscale = 1,
2579	fail_on_NaN = True,
2580	NaN_filler = 0,
2581	elevation = None,
2582	inverted_bathymetry = True
2583	): #FIXME: Bathymetry should be obtained
2584	#from MOST somehow.
2585	#Alternatively from elsewhere
2586	#or, as a last resort,
2587	#specified here.
2588	#The value of -100 will work
2589	#for the Wollongong tsunami
2590	#scenario but is very hacky
2591	"""Convert MOST and 'Ferret' NetCDF format for wave propagation to
2592	sww format native to abstract_2d_finite_volumes.
2593
2594	Specify only basename_in and read files of the form
2595	basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
2596	relative height, x-velocity and y-velocity, respectively.
2597
2598	Also convert latitude and longitude to UTM. All coordinates are
2599	assumed to be given in the GDA94 datum.
2600
2601	min's and max's: If omitted - full extend is used.
2602	To include a value min may equal it, while max must exceed it.
2603	Lat and lon are assuemd to be in decimal degrees
2604
2605	origin is a 3-tuple with geo referenced
2606	UTM coordinates (zone, easting, northing)
2607
2608	nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
2609	which means that longitude is the fastest
2610	varying dimension (row major order, so to speak)
2611
2612	ferret2sww uses grid points as vertices in a triangular grid
2613	counting vertices from lower left corner upwards, then right
2614	"""
2615
2616	import os
2617	from Scientific.IO.NetCDF import NetCDFFile
2618	from Numeric import Float, Int, Int32, searchsorted, zeros, array
2619	from Numeric import allclose, around
2620
2621	precision = Float
2622
2623	msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
2624
2625	if minlat != None:
2626	assert -90 < minlat < 90 , msg
2627	if maxlat != None:
2628	assert -90 < maxlat < 90 , msg
2629	if minlat != None:
2630	assert maxlat > minlat
2631	if minlon != None:
2632	assert -180 < minlon < 180 , msg
2633	if maxlon != None:
2634	assert -180 < maxlon < 180 , msg
2635	if minlon != None:
2636	assert maxlon > minlon
2637
2638
2639
2640	#Get NetCDF data
2641	if verbose: print 'Reading files %s_*.nc' %basename_in
2642	#print "basename_in + '_ha.nc'",basename_in + '_ha.nc'
2643	file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
2644	file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
2645	file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
2646	file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)
2647
2648	if basename_out is None:
2649	swwname = basename_in + '.sww'
2650	else:
2651	swwname = basename_out + '.sww'
2652
2653	# Get dimensions of file_h
2654	for dimension in file_h.dimensions.keys():
2655	if dimension[:3] == 'LON':
2656	dim_h_longitude = dimension
2657	if dimension[:3] == 'LAT':
2658	dim_h_latitude = dimension
2659	if dimension[:4] == 'TIME':
2660	dim_h_time = dimension
2661
2662	# print 'long:', dim_h_longitude
2663	# print 'lats:', dim_h_latitude
2664	# print 'times:', dim_h_time
2665
2666	times = file_h.variables[dim_h_time]
2667	latitudes = file_h.variables[dim_h_latitude]
2668	longitudes = file_h.variables[dim_h_longitude]
2669
2670	# get dimensions for file_e
2671	for dimension in file_e.dimensions.keys():
2672	if dimension[:3] == 'LON':
2673	dim_e_longitude = dimension
2674	if dimension[:3] == 'LAT':
2675	dim_e_latitude = dimension
2676
2677	# get dimensions for file_u
2678	for dimension in file_u.dimensions.keys():
2679	if dimension[:3] == 'LON':
2680	dim_u_longitude = dimension
2681	if dimension[:3] == 'LAT':
2682	dim_u_latitude = dimension
2683	if dimension[:4] == 'TIME':
2684	dim_u_time = dimension
2685
2686	# get dimensions for file_v
2687	for dimension in file_v.dimensions.keys():
2688	if dimension[:3] == 'LON':
2689	dim_v_longitude = dimension
2690	if dimension[:3] == 'LAT':
2691	dim_v_latitude = dimension
2692	if dimension[:4] == 'TIME':
2693	dim_v_time = dimension
2694
2695
2696	# Precision used by most for lat/lon is 4 or 5 decimals
2697	e_lat = around(file_e.variables[dim_e_latitude][:], 5)
2698	e_lon = around(file_e.variables[dim_e_longitude][:], 5)
2699
2700	# Check that files are compatible
2701	assert allclose(latitudes, file_u.variables[dim_u_latitude])
2702	assert allclose(latitudes, file_v.variables[dim_v_latitude])
2703	assert allclose(latitudes, e_lat)
2704
2705	assert allclose(longitudes, file_u.variables[dim_u_longitude])
2706	assert allclose(longitudes, file_v.variables[dim_v_longitude])
2707	assert allclose(longitudes, e_lon)
2708
2709	if mint is None:
2710	jmin = 0
2711	mint = times[0]
2712	else:
2713	jmin = searchsorted(times, mint)
2714
2715	if maxt is None:
2716	jmax = len(times)
2717	maxt = times[-1]
2718	else:
2719	jmax = searchsorted(times, maxt)
2720
2721	#print "latitudes[:]",latitudes[:]
2722	#print "longitudes[:]",longitudes [:]
2723	kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
2724	longitudes[:],
2725	minlat, maxlat,
2726	minlon, maxlon)
2727
2728
2729	times = times[jmin:jmax]
2730	latitudes = latitudes[kmin:kmax]
2731	longitudes = longitudes[lmin:lmax]
2732
2733	#print "latitudes[:]",latitudes[:]
2734	#print "longitudes[:]",longitudes [:]
2735
2736	if verbose: print 'cropping'
2737	zname = 'ELEVATION'
2738
2739	amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
2740	uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
2741	vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
2742	elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
2743
2744	# if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
2745	# elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
2746	# elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
2747	# from Numeric import asarray
2748	# elevations=elevations.tolist()
2749	# elevations.reverse()
2750	# elevations=asarray(elevations)
2751	# else:
2752	# from Numeric import asarray
2753	# elevations=elevations.tolist()
2754	# elevations.reverse()
2755	# elevations=asarray(elevations)
2756	# 'print hmmm'
2757
2758
2759
2760	#Get missing values
2761	nan_ha = file_h.variables['HA'].missing_value[0]
2762	nan_ua = file_u.variables['UA'].missing_value[0]
2763	nan_va = file_v.variables['VA'].missing_value[0]
2764	if hasattr(file_e.variables[zname],'missing_value'):
2765	nan_e = file_e.variables[zname].missing_value[0]
2766	else:
2767	nan_e = None
2768
2769	#Cleanup
2770	from Numeric import sometrue
2771
2772	missing = (amplitudes == nan_ha)
2773	if sometrue (missing):
2774	if fail_on_NaN:
2775	msg = 'NetCDFFile %s contains missing values'\
2776	%(basename_in+'_ha.nc')
2777	raise DataMissingValuesError, msg
2778	else:
2779	amplitudes = amplitudes(missing==0) + missingNaN_filler
2780
2781	missing = (uspeed == nan_ua)
2782	if sometrue (missing):
2783	if fail_on_NaN:
2784	msg = 'NetCDFFile %s contains missing values'\
2785	%(basename_in+'_ua.nc')
2786	raise DataMissingValuesError, msg
2787	else:
2788	uspeed = uspeed(missing==0) + missingNaN_filler
2789
2790	missing = (vspeed == nan_va)
2791	if sometrue (missing):
2792	if fail_on_NaN:
2793	msg = 'NetCDFFile %s contains missing values'\
2794	%(basename_in+'_va.nc')
2795	raise DataMissingValuesError, msg
2796	else:
2797	vspeed = vspeed(missing==0) + missingNaN_filler
2798
2799
2800	missing = (elevations == nan_e)
2801	if sometrue (missing):
2802	if fail_on_NaN:
2803	msg = 'NetCDFFile %s contains missing values'\
2804	%(basename_in+'_e.nc')
2805	raise DataMissingValuesError, msg
2806	else:
2807	elevations = elevations(missing==0) + missingNaN_filler
2808
2809	#######
2810
2811
2812
2813	number_of_times = times.shape[0]
2814	number_of_latitudes = latitudes.shape[0]
2815	number_of_longitudes = longitudes.shape[0]
2816
2817	assert amplitudes.shape[0] == number_of_times
2818	assert amplitudes.shape[1] == number_of_latitudes
2819	assert amplitudes.shape[2] == number_of_longitudes
2820
2821	if verbose:
2822	print '------------------------------------------------'
2823	print 'Statistics:'
2824	print ' Extent (lat/lon):'
2825	print ' lat in [%f, %f], len(lat) == %d'\
2826	%(min(latitudes.flat), max(latitudes.flat),
2827	len(latitudes.flat))
2828	print ' lon in [%f, %f], len(lon) == %d'\
2829	%(min(longitudes.flat), max(longitudes.flat),
2830	len(longitudes.flat))
2831	print ' t in [%f, %f], len(t) == %d'\
2832	%(min(times.flat), max(times.flat), len(times.flat))
2833
2834	q = amplitudes.flat
2835	name = 'Amplitudes (ha) [cm]'
2836	print ' %s in [%f, %f]' %(name, min(q), max(q))
2837
2838	q = uspeed.flat
2839	name = 'Speeds (ua) [cm/s]'
2840	print ' %s in [%f, %f]' %(name, min(q), max(q))
2841
2842	q = vspeed.flat
2843	name = 'Speeds (va) [cm/s]'
2844	print ' %s in [%f, %f]' %(name, min(q), max(q))
2845
2846	q = elevations.flat
2847	name = 'Elevations (e) [m]'
2848	print ' %s in [%f, %f]' %(name, min(q), max(q))
2849
2850
2851	# print number_of_latitudes, number_of_longitudes
2852	number_of_points = number_of_latitudes*number_of_longitudes
2853	number_of_volumes = (number_of_latitudes-1)(number_of_longitudes-1)2
2854
2855
2856	file_h.close()
2857	file_u.close()
2858	file_v.close()
2859	file_e.close()
2860
2861
2862	# NetCDF file definition
2863	outfile = NetCDFFile(swwname, 'w')
2864
2865	description = 'Converted from Ferret files: %s, %s, %s, %s'\
2866	%(basename_in + '_ha.nc',
2867	basename_in + '_ua.nc',
2868	basename_in + '_va.nc',
2869	basename_in + '_e.nc')
2870
2871	# Create new file
2872	starttime = times[0]
2873	sww = Write_sww()
2874	sww.store_header(outfile, times, number_of_volumes,
2875	number_of_points, description=description,
2876	verbose=verbose)
2877
2878	# Store
2879	from anuga.coordinate_transforms.redfearn import redfearn
2880	x = zeros(number_of_points, Float) #Easting
2881	y = zeros(number_of_points, Float) #Northing
2882
2883
2884	if verbose: print 'Making triangular grid'
2885
2886	# Check zone boundaries
2887	refzone, _, _ = redfearn(latitudes[0],longitudes[0])
2888
2889	vertices = {}
2890	i = 0
2891	for k, lat in enumerate(latitudes): #Y direction
2892	for l, lon in enumerate(longitudes): #X direction
2893
2894	vertices[l,k] = i
2895
2896	zone, easting, northing = redfearn(lat,lon)
2897
2898	msg = 'Zone boundary crossed at longitude =', lon
2899	#assert zone == refzone, msg
2900	#print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
2901	x[i] = easting
2902	y[i] = northing
2903	i += 1
2904
2905	#Construct 2 triangles per 'rectangular' element
2906	volumes = []
2907	for l in range(number_of_longitudes-1): #X direction
2908	for k in range(number_of_latitudes-1): #Y direction
2909	v1 = vertices[l,k+1]
2910	v2 = vertices[l,k]
2911	v3 = vertices[l+1,k+1]
2912	v4 = vertices[l+1,k]
2913
2914	volumes.append([v1,v2,v3]) #Upper element
2915	volumes.append([v4,v3,v2]) #Lower element
2916
2917	volumes = array(volumes)
2918
2919	if origin is None:
2920	origin = Geo_reference(refzone,min(x),min(y))
2921	geo_ref = write_NetCDF_georeference(origin, outfile)
2922
2923	if elevation is not None:
2924	z = elevation
2925	else:
2926	if inverted_bathymetry:
2927	z = -1*elevations
2928	else:
2929	z = elevations
2930	#FIXME: z should be obtained from MOST and passed in here
2931
2932	from Numeric import resize
2933	z = resize(z,outfile.variables['z'][:].shape)
2934	outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
2935	outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
2936	outfile.variables['z'][:] = z #FIXME HACK for bacwards compat.
2937	outfile.variables['elevation'][:] = z
2938	outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
2939
2940
2941
2942	#Time stepping
2943	stage = outfile.variables['stage']
2944	xmomentum = outfile.variables['xmomentum']
2945	ymomentum = outfile.variables['ymomentum']
2946
2947	if verbose: print 'Converting quantities'
2948	n = len(times)
2949	for j in range(n):
2950	if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n)
2951	i = 0
2952	for k in range(number_of_latitudes): #Y direction
2953	for l in range(number_of_longitudes): #X direction
2954	w = zscale*amplitudes[j,k,l]/100 + mean_stage
2955	stage[j,i] = w
2956	h = w - z[i]
2957	xmomentum[j,i] = uspeed[j,k,l]/100*h
2958	ymomentum[j,i] = vspeed[j,k,l]/100*h
2959	i += 1
2960
2961	#outfile.close()
2962
2963	#FIXME: Refactor using code from file_function.statistics
2964	#Something like print swwstats(swwname)
2965	if verbose:
2966	x = outfile.variables['x'][:]
2967	y = outfile.variables['y'][:]
2968	print '------------------------------------------------'
2969	print 'Statistics of output file:'
2970	print ' Name: %s' %swwname
2971	print ' Reference:'
2972	print ' Lower left corner: [%f, %f]'\
2973	%(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
2974	print ' Start time: %f' %starttime
2975	print ' Min time: %f' %mint
2976	print ' Max time: %f' %maxt
2977	print ' Extent:'
2978	print ' x [m] in [%f, %f], len(x) == %d'\
2979	%(min(x.flat), max(x.flat), len(x.flat))
2980	print ' y [m] in [%f, %f], len(y) == %d'\
2981	%(min(y.flat), max(y.flat), len(y.flat))
2982	print ' t [s] in [%f, %f], len(t) == %d'\
2983	%(min(times), max(times), len(times))
2984	print ' Quantities [SI units]:'
2985	for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2986	q = outfile.variables[name][:].flat
2987	print ' %s in [%f, %f]' %(name, min(q), max(q))
2988
2989
2990
2991	outfile.close()
2992
2993
2994
2995
2996
2997	def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
2998	"""Template for converting typical text files with time series to
2999	NetCDF tms file.
3000
3001
3002	The file format is assumed to be either two fields separated by a comma:
3003
3004	time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
3005
3006	E.g
3007
3008	31/08/04 00:00:00, 1.328223 0 0
3009	31/08/04 00:15:00, 1.292912 0 0
3010
3011	or time (seconds), value0 value1 value2 ...
3012
3013	0.0, 1.328223 0 0
3014	0.1, 1.292912 0 0
3015
3016	will provide a time dependent function f(t) with three attributes
3017
3018	filename is assumed to be the rootname with extenisons .txt and .sww
3019	"""
3020
3021	import time, calendar
3022	from Numeric import array
3023	from anuga.config import time_format
3024	from anuga.utilities.numerical_tools import ensure_numeric
3025
3026
3027
3028	fid = open(filename + '.txt')
3029	line = fid.readline()
3030	fid.close()
3031
3032	fields = line.split(',')
3033	msg = 'File %s must have the format date, value0 value1 value2 ...'
3034	assert len(fields) == 2, msg
3035
3036	if not time_as_seconds:
3037	try:
3038	starttime = calendar.timegm(time.strptime(fields[0], time_format))
3039	except ValueError:
3040	msg = 'First field in file %s must be' %filename
3041	msg += ' date-time with format %s.\n' %time_format
3042	msg += 'I got %s instead.' %fields[0]
3043	raise DataTimeError, msg
3044	else:
3045	try:
3046	starttime = float(fields[0])
3047	except Error:
3048	msg = "Bad time format"
3049	raise DataTimeError, msg
3050
3051
3052	#Split values
3053	values = []
3054	for value in fields[1].split():
3055	values.append(float(value))
3056
3057	q = ensure_numeric(values)
3058
3059	msg = 'ERROR: File must contain at least one independent value'
3060	assert len(q.shape) == 1, msg
3061
3062
3063
3064	#Read times proper
3065	from Numeric import zeros, Float, alltrue
3066	from anuga.config import time_format
3067	import time, calendar
3068
3069	fid = open(filename + '.txt')
3070	lines = fid.readlines()
3071	fid.close()
3072
3073	N = len(lines)
3074	d = len(q)
3075
3076	T = zeros(N, Float) #Time
3077	Q = zeros((N, d), Float) #Values
3078
3079	for i, line in enumerate(lines):
3080	fields = line.split(',')
3081	if not time_as_seconds:
3082	realtime = calendar.timegm(time.strptime(fields[0], time_format))
3083	else:
3084	realtime = float(fields[0])
3085	T[i] = realtime - starttime
3086
3087	for j, value in enumerate(fields[1].split()):
3088	Q[i, j] = float(value)
3089
3090	msg = 'File %s must list time as a monotonuosly ' %filename
3091	msg += 'increasing sequence'
3092	assert alltrue( T[1:] - T[:-1] > 0 ), msg
3093
3094	#Create NetCDF file
3095	from Scientific.IO.NetCDF import NetCDFFile
3096
3097	fid = NetCDFFile(filename + '.tms', 'w')
3098
3099
3100	fid.institution = 'Geoscience Australia'
3101	fid.description = 'Time series'
3102
3103
3104	#Reference point
3105	#Start time in seconds since the epoch (midnight 1/1/1970)
3106	#FIXME: Use Georef
3107	fid.starttime = starttime
3108
3109	# dimension definitions
3110	#fid.createDimension('number_of_volumes', self.number_of_volumes)
3111	#fid.createDimension('number_of_vertices', 3)
3112
3113
3114	fid.createDimension('number_of_timesteps', len(T))
3115
3116	fid.createVariable('time', Float, ('number_of_timesteps',))
3117
3118	fid.variables['time'][:] = T
3119
3120	for i in range(Q.shape[1]):
3121	try:
3122	name = quantity_names[i]
3123	except:
3124	name = 'Attribute%d'%i
3125
3126	fid.createVariable(name, Float, ('number_of_timesteps',))
3127	fid.variables[name][:] = Q[:,i]
3128
3129	fid.close()
3130
3131
3132	def extent_sww(file_name):
3133	"""
3134	Read in an sww file.
3135
3136	Input;
3137	file_name - the sww file
3138
3139	Output;
3140	z - Vector of bed elevation
3141	volumes - Array. Each row has 3 values, representing
3142	the vertices that define the volume
3143	time - Vector of the times where there is stage information
3144	stage - array with respect to time and vertices (x,y)
3145	"""
3146
3147
3148	from Scientific.IO.NetCDF import NetCDFFile
3149
3150	#Check contents
3151	#Get NetCDF
3152	fid = NetCDFFile(file_name, 'r')
3153
3154	# Get the variables
3155	x = fid.variables['x'][:]
3156	y = fid.variables['y'][:]
3157	stage = fid.variables['stage'][:]
3158	#print "stage",stage
3159	#print "stage.shap",stage.shape
3160	#print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
3161	#print "min(stage)",min(stage)
3162
3163	fid.close()
3164
3165	return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
3166
3167
3168	def sww2domain(filename,boundary=None,t=None,\
3169	fail_if_NaN=True,NaN_filler=0\
3170	,verbose = False,very_verbose = False):
3171	"""
3172	Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
3173
3174	Boundary is not recommended if domain.smooth is not selected, as it
3175	uses unique coordinates, but not unique boundaries. This means that
3176	the boundary file will not be compatable with the coordinates, and will
3177	give a different final boundary, or crash.
3178	"""
3179	NaN=9.969209968386869e+036
3180	#initialise NaN.
3181
3182	from Scientific.IO.NetCDF import NetCDFFile
3183	from shallow_water import Domain
3184	from Numeric import asarray, transpose, resize
3185
3186	if verbose: print 'Reading from ', filename
3187	fid = NetCDFFile(filename, 'r') #Open existing file for read
3188	time = fid.variables['time'] #Timesteps
3189	if t is None:
3190	t = time[-1]
3191	time_interp = get_time_interp(time,t)
3192
3193	# Get the variables as Numeric arrays
3194	x = fid.variables['x'][:] #x-coordinates of vertices
3195	y = fid.variables['y'][:] #y-coordinates of vertices
3196	elevation = fid.variables['elevation'] #Elevation
3197	stage = fid.variables['stage'] #Water level
3198	xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction
3199	ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction
3200
3201	starttime = fid.starttime[0]
3202	volumes = fid.variables['volumes'][:] #Connectivity
3203	coordinates=transpose(asarray([x.tolist(),y.tolist()]))
3204
3205	conserved_quantities = []
3206	interpolated_quantities = {}
3207	other_quantities = []
3208
3209	# get geo_reference
3210	#sww files don't have to have a geo_ref
3211	try:
3212	geo_reference = Geo_reference(NetCDFObject=fid)
3213	except: #AttributeError, e:
3214	geo_reference = None
3215
3216	if verbose: print ' getting quantities'
3217	for quantity in fid.variables.keys():
3218	dimensions = fid.variables[quantity].dimensions
3219	if 'number_of_timesteps' in dimensions:
3220	conserved_quantities.append(quantity)
3221	interpolated_quantities[quantity]=\
3222	interpolated_quantity(fid.variables[quantity][:],time_interp)
3223	else: other_quantities.append(quantity)
3224
3225	other_quantities.remove('x')
3226	other_quantities.remove('y')
3227	other_quantities.remove('z')
3228	other_quantities.remove('volumes')
3229	try:
3230	other_quantities.remove('stage_range')
3231	other_quantities.remove('xmomentum_range')
3232	other_quantities.remove('ymomentum_range')
3233	other_quantities.remove('elevation_range')
3234	except:
3235	pass
3236
3237
3238	conserved_quantities.remove('time')
3239
3240	if verbose: print ' building domain'
3241	# From domain.Domain:
3242	# domain = Domain(coordinates, volumes,\
3243	# conserved_quantities = conserved_quantities,\
3244	# other_quantities = other_quantities,zone=zone,\
3245	# xllcorner=xllcorner, yllcorner=yllcorner)
3246
3247	# From shallow_water.Domain:
3248	coordinates=coordinates.tolist()
3249	volumes=volumes.tolist()
3250	#FIXME:should this be in mesh?(peter row)
3251	if fid.smoothing == 'Yes': unique = False
3252	else: unique = True
3253	if unique:
3254	coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
3255
3256
3257	try:
3258	domain = Domain(coordinates, volumes, boundary)
3259	except AssertionError, e:
3260	fid.close()
3261	msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
3262	raise DataDomainError, msg
3263
3264	if not boundary is None:
3265	domain.boundary = boundary
3266
3267	domain.geo_reference = geo_reference
3268
3269	domain.starttime=float(starttime)+float(t)
3270	domain.time=0.0
3271
3272	for quantity in other_quantities:
3273	try:
3274	NaN = fid.variables[quantity].missing_value
3275	except:
3276	pass #quantity has no missing_value number
3277	X = fid.variables[quantity][:]
3278	if very_verbose:
3279	print ' ',quantity
3280	print ' NaN =',NaN
3281	print ' max(X)'
3282	print ' ',max(X)
3283	print ' max(X)==NaN'
3284	print ' ',max(X)==NaN
3285	print ''
3286	if (max(X)==NaN) or (min(X)==NaN):
3287	if fail_if_NaN:
3288	msg = 'quantity "%s" contains no_data entry'%quantity
3289	raise DataMissingValuesError, msg
3290	else:
3291	data = (X<>NaN)
3292	X = (Xdata)+(data==0)NaN_filler
3293	if unique:
3294	X = resize(X,(len(X)/3,3))
3295	domain.set_quantity(quantity,X)
3296	#
3297	for quantity in conserved_quantities:
3298	try:
3299	NaN = fid.variables[quantity].missing_value
3300	except:
3301	pass #quantity has no missing_value number
3302	X = interpolated_quantities[quantity]
3303	if very_verbose:
3304	print ' ',quantity
3305	print ' NaN =',NaN
3306	print ' max(X)'
3307	print ' ',max(X)
3308	print ' max(X)==NaN'
3309	print ' ',max(X)==NaN
3310	print ''
3311	if (max(X)==NaN) or (min(X)==NaN):
3312	if fail_if_NaN:
3313	msg = 'quantity "%s" contains no_data entry'%quantity
3314	raise DataMissingValuesError, msg
3315	else:
3316	data = (X<>NaN)
3317	X = (Xdata)+(data==0)NaN_filler
3318	if unique:
3319	X = resize(X,(X.shape[0]/3,3))
3320	domain.set_quantity(quantity,X)
3321
3322	fid.close()
3323	return domain
3324
3325	def interpolated_quantity(saved_quantity,time_interp):
3326
3327	#given an index and ratio, interpolate quantity with respect to time.
3328	index,ratio = time_interp
3329	Q = saved_quantity
3330	if ratio > 0:
3331	q = (1-ratio)Q[index]+ ratioQ[index+1]
3332	else:
3333	q = Q[index]
3334	#Return vector of interpolated values
3335	return q
3336
3337	def get_time_interp(time,t=None):
3338	#Finds the ratio and index for time interpolation.
3339	#It is borrowed from previous abstract_2d_finite_volumes code.
3340	if t is None:
3341	t=time[-1]
3342	index = -1
3343	ratio = 0.
3344	else:
3345	T = time
3346	tau = t
3347	index=0
3348	msg = 'Time interval derived from file %s [%s:%s]'\
3349	%('FIXMEfilename', T[0], T[-1])
3350	msg += ' does not match model time: %s' %tau
3351	if tau < time[0]: raise DataTimeError, msg
3352	if tau > time[-1]: raise DataTimeError, msg
3353	while tau > time[index]: index += 1
3354	while tau < time[index]: index -= 1
3355	if tau == time[index]:
3356	#Protect against case where tau == time[-1] (last time)
3357	# - also works in general when tau == time[i]
3358	ratio = 0
3359	else:
3360	#t is now between index and index+1
3361	ratio = (tau - time[index])/(time[index+1] - time[index])
3362	return (index,ratio)
3363
3364
3365	def weed(coordinates,volumes,boundary = None):
3366	if type(coordinates)==ArrayType:
3367	coordinates = coordinates.tolist()
3368	if type(volumes)==ArrayType:
3369	volumes = volumes.tolist()
3370
3371	unique = False
3372	point_dict = {}
3373	same_point = {}
3374	for i in range(len(coordinates)):
3375	point = tuple(coordinates[i])
3376	if point_dict.has_key(point):
3377	unique = True
3378	same_point[i]=point
3379	#to change all point i references to point j
3380	else:
3381	point_dict[point]=i
3382	same_point[i]=point
3383
3384	coordinates = []
3385	i = 0
3386	for point in point_dict.keys():
3387	point = tuple(point)
3388	coordinates.append(list(point))
3389	point_dict[point]=i
3390	i+=1
3391
3392
3393	for volume in volumes:
3394	for i in range(len(volume)):
3395	index = volume[i]
3396	if index>-1:
3397	volume[i]=point_dict[same_point[index]]
3398
3399	new_boundary = {}
3400	if not boundary is None:
3401	for segment in boundary.keys():
3402	point0 = point_dict[same_point[segment[0]]]
3403	point1 = point_dict[same_point[segment[1]]]
3404	label = boundary[segment]
3405	#FIXME should the bounday attributes be concaterated
3406	#('exterior, pond') or replaced ('pond')(peter row)
3407
3408	if new_boundary.has_key((point0,point1)):
3409	new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
3410	#+','+label
3411
3412	elif new_boundary.has_key((point1,point0)):
3413	new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
3414	#+','+label
3415	else: new_boundary[(point0,point1)]=label
3416
3417	boundary = new_boundary
3418
3419	return coordinates,volumes,boundary
3420
3421
3422	def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
3423	verbose=False):
3424	"""Read Digitial Elevation model from the following NetCDF format (.dem)
3425
3426	Example:
3427
3428	ncols 3121
3429	nrows 1800
3430	xllcorner 722000
3431	yllcorner 5893000
3432	cellsize 25
3433	NODATA_value -9999
3434	138.3698 137.4194 136.5062 135.5558 ..........
3435
3436	Decimate data to cellsize_new using stencil and write to NetCDF dem format.
3437	"""
3438
3439	import os
3440	from Scientific.IO.NetCDF import NetCDFFile
3441	from Numeric import Float, zeros, sum, reshape, equal
3442
3443	root = basename_in
3444	inname = root + '.dem'
3445
3446	#Open existing netcdf file to read
3447	infile = NetCDFFile(inname, 'r')
3448	if verbose: print 'Reading DEM from %s' %inname
3449
3450	#Read metadata
3451	ncols = infile.ncols[0]
3452	nrows = infile.nrows[0]
3453	xllcorner = infile.xllcorner[0]
3454	yllcorner = infile.yllcorner[0]
3455	cellsize = infile.cellsize[0]
3456	NODATA_value = infile.NODATA_value[0]
3457	zone = infile.zone[0]
3458	false_easting = infile.false_easting[0]
3459	false_northing = infile.false_northing[0]
3460	projection = infile.projection
3461	datum = infile.datum
3462	units = infile.units
3463
3464	dem_elevation = infile.variables['elevation']
3465
3466	#Get output file name
3467	if basename_out == None:
3468	outname = root + '_' + repr(cellsize_new) + '.dem'
3469	else:
3470	outname = basename_out + '.dem'
3471
3472	if verbose: print 'Write decimated NetCDF file to %s' %outname
3473
3474	#Determine some dimensions for decimated grid
3475	(nrows_stencil, ncols_stencil) = stencil.shape
3476	x_offset = ncols_stencil / 2
3477	y_offset = nrows_stencil / 2
3478	cellsize_ratio = int(cellsize_new / cellsize)
3479	ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
3480	nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
3481
3482	#Open netcdf file for output
3483	outfile = NetCDFFile(outname, 'w')
3484
3485	#Create new file
3486	outfile.institution = 'Geoscience Australia'
3487	outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
3488	'of spatial point data'
3489	#Georeferencing
3490	outfile.zone = zone
3491	outfile.projection = projection
3492	outfile.datum = datum
3493	outfile.units = units
3494
3495	outfile.cellsize = cellsize_new
3496	outfile.NODATA_value = NODATA_value
3497	outfile.false_easting = false_easting
3498	outfile.false_northing = false_northing
3499
3500	outfile.xllcorner = xllcorner + (x_offset * cellsize)
3501	outfile.yllcorner = yllcorner + (y_offset * cellsize)
3502	outfile.ncols = ncols_new
3503	outfile.nrows = nrows_new
3504
3505	# dimension definition
3506	outfile.createDimension('number_of_points', nrows_new*ncols_new)
3507
3508	# variable definition
3509	outfile.createVariable('elevation', Float, ('number_of_points',))
3510
3511	# Get handle to the variable
3512	elevation = outfile.variables['elevation']
3513
3514	dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
3515
3516	#Store data
3517	global_index = 0
3518	for i in range(nrows_new):
3519	if verbose: print 'Processing row %d of %d' %(i, nrows_new)
3520	lower_index = global_index
3521	telev = zeros(ncols_new, Float)
3522	local_index = 0
3523	trow = i * cellsize_ratio
3524
3525	for j in range(ncols_new):
3526	tcol = j * cellsize_ratio
3527	tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
3528
3529	#if dem contains 1 or more NODATA_values set value in
3530	#decimated dem to NODATA_value, else compute decimated
3531	#value using stencil
3532	if sum(sum(equal(tmp, NODATA_value))) > 0:
3533	telev[local_index] = NODATA_value
3534	else:
3535	telev[local_index] = sum(sum(tmp * stencil))
3536
3537	global_index += 1
3538	local_index += 1
3539
3540	upper_index = global_index
3541
3542	elevation[lower_index:upper_index] = telev
3543
3544	assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
3545
3546	infile.close()
3547	outfile.close()
3548
3549
3550
3551
3552	def tsh2sww(filename, verbose=False): #test_tsh2sww
3553	"""
3554	to check if a tsh/msh file 'looks' good.
3555	"""
3556
3557
3558	if verbose == True:print 'Creating domain from', filename
3559	domain = pmesh_to_domain_instance(filename, Domain)
3560	if verbose == True:print "Number of triangles = ", len(domain)
3561
3562	domain.smooth = True
3563	domain.format = 'sww' #Native netcdf visualisation format
3564	file_path, filename = path.split(filename)
3565	filename, ext = path.splitext(filename)
3566	domain.set_name(filename)
3567	domain.reduction = mean
3568	if verbose == True:print "file_path",file_path
3569	if file_path == "":file_path = "."
3570	domain.set_datadir(file_path)
3571
3572	if verbose == True:
3573	print "Output written to " + domain.get_datadir() + sep + \
3574	domain.get_name() + "." + domain.format
3575	sww = get_dataobject(domain)
3576	sww.store_connectivity()
3577	sww.store_timestep('stage')
3578
3579
3580	def asc_csiro2sww(bath_dir,
3581	elevation_dir,
3582	ucur_dir,
3583	vcur_dir,
3584	sww_file,
3585	minlat = None, maxlat = None,
3586	minlon = None, maxlon = None,
3587	zscale=1,
3588	mean_stage = 0,
3589	fail_on_NaN = True,
3590	elevation_NaN_filler = 0,
3591	bath_prefix='ba',
3592	elevation_prefix='el',
3593	verbose=False):
3594	"""
3595	Produce an sww boundary file, from esri ascii data from CSIRO.
3596
3597	Also convert latitude and longitude to UTM. All coordinates are
3598	assumed to be given in the GDA94 datum.
3599
3600	assume:
3601	All files are in esri ascii format
3602
3603	4 types of information
3604	bathymetry
3605	elevation
3606	u velocity
3607	v velocity
3608
3609	Assumptions
3610	The metadata of all the files is the same
3611	Each type is in a seperate directory
3612	One bath file with extention .000
3613	The time period is less than 24hrs and uniform.
3614	"""
3615	from Scientific.IO.NetCDF import NetCDFFile
3616
3617	from anuga.coordinate_transforms.redfearn import redfearn
3618
3619	precision = Float # So if we want to change the precision its done here
3620
3621	# go in to the bath dir and load the only file,
3622	bath_files = os.listdir(bath_dir)
3623
3624	bath_file = bath_files[0]
3625	bath_dir_file = bath_dir + os.sep + bath_file
3626	bath_metadata,bath_grid = _read_asc(bath_dir_file)
3627
3628	#Use the date.time of the bath file as a basis for
3629	#the start time for other files
3630	base_start = bath_file[-12:]
3631
3632	#go into the elevation dir and load the 000 file
3633	elevation_dir_file = elevation_dir + os.sep + elevation_prefix \
3634	+ base_start
3635
3636	elevation_files = os.listdir(elevation_dir)
3637	ucur_files = os.listdir(ucur_dir)
3638	vcur_files = os.listdir(vcur_dir)
3639	elevation_files.sort()
3640	# the first elevation file should be the
3641	# file with the same base name as the bath data
3642	assert elevation_files[0] == 'el' + base_start
3643
3644	number_of_latitudes = bath_grid.shape[0]
3645	number_of_longitudes = bath_grid.shape[1]
3646	number_of_volumes = (number_of_latitudes-1)(number_of_longitudes-1)2
3647
3648	longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
3649	for x in range(number_of_longitudes)]
3650	latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
3651	for y in range(number_of_latitudes)]
3652
3653	# reverse order of lat, so the fist lat represents the first grid row
3654	latitudes.reverse()
3655
3656	kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:],
3657	minlat=minlat, maxlat=maxlat,
3658	minlon=minlon, maxlon=maxlon)
3659
3660
3661	bath_grid = bath_grid[kmin:kmax,lmin:lmax]
3662	latitudes = latitudes[kmin:kmax]
3663	longitudes = longitudes[lmin:lmax]
3664	number_of_latitudes = len(latitudes)
3665	number_of_longitudes = len(longitudes)
3666	number_of_times = len(os.listdir(elevation_dir))
3667	number_of_points = number_of_latitudes*number_of_longitudes
3668	number_of_volumes = (number_of_latitudes-1)(number_of_longitudes-1)2
3669
3670	#Work out the times
3671	if len(elevation_files) > 1:
3672	# Assume: The time period is less than 24hrs.
3673	time_period = (int(elevation_files[1][-3:]) - \
3674	int(elevation_files[0][-3:]))6060
3675	times = [x*time_period for x in range(len(elevation_files))]
3676	else:
3677	times = [0.0]
3678
3679
3680	if verbose:
3681	print '------------------------------------------------'
3682	print 'Statistics:'
3683	print ' Extent (lat/lon):'
3684	print ' lat in [%f, %f], len(lat) == %d'\
3685	%(min(latitudes), max(latitudes),
3686	len(latitudes))
3687	print ' lon in [%f, %f], len(lon) == %d'\
3688	%(min(longitudes), max(longitudes),
3689	len(longitudes))
3690	print ' t in [%f, %f], len(t) == %d'\
3691	%(min(times), max(times), len(times))
3692
3693	######### WRITE THE SWW FILE #############
3694	# NetCDF file definition
3695	outfile = NetCDFFile(sww_file, 'w')
3696
3697	#Create new file
3698	outfile.institution = 'Geoscience Australia'
3699	outfile.description = 'Converted from XXX'
3700
3701
3702	#For sww compatibility
3703	outfile.smoothing = 'Yes'
3704	outfile.order = 1
3705
3706	#Start time in seconds since the epoch (midnight 1/1/1970)
3707	outfile.starttime = starttime = times[0]
3708
3709
3710	# dimension definitions
3711	outfile.createDimension('number_of_volumes', number_of_volumes)
3712
3713	outfile.createDimension('number_of_vertices', 3)
3714	outfile.createDimension('number_of_points', number_of_points)
3715	outfile.createDimension('number_of_timesteps', number_of_times)
3716
3717	# variable definitions
3718	outfile.createVariable('x', precision, ('number_of_points',))
3719	outfile.createVariable('y', precision, ('number_of_points',))
3720	outfile.createVariable('elevation', precision, ('number_of_points',))
3721
3722	#FIXME: Backwards compatibility
3723	outfile.createVariable('z', precision, ('number_of_points',))
3724	#################################
3725
3726	outfile.createVariable('volumes', Int, ('number_of_volumes',
3727	'number_of_vertices'))
3728
3729	outfile.createVariable('time', precision,
3730	('number_of_timesteps',))
3731
3732	outfile.createVariable('stage', precision,
3733	('number_of_timesteps',
3734	'number_of_points'))
3735
3736	outfile.createVariable('xmomentum', precision,
3737	('number_of_timesteps',
3738	'number_of_points'))
3739
3740	outfile.createVariable('ymomentum', precision,
3741	('number_of_timesteps',
3742	'number_of_points'))
3743
3744	#Store
3745	from anuga.coordinate_transforms.redfearn import redfearn
3746	x = zeros(number_of_points, Float) #Easting
3747	y = zeros(number_of_points, Float) #Northing
3748
3749	if verbose: print 'Making triangular grid'
3750	#Get zone of 1st point.
3751	refzone, _, _ = redfearn(latitudes[0],longitudes[0])
3752
3753	vertices = {}
3754	i = 0
3755	for k, lat in enumerate(latitudes):
3756	for l, lon in enumerate(longitudes):
3757
3758	vertices[l,k] = i
3759
3760	zone, easting, northing = redfearn(lat,lon)
3761
3762	msg = 'Zone boundary crossed at longitude =', lon
3763	#assert zone == refzone, msg
3764	#print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
3765	x[i] = easting
3766	y[i] = northing
3767	i += 1
3768
3769
3770	#Construct 2 triangles per 'rectangular' element
3771	volumes = []
3772	for l in range(number_of_longitudes-1): #X direction
3773	for k in range(number_of_latitudes-1): #Y direction
3774	v1 = vertices[l,k+1]
3775	v2 = vertices[l,k]
3776	v3 = vertices[l+1,k+1]
3777	v4 = vertices[l+1,k]
3778
3779	#Note, this is different to the ferrit2sww code
3780	#since the order of the lats is reversed.
3781	volumes.append([v1,v3,v2]) #Upper element
3782	volumes.append([v4,v2,v3]) #Lower element
3783
3784	volumes = array(volumes)
3785
3786	geo_ref = Geo_reference(refzone,min(x),min(y))
3787	geo_ref.write_NetCDF(outfile)
3788
3789	# This will put the geo ref in the middle
3790	#geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
3791
3792
3793	if verbose:
3794	print '------------------------------------------------'
3795	print 'More Statistics:'
3796	print ' Extent (/lon):'
3797	print ' x in [%f, %f], len(lat) == %d'\
3798	%(min(x), max(x),
3799	len(x))
3800	print ' y in [%f, %f], len(lon) == %d'\
3801	%(min(y), max(y),
3802	len(y))
3803	print 'geo_ref: ',geo_ref
3804
3805	z = resize(bath_grid,outfile.variables['z'][:].shape)
3806	outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
3807	outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
3808	outfile.variables['z'][:] = z
3809	outfile.variables['elevation'][:] = z #FIXME HACK
3810	outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
3811
3812	stage = outfile.variables['stage']
3813	xmomentum = outfile.variables['xmomentum']
3814	ymomentum = outfile.variables['ymomentum']
3815
3816	outfile.variables['time'][:] = times #Store time relative
3817
3818	if verbose: print 'Converting quantities'
3819	n = number_of_times
3820	for j in range(number_of_times):
3821	# load in files
3822	elevation_meta, elevation_grid = \
3823	_read_asc(elevation_dir + os.sep + elevation_files[j])
3824
3825	_, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j])
3826	_, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j])
3827
3828	#cut matrix to desired size
3829	elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
3830	u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
3831	v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
3832
3833	# handle missing values
3834	missing = (elevation_grid == elevation_meta['NODATA_value'])
3835	if sometrue (missing):
3836	if fail_on_NaN:
3837	msg = 'File %s contains missing values'\
3838	%(elevation_files[j])
3839	raise DataMissingValuesError, msg
3840	else:
3841	elevation_grid = elevation_grid*(missing==0) + \
3842	missing*elevation_NaN_filler
3843
3844
3845	if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n)
3846	i = 0
3847	for k in range(number_of_latitudes): #Y direction
3848	for l in range(number_of_longitudes): #X direction
3849	w = zscale*elevation_grid[k,l] + mean_stage
3850	stage[j,i] = w
3851	h = w - z[i]
3852	xmomentum[j,i] = u_momentum_grid[k,l]*h
3853	ymomentum[j,i] = v_momentum_grid[k,l]*h
3854	i += 1
3855	outfile.close()
3856
3857	def _get_min_max_indexes(latitudes_ref,longitudes_ref,
3858	minlat=None, maxlat=None,
3859	minlon=None, maxlon=None):
3860	"""
3861	return max, min indexes (for slicing) of the lat and long arrays to cover the area
3862	specified with min/max lat/long
3863
3864	Think of the latitudes and longitudes describing a 2d surface.
3865	The area returned is, if possible, just big enough to cover the
3866	inputed max/min area. (This will not be possible if the max/min area
3867	has a section outside of the latitudes/longitudes area.)
3868
3869	asset longitudes are sorted,
3870	long - from low to high (west to east, eg 148 - 151)
3871	assert latitudes are sorted, ascending or decending
3872	"""
3873	latitudes = latitudes_ref[:]
3874	longitudes = longitudes_ref[:]
3875
3876	latitudes = ensure_numeric(latitudes)
3877	longitudes = ensure_numeric(longitudes)
3878
3879	assert allclose(sort(longitudes), longitudes)
3880
3881	lat_ascending = True
3882	if not allclose(sort(latitudes), latitudes):
3883	lat_ascending = False
3884	# reverse order of lat, so it's in ascending order
3885	latitudes = latitudes[::-1]
3886	assert allclose(sort(latitudes), latitudes)
3887	#print "latitudes in funct", latitudes
3888
3889	largest_lat_index = len(latitudes)-1
3890	#Cut out a smaller extent.
3891	if minlat == None:
3892	lat_min_index = 0
3893	else:
3894	lat_min_index = searchsorted(latitudes, minlat)-1
3895	if lat_min_index <0:
3896	lat_min_index = 0
3897
3898
3899	if maxlat == None:
3900	lat_max_index = largest_lat_index #len(latitudes)
3901	else:
3902	lat_max_index = searchsorted(latitudes, maxlat)
3903	if lat_max_index > largest_lat_index:
3904	lat_max_index = largest_lat_index
3905
3906	if minlon == None:
3907	lon_min_index = 0
3908	else:
3909	lon_min_index = searchsorted(longitudes, minlon)-1
3910	if lon_min_index <0:
3911	lon_min_index = 0
3912
3913	if maxlon == None:
3914	lon_max_index = len(longitudes)
3915	else:
3916	lon_max_index = searchsorted(longitudes, maxlon)
3917
3918	# Reversing the indexes, if the lat array is decending
3919	if lat_ascending is False:
3920	lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
3921	largest_lat_index - lat_min_index
3922	lat_max_index = lat_max_index + 1 # taking into account how slicing works
3923	lon_max_index = lon_max_index + 1 # taking into account how slicing works
3924
3925	return lat_min_index, lat_max_index, lon_min_index, lon_max_index
3926
3927
3928	def _read_asc(filename, verbose=False):
3929	"""Read esri file from the following ASCII format (.asc)
3930
3931	Example:
3932
3933	ncols 3121
3934	nrows 1800
3935	xllcorner 722000
3936	yllcorner 5893000
3937	cellsize 25
3938	NODATA_value -9999
3939	138.3698 137.4194 136.5062 135.5558 ..........
3940
3941	"""
3942
3943	datafile = open(filename)
3944
3945	if verbose: print 'Reading DEM from %s' %(filename)
3946	lines = datafile.readlines()
3947	datafile.close()
3948
3949	if verbose: print 'Got', len(lines), ' lines'
3950
3951	ncols = int(lines.pop(0).split()[1].strip())
3952	nrows = int(lines.pop(0).split()[1].strip())
3953	xllcorner = float(lines.pop(0).split()[1].strip())
3954	yllcorner = float(lines.pop(0).split()[1].strip())
3955	cellsize = float(lines.pop(0).split()[1].strip())
3956	NODATA_value = float(lines.pop(0).split()[1].strip())
3957
3958	assert len(lines) == nrows
3959
3960	#Store data
3961	grid = []
3962
3963	n = len(lines)
3964	for i, line in enumerate(lines):
3965	cells = line.split()
3966	assert len(cells) == ncols
3967	grid.append(array([float(x) for x in cells]))
3968	grid = array(grid)
3969
3970	return {'xllcorner':xllcorner,
3971	'yllcorner':yllcorner,
3972	'cellsize':cellsize,
3973	'NODATA_value':NODATA_value}, grid
3974
3975
3976
3977	#### URS 2 SWW ###
3978
3979	lon_name = 'LON'
3980	lat_name = 'LAT'
3981	time_name = 'TIME'
3982	precision = Float # So if we want to change the precision its done here
3983	class Write_nc:
3984	"""
3985	Write an nc file.
3986
3987	Note, this should be checked to meet cdc netcdf conventions for gridded
3988	data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
3989
3990	"""
3991	def __init__(self,
3992	quantity_name,
3993	file_name,
3994	time_step_count,
3995	time_step,
3996	lon,
3997	lat):
3998	"""
3999	time_step_count is the number of time steps.
4000	time_step is the time step size
4001
4002	pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
4003	"""
4004	self.quantity_name = quantity_name
4005	quantity_units = {'HA':'CENTIMETERS',
4006	'UA':'CENTIMETERS/SECOND',
4007	'VA':'CENTIMETERS/SECOND'}
4008
4009	multiplier_dic = {'HA':100.0, # To convert from m to cm
4010	'UA':100.0, # m/s to cm/sec
4011	'VA':-100.0} # MUX files have positve x in the
4012	# Southern direction. This corrects for it, when writing nc files.
4013
4014	self.quantity_multiplier = multiplier_dic[self.quantity_name]
4015
4016	#self.file_name = file_name
4017	self.time_step_count = time_step_count
4018	self.time_step = time_step
4019
4020	# NetCDF file definition
4021	self.outfile = NetCDFFile(file_name, 'w')
4022	outfile = self.outfile
4023
4024	#Create new file
4025	nc_lon_lat_header(outfile, lon, lat)
4026
4027	# TIME
4028	outfile.createDimension(time_name, None)
4029	outfile.createVariable(time_name, precision, (time_name,))
4030
4031	#QUANTITY
4032	outfile.createVariable(self.quantity_name, precision,
4033	(time_name, lat_name, lon_name))
4034	outfile.variables[self.quantity_name].missing_value=-1.e+034
4035	outfile.variables[self.quantity_name].units= \
4036	quantity_units[self.quantity_name]
4037	outfile.variables[lon_name][:]= ensure_numeric(lon)
4038	outfile.variables[lat_name][:]= ensure_numeric(lat)
4039
4040	#Assume no one will be wanting to read this, while we are writing
4041	#outfile.close()
4042
4043	def store_timestep(self, quantity_slice):
4044	"""
4045	Write a time slice of quantity info
4046	quantity_slice is the data to be stored at this time step
4047	"""
4048
4049	outfile = self.outfile
4050
4051	# Get the variables
4052	time = outfile.variables[time_name]
4053	quantity = outfile.variables[self.quantity_name]
4054
4055	i = len(time)
4056
4057	#Store time
4058	time[i] = i*self.time_step #self.domain.time
4059	quantity[i,:] = quantity_slice* self.quantity_multiplier
4060
4061	def close(self):
4062	self.outfile.close()
4063
4064	def urs2sww(basename_in='o', basename_out=None, verbose=False,
4065	remove_nc_files=True,
4066	minlat=None, maxlat=None,
4067	minlon= None, maxlon=None,
4068	mint=None, maxt=None,
4069	mean_stage=0,
4070	origin = None,
4071	zscale=1,
4072	fail_on_NaN=True,
4073	NaN_filler=0,
4074	elevation=None):
4075	"""
4076	Convert URS C binary format for wave propagation to
4077	sww format native to abstract_2d_finite_volumes.
4078
4079	Specify only basename_in and read files of the form
4080	basefilename_velocity-z-mux, basefilename_velocity-e-mux and
4081	basefilename_waveheight-n-mux containing relative height,
4082	x-velocity and y-velocity, respectively.
4083
4084	Also convert latitude and longitude to UTM. All coordinates are
4085	assumed to be given in the GDA94 datum. The latitude and longitude
4086	information is for a grid.
4087
4088	min's and max's: If omitted - full extend is used.
4089	To include a value min may equal it, while max must exceed it.
4090	Lat and lon are assumed to be in decimal degrees.
4091	NOTE: minlon is the most east boundary.
4092
4093	origin is a 3-tuple with geo referenced
4094	UTM coordinates (zone, easting, northing)
4095	It will be the origin of the sww file. This shouldn't be used,
4096	since all of anuga should be able to handle an arbitary origin.
4097
4098
4099	URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
4100	which means that latitude is the fastest
4101	varying dimension (row major order, so to speak)
4102
4103	In URS C binary the latitudes and longitudes are in assending order.
4104	"""
4105	if basename_out == None:
4106	basename_out = basename_in
4107	files_out = urs2nc(basename_in, basename_out)
4108	ferret2sww(basename_out,
4109	minlat=minlat,
4110	maxlat=maxlat,
4111	minlon=minlon,
4112	maxlon=maxlon,
4113	mint=mint,
4114	maxt=maxt,
4115	mean_stage=mean_stage,
4116	origin=origin,
4117	zscale=zscale,
4118	fail_on_NaN=fail_on_NaN,
4119	NaN_filler=NaN_filler,
4120	inverted_bathymetry=True,
4121	verbose=verbose)
4122	#print "files_out",files_out
4123	if remove_nc_files:
4124	for file_out in files_out:
4125	os.remove(file_out)
4126
4127	def urs2nc(basename_in = 'o', basename_out = 'urs'):
4128	"""
4129	Convert the 3 urs files to 4 nc files.
4130
4131	The name of the urs file names must be;
4132	[basename_in]_velocity-z-mux
4133	[basename_in]_velocity-e-mux
4134	[basename_in]_waveheight-n-mux
4135
4136	"""
4137
4138	files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4139	basename_in + EAST_VELOCITY_LABEL,
4140	basename_in + NORTH_VELOCITY_LABEL]
4141	files_out = [basename_out+'_ha.nc',
4142	basename_out+'_ua.nc',
4143	basename_out+'_va.nc']
4144	quantities = ['HA','UA','VA']
4145
4146	for file_name in files_in:
4147	if os.access(file_name, os.F_OK) == 0 :
4148	msg = 'File %s does not exist or is not accessible' %file_name
4149	raise IOError, msg
4150
4151	hashed_elevation = None
4152	for file_in, file_out, quantity in map(None, files_in,
4153	files_out,
4154	quantities):
4155	lonlatdep, lon, lat, depth = _binary_c2nc(file_in,
4156	file_out,
4157	quantity)
4158	#print "lonlatdep", lonlatdep
4159	if hashed_elevation == None:
4160	elevation_file = basename_out+'_e.nc'
4161	write_elevation_nc(elevation_file,
4162	lon,
4163	lat,
4164	depth)
4165	hashed_elevation = myhash(lonlatdep)
4166	else:
4167	msg = "The elevation information in the mux files is inconsistent"
4168	assert hashed_elevation == myhash(lonlatdep), msg
4169	files_out.append(elevation_file)
4170	return files_out
4171
4172	def _binary_c2nc(file_in, file_out, quantity):
4173	"""
4174	Reads in a quantity urs file and writes a quantity nc file.
4175	additionally, returns the depth and lat, long info,
4176	so it can be written to a file.
4177	"""
4178	columns = 3 # long, lat , depth
4179	mux_file = open(file_in, 'rb')
4180
4181	# Number of points/stations
4182	(points_num,)= unpack('i',mux_file.read(4))
4183
4184	# nt, int - Number of time steps
4185	(time_step_count,)= unpack('i',mux_file.read(4))
4186
4187	#dt, float - time step, seconds
4188	(time_step,) = unpack('f', mux_file.read(4))
4189
4190	msg = "Bad data in the mux file."
4191	if points_num < 0:
4192	mux_file.close()
4193	raise ANUGAError, msg
4194	if time_step_count < 0:
4195	mux_file.close()
4196	raise ANUGAError, msg
4197	if time_step < 0:
4198	mux_file.close()
4199	raise ANUGAError, msg
4200
4201	lonlatdep = p_array.array('f')
4202	lonlatdep.read(mux_file, columns * points_num)
4203	lonlatdep = array(lonlatdep, typecode=Float)
4204	lonlatdep = reshape(lonlatdep, (points_num, columns))
4205
4206	lon, lat, depth = lon_lat2grid(lonlatdep)
4207	lon_sorted = list(lon)
4208	lon_sorted.sort()
4209
4210	if not lon == lon_sorted:
4211	msg = "Longitudes in mux file are not in ascending order"
4212	raise IOError, msg
4213	lat_sorted = list(lat)
4214	lat_sorted.sort()
4215
4216	if not lat == lat_sorted:
4217	msg = "Latitudes in mux file are not in ascending order"
4218
4219	nc_file = Write_nc(quantity,
4220	file_out,
4221	time_step_count,
4222	time_step,
4223	lon,
4224	lat)
4225
4226	for i in range(time_step_count):
4227	#Read in a time slice from mux file
4228	hz_p_array = p_array.array('f')
4229	hz_p_array.read(mux_file, points_num)
4230	hz_p = array(hz_p_array, typecode=Float)
4231	hz_p = reshape(hz_p, (len(lon), len(lat)))
4232	hz_p = transpose(hz_p) #mux has lat varying fastest, nc has long v.f.
4233
4234	#write time slice to nc file
4235	nc_file.store_timestep(hz_p)
4236	mux_file.close()
4237	nc_file.close()
4238
4239	return lonlatdep, lon, lat, depth
4240
4241
4242	def write_elevation_nc(file_out, lon, lat, depth_vector):
4243	"""
4244	Write an nc elevation file.
4245	"""
4246
4247	# NetCDF file definition
4248	outfile = NetCDFFile(file_out, 'w')
4249
4250	#Create new file
4251	nc_lon_lat_header(outfile, lon, lat)
4252
4253	# ELEVATION
4254	zname = 'ELEVATION'
4255	outfile.createVariable(zname, precision, (lat_name, lon_name))
4256	outfile.variables[zname].units='CENTIMETERS'
4257	outfile.variables[zname].missing_value=-1.e+034
4258
4259	outfile.variables[lon_name][:]= ensure_numeric(lon)
4260	outfile.variables[lat_name][:]= ensure_numeric(lat)
4261
4262	depth = reshape(depth_vector, ( len(lat), len(lon)))
4263	outfile.variables[zname][:]= depth
4264
4265	outfile.close()
4266
4267	def nc_lon_lat_header(outfile, lon, lat):
4268	"""
4269	outfile is the netcdf file handle.
4270	lon - a list/array of the longitudes
4271	lat - a list/array of the latitudes
4272	"""
4273
4274	outfile.institution = 'Geoscience Australia'
4275	outfile.description = 'Converted from URS binary C'
4276
4277	# Longitude
4278	outfile.createDimension(lon_name, len(lon))
4279	outfile.createVariable(lon_name, precision, (lon_name,))
4280	outfile.variables[lon_name].point_spacing='uneven'
4281	outfile.variables[lon_name].units='degrees_east'
4282	outfile.variables[lon_name].assignValue(lon)
4283
4284
4285	# Latitude
4286	outfile.createDimension(lat_name, len(lat))
4287	outfile.createVariable(lat_name, precision, (lat_name,))
4288	outfile.variables[lat_name].point_spacing='uneven'
4289	outfile.variables[lat_name].units='degrees_north'
4290	outfile.variables[lat_name].assignValue(lat)
4291
4292
4293
4294	def lon_lat2grid(long_lat_dep):
4295	"""
4296	given a list of points that are assumed to be an a grid,
4297	return the long's and lat's of the grid.
4298	long_lat_dep is an array where each row is a position.
4299	The first column is longitudes.
4300	The second column is latitudes.
4301
4302	The latitude is the fastest varying dimension - in mux files
4303	"""
4304	LONG = 0
4305	LAT = 1
4306	QUANTITY = 2
4307
4308	long_lat_dep = ensure_numeric(long_lat_dep, Float)
4309
4310	num_points = long_lat_dep.shape[0]
4311	this_rows_long = long_lat_dep[0,LONG]
4312
4313	# Count the length of unique latitudes
4314	i = 0
4315	while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
4316	i += 1
4317	# determine the lats and longsfrom the grid
4318	lat = long_lat_dep[:i, LAT]
4319	long = long_lat_dep[::i, LONG]
4320
4321	lenlong = len(long)
4322	lenlat = len(lat)
4323	#print 'len lat', lat, len(lat)
4324	#print 'len long', long, len(long)
4325
4326	msg = 'Input data is not gridded'
4327	assert num_points % lenlat == 0, msg
4328	assert num_points % lenlong == 0, msg
4329
4330	# Test that data is gridded
4331	for i in range(lenlong):
4332	msg = 'Data is not gridded. It must be for this operation'
4333	first = i*lenlat
4334	last = first + lenlat
4335
4336	assert allclose(long_lat_dep[first:last,LAT], lat), msg
4337	assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
4338
4339
4340	# print 'range long', min(long), max(long)
4341	# print 'range lat', min(lat), max(lat)
4342	# print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0])
4343	# print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1])
4344
4345
4346
4347	msg = 'Out of range latitudes/longitudes'
4348	for l in lat:assert -90 < l < 90 , msg
4349	for l in long:assert -180 < l < 180 , msg
4350
4351	# Changing quantity from lat being the fastest varying dimension to
4352	# long being the fastest varying dimension
4353	# FIXME - make this faster/do this a better way
4354	# use numeric transpose, after reshaping the quantity vector
4355	# quantity = zeros(len(long_lat_dep), Float)
4356	quantity = zeros(num_points, Float)
4357
4358	# print 'num',num_points
4359	for lat_i, _ in enumerate(lat):
4360	for long_i, _ in enumerate(long):
4361	q_index = lat_i*lenlong+long_i
4362	lld_index = long_i*lenlat+lat_i
4363	# print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index
4364	temp = long_lat_dep[lld_index, QUANTITY]
4365	quantity[q_index] = temp
4366
4367	return long, lat, quantity
4368
4369	#### END URS 2 SWW ###
4370
4371	#### URS UNGRIDDED 2 SWW ###
4372
4373	### PRODUCING THE POINTS NEEDED FILE ###
4374	LL_LAT = -50.0
4375	LL_LONG = 80.0
4376	GRID_SPACING = 1.0/60.0
4377	LAT_AMOUNT = 4800
4378	LONG_AMOUNT = 3600
4379	def URS_points_needed_to_file(file_name, boundary_polygon, zone,
4380	ll_lat=LL_LAT, ll_long=LL_LONG,
4381	grid_spacing=GRID_SPACING,
4382	lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4383	export_csv=False, use_cache=False,
4384	verbose=False):
4385	"""
4386	file_name - name of the urs file produced for David.
4387	boundary_polygon - a list of points that describes a polygon.
4388	The last point is assumed ot join the first point.
4389	This is in UTM (lat long would be better though)
4390
4391	ll_lat - lower left latitude, in decimal degrees
4392	ll-long - lower left longitude, in decimal degrees
4393	grid_spacing - in deciamal degrees
4394
4395
4396	Don't add the file extension. It will be added.
4397	"""
4398	geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long,
4399	grid_spacing,
4400	lat_amount, long_amount,use_cache, verbose)
4401	if not file_name[-4:] == ".urs":
4402	file_name += ".urs"
4403	geo.export_points_file(file_name)
4404	if export_csv:
4405	if file_name[-4:] == ".urs":
4406	file_name = file_name[:-4] + ".csv"
4407	geo.export_points_file(file_name)
4408
4409	def URS_points_needed(boundary_polygon, zone, ll_lat=LL_LAT,
4410	ll_long=LL_LONG, grid_spacing=GRID_SPACING,
4411	lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4412	use_cache=False, verbose=False):
4413	args = (boundary_polygon,
4414	zone)
4415	kwargs = {'ll_lat': ll_lat,
4416	'll_long': ll_long,
4417	'grid_spacing': grid_spacing,
4418	'lat_amount': lat_amount,
4419	'long_amount': long_amount}
4420	if use_cache is True:
4421	try:
4422	from anuga.caching import cache
4423	except:
4424	msg = 'Caching was requested, but caching module'+\
4425	'could not be imported'
4426	raise msg
4427
4428
4429	geo = cache(_URS_points_needed,
4430	args, kwargs,
4431	verbose=verbose,
4432	compression=False)
4433	else:
4434	#I was getting 'got multiple values for keyword argument' errors
4435	#geo = apply(_URS_points_needed, args, kwargs)
4436	geo = _URS_points_needed(boundary_polygon,
4437	zone, ll_lat,
4438	ll_long, grid_spacing,
4439	lat_amount, long_amount)
4440
4441	return geo
4442	def _URS_points_needed(boundary_polygon,
4443	zone, ll_lat=LL_LAT,
4444	ll_long=LL_LONG, grid_spacing=GRID_SPACING,
4445	lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT):
4446	"""
4447
4448	boundary_polygon - a list of points that describes a polygon.
4449	The last point is assumed ot join the first point.
4450	This is in UTM (lat long would be better though)
4451
4452	ll_lat - lower left latitude, in decimal degrees
4453	ll-long - lower left longitude, in decimal degrees
4454	grid_spacing - in deciamal degrees
4455
4456	"""
4457
4458	from sets import ImmutableSet
4459
4460	msg = "grid_spacing can not be zero"
4461	assert not grid_spacing ==0, msg
4462	a = boundary_polygon
4463	# List of segments. Each segment is two points.
4464	segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
4465
4466	# convert the segs to Lat's and longs.
4467
4468	# Don't assume the zone of the segments is the same as the lower left
4469	# corner of the lat long data!! They can easily be in different zones
4470
4471	lat_long_set = ImmutableSet()
4472	for seg in segs:
4473	points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
4474	lat_amount, long_amount, zone)
4475	lat_long_set \|= ImmutableSet(points_lat_long)
4476	#print "lat_long_set",lat_long_set
4477	geo = Geospatial_data(data_points=list(lat_long_set),
4478	points_are_lats_longs=True)
4479	return geo
4480
4481	def points_needed(seg, ll_lat, ll_long, grid_spacing,
4482	lat_amount, long_amount, zone):
4483	"""
4484	return a list of the points, in lats and longs that are needed to
4485	interpolate any point on the segment.
4486	"""
4487	from math import sqrt
4488	#print "zone",zone
4489	geo_reference = Geo_reference(zone=zone)
4490	#print "seg",seg
4491	geo = Geospatial_data(seg,geo_reference=geo_reference)
4492	seg_lat_long = geo.get_data_points(as_lat_long=True)
4493	#print "seg_lat_long", seg_lat_long
4494	# 1.415 = 2^0.5, rounded up....
4495	sqrt_2_rounded_up = 1.415
4496	buffer = sqrt_2_rounded_up * grid_spacing
4497
4498	max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
4499	max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
4500	min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
4501	min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
4502
4503	first_row = (min_long - ll_long)/grid_spacing
4504	# To round up
4505	first_row_long = int(round(first_row + 0.5))
4506	#print "first_row", first_row_long
4507
4508	last_row = (max_long - ll_long)/grid_spacing # round down
4509	last_row_long = int(round(last_row))
4510	#print "last_row",last_row _long
4511
4512	first_row = (min_lat - ll_lat)/grid_spacing
4513	# To round up
4514	first_row_lat = int(round(first_row + 0.5))
4515	#print "first_row", first_row_lat
4516
4517	last_row = (max_lat - ll_lat)/grid_spacing # round down
4518	last_row_lat = int(round(last_row))
4519	#print "last_row",last_row_lat
4520
4521	# to work out the max distance -
4522	# 111120 - horizontal distance between 1 deg latitude.
4523	#max_distance = sqrt_2_rounded_up * 111120 * grid_spacing
4524	max_distance = 157147.4112 * grid_spacing
4525	#print "max_distance", max_distance #2619.12 m for 1 minute
4526	points_lat_long = []
4527	# Create a list of the lat long points to include.
4528	for index_lat in range(first_row_lat, last_row_lat + 1):
4529	for index_long in range(first_row_long, last_row_long + 1):
4530	lat = ll_lat + index_lat*grid_spacing
4531	long = ll_long + index_long*grid_spacing
4532
4533	#filter here to keep good points
4534	if keep_point(lat, long, seg, max_distance):
4535	points_lat_long.append((lat, long)) #must be hashable
4536	#print "points_lat_long", points_lat_long
4537
4538	# Now that we have these points, lets throw ones out that are too far away
4539	return points_lat_long
4540
4541	def keep_point(lat, long, seg, max_distance):
4542	"""
4543	seg is two points, UTM
4544	"""
4545	from math import sqrt
4546	_ , x0, y0 = redfearn(lat, long)
4547	x1 = seg[0][0]
4548	y1 = seg[0][1]
4549	x2 = seg[1][0]
4550	y2 = seg[1][1]
4551
4552	x2_1 = x2-x1
4553	y2_1 = y2-y1
4554	d = abs((x2_1)(y1-y0)-(x1-x0)(y2_1))/sqrt((x2_1)(x2_1)+(y2_1)(y2_1))
4555	if d <= max_distance:
4556	return True
4557	else:
4558	return False
4559
4560	#### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE ####
4561
4562	WAVEHEIGHT_MUX_LABEL = '_waveheight-z-mux'
4563	EAST_VELOCITY_LABEL = '_velocity-e-mux'
4564	NORTH_VELOCITY_LABEL = '_velocity-n-mux'
4565	def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
4566	mint=None, maxt=None,
4567	mean_stage=0,
4568	origin=None,
4569	hole_points_UTM=None,
4570	zscale=1):
4571	"""
4572	Convert URS C binary format for wave propagation to
4573	sww format native to abstract_2d_finite_volumes.
4574
4575
4576	Specify only basename_in and read files of the form
4577	basefilename_velocity-z-mux, basefilename_velocity-e-mux and
4578	basefilename_waveheight-n-mux containing relative height,
4579	x-velocity and y-velocity, respectively.
4580
4581	Also convert latitude and longitude to UTM. All coordinates are
4582	assumed to be given in the GDA94 datum. The latitude and longitude
4583	information is assumed ungridded grid.
4584
4585	min's and max's: If omitted - full extend is used.
4586	To include a value min ans max may equal it.
4587	Lat and lon are assumed to be in decimal degrees.
4588
4589	origin is a 3-tuple with geo referenced
4590	UTM coordinates (zone, easting, northing)
4591	It will be the origin of the sww file. This shouldn't be used,
4592	since all of anuga should be able to handle an arbitary origin.
4593	The mux point info is NOT relative to this origin.
4594
4595
4596	URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
4597	which means that latitude is the fastest
4598	varying dimension (row major order, so to speak)
4599
4600	In URS C binary the latitudes and longitudes are in assending order.
4601
4602	Note, interpolations of the resulting sww file will be different
4603	from results of urs2sww. This is due to the interpolation
4604	function used, and the different grid structure between urs2sww
4605	and this function.
4606
4607	Interpolating data that has an underlying gridded source can
4608	easily end up with different values, depending on the underlying
4609	mesh.
4610
4611	consider these 4 points
4612	50 -50
4613
4614	0 0
4615
4616	The grid can be
4617	-
4618	\|\\| A
4619	-
4620	or;
4621	-
4622	\|/\| B
4623	-
4624	If a point is just below the center of the midpoint, it will have a
4625	+ve value in grid A and a -ve value in grid B.
4626	"""
4627	from anuga.pmesh.mesh import Mesh, NoTrianglesError
4628
4629	files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4630	basename_in + EAST_VELOCITY_LABEL,
4631	basename_in + NORTH_VELOCITY_LABEL]
4632	quantities = ['HA','UA','VA']
4633
4634	# instanciate urs_points of the three mux files.
4635	mux = {}
4636	for quantity, file in map(None, quantities, files_in):
4637	mux[quantity] = Urs_points(file)
4638
4639	# Could check that the depth is the same. (hashing)
4640
4641	# handle to a mux file to do depth stuff
4642	a_mux = mux[quantities[0]]
4643
4644	# Convert to utm
4645	lat = a_mux.lonlatdep[:,1]
4646	long = a_mux.lonlatdep[:,0]
4647	points_utm, zone = convert_from_latlon_to_utm( \
4648	latitudes=lat, longitudes=long)
4649	#print "points_utm", points_utm
4650	#print "zone", zone
4651
4652	elevation = a_mux.lonlatdep[:,2] * -1 #
4653
4654	# grid ( create a mesh from the selected points)
4655	# This mesh has a problem. Triangles are streched over ungridded areas.
4656	# If these areas could be described as holes in pmesh, that would be great
4657
4658	# I can't just get the user to selection a point in the middle.
4659	# A boundary is needed around these points.
4660	# But if the zone of points is obvious enough auto-segment should do
4661	# a good boundary.
4662	mesh = Mesh()
4663	mesh.add_vertices(points_utm)
4664	mesh.auto_segment(smooth_indents=True, expand_pinch=True)
4665	# To try and avoid alpha shape 'hugging' too much
4666	mesh.auto_segment( mesh.shape.get_alpha()*1.1 )
4667	if hole_points_UTM is not None:
4668	point = ensure_absolute(hole_points_UTM)
4669	mesh.add_hole(point[0], point[1])
4670	try:
4671	mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4672	except NoTrianglesError:
4673	# This is a bit of a hack, going in and changing the
4674	# data structure.
4675	mesh.holes = []
4676	mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4677	mesh_dic = mesh.Mesh2MeshList()
4678
4679	#mesh.export_mesh_file(basename_in + '_168.tsh')
4680	#import sys; sys.exit()
4681	# These are the times of the mux file
4682	mux_times = []
4683	for i in range(a_mux.time_step_count):
4684	mux_times.append(a_mux.time_step * i)
4685	mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt)
4686	times = mux_times[mux_times_start_i:mux_times_fin_i]
4687
4688	if mux_times_start_i == mux_times_fin_i:
4689	# Close the mux files
4690	for quantity, file in map(None, quantities, files_in):
4691	mux[quantity].close()
4692	msg="Due to mint and maxt there's no time info in the boundary SWW."
4693	raise Exception, msg
4694
4695	# If this raise is removed there is currently no downstream errors
4696
4697	points_utm=ensure_numeric(points_utm)
4698	#print "mesh_dic['generatedpointlist']", mesh_dic['generatedpointlist']
4699	#print "points_utm", points_utm
4700	assert ensure_numeric(mesh_dic['generatedpointlist']) == \
4701	ensure_numeric(points_utm)
4702
4703	volumes = mesh_dic['generatedtrianglelist']
4704
4705	# write sww intro and grid stuff.
4706	if basename_out is None:
4707	swwname = basename_in + '.sww'
4708	else:
4709	swwname = basename_out + '.sww'
4710
4711	if verbose: print 'Output to ', swwname
4712	outfile = NetCDFFile(swwname, 'w')
4713	# For a different way of doing this, check out tsh2sww
4714	# work out sww_times and the index range this covers
4715	sww = Write_sww()
4716	sww.store_header(outfile, times, len(volumes), len(points_utm),
4717	verbose=verbose)
4718	outfile.mean_stage = mean_stage
4719	outfile.zscale = zscale
4720
4721	sww.store_triangulation(outfile, points_utm, volumes,
4722	elevation, zone, new_origin=origin,
4723	verbose=verbose)
4724
4725	if verbose: print 'Converting quantities'
4726	j = 0
4727	# Read in a time slice from each mux file and write it to the sww file
4728	for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
4729	if j >= mux_times_start_i and j < mux_times_fin_i:
4730	stage = zscale*ha + mean_stage
4731	h = stage - elevation
4732	xmomentum = ua*h
4733	ymomentum = -1vah # -1 since in mux files south is positive.
4734	sww.store_quantities(outfile,
4735	slice_index=j - mux_times_start_i,
4736	verbose=verbose,
4737	stage=stage,
4738	xmomentum=xmomentum,
4739	ymomentum=ymomentum)
4740	j += 1
4741	if verbose: sww.verbose_quantities(outfile)
4742	outfile.close()
4743	#
4744	# Do some conversions while writing the sww file
4745
4746
4747	def mux2sww_time(mux_times, mint, maxt):
4748	"""
4749	"""
4750
4751	if mint == None:
4752	mux_times_start_i = 0
4753	else:
4754	mux_times_start_i = searchsorted(mux_times, mint)
4755
4756	if maxt == None:
4757	mux_times_fin_i = len(mux_times)
4758	else:
4759	maxt += 0.5 # so if you specify a time where there is
4760	# data that time will be included
4761	mux_times_fin_i = searchsorted(mux_times, maxt)
4762
4763	return mux_times_start_i, mux_times_fin_i
4764
4765
4766	class Write_sww:
4767	from anuga.shallow_water.shallow_water_domain import Domain
4768
4769	# FIXME (Ole): Hardwiring the conserved quantities like
4770	# this could be a problem. I would prefer taking them from
4771	# the instantiation of Domain.
4772	sww_quantities = Domain.conserved_quantities
4773
4774
4775	RANGE = '_range'
4776	EXTREMA = ':extrema'
4777
4778	def __init__(self):
4779	pass
4780
4781	def store_header(self,
4782	outfile,
4783	times,
4784	number_of_volumes,
4785	number_of_points,
4786	description='Converted from XXX',
4787	smoothing=True,
4788	order=1, verbose=False):
4789	"""
4790	outfile - the name of the file that will be written
4791	times - A list of the time slice times OR a start time
4792	Note, if a list is given the info will be made relative.
4793	number_of_volumes - the number of triangles
4794	"""
4795
4796	outfile.institution = 'Geoscience Australia'
4797	outfile.description = description
4798
4799	# For sww compatibility
4800	if smoothing is True:
4801	# Smoothing to be depreciated
4802	outfile.smoothing = 'Yes'
4803	outfile.vertices_are_stored_uniquely = 'False'
4804	else:
4805	# Smoothing to be depreciated
4806	outfile.smoothing = 'No'
4807	outfile.vertices_are_stored_uniquely = 'True'
4808	outfile.order = order
4809
4810	try:
4811	revision_number = get_revision_number()
4812	except:
4813	revision_number = None
4814	# Allow None to be stored as a string
4815	outfile.revision_number = str(revision_number)
4816
4817
4818
4819	# times - A list or array of the time slice times OR a start time
4820	# times = ensure_numeric(times)
4821	# Start time in seconds since the epoch (midnight 1/1/1970)
4822
4823	# This is being used to seperate one number from a list.
4824	# what it is actually doing is sorting lists from numeric arrays.
4825	if type(times) is list or type(times) is ArrayType:
4826	number_of_times = len(times)
4827	times = ensure_numeric(times)
4828	if number_of_times == 0:
4829	starttime = 0
4830	else:
4831	starttime = times[0]
4832	times = times - starttime #Store relative times
4833	else:
4834	number_of_times = 0
4835	starttime = times
4836	#times = ensure_numeric([])
4837	outfile.starttime = starttime
4838	# dimension definitions
4839	outfile.createDimension('number_of_volumes', number_of_volumes)
4840	outfile.createDimension('number_of_vertices', 3)
4841	outfile.createDimension('numbers_in_range', 2)
4842
4843	if smoothing is True:
4844	outfile.createDimension('number_of_points', number_of_points)
4845
4846	# FIXME(Ole): This will cause sww files for paralle domains to
4847	# have ghost nodes stored (but not used by triangles).
4848	# To clean this up, we have to change get_vertex_values and
4849	# friends in quantity.py (but I can't be bothered right now)
4850	else:
4851	outfile.createDimension('number_of_points', 3*number_of_volumes)
4852	outfile.createDimension('number_of_timesteps', number_of_times)
4853
4854	# variable definitions
4855	outfile.createVariable('x', precision, ('number_of_points',))
4856	outfile.createVariable('y', precision, ('number_of_points',))
4857	outfile.createVariable('elevation', precision, ('number_of_points',))
4858	q = 'elevation'
4859	outfile.createVariable(q+Write_sww.RANGE, precision,
4860	('numbers_in_range',))
4861
4862
4863	# Initialise ranges with small and large sentinels.
4864	# If this was in pure Python we could have used None sensibly
4865	outfile.variables[q+Write_sww.RANGE][0] = max_float # Min
4866	outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
4867
4868	# FIXME: Backwards compatibility
4869	outfile.createVariable('z', precision, ('number_of_points',))
4870	#################################
4871
4872	outfile.createVariable('volumes', Int, ('number_of_volumes',
4873	'number_of_vertices'))
4874
4875	outfile.createVariable('time', precision,
4876	('number_of_timesteps',))
4877
4878	for q in Write_sww.sww_quantities:
4879	outfile.createVariable(q, precision,
4880	('number_of_timesteps',
4881	'number_of_points'))
4882	outfile.createVariable(q+Write_sww.RANGE, precision,
4883	('numbers_in_range',))
4884
4885	# Initialise ranges with small and large sentinels.
4886	# If this was in pure Python we could have used None sensibly
4887	outfile.variables[q+Write_sww.RANGE][0] = max_float # Min
4888	outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
4889
4890	if type(times) is list or type(times) is ArrayType:
4891	outfile.variables['time'][:] = times #Store time relative
4892
4893	if verbose:
4894	print '------------------------------------------------'
4895	print 'Statistics:'
4896	print ' t in [%f, %f], len(t) == %d'\
4897	%(min(times.flat), max(times.flat), len(times.flat))
4898
4899
4900	def store_triangulation(self,
4901	outfile,
4902	points_utm,
4903	volumes,
4904	elevation, zone=None, new_origin=None,
4905	points_georeference=None, verbose=False):
4906	"""
4907
4908	new_origin - qa georeference that the points can be set to. (Maybe
4909	do this before calling this function.)
4910
4911	points_utm - currently a list or array of the points in UTM.
4912	points_georeference - the georeference of the points_utm
4913
4914	How about passing new_origin and current_origin.
4915	If you get both, do a convertion from the old to the new.
4916
4917	If you only get new_origin, the points are absolute,
4918	convert to relative
4919
4920	if you only get the current_origin the points are relative, store
4921	as relative.
4922
4923	if you get no georefs create a new georef based on the minimums of
4924	points_utm. (Another option would be to default to absolute)
4925
4926	Yes, and this is done in another part of the code.
4927	Probably geospatial.
4928
4929	If you don't supply either geo_refs, then supply a zone. If not
4930	the default zone will be used.
4931
4932
4933	precon
4934
4935	header has been called.
4936	"""
4937
4938	number_of_points = len(points_utm)
4939	volumes = array(volumes)
4940	points_utm = array(points_utm)
4941
4942	# given the two geo_refs and the points, do the stuff
4943	# described in the method header
4944	# if this is needed else where, pull out as a function
4945	points_georeference = ensure_geo_reference(points_georeference)
4946	new_origin = ensure_geo_reference(new_origin)
4947	if new_origin is None and points_georeference is not None:
4948	points = points_utm
4949	geo_ref = points_georeference
4950	else:
4951	if new_origin is None:
4952	new_origin = Geo_reference(zone,min(points_utm[:,0]),
4953	min(points_utm[:,1]))
4954	points = new_origin.change_points_geo_ref(points_utm,
4955	points_georeference)
4956	geo_ref = new_origin
4957
4958	# At this stage I need a georef and points
4959	# the points are relative to the georef
4960	geo_ref.write_NetCDF(outfile)
4961
4962	# This will put the geo ref in the middle
4963	#geo_ref=Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
4964
4965	x = points[:,0]
4966	y = points[:,1]
4967	z = outfile.variables['z'][:]
4968
4969	if verbose:
4970	print '------------------------------------------------'
4971	print 'More Statistics:'
4972	print ' Extent (/lon):'
4973	print ' x in [%f, %f], len(lat) == %d'\
4974	%(min(x), max(x),
4975	len(x))
4976	print ' y in [%f, %f], len(lon) == %d'\
4977	%(min(y), max(y),
4978	len(y))
4979	print ' z in [%f, %f], len(z) == %d'\
4980	%(min(elevation), max(elevation),
4981	len(elevation))
4982	print 'geo_ref: ',geo_ref
4983	print '------------------------------------------------'
4984
4985	#z = resize(bath_grid,outfile.variables['z'][:].shape)
4986	outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
4987	outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
4988	outfile.variables['z'][:] = elevation
4989	outfile.variables['elevation'][:] = elevation #FIXME HACK
4990	outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
4991
4992	q = 'elevation'
4993	# This updates the _range values
4994	outfile.variables[q+Write_sww.RANGE][0] = min(elevation)
4995	outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
4996
4997
4998	def store_quantities(self, outfile, precision=Float,
4999	slice_index=None, time=None,
5000	verbose=False, **quant):
5001	"""
5002	Write the quantity info.
5003
5004	**quant is extra keyword arguments passed in. These must be
5005	the sww quantities, currently; stage, xmomentum, ymomentum.
5006
5007	if the time array is already been built, use the slice_index
5008	to specify the index.
5009
5010	Otherwise, use time to increase the time dimension
5011
5012	Maybe make this general, but the viewer assumes these quantities,
5013	so maybe we don't want it general - unless the viewer is general
5014
5015	precon
5016	triangulation and
5017	header have been called.
5018	"""
5019
5020	if time is not None:
5021	file_time = outfile.variables['time']
5022	slice_index = len(file_time)
5023	file_time[slice_index] = time
5024
5025	# write the conserved quantities from Domain.
5026	# Typically stage, xmomentum, ymomentum
5027	# other quantities will be ignored, silently.
5028	for q in Write_sww.sww_quantities:
5029	if not quant.has_key(q):
5030	msg = 'SWW file can not write quantity %s' %q
5031	raise NewQuantity, msg
5032	else:
5033	q_values = quant[q]
5034	outfile.variables[q][slice_index] = q_values.astype(precision)
5035
5036	# This updates the _range values
5037	q_range = outfile.variables[q+Write_sww.RANGE][:]
5038	q_values_min = min(q_values)
5039	if q_values_min < q_range[0]:
5040	outfile.variables[q+Write_sww.RANGE][0] = q_values_min
5041	q_values_max = max(q_values)
5042	if q_values_max > q_range[1]:
5043	outfile.variables[q+Write_sww.RANGE][1] = q_values_max
5044
5045	def verbose_quantities(self, outfile):
5046	print '------------------------------------------------'
5047	print 'More Statistics:'
5048	for q in Write_sww.sww_quantities:
5049	print ' %s in [%f, %f]' %(q,
5050	outfile.variables[q+Write_sww.RANGE][0],
5051	outfile.variables[q+Write_sww.RANGE][1])
5052	print '------------------------------------------------'
5053
5054
5055
5056	def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation,
5057	mean_stage=0, zscale=1,
5058	verbose=False):
5059	#Time stepping
5060	stage = outfile.variables['stage']
5061	xmomentum = outfile.variables['xmomentum']
5062	ymomentum = outfile.variables['ymomentum']
5063
5064	n = len(has)
5065	j=0
5066	for ha, ua, va in map(None, has, uas, vas):
5067	if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n)
5068	w = zscale*ha + mean_stage
5069	stage[j] = w
5070	h = w - elevation
5071	xmomentum[j] = ua*h
5072	ymomentum[j] = -1vah # -1 since in mux files south is positive.
5073	j += 1
5074
5075	def urs2txt(basename_in, location_index=None):
5076	"""
5077	Not finished or tested
5078	"""
5079
5080	files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
5081	basename_in + EAST_VELOCITY_LABEL,
5082	basename_in + NORTH_VELOCITY_LABEL]
5083	quantities = ['HA','UA','VA']
5084
5085	d = ","
5086
5087	# instanciate urs_points of the three mux files.
5088	mux = {}
5089	for quantity, file in map(None, quantities, files_in):
5090	mux[quantity] = Urs_points(file)
5091
5092	# Could check that the depth is the same. (hashing)
5093
5094	# handle to a mux file to do depth stuff
5095	a_mux = mux[quantities[0]]
5096
5097	# Convert to utm
5098	latitudes = a_mux.lonlatdep[:,1]
5099	longitudes = a_mux.lonlatdep[:,0]
5100	points_utm, zone = convert_from_latlon_to_utm( \
5101	latitudes=latitudes, longitudes=longitudes)
5102	#print "points_utm", points_utm
5103	#print "zone", zone
5104	depths = a_mux.lonlatdep[:,2] #
5105
5106	fid = open(basename_in + '.txt', 'w')
5107
5108	fid.write("zone: " + str(zone) + "\n")
5109
5110	if location_index is not None:
5111	#Title
5112	li = location_index
5113	fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5114	'Northing' + "\n")
5115	fid.write(str(li) +d+ str(latitudes[li])+d+ \
5116	str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \
5117	str(points_utm[li][01]) + "\n")
5118
5119	# the non-time dependent stuff
5120	#Title
5121	fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5122	'Northing' +d+ 'depth m' + "\n")
5123	i = 0
5124	for depth, point_utm, lat, long in map(None, depths,
5125	points_utm, latitudes,
5126	longitudes):
5127
5128	fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \
5129	str(point_utm[01]) +d+ str(depth) + "\n")
5130	i +=1
5131	#Time dependent
5132	if location_index is not None:
5133	time_step = a_mux.time_step
5134	i = 0
5135	#Title
5136	fid.write('time' +d+ 'HA depth m'+d+ \
5137	'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \
5138	+ "\n")
5139	for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
5140	fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \
5141	str(UA[location_index]) +d+ str(VA[location_index]) \
5142	+ "\n")
5143
5144	i +=1
5145
5146	class Urs_points:
5147	"""
5148	Read the info in URS mux files.
5149
5150	for the quantities heres a correlation between the file names and
5151	what they mean;
5152	z-mux is height above sea level, m
5153	e-mux is velocity is Eastern direction, m/s
5154	n-mux is velocity is Northern direction, m/s
5155	"""
5156	def __init__(self,urs_file):
5157	self.iterated = False
5158	columns = 3 # long, lat , depth
5159	mux_file = open(urs_file, 'rb')
5160
5161	# Number of points/stations
5162	(self.points_num,)= unpack('i',mux_file.read(4))
5163
5164	# nt, int - Number of time steps
5165	(self.time_step_count,)= unpack('i',mux_file.read(4))
5166	#print "self.time_step_count", self.time_step_count
5167	#dt, float - time step, seconds
5168	(self.time_step,) = unpack('f', mux_file.read(4))
5169	#print "self.time_step", self.time_step
5170	msg = "Bad data in the urs file."
5171	if self.points_num < 0:
5172	mux_file.close()
5173	raise ANUGAError, msg
5174	if self.time_step_count < 0:
5175	mux_file.close()
5176	raise ANUGAError, msg
5177	if self.time_step < 0:
5178	mux_file.close()
5179	raise ANUGAError, msg
5180
5181	# the depth is in meters, and it is the distance from the ocean
5182	# to the sea bottom.
5183	lonlatdep = p_array.array('f')
5184	lonlatdep.read(mux_file, columns * self.points_num)
5185	lonlatdep = array(lonlatdep, typecode=Float)
5186	lonlatdep = reshape(lonlatdep, (self.points_num, columns))
5187	#print 'lonlatdep',lonlatdep
5188	self.lonlatdep = lonlatdep
5189
5190	self.mux_file = mux_file
5191	# check this array
5192
5193	def __iter__(self):
5194	"""
5195	iterate over quantity data which is with respect to time.
5196
5197	Note: You can only interate once over an object
5198
5199	returns quantity infomation for each time slice
5200	"""
5201	msg = "You can only interate once over a urs file."
5202	assert not self.iterated, msg
5203	self.iter_time_step = 0
5204	self.iterated = True
5205	return self
5206
5207	def next(self):
5208	if self.time_step_count == self.iter_time_step:
5209	self.close()
5210	raise StopIteration
5211	#Read in a time slice from mux file
5212	hz_p_array = p_array.array('f')
5213	hz_p_array.read(self.mux_file, self.points_num)
5214	hz_p = array(hz_p_array, typecode=Float)
5215	self.iter_time_step += 1
5216
5217	return hz_p
5218
5219	def close(self):
5220	self.mux_file.close()
5221
5222	#### END URS UNGRIDDED 2 SWW ###
5223
5224
5225	def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
5226	print_to_screen=False, verbose=False):
5227	"""
5228	Used to store screen output and errors to file, if run on multiple
5229	processes eachprocessor will have its own output and error file.
5230
5231	extra_info - is used as a string that can identify outputs with another
5232	string eg. '_other'
5233	"""
5234	import sys
5235	dir_name = dir_name
5236	if access(dir_name,W_OK) == 0:
5237	if verbose: print 'Make directory %s' %dir_name
5238	if verbose: print "myid", myid
5239	mkdir (dir_name,0777)
5240	if myid <>'':
5241	myid = '_'+str(myid)
5242	if numprocs <>'':
5243	numprocs = '_'+str(numprocs)
5244	if extra_info <>'':
5245	extra_info = '_'+str(extra_info)
5246	screen_output_name = dir_name + "screen_output%s%s%s.txt" %(myid,
5247	numprocs,
5248	extra_info)
5249	screen_error_name = dir_name + "screen_error%s%s%s.txt" %(myid,
5250	numprocs,
5251	extra_info)
5252	print screen_output_name
5253	#used to catch screen output to file
5254	sys.stdout = Screen_Catcher(screen_output_name)
5255	sys.stderr = Screen_Catcher(screen_error_name)
5256
5257	class Screen_Catcher:
5258	"""this simply catches the screen output and stores it to file defined by
5259	start_screen_catcher (above)
5260	"""
5261
5262	def __init__(self, filename):
5263	self.filename = filename
5264
5265	if exists(self.filename)is True:
5266	print'Old existing file "%s" has been deleted' %(self.filename)
5267	remove(self.filename)
5268
5269	def write(self, stuff):
5270	fid = open(self.filename, 'a')
5271	fid.write(stuff)
5272	# if print_to_screen: print stuff
5273
5274	def copy_code_files(dir_name, filename1, filename2):
5275	"""Copies "filename1" and "filename2" to "dir_name". Very useful for
5276	information management
5277	filename1 and filename2 are both absolute pathnames
5278	"""
5279
5280	if access(dir_name,F_OK) == 0:
5281	print 'Make directory %s' %dir_name
5282	mkdir (dir_name,0777)
5283	shutil.copy(filename1, dir_name + sep + basename(filename1))
5284	shutil.copy(filename2, dir_name + sep + basename(filename2))
5285	# copy (__file__, project.output_run_time_dir + basename(__file__))
5286	print 'Files %s and %s copied' %(filename1, filename2)
5287
5288	def get_data_from_file(filename,separator_value = ','):
5289	"""
5290	Read in data information from file and
5291
5292	Returns:
5293	header_fields, a string? of the first line separated
5294	by the 'separator_value'
5295
5296	data, a array (N data columns X M lines) in the file
5297	excluding the header
5298
5299	NOTE: wont deal with columns with different lenghts and there must be
5300	no blank lines at the end.
5301	"""
5302
5303	fid = open(filename)
5304	lines = fid.readlines()
5305
5306	fid.close()
5307
5308	header_line = lines[0]
5309	header_fields = header_line.split(separator_value)
5310
5311	#array to store data, number in there is to allow float...
5312	#i'm sure there is a better way!
5313	data=array([],typecode=Float)
5314	data=resize(data,((len(lines)-1),len(header_fields)))
5315	# print 'number of fields',range(len(header_fields))
5316	# print 'number of lines',len(lines), shape(data)
5317	# print'data',data[1,1],header_line
5318
5319	array_number = 0
5320	line_number = 1
5321	while line_number < (len(lines)):
5322	for i in range(len(header_fields)):
5323	#this get line below the header, explaining the +1
5324	#and also the line_number can be used as the array index
5325	fields = lines[line_number].split(separator_value)
5326	#assign to array
5327	data[array_number,i] = float(fields[i])
5328
5329	line_number = line_number +1
5330	array_number = array_number +1
5331
5332	return header_fields, data
5333
5334	def store_parameters(verbose=False,**kwargs):
5335	"""
5336	Store "kwargs" into a temp csv file, if "completed" is a kwargs csv file is
5337	kwargs[file_name] else it is kwargs[output_dir] + details_temp.csv
5338
5339	Must have a file_name keyword arg, this is what is writing to.
5340	might be a better way to do this using CSV module Writer and writeDict
5341
5342	writes file to "output_dir" unless "completed" is in kwargs, then
5343	it writes to "file_name" kwargs
5344
5345	"""
5346	import types
5347	# import os
5348
5349	# Check that kwargs is a dictionary
5350	if type(kwargs) != types.DictType:
5351	raise TypeError
5352
5353	#is completed is kwargs?
5354	try:
5355	kwargs['completed']
5356	completed=True
5357	except:
5358	completed=False
5359
5360	#get file name and removes from dict and assert that a file_name exists
5361	if completed:
5362	try:
5363	file = str(kwargs['file_name'])
5364	except:
5365	raise 'kwargs must have file_name'
5366	else:
5367	#write temp file in output directory
5368	try:
5369	file = str(kwargs['output_dir'])+'detail_temp.csv'
5370	except:
5371	raise 'kwargs must have output_dir'
5372
5373	#extracts the header info and the new line info
5374	line=''
5375	header=''
5376	count=0
5377	keys = kwargs.keys()
5378	keys.sort()
5379
5380	#used the sorted keys to create the header and line data
5381	for k in keys:
5382	# print "%s = %s" %(k, kwargs[k])
5383	header = header+str(k)
5384	line = line+str(kwargs[k])
5385	count+=1
5386	if count <len(kwargs):
5387	header = header+','
5388	line = line+','
5389	header+='\n'
5390	line+='\n'
5391
5392	# checks the header info, if the same, then write, if not create a new file
5393	#try to open!
5394	try:
5395	fid = open(file,"r")
5396	file_header=fid.readline()
5397	fid.close()
5398	if verbose: print 'read file header %s' %file_header
5399
5400	except:
5401	msg = 'try to create new file',file
5402	if verbose: print msg
5403	#tries to open file, maybe directory is bad
5404	try:
5405	fid = open(file,"w")
5406	fid.write(header)
5407	fid.close()
5408	file_header=header
5409	except:
5410	msg = 'cannot create new file',file
5411	raise msg
5412
5413	#if header is same or this is a new file
5414	if file_header==str(header):
5415	fid=open(file,"a")
5416	#write new line
5417	fid.write(line)
5418	fid.close()
5419	else:
5420	#backup plan,
5421	# if header is different and has completed will append info to
5422	#end of details_temp.cvs file in output directory
5423	file = str(kwargs['output_dir'])+'detail_temp.csv'
5424	fid=open(file,"a")
5425	fid.write(header)
5426	fid.write(line)
5427	fid.close()
5428	if verbose: print 'file',file_header.strip('\n')
5429	if verbose: print 'head',header.strip('\n')
5430	if file_header.strip('\n')==str(header): print 'they equal'
5431	msg = 'WARNING: File header does not match input info, the input variables have changed, suggest to change file name'
5432	print msg
5433
5434
5435
5436
5437
5438	def get_maximum_inundation_elevation(filename,
5439	polygon=None,
5440	time_interval=None,
5441	verbose=False):
5442
5443	"""Return highest elevation where depth > 0
5444
5445	Usage:
5446	max_runup = get_maximum_inundation_elevation(filename,
5447	polygon=None,
5448	time_interval=None,
5449	verbose=False)
5450
5451	filename is a NetCDF sww file containing ANUGA model output.
5452	Optional arguments polygon and time_interval restricts the maximum
5453	runup calculation
5454	to a points that lie within the specified polygon and time interval.
5455
5456	If no inundation is found within polygon and time_interval the return value
5457	is None signifying "No Runup" or "Everything is dry".
5458
5459	See general function get_maximum_inundation_data for details.
5460
5461	"""
5462
5463	runup, _ = get_maximum_inundation_data(filename,
5464	polygon=polygon,
5465	time_interval=time_interval,
5466	verbose=verbose)
5467	return runup
5468
5469
5470
5471
5472	def get_maximum_inundation_location(filename,
5473	polygon=None,
5474	time_interval=None,
5475	verbose=False):
5476	"""Return location of highest elevation where h > 0
5477
5478
5479	Usage:
5480	max_runup_location = get_maximum_inundation_location(filename,
5481	polygon=None,
5482	time_interval=None,
5483	verbose=False)
5484
5485	filename is a NetCDF sww file containing ANUGA model output.
5486	Optional arguments polygon and time_interval restricts the maximum
5487	runup calculation
5488	to a points that lie within the specified polygon and time interval.
5489
5490	If no inundation is found within polygon and time_interval the return value
5491	is None signifying "No Runup" or "Everything is dry".
5492
5493	See general function get_maximum_inundation_data for details.
5494	"""
5495
5496	_, max_loc = get_maximum_inundation_data(filename,
5497	polygon=polygon,
5498	time_interval=time_interval,
5499	verbose=verbose)
5500	return max_loc
5501
5502
5503
5504	def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
5505	use_centroid_values=False,
5506	verbose=False):
5507	"""Compute maximum run up height from sww file.
5508
5509
5510	Usage:
5511	runup, location = get_maximum_inundation_data(filename,
5512	polygon=None,
5513	time_interval=None,
5514	verbose=False)
5515
5516
5517	Algorithm is as in get_maximum_inundation_elevation from
5518	shallow_water_domain
5519	except that this function works with the sww file and computes the maximal
5520	runup height over multiple timesteps.
5521
5522	Optional arguments polygon and time_interval restricts the
5523	maximum runup calculation
5524	to a points that lie within the specified polygon and time interval.
5525	Polygon is
5526	assumed to be in (absolute) UTM coordinates in the same zone as domain.
5527
5528	If no inundation is found within polygon and time_interval the return value
5529	is None signifying "No Runup" or "Everything is dry".
5530	"""
5531
5532	# We are using nodal values here as that is what is stored in sww files.
5533
5534	# Water depth below which it is considered to be 0 in the model
5535	# FIXME (Ole): Allow this to be specified as a keyword argument as well
5536
5537	from anuga.utilities.polygon import inside_polygon
5538	from anuga.config import minimum_allowed_height
5539	from Scientific.IO.NetCDF import NetCDFFile
5540
5541	dir, base = os.path.split(filename)
5542
5543	iterate_over = get_all_swwfiles(dir,base)
5544
5545	# Read sww file
5546	if verbose:
5547	print 'Reading from %s' %filename
5548	# FIXME: Use general swwstats (when done)
5549
5550	maximal_runup = None
5551	maximal_runup_location = None
5552
5553	for file, swwfile in enumerate (iterate_over):
5554
5555	# Read sww file
5556	filename = join(dir,swwfile+'.sww')
5557
5558	if verbose:
5559	print 'Reading from %s' %filename
5560	# FIXME: Use general swwstats (when done)
5561
5562	fid = NetCDFFile(filename)
5563
5564	# Get geo_reference
5565	# sww files don't have to have a geo_ref
5566	try:
5567	geo_reference = Geo_reference(NetCDFObject=fid)
5568	except AttributeError, e:
5569	geo_reference = Geo_reference() # Default georef object
5570
5571	xllcorner = geo_reference.get_xllcorner()
5572	yllcorner = geo_reference.get_yllcorner()
5573	zone = geo_reference.get_zone()
5574
5575	# Get extent
5576	volumes = fid.variables['volumes'][:]
5577	x = fid.variables['x'][:] + xllcorner
5578	y = fid.variables['y'][:] + yllcorner
5579
5580
5581	# Get the relevant quantities
5582	elevation = fid.variables['elevation'][:]
5583	stage = fid.variables['stage'][:]
5584
5585
5586	# Here's where one could convert nodal information to centroid
5587	# information
5588	# but is probably something we need to write in C.
5589	# Here's a Python thought which is NOT finished!!!
5590	if use_centroid_values is True:
5591	x = get_centroid_values(x, volumes)
5592	y = get_centroid_values(y, volumes)
5593	elevation = get_centroid_values(elevation, volumes)
5594
5595
5596	# Spatial restriction
5597	if polygon is not None:
5598	msg = 'polygon must be a sequence of points.'
5599	assert len(polygon[0]) == 2, msg
5600	# FIXME (Ole): Make a generic polygon input check in polygon.py
5601	# and call it here
5602
5603	points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1)
5604
5605	point_indices = inside_polygon(points, polygon)
5606
5607	# Restrict quantities to polygon
5608	elevation = take(elevation, point_indices)
5609	stage = take(stage, point_indices, axis=1)
5610
5611	# Get info for location of maximal runup
5612	points_in_polygon = take(points, point_indices)
5613	x = points_in_polygon[:,0]
5614	y = points_in_polygon[:,1]
5615	else:
5616	# Take all points
5617	point_indices = arange(len(x))
5618
5619
5620	# Temporal restriction
5621	time = fid.variables['time'][:]
5622	all_timeindices = arange(len(time))
5623	if time_interval is not None:
5624
5625	msg = 'time_interval must be a sequence of length 2.'
5626	assert len(time_interval) == 2, msg
5627	msg = 'time_interval %s must not be decreasing.' %(time_interval)
5628	assert time_interval[1] >= time_interval[0], msg
5629
5630	msg = 'Specified time interval [%.8f:%.8f]' %tuple(time_interval)
5631	msg += ' must does not match model time interval: [%.8f, %.8f]\n'\
5632	%(time[0], time[-1])
5633	if time_interval[1] < time[0]: raise ValueError(msg)
5634	if time_interval[0] > time[-1]: raise ValueError(msg)
5635
5636	# Take time indices corresponding to interval (& is bitwise AND)
5637	timesteps = compress((time_interval[0] <= time) & (time <= time_interval[1]),
5638	all_timeindices)
5639
5640
5641	msg = 'time_interval %s did not include any model timesteps.' %(time_interval)
5642	assert not alltrue(timesteps == 0), msg
5643
5644
5645	else:
5646	# Take them all
5647	timesteps = all_timeindices
5648
5649
5650	fid.close()
5651
5652	# Compute maximal runup for each timestep
5653	#maximal_runup = None
5654	#maximal_runup_location = None
5655	#maximal_runups = [None]
5656	#maximal_runup_locations = [None]
5657
5658	for i in timesteps:
5659	if use_centroid_values is True:
5660	stage_i = get_centroid_values(stage[i,:], volumes)
5661	else:
5662	stage_i = stage[i,:]
5663
5664	depth = stage_i - elevation
5665
5666	# Get wet nodes i.e. nodes with depth>0 within given region and timesteps
5667	wet_nodes = compress(depth > minimum_allowed_height, arange(len(depth)))
5668
5669	if alltrue(wet_nodes == 0):
5670	runup = None
5671	else:
5672	# Find maximum elevation among wet nodes
5673	wet_elevation = take(elevation, wet_nodes)
5674
5675	runup_index = argmax(wet_elevation)
5676	runup = max(wet_elevation)
5677	assert wet_elevation[runup_index] == runup # Must always be True
5678	if runup > maximal_runup:
5679	maximal_runup = runup # This works even if maximal_runups is None
5680	#print "NEW RUNUP",runup
5681
5682	# Record location
5683	wet_x = take(x, wet_nodes)
5684	wet_y = take(y, wet_nodes)
5685	maximal_runup_location = [wet_x[runup_index], wet_y[runup_index]]
5686
5687	#print 'maximal_runup, maximal_runup_location',maximal_runup, maximal_runup_location
5688	return maximal_runup, maximal_runup_location
5689
5690	def get_all_swwfiles(look_in_dir='',base_name='',verbose=False):
5691	'''
5692	Finds all the sww files in a "look_in_dir" which contains a "base_name".
5693	will accept base_name with or without the extension ".sww"
5694
5695	Returns: a list of strings
5696
5697	Usage: iterate_over = get_all_swwfiles(dir, name)
5698	then
5699	for swwfile in iterate_over:
5700	do stuff
5701
5702	Check "export_grids" and "get_maximum_inundation_data" for examples
5703	'''
5704
5705	#plus tests the extension
5706	name, extension = os.path.splitext(base_name)
5707
5708	if extension <>'' and extension <> '.sww':
5709	msg = msg='file %s %s must be an NetCDF sww file!'%(base_name,extension)
5710	raise IOError, msg
5711
5712	if look_in_dir == "":
5713	look_in_dir = "." # Unix compatibility
5714
5715	dir_ls = os.listdir(look_in_dir)
5716	#print 'dir_ls',dir_ls, base
5717	iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == '.sww']
5718	if len(iterate_over) == 0:
5719	msg = 'No files of the base name %s'\
5720	%(name)
5721	raise IOError, msg
5722	if verbose: print 'iterate over %s' %(iterate_over)
5723
5724	return iterate_over
5725
5726	def get_all_files_with_extension(look_in_dir='',base_name='',extension='.sww',verbose=False):
5727	'''
5728	Finds all the sww files in a "look_in_dir" which contains a "base_name".
5729
5730
5731	Returns: a list of strings
5732
5733	Usage: iterate_over = get_all_swwfiles(dir, name)
5734	then
5735	for swwfile in iterate_over:
5736	do stuff
5737
5738	Check "export_grids" and "get_maximum_inundation_data" for examples
5739	'''
5740
5741	#plus tests the extension
5742	name, ext = os.path.splitext(base_name)
5743	# print 'look_in_dir',look_in_dir
5744
5745	if ext <>'' and ext <> extension:
5746	msg = msg='base_name %s must be an file with %s extension!'%(base_name,extension)
5747	raise IOError, msg
5748
5749	if look_in_dir == "":
5750	look_in_dir = "." # Unix compatibility
5751	# print 'look_in_dir',look_in_dir, getcwd()
5752	dir_ls = os.listdir(look_in_dir)
5753	#print 'dir_ls',dir_ls, base_name
5754	iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == extension]
5755	if len(iterate_over) == 0:
5756	msg = 'No files of the base name %s in %s'\
5757	%(name, look_in_dir)
5758	raise IOError, msg
5759	if verbose: print 'iterate over %s' %(iterate_over)
5760
5761	return iterate_over
5762
5763	def get_all_directories_with_name(look_in_dir='',base_name='',verbose=False):
5764	'''
5765	Finds all the sww files in a "look_in_dir" which contains a "base_name".
5766
5767
5768	Returns: a list of strings
5769
5770	Usage: iterate_over = get_all_swwfiles(dir, name)
5771	then
5772	for swwfile in iterate_over:
5773	do stuff
5774
5775	Check "export_grids" and "get_maximum_inundation_data" for examples
5776	'''
5777
5778	#plus tests the extension
5779
5780	if look_in_dir == "":
5781	look_in_dir = "." # Unix compatibility
5782	# print 'look_in_dir',look_in_dir
5783	dir_ls = os.listdir(look_in_dir)
5784	# print 'dir_ls',dir_ls
5785	iterate_over = [x for x in dir_ls if base_name in x]
5786	if len(iterate_over) == 0:
5787	msg = 'No files of the base name %s'\
5788	%(name)
5789	raise IOError, msg
5790	if verbose: print 'iterate over %s' %(iterate_over)
5791
5792	return iterate_over
5793
5794
5795	#-------------------------------------------------------------
5796	if __name__ == "__main__":
5797	#setting umask from config to force permissions for all files and directories
5798	# created to the same. (it was noticed the "mpirun" doesn't honour the umask
5799	# set in your .bashrc etc file)
5800	from config import umask
5801	import os
5802	os.umask(umask)
5803

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

Download in other formats: