Context Navigation

source: trunk/anuga_core/source/anuga/load_mesh/loadASCII.py @ 8050

Last change on this file since 8050 was 7317, checked in by rwilson, 15 years ago
Replaced 'print' statements with log.critical() calls.
File size: 40.1 KB

Line
1	"""
2	The format for a Points dictionary is:
3
4	['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
5	['attributelist'], a dictionary of 1D arrays, representing attribute values
6	at the point. The dictionary key is the attribute header.
7	['geo_reference'] a Geo_refernece object. Use if the point information
8	is relative. This is optional.
9	eg
10	dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
11	dic['attributelist']['elevation'] = [[7.0,5.0]
12
13
14	The dict format for IO with mesh files is;
15	(the triangulation)
16	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
17	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
18	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
19	segments: [[v1,v2],[v3,v4],...] (lists of integers)
20	segment_tags : [tag,tag,...] list of strings
21	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
22	triangle_tags: [s1,s2,...] A list of strings
23	triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
24
25	(the outline)
26	points: [[x1,y1],[x2,y2],...] (lists of doubles)
27	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
28	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
29	outline_segment_tags : [tag1,tag2,...] list of strings
30	holes : [[x1,y1],...](List of doubles, one inside each hole region)
31	regions : [ [x1,y1],...] (List of 4 doubles)
32	region_tags : [tag1,tag2,...] (list of strings)
33	region_max_areas: [ma1,ma2,...] (A list of doubles)
34	{Convension: A -ve max area means no max area}
35	geo_reference: a Geo_refernece object. Use if the point information
36	is relative. This is optional.
37
38	Points files are .csv for ascii and .pts for NetCDF
39	Mesh files are .tsh for ascii and .msh for NetCDF
40
41	Note: point att's have no titles - that's currently ok, since least
42	squares adds the point info to vertices, and they have titles
43
44	The format for a .tsh file is (FIXME update this)
45
46	First line: <# of vertices> <# of attributes>
47	Following lines: <vertex #> <x> <y> [attributes]
48	One line: <# of triangles>
49	Following lines: <triangle #> <vertex #> <vertex #> <vertex #>
50	<neigbouring triangle #> <neigbouring triangle #>
51	<neigbouring triangle #> [attribute of region]
52	One line: <# of segments>
53	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
54	"""
55
56	##FIXME (DSG-DSG) Is the dict format mentioned above a list of a numeric array?
57	# Needs to be defined
58
59
60	from string import find, rfind
61	from os.path import splitext
62	import exceptions
63
64	from anuga.coordinate_transforms.geo_reference import Geo_reference, TITLE, \
65	TitleError
66	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
67	from anuga.config import netcdf_float, netcdf_char, netcdf_int
68	from anuga.utilities.system_tools import *
69	import anuga.utilities.log as log
70
71	from Scientific.IO.NetCDF import NetCDFFile
72
73	import numpy as num
74
75
76	class TitleAmountError(exceptions.Exception): pass
77
78
79	NOMAXAREA=-999
80
81
82	##
83	# @brief Read a mesh file (.tsh or .msh) and return a dictionary.
84	# @param ofile Path to the file to read.
85	# @return A dictionary of data from the input file.
86	# @note Throws IOError if can't read file.
87	# @note This is used by pmesh.
88	def import_mesh_file(ofile):
89	"""Read a mesh file, either .tsh or .msh
90
91	Note: will throw an IOError if it can't load the file.
92	"""
93
94	try:
95	if ofile[-4:] == ".tsh":
96	dict = _read_tsh_file(ofile)
97	elif ofile[-4:] == ".msh":
98	dict = _read_msh_file(ofile)
99	else:
100	msg = 'Extension .%s is unknown' % ofile[-4:]
101	raise IOError, msg
102	#FIXME No test for ValueError
103	except (TitleError, SyntaxError, IndexError, ValueError):
104	msg = 'File could not be opened'
105	raise IOError, msg
106	return dict
107
108
109	##
110	# @brief Write a mesh file from a dictionary.
111	# @param ofile The path of the mesh file to write.
112	# @param mesh_dict Dictionary containing data to write to output mesh file.
113	# @note Raises IOError if file extension is unrecognized.
114	def export_mesh_file(ofile, mesh_dict):
115	"""Write a mesh file given a dictionary.
116
117	First line: <# of vertices> <# of attributes>
118	Following lines: <vertex #> <x> <y> [attributes]
119	One line: <# of triangles>
120	Following lines: <triangle #> <vertex #> <vertex #> <vertex #>
121	<neigbouring triangle #> <neigbouring triangle #>
122	<neigbouring triangle #> [attribute of region]
123	One line: <# of segments>
124	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
125	"""
126
127	# Ensure that if required key isn't present, we add it with [] value.
128	# .setdefault() for dictionaries was added in python 2.0.
129	# The only other neat way to do this is to override the dictionary
130	# .__getitem__() method to return [] if key not found.
131	reqd_keys = ['points', 'point_attributes', 'outline_segments',
132	'outline_segment_tags', 'holes', 'regions', 'region_tags',
133	'region_max_areas', 'vertices', 'vertex_attributes',
134	'vertex_attribute_titles', 'segments', 'segment_tags',
135	'triangles', 'triangle_tags', 'triangle_neighbors']
136
137	for k in reqd_keys:
138	mesh_dict.setdefault(k, [])
139
140	# hand-off to appropriate function
141	if (ofile[-4:] == ".tsh"):
142	_write_tsh_file(ofile, mesh_dict)
143	elif (ofile[-4:] == ".msh"):
144	_write_msh_file(ofile, mesh_dict)
145	else:
146	msg = 'Unknown file type %s ' % ofile
147	raise IOError, msg
148
149
150	##
151	# @brief Read a mesh file into a dictionary.
152	# @param ofile Path of the file to read.
153	# @return Points dictionary from the mesh (.tsh) file.
154	def _read_tsh_file(ofile):
155	"""Read the text file format for meshes"""
156
157	fd = open(ofile, 'r')
158	dict = _read_triangulation(fd)
159	dict_mesh = _read_outline(fd)
160	for element in dict_mesh.keys():
161	dict[element] = dict_mesh[element]
162	fd.close()
163
164	return dict
165
166
167	##
168	# @brief Read triangulation data from a file, leave file open.
169	# @param fd An open descriptor of the file to read.
170	# @return A dictionary ...
171	def _read_triangulation(fd):
172	"""Read the generated triangulation, NOT the outline."""
173
174	delimiter = " "
175
176	# loading the point info
177	line = fd.readline()
178	fragments = line.split()
179	if fragments == []:
180	NumOfVertices = 0
181	NumOfVertAttributes = 0
182	else:
183	NumOfVertices = fragments[0]
184	NumOfVertAttributes = fragments[1]
185	points = []
186	pointattributes = []
187	for index in range(int(NumOfVertices)):
188	fragments = fd.readline().split()
189	fragments.pop(0) # pop off the index
190
191	# pop the x & y off so we're left with a list of attributes
192	vert = [float(fragments.pop(0)), float(fragments.pop(0))]
193	points.append(vert)
194	apointattributes = []
195	for fragment in fragments:
196	apointattributes.append(float(fragment))
197	if apointattributes != []:
198	pointattributes.append(apointattributes)
199
200	# loading the point title info
201	line = fd.readline()
202	vertTitle = []
203	for index in range(int(NumOfVertAttributes)):
204	fragments = fd.readline().strip()
205	vertTitle.append(fragments)
206
207	# loading the triangle info
208	line = fd.readline()
209	fragments = line.split()
210	NumOfTriangles = fragments[0]
211	triangles = []
212	triangleattributes = []
213	triangleneighbors = []
214	for index in range(int(NumOfTriangles)):
215	line = fd.readline()
216	line.strip() # so we can get the region string
217	fragments = line.split()
218	fragments.pop(0) # pop off the index
219
220	tri = [int(fragments[0]), int(fragments[1]), int(fragments[2])]
221	triangles.append(tri)
222	neighbors = [int(fragments[3]), int(fragments[4]), int(fragments[5])]
223	triangleneighbors.append(neighbors)
224	for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
225	line = line[find(line, delimiter):] # remove index
226	line = line.lstrip()
227	stringtag = line.strip()
228	triangleattributes.append(stringtag)
229
230	# loading the segment info
231	line = fd.readline()
232	fragments = line.split()
233	NumOfSegments = fragments[0]
234	segments = []
235	segmenttags = []
236	for index in range(int(NumOfSegments)):
237	line = fd.readline()
238	line.strip() # to get the segment string
239	fragments = line.split()
240	fragments.pop(0) #pop off the index
241	seg = [int(fragments[0]), int(fragments[1])]
242	segments.append(seg)
243	line = line[find(line, delimiter):] # remove index
244	line = line.lstrip()
245	line = line[find(line, delimiter):] # remove x
246	line = line.lstrip()
247	line = line[find(line, delimiter):] # remove y
248	stringtag = line.strip()
249	segmenttags.append(stringtag)
250
251	meshDict = {}
252	meshDict['vertices'] = points
253	if pointattributes == []:
254	meshDict['vertex_attributes'] = None
255	else:
256	meshDict['vertex_attributes'] = pointattributes
257	meshDict['triangles'] = triangles
258	meshDict['triangle_tags'] = triangleattributes
259	meshDict['triangle_neighbors'] = triangleneighbors
260	meshDict['segments'] = segments
261	meshDict['segment_tags'] = segmenttags
262	meshDict['vertex_attribute_titles'] = vertTitle
263
264	return meshDict
265
266
267	##
268	# @brief Read a mesh outline file.
269	# @param fd An open descriptor of the file to read.
270	# @return A Points dictionary.
271	# @note If file has no mesh info, an empty dict will be returned.
272	def _read_outline(fd):
273	"""Read a mesh file outline.
274
275	Note, if a file has no mesh info, it can still be read -
276	the meshdic returned will be 'empty'.
277	"""
278
279	delimiter = " " # warning: split() calls are using default whitespace
280
281	# loading the point info
282	line = fd.readline()
283	fragments = line.split()
284	if fragments == []:
285	NumOfVertices = 0
286	NumOfVertAttributes = 0
287	else:
288	NumOfVertices = fragments[0]
289	NumOfVertAttributes = fragments[1]
290	points = []
291	pointattributes = []
292	for index in range(int(NumOfVertices)):
293	fragments = fd.readline().split()
294	fragments.pop(0) # pop off the index
295	# pop the x & y off so we're left with a list of attributes
296	vert = [float(fragments.pop(0)), float(fragments.pop(0))]
297	points.append(vert)
298	apointattributes = []
299	for fragment in fragments:
300	apointattributes.append(float(fragment))
301	pointattributes.append(apointattributes)
302
303	# loading the segment info
304	line = fd.readline()
305	fragments = line.split()
306	if fragments == []:
307	NumOfSegments = 0
308	else:
309	NumOfSegments = fragments[0]
310	segments = []
311	segmenttags = []
312	for index in range(int(NumOfSegments)):
313	line = fd.readline()
314	fragments = line.split()
315	fragments.pop(0) # pop off the index
316	seg = [int(fragments[0]), int(fragments[1])]
317	segments.append(seg)
318	line = line[find(line, delimiter):] # remove index
319	line = line.lstrip()
320	line = line[find(line, delimiter):] # remove x
321	line = line.lstrip()
322	line = line[find(line, delimiter):] # remove y
323	stringtag = line.strip()
324	segmenttags.append(stringtag)
325
326	# loading the hole info
327	line = fd.readline()
328	fragments = line.split()
329	if fragments == []:
330	numOfHoles = 0
331	else:
332	numOfHoles = fragments[0]
333	holes = []
334	for index in range(int(numOfHoles)):
335	fragments = fd.readline().split()
336	fragments.pop(0) #pop off the index
337	hole = [float(fragments[0]), float(fragments[1])]
338	holes.append(hole)
339
340	# loading the region info
341	line = fd.readline()
342	fragments = line.split()
343	if fragments == []:
344	numOfRegions = 0
345	else:
346	numOfRegions = fragments[0]
347	regions = []
348	regionattributes = []
349	for index in range(int(numOfRegions)):
350	line = fd.readline()
351	fragments = line.split()
352	fragments.pop(0) # pop off the index
353	region = [float(fragments[0]), float(fragments[1])]
354	regions.append(region)
355
356	line = line[find(line, delimiter):] # remove index
357	line = line.lstrip()
358	line = line[find(line, delimiter):] # remove x
359	line = line.lstrip()
360	line = line[find(line, delimiter):] # remove y
361	stringtag = line.strip()
362	regionattributes.append(stringtag)
363
364	regionmaxareas = []
365	line = fd.readline()
366	for index in range(int(numOfRegions)): # Read in the Max area info
367	line = fd.readline()
368	fragments = line.split()
369	# The try is here for format compatibility
370	try:
371	fragments.pop(0) # pop off the index
372	if len(fragments) == 0: # no max area
373	regionmaxareas.append(None)
374	else:
375	regionmaxareas.append(float(fragments[0]))
376	except IndexError, e:
377	regionmaxareas.append(None)
378
379	try:
380	geo_reference = Geo_reference(ASCIIFile=fd)
381	except:
382	#geo_ref not compulsory
383	geo_reference = None
384
385	meshDict = {}
386	meshDict['points'] = points
387	meshDict['point_attributes'] = pointattributes
388	meshDict['outline_segments'] = segments
389	meshDict['outline_segment_tags'] = segmenttags
390	meshDict['holes'] = holes
391	meshDict['regions'] = regions
392	meshDict['region_tags'] = regionattributes
393	meshDict['region_max_areas'] = regionmaxareas
394	meshDict['geo_reference'] = geo_reference
395
396	return meshDict
397
398
399	##
400	# @brief Write an ASCII triangulation file.
401	# @param fd An open descriptor to the file to write.
402	# @param gen_dict Dictionary containing data to write.
403	def _write_ASCII_triangulation(fd, gen_dict):
404	vertices = gen_dict['vertices']
405	vertices_attributes = gen_dict['vertex_attributes']
406	try:
407	vertices_attribute_titles = gen_dict['vertex_attribute_titles']
408	except KeyError, e:
409	#FIXME is this the best way?
410	if vertices_attributes == [] or vertices_attributes[0] == []:
411	vertices_attribute_titles = []
412	else:
413	raise KeyError, e
414
415	triangles = gen_dict['triangles']
416	triangles_attributes = gen_dict['triangle_tags']
417	triangle_neighbors = gen_dict['triangle_neighbors']
418
419	segments = gen_dict['segments']
420	segment_tags = gen_dict['segment_tags']
421
422	numVert = str(len(vertices))
423	# Don't understand why we have to do vertices_attributes[0] is None,
424	# but it doesn't work otherwise...
425	if (vertices_attributes == None or
426	numVert == "0" or
427	len(vertices_attributes) == 0):
428	numVertAttrib = "0"
429	else:
430	numVertAttrib = str(len(vertices_attributes[0]))
431
432	fd.write(numVert + " " + numVertAttrib +
433	" # <# of verts> <# of vert attributes>, next lines <vertex #> "
434	"<x> <y> [attributes] ...Triangulation Vertices...\n")
435
436	#<vertex #> <x> <y> [attributes]
437	index = 0
438	for vert in vertices:
439	attlist = ""
440
441	if vertices_attributes == None or vertices_attributes == []:
442	attlist = ""
443	else:
444	for att in vertices_attributes[index]:
445	attlist = attlist + str(att) + " "
446	attlist.strip()
447	fd.write(str(index) + " " + str(vert[0]) + " " + str(vert[1])
448	+ " " + attlist + "\n")
449	index += 1
450
451	# write comments for title
452	fd.write("# attribute column titles ...Triangulation Vertex Titles...\n")
453	for title in vertices_attribute_titles:
454	fd.write(title + "\n")
455
456	# <# of triangles>
457	n = len(triangles)
458	fd.write(str(n)
459	+ " # <# of triangles>, next lines <triangle #> [<vertex #>] "
460	"[<neigbouring triangle #>] [attribute of region] ..."
461	"Triangulation Triangles...\n")
462
463	# <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #>
464	# <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
465	for index in range(n):
466	neighbors = ""
467	tri = triangles[index]
468	if triangle_neighbors == []:
469	neighbors = "-1 -1 -1 "
470	else:
471	for neighbor in triangle_neighbors[index]:
472	if neighbor:
473	neighbors += str(neighbor) + " "
474	else:
475	if neighbor == 0:
476	neighbors += "0 "
477	else:
478	neighbors += "-1 "
479	# Warning even though a list is passed, only the first value
480	# is written. There's an assumption that the list only
481	# contains one item. This assumption is made since the
482	# dict that's being passed around is also be used to communicate
483	# with triangle, and it seems to have the option of returning
484	# more than one value for triangle attributex
485	if (triangles_attributes == None or
486	triangles_attributes == [] or
487	triangles_attributes[index] == ['']):
488	att = ""
489	else:
490	att = str(triangles_attributes[index])
491	fd.write(str(index) + " "
492	+ str(tri[0]) + " "
493	+ str(tri[1]) + " "
494	+ str(tri[2]) + " "
495	+ neighbors + " "
496	+ att + "\n")
497
498	# One line: <# of segments>
499	fd.write(str(len(segments)) +
500	" # <# of segments>, next lines <segment #> <vertex #> "
501	"<vertex #> [boundary tag] ...Triangulation Segments...\n")
502
503	# Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
504	for i in range(len(segments)):
505	seg = segments[i]
506	fd.write(str(i) + " "
507	+ str(seg[0]) + " "
508	+ str(seg[1]) + " "
509	+ str(segment_tags[i]) + "\n")
510
511
512	##
513	# @brief Write triangulation and outline to a .tsh file.
514	# @param ofile Path of the file to write.
515	# @mesh_dict Dictionary of data to write.
516	def _write_tsh_file(ofile, mesh_dict):
517	fd = open(ofile, 'w')
518	_write_ASCII_triangulation(fd, mesh_dict)
519	_write_ASCII_outline(fd, mesh_dict)
520	fd.close()
521
522
523	##
524	# @brief Write an ASCII outline to a file.
525	# @param fd An open descriptor of the file to write.
526	# @param dict Dictionary holding data to write.
527	def _write_ASCII_outline(fd, dict):
528	points = dict['points']
529	point_attributes = dict['point_attributes']
530	segments = dict['outline_segments']
531	segment_tags = dict['outline_segment_tags']
532	holes = dict['holes']
533	regions = dict['regions']
534	region_tags = dict['region_tags']
535	region_max_areas = dict['region_max_areas']
536
537	num_points = str(len(points))
538	if (num_points == '0'):
539	num_point_atts = '0'
540	else:
541	num_point_atts = str(len(point_attributes[0]))
542
543	fd.write(num_points + ' ' + num_point_atts +
544	' # <# of verts> <# of vert attributes>, next lines <vertex #> '
545	'<x> <y> [attributes] ...Mesh Vertices...\n')
546
547	# <x> <y> [attributes]
548	for i, point in enumerate(points):
549	attlist = ''
550	for att in point_attributes[i]:
551	attlist = attlist + str(att) + ' '
552	attlist.strip()
553	fd.write(str(i) + ' ' + str(point[0]) + ' ' + str(point[1]) + ' ' +
554	attlist + '\n')
555
556	# One line: <# of segments>
557	fd.write(str(len(segments)) +
558	' # <# of segments>, next lines <segment #> <vertex #> '
559	'<vertex #> [boundary tag] ...Mesh Segments...\n')
560
561	# Following lines: <vertex #> <vertex #> [boundary tag]
562	for i,seg in enumerate(segments):
563	fd.write(str(i) + ' ' + str(seg[0]) + ' ' + str(seg[1]) + ' ' +
564	str(segment_tags[i]) + '\n')
565
566	# One line: <# of holes>
567	fd.write(str(len(holes)) +
568	' # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...\n')
569	# <x> <y>
570	for i,h in enumerate(holes):
571	fd.write(str(i) + ' ' + str(h[0]) + ' ' + str(h[1]) + '\n')
572
573	# One line: <# of regions>
574	fd.write(str(len(regions)) +
575	' # <# of regions>, next lines <Region #> <x> <y> <tag>'
576	'...Mesh Regions...\n')
577
578	# <index> <x> <y> <tag>
579	for i,r in enumerate(regions):
580	fd.write(str(i) + ' ' + str(r[0]) + ' ' + str(r[1])+ ' ' +
581	str(region_tags[i]) + '\n')
582
583	# <index> [<MaxArea>\|'']
584
585	# One line: <# of regions>
586	fd.write(str(len(regions)) +
587	' # <# of regions>, next lines <Region #> [Max Area] '
588	'...Mesh Regions...\n')
589	for i,r in enumerate(regions):
590	area = str(region_max_areas[i])
591
592	fd.write(str(i) + ' ' + area + '\n')
593
594	# geo_reference info
595	if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
596	dict['geo_reference'].write_ASCII(fd)
597
598
599	##
600	# @brief Write a .msh file.
601	# @param file Path to the file to write.
602	# @param mesh Dictionary holding data to write.
603	def _write_msh_file(file_name, mesh):
604	"""Write .msh NetCDF file
605
606	WARNING: This function mangles the mesh data structure
607	"""
608
609	# FIXME(Ole and John): We ran into a problem on Bogong (64 bit)
610	# where integers appeared as arrays. This may be similar to
611	# problem seen by Steve in changeset:2778 where he had to wrap
612	# them in int. Now we are trying with the native Integer format
613	# (Int == 'l' == Int64). However, that caused casting errors, when
614	# 64bit arrays are to be assigned to their NetCDF counterparts. It
615	# seems that the NetCDF arrays are 32bit even though they are
616	# created with the type Int64. Need to look at the NetCDF library
617	# in more detail.
618
619	IntType = num.int32
620	#IntType = Int
621
622	#the triangulation
623	mesh['vertices'] = num.array(mesh['vertices'], num.float)
624	if mesh['vertex_attributes'] != None:
625	mesh['vertex_attributes'] = \
626	num.array(mesh['vertex_attributes'], num.float)
627	mesh['vertex_attribute_titles'] = \
628	num.array(string_to_char(mesh['vertex_attribute_titles']), num.character)
629	mesh['segments'] = num.array(mesh['segments'], IntType)
630	mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']),
631	num.character)
632	mesh['triangles'] = num.array(mesh['triangles'], IntType)
633	mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']),
634	num.character)
635	mesh['triangle_neighbors'] = \
636	num.array(mesh['triangle_neighbors'], IntType)
637
638	#the outline
639	mesh['points'] = num.array(mesh['points'], num.float)
640	mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float)
641	mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
642	mesh['outline_segment_tags'] = \
643	num.array(string_to_char(mesh['outline_segment_tags']), num.character)
644	mesh['holes'] = num.array(mesh['holes'], num.float)
645	mesh['regions'] = num.array(mesh['regions'], num.float)
646	mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character)
647	mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)
648
649	# NetCDF file definition
650	try:
651	outfile = NetCDFFile(file_name, netcdf_mode_w)
652	except IOError:
653	msg = 'File %s could not be created' % file_name
654	raise Exception, msg
655
656	#Create new file
657	outfile.institution = 'Geoscience Australia'
658	outfile.description = 'NetCDF format for compact and portable storage ' + \
659	'of spatial point data'
660
661	# dimension definitions - fixed
662	outfile.createDimension('num_of_dimensions', 2) # This is 2d data
663	outfile.createDimension('num_of_segment_ends', 2) # Segs have two points
664	outfile.createDimension('num_of_triangle_vertices', 3)
665	outfile.createDimension('num_of_triangle_faces', 3)
666	outfile.createDimension('num_of_region_max_area', 1)
667
668	# Create dimensions, variables and set the variables
669
670	# trianglulation
671	# vertices
672	if (mesh['vertices'].shape[0] > 0):
673	outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
674	outfile.createVariable('vertices', netcdf_float, ('num_of_vertices',
675	'num_of_dimensions'))
676	outfile.variables['vertices'][:] = mesh['vertices']
677	if (mesh['vertex_attributes'] != None and
678	(mesh['vertex_attributes'].shape[0] > 0 and
679	mesh['vertex_attributes'].shape[1] > 0)):
680	outfile.createDimension('num_of_vertex_attributes',
681	mesh['vertex_attributes'].shape[1])
682	outfile.createDimension('num_of_vertex_attribute_title_chars',
683	mesh['vertex_attribute_titles'].shape[1])
684	outfile.createVariable('vertex_attributes',
685	netcdf_float,
686	('num_of_vertices',
687	'num_of_vertex_attributes'))
688	outfile.createVariable('vertex_attribute_titles',
689	netcdf_char,
690	('num_of_vertex_attributes',
691	'num_of_vertex_attribute_title_chars'))
692	outfile.variables['vertex_attributes'][:] = \
693	mesh['vertex_attributes']
694	outfile.variables['vertex_attribute_titles'][:] = \
695	mesh['vertex_attribute_titles']
696
697	# segments
698	if (mesh['segments'].shape[0] > 0):
699	outfile.createDimension('num_of_segments', mesh['segments'].shape[0])
700	outfile.createVariable('segments', netcdf_int,
701	('num_of_segments', 'num_of_segment_ends'))
702	outfile.variables['segments'][:] = mesh['segments']
703	if (mesh['segment_tags'].shape[1] > 0):
704	outfile.createDimension('num_of_segment_tag_chars',
705	mesh['segment_tags'].shape[1])
706	outfile.createVariable('segment_tags',
707	netcdf_char,
708	('num_of_segments',
709	'num_of_segment_tag_chars'))
710	outfile.variables['segment_tags'][:] = mesh['segment_tags']
711
712	# triangles
713	if (mesh['triangles'].shape[0] > 0):
714	outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0])
715	outfile.createVariable('triangles', netcdf_int,
716	('num_of_triangles', 'num_of_triangle_vertices'))
717	outfile.createVariable('triangle_neighbors', netcdf_int,
718	('num_of_triangles', 'num_of_triangle_faces'))
719	outfile.variables['triangles'][:] = mesh['triangles']
720	outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
721	if (mesh['triangle_tags'] != None and
722	(mesh['triangle_tags'].shape[1] > 0)):
723	outfile.createDimension('num_of_triangle_tag_chars',
724	mesh['triangle_tags'].shape[1])
725	outfile.createVariable('triangle_tags', netcdf_char,
726	('num_of_triangles',
727	'num_of_triangle_tag_chars'))
728	outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
729
730	# outline
731	# points
732	if (mesh['points'].shape[0] > 0):
733	outfile.createDimension('num_of_points', mesh['points'].shape[0])
734	outfile.createVariable('points', netcdf_float,
735	('num_of_points', 'num_of_dimensions'))
736	outfile.variables['points'][:] = mesh['points']
737	if mesh['point_attributes'].shape[0] > 0 \
738	and mesh['point_attributes'].shape[1] > 0:
739	outfile.createDimension('num_of_point_attributes',
740	mesh['point_attributes'].shape[1])
741	outfile.createVariable('point_attributes', netcdf_float,
742	('num_of_points', 'num_of_point_attributes'))
743	outfile.variables['point_attributes'][:] = mesh['point_attributes']
744
745	# outline_segments
746	if mesh['outline_segments'].shape[0] > 0:
747	outfile.createDimension('num_of_outline_segments',
748	mesh['outline_segments'].shape[0])
749	outfile.createVariable('outline_segments', netcdf_int,
750	('num_of_outline_segments',
751	'num_of_segment_ends'))
752	outfile.variables['outline_segments'][:] = mesh['outline_segments']
753	if mesh['outline_segment_tags'].shape[1] > 0:
754	outfile.createDimension('num_of_outline_segment_tag_chars',
755	mesh['outline_segment_tags'].shape[1])
756	outfile.createVariable('outline_segment_tags', netcdf_char,
757	('num_of_outline_segments',
758	'num_of_outline_segment_tag_chars'))
759	outfile.variables['outline_segment_tags'][:] = \
760	mesh['outline_segment_tags']
761
762	# holes
763	if (mesh['holes'].shape[0] > 0):
764	outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
765	outfile.createVariable('holes', netcdf_float,
766	('num_of_holes', 'num_of_dimensions'))
767	outfile.variables['holes'][:] = mesh['holes']
768
769	# regions
770	if (mesh['regions'].shape[0] > 0):
771	outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
772	outfile.createVariable('regions', netcdf_float,
773	('num_of_regions', 'num_of_dimensions'))
774	outfile.createVariable('region_max_areas', netcdf_float,
775	('num_of_regions',))
776	outfile.variables['regions'][:] = mesh['regions']
777	outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
778	if (mesh['region_tags'].shape[1] > 0):
779	outfile.createDimension('num_of_region_tag_chars',
780	mesh['region_tags'].shape[1])
781	outfile.createVariable('region_tags', netcdf_char,
782	('num_of_regions',
783	'num_of_region_tag_chars'))
784	outfile.variables['region_tags'][:] = mesh['region_tags']
785
786	# geo_reference info
787	if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
788	mesh['geo_reference'].write_NetCDF(outfile)
789
790	outfile.close()
791
792
793	##
794	# @brief Read a mesh file.
795	# @param file_name Path to the file to read.
796	# @return A dictionary containing the .msh file data.
797	# @note Throws IOError if file not found.
798	def _read_msh_file(file_name):
799	""" Read in an msh file."""
800
801	#Check contents. Get NetCDF
802	fd = open(file_name, 'r')
803	fd.close()
804
805	# throws prints to screen if file not present
806	fid = NetCDFFile(file_name, netcdf_mode_r)
807	mesh = {}
808
809	# Get the variables - the triangulation
810	try:
811	mesh['vertices'] = fid.variables['vertices'][:]
812	except KeyError:
813	mesh['vertices'] = num.array([], num.int) #array default#
814
815	try:
816	mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
817	except KeyError:
818	mesh['vertex_attributes'] = None
819
820	mesh['vertex_attribute_titles'] = []
821	try:
822	titles = fid.variables['vertex_attribute_titles'][:]
823	mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles]
824	except KeyError:
825	pass
826
827	try:
828	mesh['segments'] = fid.variables['segments'][:]
829	except KeyError:
830	mesh['segments'] = num.array([], num.int) #array default#
831
832	mesh['segment_tags'] = []
833	try:
834	tags = fid.variables['segment_tags'][:]
835	mesh['segment_tags'] = [x.tostring().strip() for x in tags]
836	except KeyError:
837	for ob in mesh['segments']:
838	mesh['segment_tags'].append('')
839
840	try:
841	mesh['triangles'] = fid.variables['triangles'][:]
842	mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
843	except KeyError:
844	mesh['triangles'] = num.array([], num.int) #array default#
845	mesh['triangle_neighbors'] = num.array([], num.int) #array default#
846
847	mesh['triangle_tags'] = []
848	try:
849	tags = fid.variables['triangle_tags'][:]
850	mesh['triangle_tags'] = [x.tostring().strip() for x in tags]
851	except KeyError:
852	for ob in mesh['triangles']:
853	mesh['triangle_tags'].append('')
854
855	#the outline
856	try:
857	mesh['points'] = fid.variables['points'][:]
858	except KeyError:
859	mesh['points'] = []
860
861	try:
862	mesh['point_attributes'] = fid.variables['point_attributes'][:]
863	except KeyError:
864	mesh['point_attributes'] = []
865	for point in mesh['points']:
866	mesh['point_attributes'].append([])
867
868	try:
869	mesh['outline_segments'] = fid.variables['outline_segments'][:]
870	except KeyError:
871	mesh['outline_segments'] = num.array([], num.int) #array default#
872
873	mesh['outline_segment_tags'] =[]
874	try:
875	tags = fid.variables['outline_segment_tags'][:]
876	for i, tag in enumerate(tags):
877	mesh['outline_segment_tags'].append(tags[i].tostring().strip())
878	except KeyError:
879	for ob in mesh['outline_segments']:
880	mesh['outline_segment_tags'].append('')
881
882	try:
883	mesh['holes'] = fid.variables['holes'][:]
884	except KeyError:
885	mesh['holes'] = num.array([], num.int) #array default#
886
887	try:
888	mesh['regions'] = fid.variables['regions'][:]
889	except KeyError:
890	mesh['regions'] = num.array([], num.int) #array default#
891
892	mesh['region_tags'] =[]
893	try:
894	tags = fid.variables['region_tags'][:]
895	for i, tag in enumerate(tags):
896	mesh['region_tags'].append(tags[i].tostring().strip())
897	except KeyError:
898	for ob in mesh['regions']:
899	mesh['region_tags'].append('')
900
901	try:
902	mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
903	except KeyError:
904	mesh['region_max_areas'] = num.array([], num.int) #array default#
905
906	try:
907	geo_reference = Geo_reference(NetCDFObject=fid)
908	mesh['geo_reference'] = geo_reference
909	except AttributeError, e:
910	#geo_ref not compulsory
911	mesh['geo_reference'] = None
912
913	fid.close()
914
915	return mesh
916
917
918	##
919	# @brief Export a boundary file.
920	# @param file_name Path to the file to write.
921	# @param points List of point index pairs.
922	# @paran title Title string to write into file (first line).
923	# @param delimiter Delimiter string to write between points.
924	# @note Used by alpha shapes
925	def export_boundary_file(file_name, points, title, delimiter=','):
926	"""Export a boundary file.
927
928	Format:
929	First line: Title variable
930	Following lines: [point index][delimiter][point index]
931
932	file_name - the name of the new file
933	points - List of point index pairs [[p1, p2],[p3, p4]..]
934	title - info to write in the first line
935	"""
936
937	fd = open(file_name, 'w')
938
939	fd.write(title + '\n')
940
941	# [point index][delimiter][point index]
942	for point in points:
943	fd.write(str(point[0]) + delimiter + str(point[1]) + '\n')
944
945	fd.close()
946
947	################################################################################
948	# IMPORT/EXPORT POINTS FILES
949	################################################################################
950
951	##
952	# @brief
953	# @param point_atts
954	def extent_point_atts(point_atts):
955	"""Returns 4 points representing the extent
956	This loses attribute info.
957	"""
958
959	point_atts['pointlist'] = extent(point_atts['pointlist'])
960	point_atts['attributelist'] = {}
961
962	return point_atts
963
964
965	##
966	# @brief Get an extent array for a set of points.
967	# @param points A set of points.
968	# @return An extent array of form [[min_x, min_y], [max_x, min_y],
969	# [max_x, max_y], [min_x, max_y]]
970	def extent(points):
971	points = num.array(points, num.float)
972
973	max_x = min_x = points[0][0]
974	max_y = min_y = points[0][1]
975
976	for point in points[1:]:
977	x = point[0]
978	if x > max_x:
979	max_x = x
980	if x < min_x:
981	min_x = x
982
983	y = point[1]
984	if y > max_y:
985	max_y = y
986	if y < min_y:
987	min_y = y
988
989	extent = num.array([[min_x, min_y],
990	[max_x, min_y],
991	[max_x, max_y],
992	[min_x, max_y]])
993
994	return extent
995
996
997	##
998	# @brief Reduce a points file until number of points is less than limit.
999	# @param infile Path to input file to thin.
1000	# @param outfile Path to output thinned file.
1001	# @param max_points Max number of points in output file.
1002	# @param verbose True if this function is to be verbose.
1003	def reduce_pts(infile, outfile, max_points, verbose = False):
1004	"""Reduce a points file until less than given size.
1005
1006	Reduces a points file by removing every second point until the # of points
1007	is less than max_points.
1008	"""
1009
1010	# check out pts2rectangular in least squares, and the use of reduction.
1011	# Maybe it does the same sort of thing?
1012	point_atts = _read_pts_file(infile)
1013
1014	while point_atts['pointlist'].shape[0] > max_points:
1015	if verbose: log.critical("point_atts['pointlist'].shape[0]")
1016	point_atts = half_pts(point_atts)
1017
1018	export_points_file(outfile, point_atts)
1019
1020
1021	##
1022	# @brief
1023	# @param infile
1024	# @param max_points
1025	# @param delimiter
1026	# @param verbose True if this function is to be verbose.
1027	# @return A list of
1028	def produce_half_point_files(infile, max_points, delimiter, verbose=False):
1029	point_atts = _read_pts_file(infile)
1030	root, ext = splitext(infile)
1031	outfiles = []
1032
1033	if verbose: log.critical("# of points", point_atts['pointlist'].shape[0])
1034
1035	while point_atts['pointlist'].shape[0] > max_points:
1036	point_atts = half_pts(point_atts)
1037
1038	if verbose: log.critical("# of points = %s"
1039	% str(point_atts['pointlist'].shape[0]))
1040
1041	outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1042	outfiles.append(outfile)
1043	export_points_file(outfile, point_atts)
1044
1045	return outfiles
1046
1047
1048	##
1049	# @brief
1050	# @param point_atts Object with attribute of 'pointlist'.
1051	# @return
1052	def point_atts2array(point_atts):
1053	# convert attribute list to array of floats
1054	point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
1055
1056	for key in point_atts['attributelist'].keys():
1057	point_atts['attributelist'][key] = \
1058	num.array(point_atts['attributelist'][key], num.float)
1059
1060	return point_atts
1061
1062
1063	##
1064	# @brief ??
1065	# @param point_atts ??
1066	# @return ??
1067	def half_pts(point_atts):
1068	point_atts2array(point_atts)
1069	point_atts['pointlist'] = point_atts['pointlist'][::2]
1070
1071	for key in point_atts['attributelist'].keys():
1072	point_atts['attributelist'][key] = point_atts['attributelist'][key][::2]
1073
1074	return point_atts
1075
1076	##
1077	# @brief ??
1078	# @param dic ??
1079	# @return ??
1080	def concatinate_attributelist(dic):
1081	"""
1082	giving a dic[attribute title] = attribute
1083	return list of attribute titles, array of attributes
1084	"""
1085
1086	point_attributes = num.array([], num.float)
1087	keys = dic.keys()
1088	key = keys.pop(0)
1089	point_attributes = num.reshape(dic[key], (dic[key].shape[0], 1))
1090	for key in keys:
1091	reshaped = num.reshape(dic[key], (dic[key].shape[0], 1))
1092	point_attributes = num.concatenate([point_attributes, reshaped], axis=1)
1093
1094	return dic.keys(), point_attributes
1095
1096
1097	##
1098	# @brief
1099	# @param dict
1100	# @param indices_to_keep
1101	# @return
1102	# @note FIXME(dsg), turn this dict plus methods into a class?
1103	def take_points(dict, indices_to_keep):
1104	dict = point_atts2array(dict)
1105	# FIXME maybe the points data structure should become a class?
1106	dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0)
1107
1108	for key in dict['attributelist'].keys():
1109	dict['attributelist'][key] = num.take(dict['attributelist'][key],
1110	indices_to_keep, axis=0)
1111
1112	return dict
1113
1114	##
1115	# @brief ??
1116	# @param dict1 ??
1117	# @param dict2 ??
1118	# @return ??
1119	def add_point_dictionaries(dict1, dict2):
1120	"""
1121	"""
1122
1123	dict1 = point_atts2array(dict1)
1124	dict2 = point_atts2array(dict2)
1125
1126	combined = {}
1127	combined['pointlist'] = num.concatenate((dict2['pointlist'],
1128	dict1['pointlist']), axis=0)
1129
1130	atts = {}
1131	for key in dict2['attributelist'].keys():
1132	atts[key]= num.concatenate((dict2['attributelist'][key],
1133	dict1['attributelist'][key]), axis=0)
1134	combined['attributelist'] = atts
1135	combined['geo_reference'] = dict1['geo_reference']
1136
1137	return combined
1138
1139
1140	if __name__ == "__main__":
1141	pass

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

Download in other formats: