Context Navigation

source: branches/numpy/anuga/load_mesh/loadASCII.py @ 6410

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

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

Download in other formats: