Context Navigation

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

Last change on this file since 6360 was 6360, checked in by rwilson, 15 years ago
Ongoing conversion changes.
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
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	print ("Before num.array(), mesh['segment_tags']=%s"
629	% str(mesh['segment_tags']))
630	mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character)
631	print ("After num.array(), mesh['segment_tags']=%s"
632	% str(mesh['segment_tags']))
633	print "mesh['segment_tags'].shape=%s" % str(mesh['segment_tags'].shape)
634	mesh['triangles'] = num.array(mesh['triangles'], IntType)
635	mesh['triangle_tags'] = num.array(mesh['triangle_tags'])
636	mesh['triangle_neighbors'] = \
637	num.array(mesh['triangle_neighbors'], IntType)
638
639	#the outline
640	mesh['points'] = num.array(mesh['points'], num.float)
641	mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float)
642	mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
643	mesh['outline_segment_tags'] = \
644	num.array(mesh['outline_segment_tags'], num.character)
645	mesh['holes'] = num.array(mesh['holes'], num.float)
646	mesh['regions'] = num.array(mesh['regions'], num.float)
647	mesh['region_tags'] = num.array(mesh['region_tags'], num.character)
648	mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)
649
650	# NetCDF file definition
651	try:
652	outfile = NetCDFFile(file_name, netcdf_mode_w)
653	except IOError:
654	msg = 'File %s could not be created' % file_name
655	raise Exception, msg
656
657	#Create new file
658	outfile.institution = 'Geoscience Australia'
659	outfile.description = 'NetCDF format for compact and portable storage ' + \
660	'of spatial point data'
661
662	# dimension definitions - fixed
663	outfile.createDimension('num_of_dimensions', 2) # This is 2d data
664	outfile.createDimension('num_of_segment_ends', 2) # Segs have two points
665	outfile.createDimension('num_of_triangle_vertices', 3)
666	outfile.createDimension('num_of_triangle_faces', 3)
667	outfile.createDimension('num_of_region_max_area', 1)
668
669	# Create dimensions, variables and set the variables
670
671	# trianglulation
672	# vertices
673	if (mesh['vertices'].shape[0] > 0):
674	outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
675	outfile.createVariable('vertices', netcdf_float, ('num_of_vertices',
676	'num_of_dimensions'))
677	outfile.variables['vertices'][:] = mesh['vertices']
678	if (mesh['vertex_attributes'] != None and
679	(mesh['vertex_attributes'].shape[0] > 0 and
680	mesh['vertex_attributes'].shape[1] > 0)):
681	outfile.createDimension('num_of_vertex_attributes',
682	mesh['vertex_attributes'].shape[1])
683	outfile.createDimension('num_of_vertex_attribute_title_chars',
684	mesh['vertex_attribute_titles'].shape[1])
685	outfile.createVariable('vertex_attributes',
686	netcdf_float,
687	('num_of_vertices',
688	'num_of_vertex_attributes'))
689	outfile.createVariable('vertex_attribute_titles',
690	netcdf_char,
691	('num_of_vertex_attributes',
692	'num_of_vertex_attribute_title_chars'))
693	outfile.variables['vertex_attributes'][:] = \
694	mesh['vertex_attributes']
695	outfile.variables['vertex_attribute_titles'][:] = \
696	mesh['vertex_attribute_titles']
697
698	# segments
699	if (mesh['segments'].shape[0] > 0):
700	outfile.createDimension('num_of_segments', mesh['segments'].shape[0])
701	outfile.createVariable('segments', netcdf_int,
702	('num_of_segments', 'num_of_segment_ends'))
703	outfile.variables['segments'][:] = mesh['segments']
704	if mesh['segment_tags'].shape[1] > 0:
705	outfile.createDimension('num_of_segment_tag_chars',
706	mesh['segment_tags'].shape[1])
707	outfile.createVariable('segment_tags',
708	netcdf_char,
709	('num_of_segments',
710	'num_of_segment_tag_chars'))
711	outfile.variables['segment_tags'][:] = mesh['segment_tags']
712
713	# triangles
714	if (mesh['triangles'].shape[0] > 0):
715	outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0])
716	outfile.createVariable('triangles', netcdf_int,
717	('num_of_triangles', 'num_of_triangle_vertices'))
718	outfile.createVariable('triangle_neighbors', netcdf_int,
719	('num_of_triangles', 'num_of_triangle_faces'))
720	outfile.variables['triangles'][:] = mesh['triangles']
721	outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
722	if (mesh['triangle_tags'] != None and
723	(mesh['triangle_tags'].shape[1] > 0)):
724	outfile.createDimension('num_of_triangle_tag_chars',
725	mesh['triangle_tags'].shape[1])
726	outfile.createVariable('triangle_tags', netcdf_char,
727	('num_of_triangles',
728	'num_of_triangle_tag_chars'))
729	outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
730
731	# outline
732	# points
733	if (mesh['points'].shape[0] > 0):
734	outfile.createDimension('num_of_points', mesh['points'].shape[0])
735	outfile.createVariable('points', netcdf_float,
736	('num_of_points', 'num_of_dimensions'))
737	outfile.variables['points'][:] = mesh['points']
738	if mesh['point_attributes'].shape[0] > 0 \
739	and mesh['point_attributes'].shape[1] > 0:
740	outfile.createDimension('num_of_point_attributes',
741	mesh['point_attributes'].shape[1])
742	outfile.createVariable('point_attributes', netcdf_float,
743	('num_of_points', 'num_of_point_attributes'))
744	outfile.variables['point_attributes'][:] = mesh['point_attributes']
745
746	# outline_segments
747	if mesh['outline_segments'].shape[0] > 0:
748	outfile.createDimension('num_of_outline_segments',
749	mesh['outline_segments'].shape[0])
750	outfile.createVariable('outline_segments', netcdf_int,
751	('num_of_outline_segments',
752	'num_of_segment_ends'))
753	outfile.variables['outline_segments'][:] = mesh['outline_segments']
754	if (mesh['outline_segment_tags'].shape[1] > 0):
755	outfile.createDimension('num_of_outline_segment_tag_chars',
756	mesh['outline_segment_tags'].shape[1])
757	outfile.createVariable('outline_segment_tags', netcdf_char,
758	('num_of_outline_segments',
759	'num_of_outline_segment_tag_chars'))
760	outfile.variables['outline_segment_tags'][:] = \
761	mesh['outline_segment_tags']
762
763	# holes
764	if (mesh['holes'].shape[0] > 0):
765	outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
766	outfile.createVariable('holes', netcdf_float,
767	('num_of_holes', 'num_of_dimensions'))
768	outfile.variables['holes'][:] = mesh['holes']
769
770	# regions
771	if (mesh['regions'].shape[0] > 0):
772	outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
773	outfile.createVariable('regions', netcdf_float,
774	('num_of_regions', 'num_of_dimensions'))
775	outfile.createVariable('region_max_areas', netcdf_float,
776	('num_of_regions',))
777	outfile.variables['regions'][:] = mesh['regions']
778	outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
779	if (mesh['region_tags'].shape[1] > 0):
780	outfile.createDimension('num_of_region_tag_chars',
781	mesh['region_tags'].shape[1])
782	outfile.createVariable('region_tags', netcdf_char,
783	('num_of_regions',
784	'num_of_region_tag_chars'))
785	outfile.variables['region_tags'][:] = mesh['region_tags']
786
787	# geo_reference info
788	if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
789	mesh['geo_reference'].write_NetCDF(outfile)
790
791	outfile.close()
792
793
794	##
795	# @brief Read a mesh file.
796	# @param file_name Path to the file to read.
797	# @return A dictionary containing the .msh file data.
798	# @note Throws IOError if file not found.
799	def _read_msh_file(file_name):
800	""" Read in an msh file."""
801
802	#Check contents. Get NetCDF
803	fd = open(file_name, 'r')
804	fd.close()
805
806	# throws prints to screen if file not present
807	fid = NetCDFFile(file_name, netcdf_mode_r)
808	mesh = {}
809
810	# Get the variables - the triangulation
811	try:
812	mesh['vertices'] = fid.variables['vertices'][:]
813	except KeyError:
814	mesh['vertices'] = num.array([], num.int) #array default#
815
816	try:
817	mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
818	except KeyError:
819	mesh['vertex_attributes'] = None
820
821	mesh['vertex_attribute_titles'] = []
822	try:
823	titles = fid.variables['vertex_attribute_titles'][:]
824	for i, title in enumerate(titles):
825	mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
826	except KeyError:
827	pass
828
829	try:
830	mesh['segments'] = fid.variables['segments'][:]
831	except KeyError:
832	mesh['segments'] = num.array([], num.int) #array default#
833
834	mesh['segment_tags'] =[]
835	try:
836	tags = fid.variables['segment_tags'][:]
837	for i, tag in enumerate(tags):
838	mesh['segment_tags'].append(tags[i].tostring().strip())
839	except KeyError:
840	for ob in mesh['segments']:
841	mesh['segment_tags'].append('')
842
843	try:
844	mesh['triangles'] = fid.variables['triangles'][:]
845	mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
846	except KeyError:
847	mesh['triangles'] = num.array([], num.int) #array default#
848	mesh['triangle_neighbors'] = num.array([], num.int) #array default#
849
850	mesh['triangle_tags'] = []
851	try:
852	tags = fid.variables['triangle_tags'][:]
853	for i, tag in enumerate(tags):
854	mesh['triangle_tags'].append(tags[i].tostring().strip())
855	except KeyError:
856	for ob in mesh['triangles']:
857	mesh['triangle_tags'].append('')
858
859	#the outline
860	try:
861	mesh['points'] = fid.variables['points'][:]
862	except KeyError:
863	mesh['points'] = []
864
865	try:
866	mesh['point_attributes'] = fid.variables['point_attributes'][:]
867	except KeyError:
868	mesh['point_attributes'] = []
869	for point in mesh['points']:
870	mesh['point_attributes'].append([])
871
872	try:
873	mesh['outline_segments'] = fid.variables['outline_segments'][:]
874	except KeyError:
875	mesh['outline_segments'] = num.array([], num.int) #array default#
876
877	mesh['outline_segment_tags'] =[]
878	try:
879	tags = fid.variables['outline_segment_tags'][:]
880	for i, tag in enumerate(tags):
881	mesh['outline_segment_tags'].append(tags[i].tostring().strip())
882	except KeyError:
883	for ob in mesh['outline_segments']:
884	mesh['outline_segment_tags'].append('')
885
886	try:
887	mesh['holes'] = fid.variables['holes'][:]
888	except KeyError:
889	mesh['holes'] = num.array([], num.int) #array default#
890
891	try:
892	mesh['regions'] = fid.variables['regions'][:]
893	except KeyError:
894	mesh['regions'] = num.array([], num.int) #array default#
895
896	mesh['region_tags'] =[]
897	try:
898	tags = fid.variables['region_tags'][:]
899	for i, tag in enumerate(tags):
900	mesh['region_tags'].append(tags[i].tostring().strip())
901	except KeyError:
902	for ob in mesh['regions']:
903	mesh['region_tags'].append('')
904
905	try:
906	mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
907	except KeyError:
908	mesh['region_max_areas'] = num.array([], num.int) #array default#
909
910	try:
911	geo_reference = Geo_reference(NetCDFObject=fid)
912	mesh['geo_reference'] = geo_reference
913	except AttributeError, e:
914	#geo_ref not compulsory
915	mesh['geo_reference'] = None
916
917	fid.close()
918
919	return mesh
920
921
922	##
923	# @brief Export a boundary file.
924	# @param file_name Path to the file to write.
925	# @param points List of point index pairs.
926	# @paran title Title string to write into file (first line).
927	# @param delimiter Delimiter string to write between points.
928	# @note Used by alpha shapes
929	def export_boundary_file(file_name, points, title, delimiter=','):
930	"""Export a boundary file.
931
932	Format:
933	First line: Title variable
934	Following lines: [point index][delimiter][point index]
935
936	file_name - the name of the new file
937	points - List of point index pairs [[p1, p2],[p3, p4]..]
938	title - info to write in the first line
939	"""
940
941	fd = open(file_name, 'w')
942
943	fd.write(title + '\n')
944
945	# [point index][delimiter][point index]
946	for point in points:
947	fd.write(str(point[0]) + delimiter + str(point[1]) + '\n')
948
949	fd.close()
950
951	################################################################################
952	# IMPORT/EXPORT POINTS FILES
953	################################################################################
954
955	##
956	# @brief
957	# @param point_atts
958	def extent_point_atts(point_atts):
959	"""Returns 4 points representing the extent
960	This loses attribute info.
961	"""
962
963	point_atts['pointlist'] = extent(point_atts['pointlist'])
964	point_atts['attributelist'] = {}
965
966	return point_atts
967
968
969	##
970	# @brief Get an extent array for a set of points.
971	# @param points A set of points.
972	# @return An extent array of form [[min_x, min_y], [max_x, min_y],
973	# [max_x, max_y], [min_x, max_y]]
974	def extent(points):
975	points = num.array(points, num.float)
976
977	max_x = min_x = points[0][0]
978	max_y = min_y = points[0][1]
979
980	for point in points[1:]:
981	x = point[0]
982	if x > max_x:
983	max_x = x
984	if x < min_x:
985	min_x = x
986
987	y = point[1]
988	if y > max_y:
989	max_y = y
990	if y < min_y:
991	min_y = y
992
993	extent = num.array([[min_x, min_y],
994	[max_x, min_y],
995	[max_x, max_y],
996	[min_x, max_y]])
997
998	return extent
999
1000
1001	##
1002	# @brief Reduce a points file until number of points is less than limit.
1003	# @param infile Path to input file to thin.
1004	# @param outfile Path to output thinned file.
1005	# @param max_points Max number of points in output file.
1006	# @param verbose True if this function is to be verbose.
1007	def reduce_pts(infile, outfile, max_points, verbose = False):
1008	"""Reduce a points file until less than given size.
1009
1010	Reduces a points file by removing every second point until the # of points
1011	is less than max_points.
1012	"""
1013
1014	# check out pts2rectangular in least squares, and the use of reduction.
1015	# Maybe it does the same sort of thing?
1016	point_atts = _read_pts_file(infile)
1017
1018	while point_atts['pointlist'].shape[0] > max_points:
1019	if verbose: print "point_atts['pointlist'].shape[0]"
1020	point_atts = half_pts(point_atts)
1021
1022	export_points_file(outfile, point_atts)
1023
1024
1025	##
1026	# @brief
1027	# @param infile
1028	# @param max_points
1029	# @param delimiter
1030	# @param verbose True if this function is to be verbose.
1031	# @return A list of
1032	def produce_half_point_files(infile, max_points, delimiter, verbose=False):
1033	point_atts = _read_pts_file(infile)
1034	root, ext = splitext(infile)
1035	outfiles = []
1036
1037	if verbose: print "# of points", point_atts['pointlist'].shape[0]
1038
1039	while point_atts['pointlist'].shape[0] > max_points:
1040	point_atts = half_pts(point_atts)
1041
1042	if verbose: print "# of points", point_atts['pointlist'].shape[0]
1043
1044	outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1045	outfiles.append(outfile)
1046	export_points_file(outfile, point_atts)
1047
1048	return outfiles
1049
1050
1051	##
1052	# @brief
1053	# @param point_atts Object with attribute of 'pointlist'.
1054	# @return
1055	def point_atts2array(point_atts):
1056	# convert attribute list to array of floats
1057	point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
1058
1059	for key in point_atts['attributelist'].keys():
1060	point_atts['attributelist'][key] = \
1061	num.array(point_atts['attributelist'][key], num.float)
1062
1063	return point_atts
1064
1065
1066	##
1067	# @brief ??
1068	# @param point_atts ??
1069	# @return ??
1070	def half_pts(point_atts):
1071	point_atts2array(point_atts)
1072	point_atts['pointlist'] = point_atts['pointlist'][::2]
1073
1074	for key in point_atts['attributelist'].keys():
1075	point_atts['attributelist'][key] = point_atts['attributelist'][key][::2]
1076
1077	return point_atts
1078
1079	##
1080	# @brief ??
1081	# @param dic ??
1082	# @return ??
1083	def concatinate_attributelist(dic):
1084	"""
1085	giving a dic[attribute title] = attribute
1086	return list of attribute titles, array of attributes
1087	"""
1088
1089	point_attributes = num.array([], num.float)
1090	keys = dic.keys()
1091	key = keys.pop(0)
1092	point_attributes = num.reshape(dic[key], (dic[key].shape[0], 1))
1093	for key in keys:
1094	reshaped = num.reshape(dic[key], (dic[key].shape[0], 1))
1095	point_attributes = num.concatenate([point_attributes, reshaped], axis=1)
1096
1097	return dic.keys(), point_attributes
1098
1099
1100	##
1101	# @brief
1102	# @param dict
1103	# @param indices_to_keep
1104	# @return
1105	# @note FIXME(dsg), turn this dict plus methods into a class?
1106	def take_points(dict, indices_to_keep):
1107	dict = point_atts2array(dict)
1108	dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep)
1109
1110	for key in dict['attributelist'].keys():
1111	dict['attributelist'][key] = num.take(dict['attributelist'][key],
1112	indices_to_keep)
1113
1114	return dict
1115
1116	##
1117	# @brief ??
1118	# @param dict1 ??
1119	# @param dict2 ??
1120	# @return ??
1121	def add_point_dictionaries(dict1, dict2):
1122	"""
1123	"""
1124
1125	dict1 = point_atts2array(dict1)
1126	dict2 = point_atts2array(dict2)
1127
1128	combined = {}
1129	combined['pointlist'] = num.concatenate((dict2['pointlist'],
1130	dict1['pointlist']), axis=0)
1131
1132	atts = {}
1133	for key in dict2['attributelist'].keys():
1134	atts[key]= num.concatenate((dict2['attributelist'][key],
1135	dict1['attributelist'][key]), axis=0)
1136	combined['attributelist'] = atts
1137	combined['geo_reference'] = dict1['geo_reference']
1138
1139	return combined
1140
1141
1142	if __name__ == "__main__":
1143	pass

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

Download in other formats: