Context Navigation

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

Last change on this file since 7024 was 6904, checked in by rwilson, 16 years ago
Finished back-merge from Numeric trunk to numpy branch (and some code the other way).
File size: 39.9 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
70	from Scientific.IO.NetCDF import NetCDFFile
71
72	import numpy as num
73
74
75	class TitleAmountError(exceptions.Exception): pass
76
77
78	NOMAXAREA=-999
79
80
81	##
82	# @brief Read a mesh file (.tsh or .msh) and return a dictionary.
83	# @param ofile Path to the file to read.
84	# @return A dictionary of data from the input file.
85	# @note Throws IOError if can't read file.
86	# @note This is used by pmesh.
87	def import_mesh_file(ofile):
88	"""Read a mesh file, either .tsh or .msh
89
90	Note: will throw an IOError if it can't load the file.
91	"""
92
93	try:
94	if ofile[-4:] == ".tsh":
95	dict = _read_tsh_file(ofile)
96	elif ofile[-4:] == ".msh":
97	dict = _read_msh_file(ofile)
98	else:
99	msg = 'Extension .%s is unknown' % ofile[-4:]
100	raise IOError, msg
101	#FIXME No test for ValueError
102	except (TitleError, SyntaxError, IndexError, ValueError):
103	msg = 'File could not be opened'
104	raise IOError, msg
105	return dict
106
107
108	##
109	# @brief Write a mesh file from a dictionary.
110	# @param ofile The path of the mesh file to write.
111	# @param mesh_dict Dictionary containing data to write to output mesh file.
112	# @note Raises IOError if file extension is unrecognized.
113	def export_mesh_file(ofile, mesh_dict):
114	"""Write a mesh file given a dictionary.
115
116	First line: <# of vertices> <# of attributes>
117	Following lines: <vertex #> <x> <y> [attributes]
118	One line: <# of triangles>
119	Following lines: <triangle #> <vertex #> <vertex #> <vertex #>
120	<neigbouring triangle #> <neigbouring triangle #>
121	<neigbouring triangle #> [attribute of region]
122	One line: <# of segments>
123	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
124	"""
125
126	# Ensure that if required key isn't present, we add it with [] value.
127	# .setdefault() for dictionaries was added in python 2.0.
128	# The only other neat way to do this is to override the dictionary
129	# .__getitem__() method to return [] if key not found.
130	reqd_keys = ['points', 'point_attributes', 'outline_segments',
131	'outline_segment_tags', 'holes', 'regions', 'region_tags',
132	'region_max_areas', 'vertices', 'vertex_attributes',
133	'vertex_attribute_titles', 'segments', 'segment_tags',
134	'triangles', 'triangle_tags', 'triangle_neighbors']
135
136	for k in reqd_keys:
137	mesh_dict.setdefault(k, [])
138
139	# hand-off to appropriate function
140	if (ofile[-4:] == ".tsh"):
141	_write_tsh_file(ofile, mesh_dict)
142	elif (ofile[-4:] == ".msh"):
143	_write_msh_file(ofile, mesh_dict)
144	else:
145	msg = 'Unknown file type %s ' % ofile
146	raise IOError, msg
147
148
149	##
150	# @brief Read a mesh file into a dictionary.
151	# @param ofile Path of the file to read.
152	# @return Points dictionary from the mesh (.tsh) file.
153	def _read_tsh_file(ofile):
154	"""Read the text file format for meshes"""
155
156	fd = open(ofile, 'r')
157	dict = _read_triangulation(fd)
158	dict_mesh = _read_outline(fd)
159	for element in dict_mesh.keys():
160	dict[element] = dict_mesh[element]
161	fd.close()
162
163	return dict
164
165
166	##
167	# @brief Read triangulation data from a file, leave file open.
168	# @param fd An open descriptor of the file to read.
169	# @return A dictionary ...
170	def _read_triangulation(fd):
171	"""Read the generated triangulation, NOT the outline."""
172
173	delimiter = " "
174
175	# loading the point info
176	line = fd.readline()
177	fragments = line.split()
178	if fragments == []:
179	NumOfVertices = 0
180	NumOfVertAttributes = 0
181	else:
182	NumOfVertices = fragments[0]
183	NumOfVertAttributes = fragments[1]
184	points = []
185	pointattributes = []
186	for index in range(int(NumOfVertices)):
187	fragments = fd.readline().split()
188	fragments.pop(0) # pop off the index
189
190	# pop the x & y off so we're left with a list of attributes
191	vert = [float(fragments.pop(0)), float(fragments.pop(0))]
192	points.append(vert)
193	apointattributes = []
194	for fragment in fragments:
195	apointattributes.append(float(fragment))
196	if apointattributes != []:
197	pointattributes.append(apointattributes)
198
199	# loading the point title info
200	line = fd.readline()
201	vertTitle = []
202	for index in range(int(NumOfVertAttributes)):
203	fragments = fd.readline().strip()
204	vertTitle.append(fragments)
205
206	# loading the triangle info
207	line = fd.readline()
208	fragments = line.split()
209	NumOfTriangles = fragments[0]
210	triangles = []
211	triangleattributes = []
212	triangleneighbors = []
213	for index in range(int(NumOfTriangles)):
214	line = fd.readline()
215	line.strip() # so we can get the region string
216	fragments = line.split()
217	fragments.pop(0) # pop off the index
218
219	tri = [int(fragments[0]), int(fragments[1]), int(fragments[2])]
220	triangles.append(tri)
221	neighbors = [int(fragments[3]), int(fragments[4]), int(fragments[5])]
222	triangleneighbors.append(neighbors)
223	for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
224	line = line[find(line, delimiter):] # remove index
225	line = line.lstrip()
226	stringtag = line.strip()
227	triangleattributes.append(stringtag)
228
229	# loading the segment info
230	line = fd.readline()
231	fragments = line.split()
232	NumOfSegments = fragments[0]
233	segments = []
234	segmenttags = []
235	for index in range(int(NumOfSegments)):
236	line = fd.readline()
237	line.strip() # to get the segment string
238	fragments = line.split()
239	fragments.pop(0) #pop off the index
240	seg = [int(fragments[0]), int(fragments[1])]
241	segments.append(seg)
242	line = line[find(line, delimiter):] # remove index
243	line = line.lstrip()
244	line = line[find(line, delimiter):] # remove x
245	line = line.lstrip()
246	line = line[find(line, delimiter):] # remove y
247	stringtag = line.strip()
248	segmenttags.append(stringtag)
249
250	meshDict = {}
251	meshDict['vertices'] = points
252	if pointattributes == []:
253	meshDict['vertex_attributes'] = None
254	else:
255	meshDict['vertex_attributes'] = pointattributes
256	meshDict['triangles'] = triangles
257	meshDict['triangle_tags'] = triangleattributes
258	meshDict['triangle_neighbors'] = triangleneighbors
259	meshDict['segments'] = segments
260	meshDict['segment_tags'] = segmenttags
261	meshDict['vertex_attribute_titles'] = vertTitle
262
263	return meshDict
264
265
266	##
267	# @brief Read a mesh outline file.
268	# @param fd An open descriptor of the file to read.
269	# @return A Points dictionary.
270	# @note If file has no mesh info, an empty dict will be returned.
271	def _read_outline(fd):
272	"""Read a mesh file outline.
273
274	Note, if a file has no mesh info, it can still be read -
275	the meshdic returned will be 'empty'.
276	"""
277
278	delimiter = " " # warning: split() calls are using default whitespace
279
280	# loading the point info
281	line = fd.readline()
282	fragments = line.split()
283	if fragments == []:
284	NumOfVertices = 0
285	NumOfVertAttributes = 0
286	else:
287	NumOfVertices = fragments[0]
288	NumOfVertAttributes = fragments[1]
289	points = []
290	pointattributes = []
291	for index in range(int(NumOfVertices)):
292	fragments = fd.readline().split()
293	fragments.pop(0) # pop off the index
294	# pop the x & y off so we're left with a list of attributes
295	vert = [float(fragments.pop(0)), float(fragments.pop(0))]
296	points.append(vert)
297	apointattributes = []
298	for fragment in fragments:
299	apointattributes.append(float(fragment))
300	pointattributes.append(apointattributes)
301
302	# loading the segment info
303	line = fd.readline()
304	fragments = line.split()
305	if fragments == []:
306	NumOfSegments = 0
307	else:
308	NumOfSegments = fragments[0]
309	segments = []
310	segmenttags = []
311	for index in range(int(NumOfSegments)):
312	line = fd.readline()
313	fragments = line.split()
314	fragments.pop(0) # pop off the index
315	seg = [int(fragments[0]), int(fragments[1])]
316	segments.append(seg)
317	line = line[find(line, delimiter):] # remove index
318	line = line.lstrip()
319	line = line[find(line, delimiter):] # remove x
320	line = line.lstrip()
321	line = line[find(line, delimiter):] # remove y
322	stringtag = line.strip()
323	segmenttags.append(stringtag)
324
325	# loading the hole info
326	line = fd.readline()
327	fragments = line.split()
328	if fragments == []:
329	numOfHoles = 0
330	else:
331	numOfHoles = fragments[0]
332	holes = []
333	for index in range(int(numOfHoles)):
334	fragments = fd.readline().split()
335	fragments.pop(0) #pop off the index
336	hole = [float(fragments[0]), float(fragments[1])]
337	holes.append(hole)
338
339	# loading the region info
340	line = fd.readline()
341	fragments = line.split()
342	if fragments == []:
343	numOfRegions = 0
344	else:
345	numOfRegions = fragments[0]
346	regions = []
347	regionattributes = []
348	for index in range(int(numOfRegions)):
349	line = fd.readline()
350	fragments = line.split()
351	fragments.pop(0) # pop off the index
352	region = [float(fragments[0]), float(fragments[1])]
353	regions.append(region)
354
355	line = line[find(line, delimiter):] # remove index
356	line = line.lstrip()
357	line = line[find(line, delimiter):] # remove x
358	line = line.lstrip()
359	line = line[find(line, delimiter):] # remove y
360	stringtag = line.strip()
361	regionattributes.append(stringtag)
362
363	regionmaxareas = []
364	line = fd.readline()
365	for index in range(int(numOfRegions)): # Read in the Max area info
366	line = fd.readline()
367	fragments = line.split()
368	# The try is here for format compatibility
369	try:
370	fragments.pop(0) # pop off the index
371	if len(fragments) == 0: # no max area
372	regionmaxareas.append(None)
373	else:
374	regionmaxareas.append(float(fragments[0]))
375	except IndexError, e:
376	regionmaxareas.append(None)
377
378	try:
379	geo_reference = Geo_reference(ASCIIFile=fd)
380	except:
381	#geo_ref not compulsory
382	geo_reference = None
383
384	meshDict = {}
385	meshDict['points'] = points
386	meshDict['point_attributes'] = pointattributes
387	meshDict['outline_segments'] = segments
388	meshDict['outline_segment_tags'] = segmenttags
389	meshDict['holes'] = holes
390	meshDict['regions'] = regions
391	meshDict['region_tags'] = regionattributes
392	meshDict['region_max_areas'] = regionmaxareas
393	meshDict['geo_reference'] = geo_reference
394
395	return meshDict
396
397
398	##
399	# @brief Write an ASCII triangulation file.
400	# @param fd An open descriptor to the file to write.
401	# @param gen_dict Dictionary containing data to write.
402	def _write_ASCII_triangulation(fd, gen_dict):
403	vertices = gen_dict['vertices']
404	vertices_attributes = gen_dict['vertex_attributes']
405	try:
406	vertices_attribute_titles = gen_dict['vertex_attribute_titles']
407	except KeyError, e:
408	#FIXME is this the best way?
409	if vertices_attributes == [] or vertices_attributes[0] == []:
410	vertices_attribute_titles = []
411	else:
412	raise KeyError, e
413
414	triangles = gen_dict['triangles']
415	triangles_attributes = gen_dict['triangle_tags']
416	triangle_neighbors = gen_dict['triangle_neighbors']
417
418	segments = gen_dict['segments']
419	segment_tags = gen_dict['segment_tags']
420
421	numVert = str(len(vertices))
422	# Don't understand why we have to do vertices_attributes[0] is None,
423	# but it doesn't work otherwise...
424	if (vertices_attributes == None or
425	numVert == "0" or
426	len(vertices_attributes) == 0):
427	numVertAttrib = "0"
428	else:
429	numVertAttrib = str(len(vertices_attributes[0]))
430
431	fd.write(numVert + " " + numVertAttrib +
432	" # <# of verts> <# of vert attributes>, next lines <vertex #> "
433	"<x> <y> [attributes] ...Triangulation Vertices...\n")
434
435	#<vertex #> <x> <y> [attributes]
436	index = 0
437	for vert in vertices:
438	attlist = ""
439
440	if vertices_attributes == None or vertices_attributes == []:
441	attlist = ""
442	else:
443	for att in vertices_attributes[index]:
444	attlist = attlist + str(att) + " "
445	attlist.strip()
446	fd.write(str(index) + " " + str(vert[0]) + " " + str(vert[1])
447	+ " " + attlist + "\n")
448	index += 1
449
450	# write comments for title
451	fd.write("# attribute column titles ...Triangulation Vertex Titles...\n")
452	for title in vertices_attribute_titles:
453	fd.write(title + "\n")
454
455	# <# of triangles>
456	n = len(triangles)
457	fd.write(str(n)
458	+ " # <# of triangles>, next lines <triangle #> [<vertex #>] "
459	"[<neigbouring triangle #>] [attribute of region] ..."
460	"Triangulation Triangles...\n")
461
462	# <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #>
463	# <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
464	for index in range(n):
465	neighbors = ""
466	tri = triangles[index]
467	if triangle_neighbors == []:
468	neighbors = "-1 -1 -1 "
469	else:
470	for neighbor in triangle_neighbors[index]:
471	if neighbor:
472	neighbors += str(neighbor) + " "
473	else:
474	if neighbor == 0:
475	neighbors += "0 "
476	else:
477	neighbors += "-1 "
478	# Warning even though a list is passed, only the first value
479	# is written. There's an assumption that the list only
480	# contains one item. This assumption is made since the
481	# dict that's being passed around is also be used to communicate
482	# with triangle, and it seems to have the option of returning
483	# more than one value for triangle attributex
484	if (triangles_attributes == None or
485	triangles_attributes == [] or
486	triangles_attributes[index] == ['']):
487	att = ""
488	else:
489	att = str(triangles_attributes[index])
490	fd.write(str(index) + " "
491	+ str(tri[0]) + " "
492	+ str(tri[1]) + " "
493	+ str(tri[2]) + " "
494	+ neighbors + " "
495	+ att + "\n")
496
497	# One line: <# of segments>
498	fd.write(str(len(segments)) +
499	" # <# of segments>, next lines <segment #> <vertex #> "
500	"<vertex #> [boundary tag] ...Triangulation Segments...\n")
501
502	# Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
503	for i in range(len(segments)):
504	seg = segments[i]
505	fd.write(str(i) + " "
506	+ str(seg[0]) + " "
507	+ str(seg[1]) + " "
508	+ str(segment_tags[i]) + "\n")
509
510
511	##
512	# @brief Write triangulation and outline to a .tsh file.
513	# @param ofile Path of the file to write.
514	# @mesh_dict Dictionary of data to write.
515	def _write_tsh_file(ofile, mesh_dict):
516	fd = open(ofile, 'w')
517	_write_ASCII_triangulation(fd, mesh_dict)
518	_write_ASCII_outline(fd, mesh_dict)
519	fd.close()
520
521
522	##
523	# @brief Write an ASCII outline to a file.
524	# @param fd An open descriptor of the file to write.
525	# @param dict Dictionary holding data to write.
526	def _write_ASCII_outline(fd, dict):
527	points = dict['points']
528	point_attributes = dict['point_attributes']
529	segments = dict['outline_segments']
530	segment_tags = dict['outline_segment_tags']
531	holes = dict['holes']
532	regions = dict['regions']
533	region_tags = dict['region_tags']
534	region_max_areas = dict['region_max_areas']
535
536	num_points = str(len(points))
537	if (num_points == '0'):
538	num_point_atts = '0'
539	else:
540	num_point_atts = str(len(point_attributes[0]))
541
542	fd.write(num_points + ' ' + num_point_atts +
543	' # <# of verts> <# of vert attributes>, next lines <vertex #> '
544	'<x> <y> [attributes] ...Mesh Vertices...\n')
545
546	# <x> <y> [attributes]
547	for i, point in enumerate(points):
548	attlist = ''
549	for att in point_attributes[i]:
550	attlist = attlist + str(att) + ' '
551	attlist.strip()
552	fd.write(str(i) + ' ' + str(point[0]) + ' ' + str(point[1]) + ' ' +
553	attlist + '\n')
554
555	# One line: <# of segments>
556	fd.write(str(len(segments)) +
557	' # <# of segments>, next lines <segment #> <vertex #> '
558	'<vertex #> [boundary tag] ...Mesh Segments...\n')
559
560	# Following lines: <vertex #> <vertex #> [boundary tag]
561	for i,seg in enumerate(segments):
562	fd.write(str(i) + ' ' + str(seg[0]) + ' ' + str(seg[1]) + ' ' +
563	str(segment_tags[i]) + '\n')
564
565	# One line: <# of holes>
566	fd.write(str(len(holes)) +
567	' # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...\n')
568	# <x> <y>
569	for i,h in enumerate(holes):
570	fd.write(str(i) + ' ' + str(h[0]) + ' ' + str(h[1]) + '\n')
571
572	# One line: <# of regions>
573	fd.write(str(len(regions)) +
574	' # <# of regions>, next lines <Region #> <x> <y> <tag>'
575	'...Mesh Regions...\n')
576
577	# <index> <x> <y> <tag>
578	for i,r in enumerate(regions):
579	fd.write(str(i) + ' ' + str(r[0]) + ' ' + str(r[1])+ ' ' +
580	str(region_tags[i]) + '\n')
581
582	# <index> [<MaxArea>\|'']
583
584	# One line: <# of regions>
585	fd.write(str(len(regions)) +
586	' # <# of regions>, next lines <Region #> [Max Area] '
587	'...Mesh Regions...\n')
588	for i,r in enumerate(regions):
589	area = str(region_max_areas[i])
590
591	fd.write(str(i) + ' ' + area + '\n')
592
593	# geo_reference info
594	if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
595	dict['geo_reference'].write_ASCII(fd)
596
597
598	##
599	# @brief Write a .msh file.
600	# @param file Path to the file to write.
601	# @param mesh Dictionary holding data to write.
602	def _write_msh_file(file_name, mesh):
603	"""Write .msh NetCDF file
604
605	WARNING: This function mangles the mesh data structure
606	"""
607
608	# FIXME(Ole and John): We ran into a problem on Bogong (64 bit)
609	# where integers appeared as arrays. This may be similar to
610	# problem seen by Steve in changeset:2778 where he had to wrap
611	# them in int. Now we are trying with the native Integer format
612	# (Int == 'l' == Int64). However, that caused casting errors, when
613	# 64bit arrays are to be assigned to their NetCDF counterparts. It
614	# seems that the NetCDF arrays are 32bit even though they are
615	# created with the type Int64. Need to look at the NetCDF library
616	# in more detail.
617
618	IntType = num.int32
619	#IntType = Int
620
621	#the triangulation
622	mesh['vertices'] = num.array(mesh['vertices'], num.float)
623	if mesh['vertex_attributes'] != None:
624	mesh['vertex_attributes'] = \
625	num.array(mesh['vertex_attributes'], num.float)
626	mesh['vertex_attribute_titles'] = \
627	num.array(string_to_char(mesh['vertex_attribute_titles']), num.character)
628	mesh['segments'] = num.array(mesh['segments'], IntType)
629	mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']),
630	num.character)
631	mesh['triangles'] = num.array(mesh['triangles'], IntType)
632	mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']),
633	num.character)
634	mesh['triangle_neighbors'] = \
635	num.array(mesh['triangle_neighbors'], IntType)
636
637	#the outline
638	mesh['points'] = num.array(mesh['points'], num.float)
639	mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float)
640	mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
641	mesh['outline_segment_tags'] = \
642	num.array(string_to_char(mesh['outline_segment_tags']), num.character)
643	mesh['holes'] = num.array(mesh['holes'], num.float)
644	mesh['regions'] = num.array(mesh['regions'], num.float)
645	mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character)
646	mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)
647
648	# NetCDF file definition
649	try:
650	outfile = NetCDFFile(file_name, netcdf_mode_w)
651	except IOError:
652	msg = 'File %s could not be created' % file_name
653	raise Exception, msg
654
655	#Create new file
656	outfile.institution = 'Geoscience Australia'
657	outfile.description = 'NetCDF format for compact and portable storage ' + \
658	'of spatial point data'
659
660	# dimension definitions - fixed
661	outfile.createDimension('num_of_dimensions', 2) # This is 2d data
662	outfile.createDimension('num_of_segment_ends', 2) # Segs have two points
663	outfile.createDimension('num_of_triangle_vertices', 3)
664	outfile.createDimension('num_of_triangle_faces', 3)
665	outfile.createDimension('num_of_region_max_area', 1)
666
667	# Create dimensions, variables and set the variables
668
669	# trianglulation
670	# vertices
671	if (mesh['vertices'].shape[0] > 0):
672	outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
673	outfile.createVariable('vertices', netcdf_float, ('num_of_vertices',
674	'num_of_dimensions'))
675	outfile.variables['vertices'][:] = mesh['vertices']
676	if (mesh['vertex_attributes'] != None and
677	(mesh['vertex_attributes'].shape[0] > 0 and
678	mesh['vertex_attributes'].shape[1] > 0)):
679	outfile.createDimension('num_of_vertex_attributes',
680	mesh['vertex_attributes'].shape[1])
681	outfile.createDimension('num_of_vertex_attribute_title_chars',
682	mesh['vertex_attribute_titles'].shape[1])
683	outfile.createVariable('vertex_attributes',
684	netcdf_float,
685	('num_of_vertices',
686	'num_of_vertex_attributes'))
687	outfile.createVariable('vertex_attribute_titles',
688	netcdf_char,
689	('num_of_vertex_attributes',
690	'num_of_vertex_attribute_title_chars'))
691	outfile.variables['vertex_attributes'][:] = \
692	mesh['vertex_attributes']
693	outfile.variables['vertex_attribute_titles'][:] = \
694	mesh['vertex_attribute_titles']
695
696	# segments
697	if (mesh['segments'].shape[0] > 0):
698	outfile.createDimension('num_of_segments', mesh['segments'].shape[0])
699	outfile.createVariable('segments', netcdf_int,
700	('num_of_segments', 'num_of_segment_ends'))
701	outfile.variables['segments'][:] = mesh['segments']
702	if mesh['segment_tags'].shape[0] > 0:
703	outfile.createDimension('num_of_segment_tag_chars',
704	mesh['segment_tags'].shape[1])
705	outfile.createVariable('segment_tags',
706	netcdf_char,
707	('num_of_segments',
708	'num_of_segment_tag_chars'))
709	outfile.variables['segment_tags'][:] = mesh['segment_tags']
710
711	# triangles
712	if (mesh['triangles'].shape[0] > 0):
713	outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0])
714	outfile.createVariable('triangles', netcdf_int,
715	('num_of_triangles', 'num_of_triangle_vertices'))
716	outfile.createVariable('triangle_neighbors', netcdf_int,
717	('num_of_triangles', 'num_of_triangle_faces'))
718	outfile.variables['triangles'][:] = mesh['triangles']
719	outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
720	if (mesh['triangle_tags'] != None and
721	(mesh['triangle_tags'].shape[1] > 0)):
722	outfile.createDimension('num_of_triangle_tag_chars',
723	mesh['triangle_tags'].shape[1])
724	outfile.createVariable('triangle_tags', netcdf_char,
725	('num_of_triangles',
726	'num_of_triangle_tag_chars'))
727	outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
728
729	# outline
730	# points
731	if (mesh['points'].shape[0] > 0):
732	outfile.createDimension('num_of_points', mesh['points'].shape[0])
733	outfile.createVariable('points', netcdf_float,
734	('num_of_points', 'num_of_dimensions'))
735	outfile.variables['points'][:] = mesh['points']
736	if mesh['point_attributes'].shape[0] > 0 \
737	and mesh['point_attributes'].shape[1] > 0:
738	outfile.createDimension('num_of_point_attributes',
739	mesh['point_attributes'].shape[1])
740	outfile.createVariable('point_attributes', netcdf_float,
741	('num_of_points', 'num_of_point_attributes'))
742	outfile.variables['point_attributes'][:] = mesh['point_attributes']
743
744	# outline_segments
745	if mesh['outline_segments'].shape[0] > 0:
746	outfile.createDimension('num_of_outline_segments',
747	mesh['outline_segments'].shape[0])
748	outfile.createVariable('outline_segments', netcdf_int,
749	('num_of_outline_segments',
750	'num_of_segment_ends'))
751	outfile.variables['outline_segments'][:] = mesh['outline_segments']
752	if mesh['outline_segment_tags'].shape[1] > 0:
753	outfile.createDimension('num_of_outline_segment_tag_chars',
754	mesh['outline_segment_tags'].shape[1])
755	outfile.createVariable('outline_segment_tags', netcdf_char,
756	('num_of_outline_segments',
757	'num_of_outline_segment_tag_chars'))
758	outfile.variables['outline_segment_tags'][:] = \
759	mesh['outline_segment_tags']
760
761	# holes
762	if (mesh['holes'].shape[0] > 0):
763	outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
764	outfile.createVariable('holes', netcdf_float,
765	('num_of_holes', 'num_of_dimensions'))
766	outfile.variables['holes'][:] = mesh['holes']
767
768	# regions
769	if (mesh['regions'].shape[0] > 0):
770	outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
771	outfile.createVariable('regions', netcdf_float,
772	('num_of_regions', 'num_of_dimensions'))
773	outfile.createVariable('region_max_areas', netcdf_float,
774	('num_of_regions',))
775	outfile.variables['regions'][:] = mesh['regions']
776	outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
777	if (mesh['region_tags'].shape[1] > 0):
778	outfile.createDimension('num_of_region_tag_chars',
779	mesh['region_tags'].shape[1])
780	outfile.createVariable('region_tags', netcdf_char,
781	('num_of_regions',
782	'num_of_region_tag_chars'))
783	outfile.variables['region_tags'][:] = mesh['region_tags']
784
785	# geo_reference info
786	if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
787	mesh['geo_reference'].write_NetCDF(outfile)
788
789	outfile.close()
790
791
792	##
793	# @brief Read a mesh file.
794	# @param file_name Path to the file to read.
795	# @return A dictionary containing the .msh file data.
796	# @note Throws IOError if file not found.
797	def _read_msh_file(file_name):
798	""" Read in an msh file."""
799
800	#Check contents. Get NetCDF
801	fd = open(file_name, 'r')
802	fd.close()
803
804	# throws prints to screen if file not present
805	fid = NetCDFFile(file_name, netcdf_mode_r)
806	mesh = {}
807
808	# Get the variables - the triangulation
809	try:
810	mesh['vertices'] = fid.variables['vertices'][:]
811	except KeyError:
812	mesh['vertices'] = num.array([], num.int) #array default#
813
814	try:
815	mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
816	except KeyError:
817	mesh['vertex_attributes'] = None
818
819	mesh['vertex_attribute_titles'] = []
820	try:
821	titles = fid.variables['vertex_attribute_titles'][:]
822	mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles]
823	except KeyError:
824	pass
825
826	try:
827	mesh['segments'] = fid.variables['segments'][:]
828	except KeyError:
829	mesh['segments'] = num.array([], num.int) #array default#
830
831	mesh['segment_tags'] = []
832	try:
833	tags = fid.variables['segment_tags'][:]
834	mesh['segment_tags'] = [x.tostring().strip() for x in tags]
835	except KeyError:
836	for ob in mesh['segments']:
837	mesh['segment_tags'].append('')
838
839	try:
840	mesh['triangles'] = fid.variables['triangles'][:]
841	mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
842	except KeyError:
843	mesh['triangles'] = num.array([], num.int) #array default#
844	mesh['triangle_neighbors'] = num.array([], num.int) #array default#
845
846	mesh['triangle_tags'] = []
847	try:
848	tags = fid.variables['triangle_tags'][:]
849	mesh['triangle_tags'] = [x.tostring().strip() for x in tags]
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: