Context Navigation

source: inundation/ga/storm_surge/pyvolution/data_manager.py @ 1007

Last change on this file since 1007 was 995, checked in by ole, 20 years ago
Cast Integer array explictly as 32 bit for use with Netcdf (I think)
File size: 43.1 KB

Line
1	"""Functions to store and retrieve data for the Shallow Water Wave equation.
2	There are two kinds of data
3
4	1: Constant data: Vertex coordinates and field values. Stored once
5	2: Variable data: Conserved quantities. Stored once per timestep.
6
7	All data is assumed to reside at vertex locations.
8
9
10	Formats used within ANUGA:
11
12	.sww: Netcdf format for storing model output
13
14	.xya: ASCII format for storing arbitrary points and associated attributes
15	.pts: NetCDF format for storing arbitrary points and associated attributes
16
17	.asc: ASCII foramt of regular DEMs as output from ArcView
18	.dem: NetCDF representation of regular DEM data
19
20	.tsh: ASCII format for storing meshes and associated boundary and region info
21
22	.nc: Native ferret NetCDF format
23
24	A typical dataflow can be described as follows
25
26	Manually created files:
27	ASC, PRJ: Digital elevation models (gridded)
28	TSH: Triangular meshes (e.g. dreated from pmesh)
29	NC Model outputs for use as boundary conditions (e.g from MOST)
30
31
32	AUTOMATICALLY CREATED FILES:
33
34	ASC, PRJ -> DEM -> PTS: Conversion of DEM's to native pts file
35
36	NC -> SWW: Conversion of MOST bundary files to boundary sww
37
38	PTS + TSH -> TSH with elevation: Least squares fit
39
40	TSH -> SWW: Conversion of TSH to sww viewable using Swollen
41
42	TSH + Boundary SWW -> SWW: SImluation using pyvolution
43
44
45	"""
46
47	from Numeric import concatenate
48
49
50	def make_filename(s):
51	"""Transform argument string into a suitable filename
52	"""
53
54	s = s.strip()
55	s = s.replace(' ', '_')
56	s = s.replace('(', '')
57	s = s.replace(')', '')
58	s = s.replace('__', '_')
59
60	return s
61
62
63	def check_dir(path, verbose=None):
64	"""Check that specified path exists.
65	If path does not exist it will be created if possible
66
67	USAGE:
68	checkdir(path, verbose):
69
70	ARGUMENTS:
71	path -- Directory
72	verbose -- Flag verbose output (default: None)
73
74	RETURN VALUE:
75	Verified path including trailing separator
76
77	"""
78
79	import os, sys
80	import os.path
81
82	if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
83	unix = 0
84	else:
85	unix = 1
86
87
88	if path[-1] != os.sep:
89	path = path + os.sep # Add separator for directories
90
91	path = os.path.expanduser(path) # Expand ~ or ~user in pathname
92	if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
93	try:
94	exitcode=os.mkdir(path)
95
96	# Change access rights if possible
97	#
98	if unix:
99	exitcode=os.system('chmod 775 '+path)
100	else:
101	pass # FIXME: What about acces rights under Windows?
102
103	if verbose: print 'MESSAGE: Directory', path, 'created.'
104
105	except:
106	print 'WARNING: Directory', path, 'could not be created.'
107	if unix:
108	path = '/tmp/'
109	else:
110	path = 'C:'
111
112	print 'Using directory %s instead' %path
113
114	return(path)
115
116
117
118	def del_dir(path):
119	"""Recursively delete directory path and all its contents
120	"""
121
122	import os
123
124	if os.path.isdir(path):
125	for file in os.listdir(path):
126	X = os.path.join(path, file)
127
128
129	if os.path.isdir(X) and not os.path.islink(X):
130	del_dir(X)
131	else:
132	try:
133	os.remove(X)
134	except:
135	print "Could not remove file %s" %X
136
137	os.rmdir(path)
138
139
140
141	def create_filename(datadir, filename, format, size=None, time=None):
142
143	import os
144	#from config import data_dir
145
146	FN = check_dir(datadir) + filename
147
148	if size is not None:
149	FN += '_size%d' %size
150
151	if time is not None:
152	FN += '_time%.2f' %time
153
154	FN += '.' + format
155	return FN
156
157
158	def get_files(datadir, filename, format, size):
159	"""Get all file (names) with gven name, size and format
160	"""
161
162	import glob
163
164	import os
165	#from config import data_dir
166
167	dir = check_dir(datadir)
168
169	pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
170	return glob.glob(pattern)
171
172
173
174	#Generic class for storing output to e.g. visualisation or checkpointing
175	class Data_format:
176	"""Generic interface to data formats
177	"""
178
179
180	def __init__(self, domain, extension, mode = 'w'):
181	assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
182	''' 'w' (write)'''+\
183	''' 'r' (read)''' +\
184	''' 'a' (append)'''
185
186	#Create filename
187	#self.filename = create_filename(domain.get_datadir(),
188	# domain.get_name(), extension, len(domain))
189
190
191	self.filename = create_filename(domain.get_datadir(),
192	domain.get_name(), extension)
193
194	#print 'F', self.filename
195	self.timestep = 0
196	self.number_of_volumes = len(domain)
197	self.domain = domain
198
199
200	#FIXME: Should we have a general set_precision function?
201
202
203
204	#Class for storing output to e.g. visualisation
205	class Data_format_sww(Data_format):
206	"""Interface to native NetCDF format (.sww)
207	"""
208
209
210	def __init__(self, domain, mode = 'w'):
211	from Scientific.IO.NetCDF import NetCDFFile
212	from Numeric import Int, Float, Float32
213
214	self.precision = Float32 #Use single precision
215
216	Data_format.__init__(self, domain, 'sww', mode)
217
218
219	# NetCDF file definition
220	fid = NetCDFFile(self.filename, mode)
221
222	if mode == 'w':
223
224	#Create new file
225	fid.institution = 'Geoscience Australia'
226	fid.description = 'Output from pyvolution suitable for plotting'
227
228	if domain.smooth:
229	fid.smoothing = 'Yes'
230	else:
231	fid.smoothing = 'No'
232
233	fid.order = domain.default_order
234
235	#Reference point
236	#Start time in seconds since the epoch (midnight 1/1/1970)
237	fid.starttime = domain.starttime
238	fid.xllcorner = domain.xllcorner
239	fid.yllcorner = domain.yllcorner
240	fid.zone = str(domain.zone) #FIXME: ?
241
242	# dimension definitions
243	fid.createDimension('number_of_volumes', self.number_of_volumes)
244	fid.createDimension('number_of_vertices', 3)
245
246	if domain.smooth is True:
247	fid.createDimension('number_of_points', len(domain.vertexlist))
248	else:
249	fid.createDimension('number_of_points', 3*self.number_of_volumes)
250
251	fid.createDimension('number_of_timesteps', None) #extensible
252
253	# variable definitions
254	fid.createVariable('x', self.precision, ('number_of_points',))
255	fid.createVariable('y', self.precision, ('number_of_points',))
256	fid.createVariable('elevation', self.precision, ('number_of_points',))
257
258	#FIXME: Backwards compatibility
259	fid.createVariable('z', self.precision, ('number_of_points',))
260	#################################
261
262	fid.createVariable('volumes', Int, ('number_of_volumes',
263	'number_of_vertices'))
264
265	fid.createVariable('time', self.precision,
266	('number_of_timesteps',))
267
268	fid.createVariable('stage', self.precision,
269	('number_of_timesteps',
270	'number_of_points'))
271
272	fid.createVariable('xmomentum', self.precision,
273	('number_of_timesteps',
274	'number_of_points'))
275
276	fid.createVariable('ymomentum', self.precision,
277	('number_of_timesteps',
278	'number_of_points'))
279
280	#Close
281	fid.close()
282
283
284	def store_connectivity(self):
285	"""Specialisation of store_connectivity for net CDF format
286
287	Writes x,y,z coordinates of triangles constituting
288	the bed elevation.
289	"""
290
291	from Scientific.IO.NetCDF import NetCDFFile
292
293	from Numeric import concatenate, Int
294
295	domain = self.domain
296
297	#Get NetCDF
298	fid = NetCDFFile(self.filename, 'a') #Open existing file for append
299
300	# Get the variables
301	x = fid.variables['x']
302	y = fid.variables['y']
303	z = fid.variables['elevation']
304
305	volumes = fid.variables['volumes']
306
307	# Get X, Y and bed elevation Z
308	Q = domain.quantities['elevation']
309	X,Y,Z,V = Q.get_vertex_values(xy=True,
310	precision = self.precision)
311
312
313
314	x[:] = X.astype(self.precision)
315	y[:] = Y.astype(self.precision)
316	z[:] = Z.astype(self.precision)
317
318	#FIXME: Backwards compatibility
319	z = fid.variables['z']
320	z[:] = Z.astype(self.precision)
321	################################
322
323	volumes[:] = V.astype(volumes.typecode())
324
325	#Close
326	fid.close()
327
328
329	def store_timestep(self, names):
330	"""Store time and named quantities to file
331	"""
332	from Scientific.IO.NetCDF import NetCDFFile
333	import types
334	from time import sleep
335
336
337	#Get NetCDF
338	retries = 0
339	file_open = False
340	while not file_open and retries < 10:
341	try:
342	fid = NetCDFFile(self.filename, 'a') #Open existing file
343	except IOError:
344	#This could happen if someone was reading the file.
345	#In that case, wait a while and try again
346	msg = 'Warning (store_timestep): File %s could not be opened'\
347	%self.filename
348	msg += ' - trying again'
349	print msg
350	retries += 1
351	sleep(1)
352	else:
353	file_open = True
354
355	if not file_open:
356	msg = 'File %s could not be opened for append' %self.filename
357	raise msg
358
359
360	domain = self.domain
361
362	# Get the variables
363	time = fid.variables['time']
364	stage = fid.variables['stage']
365	xmomentum = fid.variables['xmomentum']
366	ymomentum = fid.variables['ymomentum']
367	i = len(time)
368
369	#Store time
370	time[i] = self.domain.time
371
372
373	if type(names) not in [types.ListType, types.TupleType]:
374	names = [names]
375
376	for name in names:
377	# Get quantity
378	Q = domain.quantities[name]
379	A,V = Q.get_vertex_values(xy=False,
380	precision = self.precision)
381
382	#FIXME: Make this general (see below)
383	if name == 'stage':
384	stage[i,:] = A.astype(self.precision)
385	elif name == 'xmomentum':
386	xmomentum[i,:] = A.astype(self.precision)
387	elif name == 'ymomentum':
388	ymomentum[i,:] = A.astype(self.precision)
389
390	#As in....
391	#eval( name + '[i,:] = A.astype(self.precision)' )
392	#FIXME: But we need a UNIT test for that before refactoring
393
394
395
396	#Flush and close
397	fid.sync()
398	fid.close()
399
400
401	#Class for handling checkpoints data
402	class Data_format_cpt(Data_format):
403	"""Interface to native NetCDF format (.cpt)
404	"""
405
406
407	def __init__(self, domain, mode = 'w'):
408	from Scientific.IO.NetCDF import NetCDFFile
409	from Numeric import Int, Float, Float
410
411	self.precision = Float #Use full precision
412
413	Data_format.__init__(self, domain, 'sww', mode)
414
415
416	# NetCDF file definition
417	fid = NetCDFFile(self.filename, mode)
418
419	if mode == 'w':
420	#Create new file
421	fid.institution = 'Geoscience Australia'
422	fid.description = 'Checkpoint data'
423	#fid.smooth = domain.smooth
424	fid.order = domain.default_order
425
426	# dimension definitions
427	fid.createDimension('number_of_volumes', self.number_of_volumes)
428	fid.createDimension('number_of_vertices', 3)
429
430	#Store info at all vertices (no smoothing)
431	fid.createDimension('number_of_points', 3*self.number_of_volumes)
432	fid.createDimension('number_of_timesteps', None) #extensible
433
434	# variable definitions
435
436	#Mesh
437	fid.createVariable('x', self.precision, ('number_of_points',))
438	fid.createVariable('y', self.precision, ('number_of_points',))
439
440
441	fid.createVariable('volumes', Int, ('number_of_volumes',
442	'number_of_vertices'))
443
444	fid.createVariable('time', self.precision,
445	('number_of_timesteps',))
446
447	#Allocate space for all quantities
448	for name in domain.quantities.keys():
449	fid.createVariable(name, self.precision,
450	('number_of_timesteps',
451	'number_of_points'))
452
453	#Close
454	fid.close()
455
456
457	def store_checkpoint(self):
458	"""
459	Write x,y coordinates of triangles.
460	Write connectivity (
461	constituting
462	the bed elevation.
463	"""
464
465	from Scientific.IO.NetCDF import NetCDFFile
466
467	from Numeric import concatenate
468
469	domain = self.domain
470
471	#Get NetCDF
472	fid = NetCDFFile(self.filename, 'a') #Open existing file for append
473
474	# Get the variables
475	x = fid.variables['x']
476	y = fid.variables['y']
477
478	volumes = fid.variables['volumes']
479
480	# Get X, Y and bed elevation Z
481	Q = domain.quantities['elevation']
482	X,Y,Z,V = Q.get_vertex_values(xy=True,
483	precision = self.precision)
484
485
486
487	x[:] = X.astype(self.precision)
488	y[:] = Y.astype(self.precision)
489	z[:] = Z.astype(self.precision)
490
491	volumes[:] = V
492
493	#Close
494	fid.close()
495
496
497	def store_timestep(self, name):
498	"""Store time and named quantity to file
499	"""
500	from Scientific.IO.NetCDF import NetCDFFile
501	from time import sleep
502
503	#Get NetCDF
504	retries = 0
505	file_open = False
506	while not file_open and retries < 10:
507	try:
508	fid = NetCDFFile(self.filename, 'a') #Open existing file
509	except IOError:
510	#This could happen if someone was reading the file.
511	#In that case, wait a while and try again
512	msg = 'Warning (store_timestep): File %s could not be opened'\
513	%self.filename
514	msg += ' - trying again'
515	print msg
516	retries += 1
517	sleep(1)
518	else:
519	file_open = True
520
521	if not file_open:
522	msg = 'File %s could not be opened for append' %self.filename
523	raise msg
524
525
526	domain = self.domain
527
528	# Get the variables
529	time = fid.variables['time']
530	stage = fid.variables['stage']
531	i = len(time)
532
533	#Store stage
534	time[i] = self.domain.time
535
536	# Get quantity
537	Q = domain.quantities[name]
538	A,V = Q.get_vertex_values(xy=False,
539	precision = self.precision)
540
541	stage[i,:] = A.astype(self.precision)
542
543	#Flush and close
544	fid.sync()
545	fid.close()
546
547
548
549
550
551
552
553	#Function for storing xya output
554	#FIXME Not done yet for this version
555	class Data_format_xya(Data_format):
556	"""Generic interface to data formats
557	"""
558
559
560	def __init__(self, domain, mode = 'w'):
561	from Scientific.IO.NetCDF import NetCDFFile
562	from Numeric import Int, Float, Float32
563
564	self.precision = Float32 #Use single precision
565
566	Data_format.__init__(self, domain, 'xya', mode)
567
568
569
570	#FIXME -This is the old xya format
571	def store_all(self):
572	"""Specialisation of store all for xya format
573
574	Writes x,y,z coordinates of triangles constituting
575	the bed elevation.
576	"""
577
578	from Numeric import concatenate
579
580	domain = self.domain
581
582	fd = open(self.filename, 'w')
583
584
585	if domain.smooth is True:
586	number_of_points = len(domain.vertexlist)
587	else:
588	number_of_points = 3*self.number_of_volumes
589
590	numVertAttrib = 3 #Three attributes is what is assumed by the xya format
591
592	fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
593	" # <vertex #> <x> <y> [attributes]" + "\n")
594
595
596	# Get X, Y, bed elevation and friction (index=0,1)
597	X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
598	indices = (0,1), precision = self.precision)
599
600	bed_eles = A[:,0]
601	fricts = A[:,1]
602
603	# Get stage (index=0)
604	B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
605	indices = (0,), precision = self.precision)
606
607	stages = B[:,0]
608
609	#<vertex #> <x> <y> [attributes]
610	for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
611	stages, fricts):
612
613	s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
614	fd.write(s)
615
616	#close
617	fd.close()
618
619
620	def store_timestep(self, t, V0, V1, V2):
621	"""Store time, water heights (and momentums) to file
622	"""
623	pass
624
625
626	#Auxiliary
627	def write_obj(filename,x,y,z):
628	"""Store x,y,z vectors into filename (obj format)
629	Vectors are assumed to have dimension (M,3) where
630	M corresponds to the number elements.
631	triangles are assumed to be disconnected
632
633	The three numbers in each vector correspond to three vertices,
634
635	e.g. the x coordinate of vertex 1 of element i is in x[i,1]
636
637	"""
638	#print 'Writing obj to %s' % filename
639
640	import os.path
641
642	root, ext = os.path.splitext(filename)
643	if ext == '.obj':
644	FN = filename
645	else:
646	FN = filename + '.obj'
647
648
649	outfile = open(FN, 'wb')
650	outfile.write("# Triangulation as an obj file\n")
651
652	M, N = x.shape
653	assert N==3 #Assuming three vertices per element
654
655	for i in range(M):
656	for j in range(N):
657	outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
658
659	for i in range(M):
660	base = i*N
661	outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
662
663	outfile.close()
664
665
666
667	#Conversion routines
668	def sww2obj(basefilename, size):
669	"""Convert netcdf based data output to obj
670	"""
671	from Scientific.IO.NetCDF import NetCDFFile
672
673	from Numeric import Float, zeros
674
675	#Get NetCDF
676	FN = create_filename('.', basefilename, 'sww', size)
677	print 'Reading from ', FN
678	fid = NetCDFFile(FN, 'r') #Open existing file for read
679
680
681	# Get the variables
682	x = fid.variables['x']
683	y = fid.variables['y']
684	z = fid.variables['elevation']
685	time = fid.variables['time']
686	stage = fid.variables['stage']
687
688	M = size #Number of lines
689	xx = zeros((M,3), Float)
690	yy = zeros((M,3), Float)
691	zz = zeros((M,3), Float)
692
693	for i in range(M):
694	for j in range(3):
695	xx[i,j] = x[i+j*M]
696	yy[i,j] = y[i+j*M]
697	zz[i,j] = z[i+j*M]
698
699	#Write obj for bathymetry
700	FN = create_filename('.', basefilename, 'obj', size)
701	write_obj(FN,xx,yy,zz)
702
703
704	#Now read all the data with variable information, combine with
705	#x,y info and store as obj
706
707	for k in range(len(time)):
708	t = time[k]
709	print 'Processing timestep %f' %t
710
711	for i in range(M):
712	for j in range(3):
713	zz[i,j] = stage[k,i+j*M]
714
715
716	#Write obj for variable data
717	#FN = create_filename(basefilename, 'obj', size, time=t)
718	FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
719	write_obj(FN,xx,yy,zz)
720
721
722	def dat2obj(basefilename):
723	"""Convert line based data output to obj
724	FIXME: Obsolete?
725	"""
726
727	import glob, os
728	from config import data_dir
729
730
731	#Get bathymetry and x,y's
732	lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
733
734	from Numeric import zeros, Float
735
736	M = len(lines) #Number of lines
737	x = zeros((M,3), Float)
738	y = zeros((M,3), Float)
739	z = zeros((M,3), Float)
740
741	##i = 0
742	for i, line in enumerate(lines):
743	tokens = line.split()
744	values = map(float,tokens)
745
746	for j in range(3):
747	x[i,j] = values[j*3]
748	y[i,j] = values[j*3+1]
749	z[i,j] = values[j*3+2]
750
751	##i += 1
752
753
754	#Write obj for bathymetry
755	write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
756
757
758	#Now read all the data files with variable information, combine with
759	#x,y info
760	#and store as obj
761
762	files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
763
764	for filename in files:
765	print 'Processing %s' % filename
766
767	lines = open(data_dir+os.sep+filename,'r').readlines()
768	assert len(lines) == M
769	root, ext = os.path.splitext(filename)
770
771	#Get time from filename
772	i0 = filename.find('_time=')
773	if i0 == -1:
774	#Skip bathymetry file
775	continue
776
777	i0 += 6 #Position where time starts
778	i1 = filename.find('.dat')
779
780	if i1 > i0:
781	t = float(filename[i0:i1])
782	else:
783	raise 'Hmmmm'
784
785
786
787	##i = 0
788	for i, line in enumerate(lines):
789	tokens = line.split()
790	values = map(float,tokens)
791
792	for j in range(3):
793	z[i,j] = values[j]
794
795	##i += 1
796
797	#Write obj for variable data
798	write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
799
800
801	def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
802	"""Read netcdf filename1, pick timesteps first:step:last and save to
803	nettcdf file filename2
804	"""
805	from Scientific.IO.NetCDF import NetCDFFile
806
807	#Get NetCDF
808	infile = NetCDFFile(filename1, 'r') #Open existing file for read
809	outfile = NetCDFFile(filename2, 'w') #Open new file
810
811
812	#Copy dimensions
813	for d in infile.dimensions:
814	outfile.createDimension(d, infile.dimensions[d])
815
816	for name in infile.variables:
817	var = infile.variables[name]
818	outfile.createVariable(name, var.typecode(), var.dimensions)
819
820
821	#Copy the static variables
822	for name in infile.variables:
823	if name == 'time' or name == 'stage':
824	pass
825	else:
826	#Copy
827	outfile.variables[name][:] = infile.variables[name][:]
828
829	#Copy selected timesteps
830	time = infile.variables['time']
831	stage = infile.variables['stage']
832
833	newtime = outfile.variables['time']
834	newstage = outfile.variables['stage']
835
836	if last is None:
837	last = len(time)
838
839	selection = range(first, last, step)
840	for i, j in enumerate(selection):
841	print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
842	newtime[i] = time[j]
843	newstage[i,:] = stage[j,:]
844
845	#Close
846	infile.close()
847	outfile.close()
848
849
850	#Get data objects
851	def get_dataobject(domain, mode='w'):
852	"""Return instance of class of given format using filename
853	"""
854
855	cls = eval('Data_format_%s' %domain.format)
856	return cls(domain, mode)
857
858	def dem2pts(basename_in, basename_out=None, verbose=False):
859	"""Read Digitial Elevation model from the following NetCDF format (.dem)
860
861	Example:
862
863	ncols 3121
864	nrows 1800
865	xllcorner 722000
866	yllcorner 5893000
867	cellsize 25
868	NODATA_value -9999
869	138.3698 137.4194 136.5062 135.5558 ..........
870
871	Convert to NetCDF pts format which is
872
873	points: (Nx2) Float array
874	elevation: N Float array
875	"""
876
877	import os
878	from Scientific.IO.NetCDF import NetCDFFile
879	from Numeric import Float, arrayrange, concatenate
880
881	root = basename_in
882
883	#Get NetCDF
884	infile = NetCDFFile(root + '.dem', 'r') #Open existing netcdf file for read
885
886	if verbose: print 'Reading DEM from %s' %(root + '.dem')
887
888	ncols = infile.ncols[0]
889	nrows = infile.nrows[0]
890	xllcorner = infile.xllcorner[0] #Easting of lower left corner
891	yllcorner = infile.yllcorner[0] #Northing of lower left corner
892	cellsize = infile.cellsize[0]
893	NODATA_value = infile.NODATA_value[0]
894	dem_elevation = infile.variables['elevation']
895
896	zone = infile.zone[0]
897	false_easting = infile.false_easting[0]
898	false_northing = infile.false_northing[0]
899
900	#Text strings
901	projection = infile.projection
902	datum = infile.datum
903	units = infile.units
904
905
906	#Get output file
907	if basename_out == None:
908	ptsname = root + '.pts'
909	else:
910	ptsname = basename_out + '.pts'
911
912	if verbose: print 'Store to NetCDF file %s' %ptsname
913	# NetCDF file definition
914	outfile = NetCDFFile(ptsname, 'w')
915
916	#Create new file
917	outfile.institution = 'Geoscience Australia'
918	outfile.description = 'NetCDF pts format for compact and portable storage ' +\
919	'of spatial point data'
920
921	#Georeferencing
922	outfile.zone = zone
923	outfile.xllcorner = xllcorner #Easting of lower left corner
924	outfile.yllcorner = yllcorner #Northing of lower left corner
925	outfile.false_easting = false_easting
926	outfile.false_northing =false_northing
927
928	outfile.projection = projection
929	outfile.datum = datum
930	outfile.units = units
931
932
933	#Grid info (FIXME: probably not going to be used, but heck)
934	outfile.ncols = ncols
935	outfile.nrows = nrows
936
937
938	# dimension definitions
939	outfile.createDimension('number_of_points', nrows*ncols)
940	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
941
942	# variable definitions
943	outfile.createVariable('points', Float, ('number_of_points',
944	'number_of_dimensions'))
945	outfile.createVariable('elevation', Float, ('number_of_points',))
946
947	# Get handles to the variables
948	points = outfile.variables['points']
949	elevation = outfile.variables['elevation']
950
951	#Store data
952	#FIXME: Could perhaps be faster using array operations
953	for i in range(nrows):
954	if verbose: print 'Processing row %d of %d' %(i, nrows)
955
956	y = (nrows-i)*cellsize
957	for j in range(ncols):
958	index = i*ncols + j
959
960	x = j*cellsize
961	points[index, :] = [x,y]
962	elevation[index] = dem_elevation[i, j]
963
964	infile.close()
965	outfile.close()
966
967
968	def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
969	verbose=False):
970	"""Read Digitial Elevation model from the following ASCII format (.asc)
971
972	Example:
973
974	ncols 3121
975	nrows 1800
976	xllcorner 722000
977	yllcorner 5893000
978	cellsize 25
979	NODATA_value -9999
980	138.3698 137.4194 136.5062 135.5558 ..........
981
982	Convert basename_in + '.asc' to NetCDF format (.dem)
983	mimicking the ASCII format closely.
984
985
986	An accompanying file with same basename_in but extension .prj must exist
987	and is used to fix the UTM zone, datum, false northings and eastings.
988
989	The prj format is assumed to be as
990
991	Projection UTM
992	Zone 56
993	Datum WGS84
994	Zunits NO
995	Units METERS
996	Spheroid WGS84
997	Xshift 0.0000000000
998	Yshift 10000000.0000000000
999	Parameters
1000	"""
1001
1002	import os
1003	from Scientific.IO.NetCDF import NetCDFFile
1004	from Numeric import Float, array
1005
1006	#root, ext = os.path.splitext(basename_in)
1007	root = basename_in
1008
1009	###########################################
1010	# Read Meta data
1011	if verbose: print 'Reading METADATA from %s' %root + '.prj'
1012	metadatafile = open(root + '.prj')
1013	metalines = metadatafile.readlines()
1014	metadatafile.close()
1015
1016	L = metalines[0].strip().split()
1017	assert L[0].strip().lower() == 'projection'
1018	projection = L[1].strip() #TEXT
1019
1020	L = metalines[1].strip().split()
1021	assert L[0].strip().lower() == 'zone'
1022	zone = int(L[1].strip())
1023
1024	L = metalines[2].strip().split()
1025	assert L[0].strip().lower() == 'datum'
1026	datum = L[1].strip() #TEXT
1027
1028	L = metalines[3].strip().split()
1029	assert L[0].strip().lower() == 'zunits' #IGNORE
1030	zunits = L[1].strip() #TEXT
1031
1032	L = metalines[4].strip().split()
1033	assert L[0].strip().lower() == 'units'
1034	units = L[1].strip() #TEXT
1035
1036	L = metalines[5].strip().split()
1037	assert L[0].strip().lower() == 'spheroid' #IGNORE
1038	spheroid = L[1].strip() #TEXT
1039
1040	L = metalines[6].strip().split()
1041	assert L[0].strip().lower() == 'xshift'
1042	false_easting = float(L[1].strip())
1043
1044	L = metalines[7].strip().split()
1045	assert L[0].strip().lower() == 'yshift'
1046	false_northing = float(L[1].strip())
1047
1048	#print false_easting, false_northing, zone, datum
1049
1050
1051	###########################################
1052	#Read DEM data
1053
1054	datafile = open(basename_in + '.asc')
1055
1056	if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
1057	lines = datafile.readlines()
1058	datafile.close()
1059
1060	if verbose: print 'Got', len(lines), ' lines'
1061
1062	ncols = int(lines[0].split()[1].strip())
1063	nrows = int(lines[1].split()[1].strip())
1064	xllcorner = float(lines[2].split()[1].strip())
1065	yllcorner = float(lines[3].split()[1].strip())
1066	cellsize = float(lines[4].split()[1].strip())
1067	NODATA_value = int(lines[5].split()[1].strip())
1068
1069	assert len(lines) == nrows + 6
1070
1071
1072	##########################################
1073
1074
1075	if basename_out == None:
1076	netcdfname = root + '.dem'
1077	else:
1078	netcdfname = basename_out + '.dem'
1079
1080	if verbose: print 'Store to NetCDF file %s' %netcdfname
1081	# NetCDF file definition
1082	fid = NetCDFFile(netcdfname, 'w')
1083
1084	#Create new file
1085	fid.institution = 'Geoscience Australia'
1086	fid.description = 'NetCDF DEM format for compact and portable storage ' +\
1087	'of spatial point data'
1088
1089	fid.ncols = ncols
1090	fid.nrows = nrows
1091	fid.xllcorner = xllcorner
1092	fid.yllcorner = yllcorner
1093	fid.cellsize = cellsize
1094	fid.NODATA_value = NODATA_value
1095
1096	fid.zone = zone
1097	fid.false_easting = false_easting
1098	fid.false_northing = false_northing
1099	fid.projection = projection
1100	fid.datum = datum
1101	fid.units = units
1102
1103
1104	# dimension definitions
1105	fid.createDimension('number_of_rows', nrows)
1106	fid.createDimension('number_of_columns', ncols)
1107
1108	# variable definitions
1109	fid.createVariable('elevation', Float, ('number_of_rows',
1110	'number_of_columns'))
1111
1112	# Get handles to the variables
1113	elevation = fid.variables['elevation']
1114
1115	#Store data
1116	for i, line in enumerate(lines[6:]):
1117	fields = line.split()
1118	if verbose: print 'Processing row %d of %d' %(i, nrows)
1119
1120	elevation[i, :] = array([float(x) for x in fields])
1121
1122	fid.close()
1123
1124
1125
1126	def ferret2sww(basename_in, basename_out = None,
1127	verbose = False,
1128	minlat = None, maxlat = None,
1129	minlon = None, maxlon = None,
1130	mint = None, maxt = None, mean_stage = 0,
1131	origin = None, zscale = 1):
1132	"""Convert 'Ferret' NetCDF format for wave propagation to
1133	sww format native to pyvolution.
1134
1135	Specify only basename_in and read files of the form
1136	basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1137	relative height, x-velocity and y-velocity, respectively.
1138
1139	Also convert latitude and longitude to UTM. All coordinates are
1140	assumed to be given in the GDA94 datum.
1141
1142	min's and max's: If omitted - full extend is used.
1143	To include a value min may equal it, while max must exceed it.
1144	Lat and lon are assuemd to be in decimal degrees
1145
1146	origin is a 3-tuple with geo referenced
1147	UTM coordinates (zone, easting, northing)
1148
1149	nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
1150	which means that longitude is the fastest
1151	varying dimension (row major order, so to speak)
1152
1153	ferret2sww uses grid points as vertices in a triangular grid
1154	counting vertices from lower left corner upwards, then right
1155	"""
1156
1157	import os
1158	from Scientific.IO.NetCDF import NetCDFFile
1159	from Numeric import Float, Int, Int32, searchsorted, zeros, array
1160	precision = Float
1161
1162
1163	#Get NetCDF data
1164	if verbose: print 'Reading files %s_*.nc' %basename_in
1165	file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
1166	file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
1167	file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
1168
1169	if basename_out is None:
1170	swwname = basename_in + '.sww'
1171	else:
1172	swwname = basename_out + '.sww'
1173
1174	times = file_h.variables['TIME']
1175	latitudes = file_h.variables['LAT']
1176	longitudes = file_h.variables['LON']
1177
1178	if mint == None:
1179	jmin = 0
1180	else:
1181	jmin = searchsorted(times, mint)
1182
1183	if maxt == None:
1184	jmax=len(times)
1185	else:
1186	jmax = searchsorted(times, maxt)
1187
1188	if minlat == None:
1189	kmin=0
1190	else:
1191	kmin = searchsorted(latitudes, minlat)
1192
1193	if maxlat == None:
1194	kmax = len(latitudes)
1195	else:
1196	kmax = searchsorted(latitudes, maxlat)
1197
1198	if minlon == None:
1199	lmin=0
1200	else:
1201	lmin = searchsorted(longitudes, minlon)
1202
1203	if maxlon == None:
1204	lmax = len(longitudes)
1205	else:
1206	lmax = searchsorted(longitudes, maxlon)
1207
1208	times = times[jmin:jmax]
1209	latitudes = latitudes[kmin:kmax]
1210	longitudes = longitudes[lmin:lmax]
1211
1212
1213	if verbose: print 'cropping'
1214	amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1215	xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax]
1216	yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]
1217
1218	number_of_times = times.shape[0]
1219	number_of_latitudes = latitudes.shape[0]
1220	number_of_longitudes = longitudes.shape[0]
1221
1222	assert amplitudes.shape[0] == number_of_times
1223	assert amplitudes.shape[1] == number_of_latitudes
1224	assert amplitudes.shape[2] == number_of_longitudes
1225
1226	#print times
1227	#print latitudes
1228	#print longitudes
1229
1230	#print 'MIN', min(min(min(amplitudes)))
1231	#print 'MAX', max(max(max(amplitudes)))
1232
1233	#print number_of_latitudes, number_of_longitudes
1234	number_of_points = number_of_latitudes*number_of_longitudes
1235	number_of_volumes = (number_of_latitudes-1)(number_of_longitudes-1)2
1236
1237	#print file_h.dimensions.keys()
1238	#print file_h.variables.keys()
1239
1240	file_h.close()
1241	file_u.close()
1242	file_v.close()
1243
1244
1245	if verbose: print 'Store to SWW file %s' %swwname
1246	# NetCDF file definition
1247	outfile = NetCDFFile(swwname, 'w')
1248
1249	#Create new file
1250	outfile.institution = 'Geoscience Australia'
1251	outfile.description = 'Converted from Ferret files: %s, %s, %s'\
1252	%(basename_in + '_ha.nc',
1253	basename_in + '_ua.nc',
1254	basename_in + '_va.nc')
1255
1256
1257	#For sww compatibility
1258	outfile.smoothing = 'Yes'
1259	outfile.order = 1
1260
1261	#Start time in seconds since the epoch (midnight 1/1/1970)
1262	outfile.starttime = times[0]
1263
1264	# dimension definitions
1265	outfile.createDimension('number_of_volumes', number_of_volumes)
1266
1267	outfile.createDimension('number_of_vertices', 3)
1268	outfile.createDimension('number_of_points', number_of_points)
1269
1270
1271	#outfile.createDimension('number_of_timesteps', len(times))
1272	outfile.createDimension('number_of_timesteps', len(times))
1273
1274	# variable definitions
1275	outfile.createVariable('x', precision, ('number_of_points',))
1276	outfile.createVariable('y', precision, ('number_of_points',))
1277	outfile.createVariable('elevation', precision, ('number_of_points',))
1278
1279	#FIXME: Backwards compatibility
1280	outfile.createVariable('z', precision, ('number_of_points',))
1281	#################################
1282
1283	outfile.createVariable('volumes', Int, ('number_of_volumes',
1284	'number_of_vertices'))
1285
1286	outfile.createVariable('time', precision,
1287	('number_of_timesteps',))
1288
1289	outfile.createVariable('stage', precision,
1290	('number_of_timesteps',
1291	'number_of_points'))
1292
1293	outfile.createVariable('xmomentum', precision,
1294	('number_of_timesteps',
1295	'number_of_points'))
1296
1297	outfile.createVariable('ymomentum', precision,
1298	('number_of_timesteps',
1299	'number_of_points'))
1300
1301
1302	#Store
1303	from coordinate_transforms.redfearn import redfearn
1304	x = zeros(number_of_points, Float) #Easting
1305	y = zeros(number_of_points, Float) #Northing
1306	#volumes = zeros(number_of_volumes, Int)
1307	i = 0
1308
1309	#Check zone boundaries
1310	refzone, _, _ = redfearn(latitudes[0],longitudes[0])
1311
1312	vertices = {}
1313	for k, lat in enumerate(latitudes):
1314	for l, lon in enumerate(longitudes):
1315
1316	vertices[l,k] = i
1317
1318	zone, easting, northing = redfearn(lat,lon)
1319
1320	msg = 'Zone boundary crossed at longitude =', lon
1321	assert zone == refzone, msg
1322	#print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
1323	x[i] = easting
1324	y[i] = northing
1325	i += 1
1326
1327
1328	#Construct 2 triangles per 'rectangular' element
1329	volumes = []
1330	for l in range(number_of_longitudes-1):
1331	for k in range(number_of_latitudes-1):
1332	v1 = vertices[l,k+1]
1333	v2 = vertices[l,k]
1334	v3 = vertices[l+1,k+1]
1335	v4 = vertices[l+1,k]
1336
1337	volumes.append([v1,v2,v3]) #Upper element
1338	volumes.append([v4,v3,v2]) #Lower element
1339
1340	volumes = array(volumes)
1341
1342	if origin == None:
1343	zone = refzone
1344	xllcorner = min(x)
1345	yllcorner = min(y)
1346	else:
1347	zone = origin[0]
1348	xllcorner = origin[1]
1349	yllcorner = origin[2]
1350
1351
1352	outfile.xllcorner = xllcorner
1353	outfile.yllcorner = yllcorner
1354	outfile.zone = zone
1355
1356	outfile.variables['x'][:] = x - xllcorner
1357	outfile.variables['y'][:] = y - yllcorner
1358	outfile.variables['z'][:] = 0.0
1359	outfile.variables['elevation'][:] = 0.0
1360	outfile.variables['time'][:] = times
1361	outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
1362
1363
1364
1365	#Time stepping
1366	stage = outfile.variables['stage']
1367	xmomentum = outfile.variables['xmomentum']
1368	ymomentum = outfile.variables['ymomentum']
1369
1370	for j in range(len(times)):
1371	i = 0
1372	for k in range(number_of_latitudes):
1373	for l in range(number_of_longitudes):
1374	h = zscale*amplitudes[j,k,l]/100 + mean_stage
1375	stage[j,i] = h
1376	xmomentum[j,i] = xspeed[j,k,l]/100*h
1377	ymomentum[j,i] = yspeed[j,k,l]/100*h
1378	i += 1
1379
1380	outfile.close()
1381
1382
1383	def extent_sww(file_name):
1384	"""
1385	Read in an sww file.
1386
1387	Input;
1388	file_name - the sww file
1389
1390	Output;
1391	z - Vector of bed elevation
1392	volumes - Array. Each row has 3 values, representing
1393	the vertices that define the volume
1394	time - Vector of the times where there is stage information
1395	stage - array with respect to time and vertices (x,y)
1396	"""
1397
1398
1399	from Scientific.IO.NetCDF import NetCDFFile
1400
1401	#Check contents
1402	#Get NetCDF
1403	fid = NetCDFFile(file_name, 'r')
1404
1405	# Get the variables
1406	x = fid.variables['x'][:]
1407	y = fid.variables['y'][:]
1408	stage = fid.variables['stage'][:]
1409	#print "stage",stage
1410	#print "stage.shap",stage.shape
1411	#print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
1412	#print "min(stage)",min(stage)
1413
1414	fid.close()
1415
1416	return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
1417
1418
1419	#OBSOLETE STUFF
1420	#Native checkpoint format.
1421	#Information needed to recreate a state is preserved
1422	#FIXME: Rethink and maybe use netcdf format
1423	def cpt_variable_writer(filename, t, v0, v1, v2):
1424	"""Store all conserved quantities to file
1425	"""
1426
1427	M, N = v0.shape
1428
1429	FN = create_filename(filename, 'cpt', M, t)
1430	#print 'Writing to %s' %FN
1431
1432	fid = open(FN, 'w')
1433	for i in range(M):
1434	for j in range(N):
1435	fid.write('%.16e ' %v0[i,j])
1436	for j in range(N):
1437	fid.write('%.16e ' %v1[i,j])
1438	for j in range(N):
1439	fid.write('%.16e ' %v2[i,j])
1440
1441	fid.write('\n')
1442	fid.close()
1443
1444
1445	def cpt_variable_reader(filename, t, v0, v1, v2):
1446	"""Store all conserved quantities to file
1447	"""
1448
1449	M, N = v0.shape
1450
1451	FN = create_filename(filename, 'cpt', M, t)
1452	#print 'Reading from %s' %FN
1453
1454	fid = open(FN)
1455
1456
1457	for i in range(M):
1458	values = fid.readline().split() #Get one line
1459
1460	for j in range(N):
1461	v0[i,j] = float(values[j])
1462	v1[i,j] = float(values[3+j])
1463	v2[i,j] = float(values[6+j])
1464
1465	fid.close()
1466
1467	def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1468	"""Writes x,y,z,z,z coordinates of triangles constituting the bed
1469	elevation.
1470	Not in use pt
1471	"""
1472
1473	M, N = v0.shape
1474
1475	print X0
1476	import sys; sys.exit()
1477	FN = create_filename(filename, 'cpt', M)
1478	print 'Writing to %s' %FN
1479
1480	fid = open(FN, 'w')
1481	for i in range(M):
1482	for j in range(2):
1483	fid.write('%.16e ' %X0[i,j]) #x, y
1484	for j in range(N):
1485	fid.write('%.16e ' %v0[i,j]) #z,z,z,
1486
1487	for j in range(2):
1488	fid.write('%.16e ' %X1[i,j]) #x, y
1489	for j in range(N):
1490	fid.write('%.16e ' %v1[i,j])
1491
1492	for j in range(2):
1493	fid.write('%.16e ' %X2[i,j]) #x, y
1494	for j in range(N):
1495	fid.write('%.16e ' %v2[i,j])
1496
1497	fid.write('\n')
1498	fid.close()
1499
1500
1501
1502	#Function for storing out to e.g. visualisation
1503	#FIXME: Do we want this?
1504	#FIXME: Not done yet for this version
1505	def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1506	"""Writes x,y,z coordinates of triangles constituting the bed elevation.
1507	"""
1508
1509	M, N = v0.shape
1510
1511	FN = create_filename(filename, 'dat', M)
1512	#print 'Writing to %s' %FN
1513
1514	fid = open(FN, 'w')
1515	for i in range(M):
1516	for j in range(2):
1517	fid.write('%f ' %X0[i,j]) #x, y
1518	fid.write('%f ' %v0[i,0]) #z
1519
1520	for j in range(2):
1521	fid.write('%f ' %X1[i,j]) #x, y
1522	fid.write('%f ' %v1[i,0]) #z
1523
1524	for j in range(2):
1525	fid.write('%f ' %X2[i,j]) #x, y
1526	fid.write('%f ' %v2[i,0]) #z
1527
1528	fid.write('\n')
1529	fid.close()
1530
1531
1532
1533	def dat_variable_writer(filename, t, v0, v1, v2):
1534	"""Store water height to file
1535	"""
1536
1537	M, N = v0.shape
1538
1539	FN = create_filename(filename, 'dat', M, t)
1540	#print 'Writing to %s' %FN
1541
1542	fid = open(FN, 'w')
1543	for i in range(M):
1544	fid.write('%.4f ' %v0[i,0])
1545	fid.write('%.4f ' %v1[i,0])
1546	fid.write('%.4f ' %v2[i,0])
1547
1548	fid.write('\n')
1549	fid.close()
1550

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

Download in other formats: