Context Navigation

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

Last change on this file since 8272 was 8150, checked in by wilsonr, 14 years ago
Removed '@brief' comments.
File size: 37.1 KB

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

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

Download in other formats: