Context Navigation

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

Last change on this file since 8818 was 8780, checked in by steve, 12 years ago
Some changes to allow netcdf4 use
File size: 37.5 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 anuga.file.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	#print 'mesh vertices',mesh['vertices'].shape
583
584
585	#the triangulation
586	mesh['vertices'] = num.array(mesh['vertices'], num.float)
587	mesh['vertex_attribute_titles'] = \
588	num.array(string_to_char(mesh['vertex_attribute_titles']), num.character)
589
590	num_attributes = len(mesh['vertex_attribute_titles'])
591	num_vertices = mesh['vertices'].shape[0]
592	#print 'num_attrib ',num_attributes
593	if mesh['vertex_attributes'] != None:
594	mesh['vertex_attributes'] = \
595	num.array(mesh['vertex_attributes'], num.float)
596
597	if num_attributes > 0 :
598	mesh['vertex_attributes'] = \
599	num.reshape(mesh['vertex_attributes'],(num_vertices,-1))
600
601
602
603	mesh['segments'] = num.array(mesh['segments'], IntType)
604	mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']),
605	num.character)
606	mesh['triangles'] = num.array(mesh['triangles'], IntType)
607	mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']),
608	num.character)
609	mesh['triangle_neighbors'] = \
610	num.array(mesh['triangle_neighbors'], IntType)
611
612	#the outline
613	mesh['points'] = num.array(mesh['points'], num.float)
614	mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float)
615	mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
616	mesh['outline_segment_tags'] = \
617	num.array(string_to_char(mesh['outline_segment_tags']), num.character)
618	mesh['holes'] = num.array(mesh['holes'], num.float)
619	mesh['regions'] = num.array(mesh['regions'], num.float)
620	mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character)
621	mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)
622
623	# NetCDF file definition
624	try:
625	outfile = NetCDFFile(file_name, netcdf_mode_w)
626	except IOError:
627	msg = 'File %s could not be created' % file_name
628	raise Exception, msg
629
630	#Create new file
631	outfile.institution = 'Geoscience Australia'
632	outfile.description = 'NetCDF format for compact and portable storage ' + \
633	'of spatial point data'
634
635	# dimension definitions - fixed
636	outfile.createDimension('num_of_dimensions', 2) # This is 2d data
637	outfile.createDimension('num_of_segment_ends', 2) # Segs have two points
638	outfile.createDimension('num_of_triangle_vertices', 3)
639	outfile.createDimension('num_of_triangle_faces', 3)
640	outfile.createDimension('num_of_region_max_area', 1)
641
642	# Create dimensions, variables and set the variables
643
644	# trianglulation
645	# vertices
646	if (mesh['vertices'].shape[0] > 0):
647	outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
648	outfile.createVariable('vertices', netcdf_float, ('num_of_vertices',
649	'num_of_dimensions'))
650	outfile.variables['vertices'][:] = mesh['vertices']
651
652	#print 'mesh vertex attributes', mesh['vertex_attributes'].shape
653
654	if (mesh['vertex_attributes'] != None and
655	(mesh['vertex_attributes'].shape[0] > 0 and
656	mesh['vertex_attributes'].shape[1] > 0)):
657	outfile.createDimension('num_of_vertex_attributes',
658	mesh['vertex_attributes'].shape[1])
659	outfile.createDimension('num_of_vertex_attribute_title_chars',
660	mesh['vertex_attribute_titles'].shape[1])
661	outfile.createVariable('vertex_attributes',
662	netcdf_float,
663	('num_of_vertices',
664	'num_of_vertex_attributes'))
665	outfile.createVariable('vertex_attribute_titles',
666	netcdf_char,
667	('num_of_vertex_attributes',
668	'num_of_vertex_attribute_title_chars'))
669	outfile.variables['vertex_attributes'][:] = \
670	mesh['vertex_attributes']
671	outfile.variables['vertex_attribute_titles'][:] = \
672	mesh['vertex_attribute_titles']
673
674	# segments
675	if (mesh['segments'].shape[0] > 0):
676	outfile.createDimension('num_of_segments', mesh['segments'].shape[0])
677	outfile.createVariable('segments', netcdf_int,
678	('num_of_segments', 'num_of_segment_ends'))
679	outfile.variables['segments'][:] = mesh['segments']
680	if (mesh['segment_tags'].shape[1] > 0):
681	outfile.createDimension('num_of_segment_tag_chars',
682	mesh['segment_tags'].shape[1])
683	outfile.createVariable('segment_tags',
684	netcdf_char,
685	('num_of_segments',
686	'num_of_segment_tag_chars'))
687	outfile.variables['segment_tags'][:] = mesh['segment_tags']
688
689	# triangles
690	if (mesh['triangles'].shape[0] > 0):
691	outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0])
692	outfile.createVariable('triangles', netcdf_int,
693	('num_of_triangles', 'num_of_triangle_vertices'))
694	outfile.createVariable('triangle_neighbors', netcdf_int,
695	('num_of_triangles', 'num_of_triangle_faces'))
696	outfile.variables['triangles'][:] = mesh['triangles']
697	outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
698	if (mesh['triangle_tags'] != None and
699	(mesh['triangle_tags'].shape[1] > 0)):
700	outfile.createDimension('num_of_triangle_tag_chars',
701	mesh['triangle_tags'].shape[1])
702	outfile.createVariable('triangle_tags', netcdf_char,
703	('num_of_triangles',
704	'num_of_triangle_tag_chars'))
705	outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
706
707	# outline
708	# points
709	if (mesh['points'].shape[0] > 0):
710	outfile.createDimension('num_of_points', mesh['points'].shape[0])
711	outfile.createVariable('points', netcdf_float,
712	('num_of_points', 'num_of_dimensions'))
713	outfile.variables['points'][:] = mesh['points']
714	if mesh['point_attributes'].shape[0] > 0 \
715	and mesh['point_attributes'].shape[1] > 0:
716	outfile.createDimension('num_of_point_attributes',
717	mesh['point_attributes'].shape[1])
718	outfile.createVariable('point_attributes', netcdf_float,
719	('num_of_points', 'num_of_point_attributes'))
720	outfile.variables['point_attributes'][:] = mesh['point_attributes']
721
722	# outline_segments
723	if mesh['outline_segments'].shape[0] > 0:
724	outfile.createDimension('num_of_outline_segments',
725	mesh['outline_segments'].shape[0])
726	outfile.createVariable('outline_segments', netcdf_int,
727	('num_of_outline_segments',
728	'num_of_segment_ends'))
729	outfile.variables['outline_segments'][:] = mesh['outline_segments']
730	if mesh['outline_segment_tags'].shape[1] > 0:
731	outfile.createDimension('num_of_outline_segment_tag_chars',
732	mesh['outline_segment_tags'].shape[1])
733	outfile.createVariable('outline_segment_tags', netcdf_char,
734	('num_of_outline_segments',
735	'num_of_outline_segment_tag_chars'))
736	outfile.variables['outline_segment_tags'][:] = \
737	mesh['outline_segment_tags']
738
739	# holes
740	if (mesh['holes'].shape[0] > 0):
741	outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
742	outfile.createVariable('holes', netcdf_float,
743	('num_of_holes', 'num_of_dimensions'))
744	outfile.variables['holes'][:] = mesh['holes']
745
746	# regions
747	if (mesh['regions'].shape[0] > 0):
748	outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
749	outfile.createVariable('regions', netcdf_float,
750	('num_of_regions', 'num_of_dimensions'))
751	outfile.createVariable('region_max_areas', netcdf_float,
752	('num_of_regions',))
753	outfile.variables['regions'][:] = mesh['regions']
754	outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
755	if (mesh['region_tags'].shape[1] > 0):
756	outfile.createDimension('num_of_region_tag_chars',
757	mesh['region_tags'].shape[1])
758	outfile.createVariable('region_tags', netcdf_char,
759	('num_of_regions',
760	'num_of_region_tag_chars'))
761	outfile.variables['region_tags'][:] = mesh['region_tags']
762
763	# geo_reference info
764	if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
765	mesh['geo_reference'].write_NetCDF(outfile)
766
767	outfile.close()
768
769
770	def _read_msh_file(file_name):
771	""" Read in an msh file."""
772
773	#Check contents. Get NetCDF
774	fd = open(file_name, 'r')
775	fd.close()
776
777	# throws prints to screen if file not present
778	fid = NetCDFFile(file_name, netcdf_mode_r)
779	mesh = {}
780
781	# Get the variables - the triangulation
782	try:
783	mesh['vertices'] = fid.variables['vertices'][:]
784	except KeyError:
785	mesh['vertices'] = num.array([], num.int) #array default#
786
787	try:
788	mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
789	except KeyError:
790	mesh['vertex_attributes'] = None
791
792	mesh['vertex_attribute_titles'] = []
793	try:
794	titles = fid.variables['vertex_attribute_titles'][:]
795	mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles]
796	except KeyError:
797	pass
798
799	try:
800	mesh['segments'] = fid.variables['segments'][:]
801	except KeyError:
802	mesh['segments'] = num.array([], num.int) #array default#
803
804	mesh['segment_tags'] = []
805	try:
806	tags = fid.variables['segment_tags'][:]
807	mesh['segment_tags'] = [x.tostring().strip() for x in tags]
808	except KeyError:
809	for ob in mesh['segments']:
810	mesh['segment_tags'].append('')
811
812	try:
813	mesh['triangles'] = fid.variables['triangles'][:]
814	mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
815	except KeyError:
816	mesh['triangles'] = num.array([], num.int) #array default#
817	mesh['triangle_neighbors'] = num.array([], num.int) #array default#
818
819	mesh['triangle_tags'] = []
820	try:
821	tags = fid.variables['triangle_tags'][:]
822	mesh['triangle_tags'] = [x.tostring().strip() for x in tags]
823	except KeyError:
824	for ob in mesh['triangles']:
825	mesh['triangle_tags'].append('')
826
827	#the outline
828	try:
829	mesh['points'] = fid.variables['points'][:]
830	except KeyError:
831	mesh['points'] = []
832
833	try:
834	mesh['point_attributes'] = fid.variables['point_attributes'][:]
835	except KeyError:
836	mesh['point_attributes'] = []
837	for point in mesh['points']:
838	mesh['point_attributes'].append([])
839
840	try:
841	mesh['outline_segments'] = fid.variables['outline_segments'][:]
842	except KeyError:
843	mesh['outline_segments'] = num.array([], num.int) #array default#
844
845	mesh['outline_segment_tags'] =[]
846	try:
847	tags = fid.variables['outline_segment_tags'][:]
848	for i, tag in enumerate(tags):
849	mesh['outline_segment_tags'].append(tags[i].tostring().strip())
850	except KeyError:
851	for ob in mesh['outline_segments']:
852	mesh['outline_segment_tags'].append('')
853
854	try:
855	mesh['holes'] = fid.variables['holes'][:]
856	except KeyError:
857	mesh['holes'] = num.array([], num.int) #array default#
858
859	try:
860	mesh['regions'] = fid.variables['regions'][:]
861	except KeyError:
862	mesh['regions'] = num.array([], num.int) #array default#
863
864	mesh['region_tags'] =[]
865	try:
866	tags = fid.variables['region_tags'][:]
867	for i, tag in enumerate(tags):
868	mesh['region_tags'].append(tags[i].tostring().strip())
869	except KeyError:
870	for ob in mesh['regions']:
871	mesh['region_tags'].append('')
872
873	try:
874	mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
875	except KeyError:
876	mesh['region_max_areas'] = num.array([], num.int) #array default#
877
878	try:
879	geo_reference = Geo_reference(NetCDFObject=fid)
880	mesh['geo_reference'] = geo_reference
881	except AttributeError, e:
882	#geo_ref not compulsory
883	mesh['geo_reference'] = None
884
885	fid.close()
886
887	return mesh
888
889
890	def export_boundary_file(file_name, points, title, delimiter=','):
891	"""Export a boundary file.
892
893	Format:
894	First line: Title variable
895	Following lines: [point index][delimiter][point index]
896
897	file_name - the name of the new file
898	points - List of point index pairs [[p1, p2],[p3, p4]..]
899	title - info to write in the first line
900	"""
901
902	fd = open(file_name, 'w')
903
904	fd.write(title + '\n')
905
906	# [point index][delimiter][point index]
907	for point in points:
908	fd.write(str(point[0]) + delimiter + str(point[1]) + '\n')
909
910	fd.close()
911
912	################################################################################
913	# IMPORT/EXPORT POINTS FILES
914	################################################################################
915
916	def extent_point_atts(point_atts):
917	"""Returns 4 points representing the extent
918	This loses attribute info.
919	"""
920
921	point_atts['pointlist'] = extent(point_atts['pointlist'])
922	point_atts['attributelist'] = {}
923
924	return point_atts
925
926
927	def extent(points):
928	points = num.array(points, num.float)
929
930	max_x = min_x = points[0][0]
931	max_y = min_y = points[0][1]
932
933	for point in points[1:]:
934	x = point[0]
935	if x > max_x:
936	max_x = x
937	if x < min_x:
938	min_x = x
939
940	y = point[1]
941	if y > max_y:
942	max_y = y
943	if y < min_y:
944	min_y = y
945
946	extent = num.array([[min_x, min_y],
947	[max_x, min_y],
948	[max_x, max_y],
949	[min_x, max_y]])
950
951	return extent
952
953
954	def reduce_pts(infile, outfile, max_points, verbose = False):
955	"""Reduce a points file until less than given size.
956
957	Reduces a points file by removing every second point until the # of points
958	is less than max_points.
959	"""
960
961	# check out pts2rectangular in least squares, and the use of reduction.
962	# Maybe it does the same sort of thing?
963	point_atts = _read_pts_file(infile)
964
965	while point_atts['pointlist'].shape[0] > max_points:
966	if verbose: log.critical("point_atts['pointlist'].shape[0]")
967	point_atts = half_pts(point_atts)
968
969	export_points_file(outfile, point_atts)
970
971
972	def produce_half_point_files(infile, max_points, delimiter, verbose=False):
973	point_atts = _read_pts_file(infile)
974	root, ext = splitext(infile)
975	outfiles = []
976
977	if verbose: log.critical("# of points", point_atts['pointlist'].shape[0])
978
979	while point_atts['pointlist'].shape[0] > max_points:
980	point_atts = half_pts(point_atts)
981
982	if verbose: log.critical("# of points = %s"
983	% str(point_atts['pointlist'].shape[0]))
984
985	outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
986	outfiles.append(outfile)
987	export_points_file(outfile, point_atts)
988
989	return outfiles
990
991
992	def point_atts2array(point_atts):
993	# convert attribute list to array of floats
994	point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
995
996	for key in point_atts['attributelist'].keys():
997	point_atts['attributelist'][key] = \
998	num.array(point_atts['attributelist'][key], num.float)
999
1000	return point_atts
1001
1002
1003	def half_pts(point_atts):
1004	point_atts2array(point_atts)
1005	point_atts['pointlist'] = point_atts['pointlist'][::2]
1006
1007	for key in point_atts['attributelist'].keys():
1008	point_atts['attributelist'][key] = point_atts['attributelist'][key][::2]
1009
1010	return point_atts
1011
1012	def concatinate_attributelist(dic):
1013	"""
1014	giving a dic[attribute title] = attribute
1015	return list of attribute titles, array of attributes
1016	"""
1017
1018	point_attributes = num.array([], num.float)
1019	keys = dic.keys()
1020	key = keys.pop(0)
1021	point_attributes = num.reshape(dic[key], (dic[key].shape[0], 1))
1022	for key in keys:
1023	reshaped = num.reshape(dic[key], (dic[key].shape[0], 1))
1024	point_attributes = num.concatenate([point_attributes, reshaped], axis=1)
1025
1026	return dic.keys(), point_attributes
1027
1028
1029	def take_points(dict, indices_to_keep):
1030	dict = point_atts2array(dict)
1031	# FIXME maybe the points data structure should become a class?
1032	dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0)
1033
1034	for key in dict['attributelist'].keys():
1035	dict['attributelist'][key] = num.take(dict['attributelist'][key],
1036	indices_to_keep, axis=0)
1037
1038	return dict
1039
1040	def add_point_dictionaries(dict1, dict2):
1041	"""
1042	"""
1043
1044	dict1 = point_atts2array(dict1)
1045	dict2 = point_atts2array(dict2)
1046
1047	combined = {}
1048	combined['pointlist'] = num.concatenate((dict2['pointlist'],
1049	dict1['pointlist']), axis=0)
1050
1051	atts = {}
1052	for key in dict2['attributelist'].keys():
1053	atts[key]= num.concatenate((dict2['attributelist'][key],
1054	dict1['attributelist'][key]), axis=0)
1055	combined['attributelist'] = atts
1056	combined['geo_reference'] = dict1['geo_reference']
1057
1058	return combined
1059
1060
1061	if __name__ == "__main__":
1062	pass

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

Download in other formats: