Context Navigation

loadASCII.py @ 3381

Last change on this file since 3381 was 3133, checked in by duncan, 19 years ago
small change
File size: 47.7 KB

Rev	Line
[2253]	1	"""
	2	The format for a .xya file is:
	3	1st line: [attribute names]
	4	other lines: x y [attributes]
	5
	6	for example:
	7	elevation, friction
	8	0.6, 0.7, 4.9, 0.3
	9	1.9, 2.8, 5, 0.3
	10	2.7, 2.4, 5.2, 0.3
	11
	12	The first two columns are always implicitly assumed to be x, y coordinates.
	13	Use the same delimiter for the attribute names and the data
	14
	15	An xya file can optionally end with
	16	#geo reference
	17	56
	18	466600.0
	19	8644444.0
	20
	21	When the 1st # is the zone,
	22	2nd # the xllcorner and
	23	3rd # the yllcorner
	24	The format for a Points dictionary is:
	25
	26	['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
	27	['attributelist'], a dictionary of 1D arrays, representing attribute values
	28	at the point. The dictionary key is the attribute header.
[2257]	29	['geo_reference'] a Geo_refernece object. Use if the point information
	30	is relative. This is optional.
[2253]	31	eg
	32	dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
	33	dic['attributelist']['elevation'] = [[7.0,5.0]
	34
	35
	36	The dict format for IO with mesh files is;
	37	(the triangulation)
	38	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
	39	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
	40	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
	41	segments: [[v1,v2],[v3,v4],...] (lists of integers)
	42	segment_tags : [tag,tag,...] list of strings
	43	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
	44	triangle_tags: [s1,s2,...] A list of list of strings (probably not neccecary. a list of string should be ok)
	45	triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
	46
	47	(the outline)
	48	points: [[x1,y1],[x2,y2],...] (lists of doubles)
	49	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
	50	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
	51	outline_segment_tags : [tag1,tag2,...] list of strings
	52	holes : [[x1,y1],...](List of doubles, one inside each hole region)
	53	regions : [ [x1,y1],...] (List of 4 doubles)
	54	region_tags : [tag1,tag2,...] (list of strings)
	55	region_max_areas: [ma1,ma2,...] (A list of doubles)
	56	{Convension: A -ve max area means no max area}
[2257]	57	geo_reference: a Geo_refernece object. Use if the point information
	58	is relative. This is optional.
[2253]	59
	60	Points files are .xya for ascii and .pts for NetCDF
	61	Mesh files are .tsh for ascii and .msh for NetCDF
	62
	63	Note: point att's have no titles - that's currently ok, since least
	64	squares adds the point info to vertices, and they have titles
	65
	66	The format for a .tsh file is (FIXME update this)
	67
	68	First line: <# of vertices> <# of attributes>
	69	Following lines: <vertex #> <x> <y> [attributes]
	70	One line: <# of triangles>
	71	Following lines: <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
	72	One line: <# of segments>
	73	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
	74	"""
	75	##FIXME (DSG-DSG) Is the dict format mentioned above a list of a numeric array?
	76	# Needs to be defined
	77	##FIXME (DSG-DSG) if the ascii file isn't .xya give a better error message.
	78
	79
	80	from string import find, rfind
	81	from Numeric import array, Float, Int16, Int32, Character,reshape, concatenate, take
	82	from os.path import splitext
	83
[2941]	84	from coordinate_transforms.geo_reference import Geo_reference,TITLE, TitleError
[2253]	85
[2993]	86	from Scientific.IO.NetCDF import NetCDFFile
	87
[2253]	88	import exceptions
	89	class TitleAmountError(exceptions.Exception): pass
	90
	91	NOMAXAREA=-999
	92
	93	# This is used by pmesh
	94	def import_mesh_file(ofile):
	95	"""
	96	read a mesh file, either .tsh or .msh
	97
	98	Note: will throw an IOError if it can't load the file.
	99	Catch these!
	100	"""
	101	try:
	102	if ofile[-4:]== ".tsh":
	103	dict = _read_tsh_file(ofile)
	104	elif ofile[-4:]== ".msh":
	105	dict = _read_msh_file(ofile)
	106	else:
[3133]	107	msg = 'Extension .%s is unknown' %ofile[-4:]
[2253]	108	raise IOError, msg
[2941]	109	except (TitleError,SyntaxError,IndexError, ValueError): #FIXME No test for ValueError
[2253]	110	msg = 'File could not be opened'
	111	raise IOError, msg
	112	return dict
	113
	114
	115	def export_mesh_file(ofile,mesh_dict):
	116	"""
	117	write a file, ofile, with the format
	118
	119	First line: <# of vertices> <# of attributes>
	120	Following lines: <vertex #> <x> <y> [attributes]
	121	One line: <# of triangles>
	122	Following lines: <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
	123	One line: <# of segments>
	124	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
	125	"""
	126	#FIXME DSG-anyone: automate
	127	if not mesh_dict.has_key('points'):
	128	mesh_dict['points'] = []
	129	if not mesh_dict.has_key('point_attributes'):
	130	mesh_dict['point_attributes'] = []
	131	if not mesh_dict.has_key('outline_segments'):
	132	mesh_dict['outline_segments'] = []
	133	if not mesh_dict.has_key('outline_segment_tags'):
	134	mesh_dict['outline_segment_tags'] = []
	135	if not mesh_dict.has_key('holes'):
	136	mesh_dict['holes'] = []
	137	if not mesh_dict.has_key('regions'):
	138	mesh_dict['regions'] = []
	139
	140	if not mesh_dict.has_key('region_tags'):
	141	mesh_dict['region_tags'] = []
	142	if not mesh_dict.has_key('region_max_areas'):
	143	mesh_dict['region_max_areas'] = []
	144	if not mesh_dict.has_key('vertices'):
	145	mesh_dict['vertices'] = []
	146	if not mesh_dict.has_key('vertex_attributes'):
	147	mesh_dict['vertex_attributes'] = []
	148	if not mesh_dict.has_key('vertex_attribute_titles'):
	149	mesh_dict['vertex_attribute_titles'] = []
	150	if not mesh_dict.has_key('segments'):
	151	mesh_dict['segments'] = []
	152	if not mesh_dict.has_key('segment_tags'):
	153	mesh_dict['segment_tags'] = []
	154	if not mesh_dict.has_key('triangles'):
	155	mesh_dict['triangles'] = []
	156	if not mesh_dict.has_key('triangle_tags'):
	157	mesh_dict['triangle_tags'] = []
	158	if not mesh_dict.has_key('triangle_neighbors'):
	159	mesh_dict['triangle_neighbors'] = []
	160	#print "DSG************"
	161	#print "mesh_dict",mesh_dict
	162	#print "DSG************"
	163
	164	if (ofile[-4:] == ".tsh"):
	165	_write_tsh_file(ofile,mesh_dict)
	166	elif (ofile[-4:] == ".msh"):
	167	_write_msh_file(ofile, mesh_dict)
	168	else:
	169	msg = 'Unknown file type %s ' %ofile
	170	raise IOError, msg
	171
	172	def _read_tsh_file(ofile):
	173	"""
	174	Read the text file format for meshes
	175	"""
	176	fd = open(ofile,'r')
	177	dict = _read_triangulation(fd)
	178	dict_mesh = _read_outline(fd)
	179	for element in dict_mesh.keys():
	180	dict[element] = dict_mesh[element]
	181	fd.close()
	182	return dict
	183
	184
	185	def _read_triangulation(fd):
	186	"""
	187	Read the generated triangulation, NOT the outline.
	188	"""
	189	delimiter = " "
	190	######### loading the point info
	191	line = fd.readline()
	192	#print line
	193	fragments = line.split()
	194	#for fragment in fragments:
	195	#print fragment
	196	if fragments ==[]:
	197	NumOfVertices = 0
	198	NumOfVertAttributes = 0
	199	else:
	200	NumOfVertices = fragments[0]
	201	NumOfVertAttributes = fragments[1]
	202	points = []
	203	pointattributes = []
	204	for index in range(int(NumOfVertices)):
	205	#print index
	206	fragments = fd.readline().split()
	207	#print fragments
	208	fragments.pop(0) #pop off the index
	209
	210	# pop the x y off so we're left with a list of attributes
	211	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
	212	points.append(vert)
	213	apointattributes = []
	214	#print fragments
	215	for fragment in fragments:
	216	apointattributes.append(float(fragment))
	217	pointattributes.append(apointattributes)
	218
	219	######### loading the point title info
	220	line = fd.readline()
	221	#print "point title comments",line
	222	vertTitle = []
	223	for index in range(int(NumOfVertAttributes)):
	224	#print index
	225	fragments = fd.readline().strip()
	226	vertTitle.append(fragments)
	227	#print vertTitle
	228
	229	######### loading the triangle info
	230	line = fd.readline()
	231	#print "triangle comments",line
	232	fragments = line.split()
	233	#for fragment in fragments:
	234	# print fragment
	235	NumOfTriangles = fragments[0]
	236	triangles = []
	237	triangleattributes = []
	238	triangleneighbors = []
	239	for index in range(int(NumOfTriangles)):
	240	#print index
	241	line = fd.readline()
	242	line.strip() # so we can get the region string
	243	fragments = line.split()
	244	#print "triangle info", fragments
	245	fragments.pop(0) #pop off the index
	246
	247	tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
	248	triangles.append(tri)
	249	neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
	250	triangleneighbors.append(neighbors)
	251	for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
	252	line = line[find(line,delimiter):] # remove index
	253	line = line.lstrip()
	254	stringtag = line.strip()
	255	triangleattributes.append(stringtag)
	256
	257	######### loading the segment info
	258	line = fd.readline()
	259	#print "seg comment line",line
	260	fragments = line.split()
	261	#for fragment in fragments:
	262	# print fragment
	263	NumOfSegments = fragments[0]
	264	segments = []
	265	segmenttags = []
	266	for index in range(int(NumOfSegments)):
	267	#print index
	268	line = fd.readline()
	269	line.strip() # to get the segment string
	270	fragments = line.split()
	271	#print fragments
	272	fragments.pop(0) #pop off the index
	273	seg = [int(fragments[0]),int(fragments[1])]
	274	segments.append(seg)
	275	line = line[find(line,delimiter):] # remove index
	276	line = line.lstrip()
	277	line = line[find(line,delimiter):] # remove x
	278	line = line.lstrip()
	279	line = line[find(line,delimiter):] # remove y
	280	stringtag = line.strip()
	281	segmenttags.append(stringtag)
	282
	283
	284	meshDict = {}
	285	meshDict['vertices'] = points
	286	meshDict['vertex_attributes'] = pointattributes
	287	meshDict['triangles'] = triangles
	288	meshDict['triangle_tags'] = triangleattributes
	289	meshDict['triangle_neighbors'] = triangleneighbors
	290	meshDict['segments'] = segments
	291	meshDict['segment_tags'] = segmenttags
	292	meshDict['vertex_attribute_titles'] = vertTitle
	293
	294	return meshDict
	295
	296	def _read_outline(fd):
	297	"""
	298	Note, if a file has no mesh info, it can still be read - the meshdic
	299	returned will be 'empty'.
	300	"""
	301	delimiter = " " # warning: split() calls are using default whitespace
	302
	303	######### loading the point info
	304	line = fd.readline()
	305	#print 'read_outline - line',line
	306	fragments = line.split()
	307	#for fragment in fragments:
	308	if fragments ==[]:
	309	NumOfVertices = 0
	310	NumOfVertAttributes = 0
	311	else:
	312	NumOfVertices = fragments[0]
	313	NumOfVertAttributes = fragments[1]
	314	points = []
	315	pointattributes = []
	316	for index in range(int(NumOfVertices)):
	317	#print index
	318	fragments = fd.readline().split()
	319	#print fragments
	320	fragments.pop(0) #pop off the index
	321	# pop the x y off so we're left with a list of attributes
	322	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
	323	points.append(vert)
	324	apointattributes = []
	325	#print fragments
	326	for fragment in fragments:
	327	apointattributes.append(float(fragment))
	328	pointattributes.append(apointattributes)
	329
	330
	331	######### loading the segment info
	332	line = fd.readline()
	333	fragments = line.split()
	334	#for fragment in fragments:
	335	# print fragment
	336	if fragments ==[]:
	337	NumOfSegments = 0
	338	else:
	339	NumOfSegments = fragments[0]
	340	segments = []
	341	segmenttags = []
	342	for index in range(int(NumOfSegments)):
	343	#print index
	344	line = fd.readline()
	345	fragments = line.split()
	346	#print fragments
	347	fragments.pop(0) #pop off the index
	348
	349	seg = [int(fragments[0]),int(fragments[1])]
	350	segments.append(seg)
	351	line = line[find(line,delimiter):] # remove index
	352	line = line.lstrip()
	353	line = line[find(line,delimiter):] # remove x
	354	line = line.lstrip()
	355	line = line[find(line,delimiter):] # remove y
	356	stringtag = line.strip()
	357	segmenttags.append(stringtag)
	358
	359	######### loading the hole info
	360	line = fd.readline()
	361	#print line
	362	fragments = line.split()
	363	#for fragment in fragments:
	364	# print fragment
	365	if fragments ==[]:
	366	numOfHoles = 0
	367	else:
	368	numOfHoles = fragments[0]
	369	holes = []
	370	for index in range(int(numOfHoles)):
	371	#print index
	372	fragments = fd.readline().split()
	373	#print fragments
	374	fragments.pop(0) #pop off the index
	375	hole = [float(fragments[0]),float(fragments[1])]
	376	holes.append(hole)
	377
	378
	379	######### loading the region info
	380	line = fd.readline()
	381	#print line
	382	fragments = line.split()
	383	#for fragment in fragments:
	384	# print fragment
	385	if fragments ==[]:
	386	numOfRegions = 0
	387	else:
	388	numOfRegions = fragments[0]
	389	regions = []
	390	regionattributes = []
	391	for index in range(int(numOfRegions)):
	392	line = fd.readline()
	393	fragments = line.split()
	394	#print fragments
	395	fragments.pop(0) #pop off the index
	396	region = [float(fragments[0]),float(fragments[1])]
	397	regions.append(region)
	398
	399	line = line[find(line,delimiter):] # remove index
	400	line = line.lstrip()
	401	line = line[find(line,delimiter):] # remove x
	402	line = line.lstrip()
	403	line = line[find(line,delimiter):] # remove y
	404	stringtag = line.strip()
	405	regionattributes.append(stringtag)
	406	regionmaxareas = []
	407	line = fd.readline()
	408	for index in range(int(numOfRegions)): # Read in the Max area info
	409	line = fd.readline()
	410	fragments = line.split()
	411	#print fragments
	412	# The try is here for format compatibility
	413	try:
	414	fragments.pop(0) #pop off the index
	415	if len(fragments) == 0: #no max area
	416	regionmaxareas.append(None)
	417	else:
	418	regionmaxareas.append(float(fragments[0]))
	419	except IndexError, e:
	420	regionmaxareas.append(None)
	421
	422	try:
	423	geo_reference = Geo_reference(ASCIIFile=fd)
	424	except:
	425	#geo_ref not compulsory
	426	geo_reference = None
	427
	428	meshDict = {}
	429	meshDict['points'] = points
	430	meshDict['point_attributes'] = pointattributes
	431	meshDict['outline_segments'] = segments
	432	meshDict['outline_segment_tags'] = segmenttags
	433	meshDict['holes'] = holes
	434	meshDict['regions'] = regions
	435	meshDict['region_tags'] = regionattributes
	436	meshDict['region_max_areas'] = regionmaxareas
	437	meshDict['geo_reference'] = geo_reference
	438
	439	#print "meshDict",meshDict
	440	return meshDict
	441
	442	def clean_line(line,delimiter):
	443	"""Remove whitespace
	444	"""
	445	#print ">%s" %line
	446	line = line.strip()
	447	#print "stripped>%s" %line
	448	numbers = line.split(delimiter)
	449	i = len(numbers) - 1
	450	while i >= 0:
	451	if numbers[i] == '':
	452	numbers.pop(i)
	453	else:
	454	numbers[i] = numbers[i].strip()
	455
	456	i += -1
	457	#for num in numbers:
	458	# print "num>%s<" %num
	459	return numbers
	460
	461	def _write_ASCII_triangulation(fd,
	462	gen_dict):
	463	vertices = gen_dict['vertices']
	464	vertices_attributes = gen_dict['vertex_attributes']
	465	try:
	466	vertices_attribute_titles = gen_dict['vertex_attribute_titles']
	467	except KeyError, e:
	468	#FIXME is this the best way?
	469	if vertices_attributes == [] or vertices_attributes[0] == []:
	470	vertices_attribute_titles = []
	471	else:
	472	raise KeyError, e
	473
	474	triangles = gen_dict['triangles']
	475	triangles_attributes = gen_dict['triangle_tags']
	476	triangle_neighbors = gen_dict['triangle_neighbors']
	477
	478	segments = gen_dict['segments']
	479	segment_tags = gen_dict['segment_tags']
	480
	481	numVert = str(len(vertices))
	482	if (numVert == "0") or len(vertices_attributes) == 0:
	483	numVertAttrib = "0"
	484	else:
	485	numVertAttrib = str(len(vertices_attributes[0]))
	486	fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
	487
	488	#<vertex #> <x> <y> [attributes]
	489	index = 0
	490	for vert in vertices:
	491	attlist = ""
	492	for att in vertices_attributes[index]:
	493	attlist = attlist + str(att)+" "
	494	attlist.strip()
	495	fd.write(str(index) + " "
	496	+ str(vert[0]) + " "
	497	+ str(vert[1]) + " "
	498	+ attlist + "\n")
	499	index += 1
	500
	501	# write comments for title
	502	fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
	503	for title in vertices_attribute_titles:
	504	fd.write(title + "\n")
	505
	506	#<# of triangles>
	507	n = len(triangles)
	508	fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
	509
	510	# <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
	511	for index in range(n):
	512	neighbors = ""
	513	tri = triangles[index]
	514	for neighbor in triangle_neighbors[index]:
	515	if neighbor:
	516	neighbors += str(neighbor) + " "
	517	else:
	518	if neighbor == 0:
	519	neighbors += "0 "
	520	else:
	521	neighbors += "-1 "
	522	#Warning even though a list is past, only the first value
	523	#is written. There's an assumption that the list only
	524	# contains one item. This assumption is made since the
	525	#dict that's being passed around is also be used to communicate
	526	# with triangle, and it seems to have the option of returning
	527	# more than one value for triangle attributex
	528	if triangles_attributes[index] == ['']:
	529	att = ""
	530	else:
	531	att = str(triangles_attributes[index])
	532	fd.write(str(index) + " "
	533	+ str(tri[0]) + " "
	534	+ str(tri[1]) + " "
	535	+ str(tri[2]) + " "
	536	+ neighbors + " "
	537	+ att + "\n")
	538
	539	#One line: <# of segments>
	540	fd.write(str(len(segments)) +
	541	" # <# of segments>, next lines <segment #> <vertex #> <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
	542
	543	#Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
	544	for i in range(len(segments)):
	545	seg = segments[i]
	546	fd.write(str(i) + " "
	547	+ str(seg[0]) + " "
	548	+ str(seg[1]) + " "
	549	+ str(segment_tags[i]) + "\n")
	550
	551
	552	def _write_tsh_file(ofile,mesh_dict):
	553	fd = open(ofile,'w')
	554	_write_ASCII_triangulation(fd,mesh_dict)
	555	_write_ASCII_outline(fd,mesh_dict)
	556	fd.close()
	557
	558	def _write_ASCII_outline(fd,
	559	dict):
	560	points = dict['points']
	561	point_attributes = dict['point_attributes']
	562	segments = dict['outline_segments']
	563	segment_tags = dict['outline_segment_tags']
	564	holes = dict['holes']
	565	regions = dict['regions']
	566	region_tags = dict['region_tags']
	567	region_max_areas = dict['region_max_areas']
	568
	569	num_points = str(len(points))
	570	if (num_points == "0"):
	571	num_point_atts = "0"
	572	else:
	573	num_point_atts = str(len(point_attributes[0]))
	574	fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
	575
	576	# <x> <y> [attributes]
	577	for i,point in enumerate(points):
	578	attlist = ""
	579	for att in point_attributes[i]:
	580	attlist = attlist + str(att)+" "
	581	attlist.strip()
	582	fd.write(str(i) + " "
	583	+str(point[0]) + " "
	584	+ str(point[1]) + " "
	585	+ attlist + "\n")
	586
	587	#One line: <# of segments>
	588	fd.write(str(len(segments)) +
	589	" # <# of segments>, next lines <segment #> <vertex #> <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
	590
	591	#Following lines: <vertex #> <vertex #> [boundary tag]
	592	for i,seg in enumerate(segments):
	593	fd.write(str(i) + " "
	594	+ str(seg[0]) + " "
	595	+ str(seg[1]) + " "
	596	+ str(segment_tags[i]) + "\n")
	597
	598	#One line: <# of holes>
	599	fd.write(str(len(holes)) +
	600	" # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")
	601	# <x> <y>
	602	for i,h in enumerate(holes):
	603	fd.write(str(i) + " "
	604	+ str(h[0]) + " "
	605	+ str(h[1]) + "\n")
	606
	607	#One line: <# of regions>
	608	fd.write(str(len(regions)) +
	609	" # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")
	610	# <index> <x> <y> <tag>
	611	#print "regions",regions
	612	for i,r in enumerate(regions):
	613	#print "r",r
	614	fd.write(str(i) + " "
	615	+ str(r[0]) + " "
	616	+ str(r[1])+ " "
	617	+ str(region_tags[i]) + "\n")
	618
	619	# <index> [<MaxArea>\|""]
	620
	621	#One line: <# of regions>
	622	fd.write(str(len(regions)) +
	623	" # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n")
	624	for i,r in enumerate(regions):
	625	area = str(region_max_areas[i])
	626
	627	fd.write(str(i) + " " + area + "\n")
	628
	629	# geo_reference info
	630	if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
	631	dict['geo_reference'].write_ASCII(fd)
	632
	633	def _write_msh_file(file_name, mesh):
	634	"""Write .msh NetCDF file
	635
	636	WARNING: This function mangles the mesh data structure
	637	"""
	638
	639
	640	IntType = Int32
	641
	642	#
	643	#the triangulation
	644	mesh['vertices'] = array(mesh['vertices']).astype(Float)
	645	mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
	646	mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character)
	647	mesh['segments'] = array(mesh['segments']).astype(IntType)
	648	mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
	649	mesh['triangles'] = array(mesh['triangles']).astype(IntType)
	650	mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
	651	mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType)
	652
	653	#the outline
	654	mesh['points'] = array(mesh['points']).astype(Float)
	655	mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
	656	mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
	657	mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
	658	mesh['holes'] = array(mesh['holes']).astype(Float)
	659	mesh['regions'] = array(mesh['regions']).astype(Float)
	660	mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
	661	mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
	662
	663	#mesh = mesh_dict2array(mesh)
	664	#print "new_mesh",new_mesh
	665	#print "mesh",mesh
	666
	667	# NetCDF file definition
	668	try:
	669	outfile = NetCDFFile(file_name, 'w')
	670	except IOError:
	671	msg = 'File %s could not be created' %file_name
	672	raise msg
	673
	674	#Create new file
	675	outfile.institution = 'Geoscience Australia'
	676	outfile.description = 'NetCDF format for compact and portable storage ' +\
	677	'of spatial point data'
	678
	679	# dimension definitions - fixed
	680	outfile.createDimension('num_of_dimensions', 2) #This is 2d data
	681	outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
	682	outfile.createDimension('num_of_triangle_vertices', 3)
	683	outfile.createDimension('num_of_triangle_faces', 3)
	684	outfile.createDimension('num_of_region_max_area', 1)
	685
	686	# Create dimensions, variables and set the variables
	687
	688	# trianglulation
	689	# vertices
	690	if (mesh['vertices'].shape[0] > 0):
	691	outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
	692	outfile.createVariable('vertices', Float, ('num_of_vertices',
	693	'num_of_dimensions'))
	694	outfile.variables['vertices'][:] = mesh['vertices']
	695	if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
	696	outfile.createDimension('num_of_vertex_attributes',
	697	mesh['vertex_attributes'].shape[1])
	698	outfile.createDimension('num_of_vertex_attribute_title_chars',
	699	mesh['vertex_attribute_titles'].shape[1])
	700	outfile.createVariable('vertex_attributes',
	701	Float,
	702	('num_of_vertices',
	703	'num_of_vertex_attributes'))
	704	outfile.createVariable('vertex_attribute_titles',
	705	Character,
	706	( 'num_of_vertex_attributes',
	707	'num_of_vertex_attribute_title_chars' ))
	708	outfile.variables['vertex_attributes'][:] = \
	709	mesh['vertex_attributes']
	710	outfile.variables['vertex_attribute_titles'][:] = \
	711	mesh['vertex_attribute_titles']
	712	# segments
	713	if (mesh['segments'].shape[0] > 0):
	714	outfile.createDimension('num_of_segments',
	715	mesh['segments'].shape[0])
	716	outfile.createVariable('segments',
	717	IntType,
	718	('num_of_segments', 'num_of_segment_ends'))
	719	outfile.variables['segments'][:] = mesh['segments']
	720	if (mesh['segment_tags'].shape[1] > 0):
	721	outfile.createDimension('num_of_segment_tag_chars',
	722	mesh['segment_tags'].shape[1])
	723	outfile.createVariable('segment_tags',
	724	Character,
	725	('num_of_segments',
	726	'num_of_segment_tag_chars'))
	727	outfile.variables['segment_tags'][:] = mesh['segment_tags']
	728	# triangles
	729	if (mesh['triangles'].shape[0] > 0):
	730	outfile.createDimension('num_of_triangles',
	731	mesh['triangles'].shape[0])
	732	outfile.createVariable('triangles',
	733	IntType,
	734	('num_of_triangles',
	735	'num_of_triangle_vertices'))
	736	outfile.createVariable('triangle_neighbors',
	737	IntType,
	738	('num_of_triangles',
	739	'num_of_triangle_faces'))
	740	outfile.variables['triangles'][:] = mesh['triangles']
	741	outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
	742	if (mesh['triangle_tags'].shape[1] > 0):
	743	outfile.createDimension('num_of_triangle_tag_chars',
	744	mesh['triangle_tags'].shape[1])
	745	outfile.createVariable('triangle_tags',
	746	Character,
	747	('num_of_triangles',
	748	'num_of_triangle_tag_chars'))
	749	outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
	750
	751
	752	# outline
	753	# points
	754	if (mesh['points'].shape[0] > 0):
	755	outfile.createDimension('num_of_points', mesh['points'].shape[0])
	756	outfile.createVariable('points', Float, ('num_of_points',
	757	'num_of_dimensions'))
	758	outfile.variables['points'][:] = mesh['points']
	759	if (mesh['point_attributes'].shape[0] > 0 and mesh['point_attributes'].shape[1] > 0):
	760	outfile.createDimension('num_of_point_attributes',
	761	mesh['point_attributes'].shape[1])
	762	outfile.createVariable('point_attributes',
	763	Float,
	764	('num_of_points',
	765	'num_of_point_attributes'))
	766	outfile.variables['point_attributes'][:] = mesh['point_attributes']
	767	# outline_segments
	768	if (mesh['outline_segments'].shape[0] > 0):
	769	outfile.createDimension('num_of_outline_segments',
	770	mesh['outline_segments'].shape[0])
	771	outfile.createVariable('outline_segments',
	772	IntType,
	773	('num_of_outline_segments',
	774	'num_of_segment_ends'))
	775	outfile.variables['outline_segments'][:] = mesh['outline_segments']
	776	if (mesh['outline_segment_tags'].shape[1] > 0):
	777	outfile.createDimension('num_of_outline_segment_tag_chars',
	778	mesh['outline_segment_tags'].shape[1])
	779	outfile.createVariable('outline_segment_tags',
	780	Character,
	781	('num_of_outline_segments',
	782	'num_of_outline_segment_tag_chars'))
	783	outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
	784	# holes
	785	if (mesh['holes'].shape[0] > 0):
	786	outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
	787	outfile.createVariable('holes', Float, ('num_of_holes',
	788	'num_of_dimensions'))
	789	outfile.variables['holes'][:] = mesh['holes']
	790	# regions
	791	if (mesh['regions'].shape[0] > 0):
	792	outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
	793	outfile.createVariable('regions', Float, ('num_of_regions',
	794	'num_of_dimensions'))
	795	outfile.createVariable('region_max_areas',
	796	Float,
	797	('num_of_regions',))
	798	outfile.variables['regions'][:] = mesh['regions']
	799	outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
	800	if (mesh['region_tags'].shape[1] > 0):
	801	outfile.createDimension('num_of_region_tag_chars',
	802	mesh['region_tags'].shape[1])
	803	outfile.createVariable('region_tags',
	804	Character,
	805	('num_of_regions',
	806	'num_of_region_tag_chars'))
	807	outfile.variables['region_tags'][:] = mesh['region_tags']
	808
	809	# geo_reference info
	810	if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
	811	mesh['geo_reference'].write_NetCDF(outfile)
	812
	813	outfile.close()
	814
	815
	816
	817	def _read_msh_file(file_name):
	818	"""
	819	Read in an msh file.
	820
	821	"""
[2993]	822
[2253]	823
	824	#Check contents
	825	#Get NetCDF
	826
	827	# see if the file is there. Throw a QUIET IO error if it isn't
	828	fd = open(file_name,'r')
	829	fd.close()
	830
	831	#throws prints to screen if file not present
	832	fid = NetCDFFile(file_name, 'r')
	833
	834	mesh = {}
	835	# Get the variables
	836	# the triangulation
	837	try:
	838	mesh['vertices'] = fid.variables['vertices'][:]
	839	except KeyError:
	840	mesh['vertices'] = array([])
	841	try:
	842	mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
	843	except KeyError:
	844	mesh['vertex_attributes'] = []
	845	for ob in mesh['vertices']:
	846	mesh['vertex_attributes'].append([])
	847	mesh['vertex_attribute_titles'] = []
	848	try:
	849	titles = fid.variables['vertex_attribute_titles'][:]
	850	for i, title in enumerate(titles):
	851	mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
	852	except KeyError:
	853	pass
	854	try:
	855	mesh['segments'] = fid.variables['segments'][:]
	856	except KeyError:
	857	mesh['segments'] = array([])
	858	mesh['segment_tags'] =[]
	859	try:
	860	tags = fid.variables['segment_tags'][:]
	861	for i, tag in enumerate(tags):
	862	mesh['segment_tags'].append(tags[i].tostring().strip())
	863	except KeyError:
	864	for ob in mesh['segments']:
	865	mesh['segment_tags'].append('')
	866	try:
	867	mesh['triangles'] = fid.variables['triangles'][:]
	868	mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
	869	except KeyError:
	870	mesh['triangles'] = array([])
	871	mesh['triangle_neighbors'] = array([])
	872	mesh['triangle_tags'] =[]
	873	try:
	874	tags = fid.variables['triangle_tags'][:]
	875	for i, tag in enumerate(tags):
	876	mesh['triangle_tags'].append(tags[i].tostring().strip())
	877	except KeyError:
	878	for ob in mesh['triangles']:
	879	mesh['triangle_tags'].append('')
	880
	881	#the outline
	882	try:
	883	mesh['points'] = fid.variables['points'][:]
	884	except KeyError:
	885	mesh['points'] = []
	886	try:
	887	mesh['point_attributes'] = fid.variables['point_attributes'][:]
	888	except KeyError:
	889	mesh['point_attributes'] = []
	890	for point in mesh['points']:
	891	mesh['point_attributes'].append([])
	892	try:
	893	mesh['outline_segments'] = fid.variables['outline_segments'][:]
	894	except KeyError:
	895	mesh['outline_segments'] = array([])
	896	mesh['outline_segment_tags'] =[]
	897	try:
	898	tags = fid.variables['outline_segment_tags'][:]
	899	for i, tag in enumerate(tags):
	900	mesh['outline_segment_tags'].append(tags[i].tostring().strip())
	901	except KeyError:
	902	for ob in mesh['outline_segments']:
	903	mesh['outline_segment_tags'].append('')
	904	try:
	905	mesh['holes'] = fid.variables['holes'][:]
	906	except KeyError:
	907	mesh['holes'] = array([])
	908	try:
	909	mesh['regions'] = fid.variables['regions'][:]
	910	except KeyError:
	911	mesh['regions'] = array([])
	912	mesh['region_tags'] =[]
	913	try:
	914	tags = fid.variables['region_tags'][:]
	915	for i, tag in enumerate(tags):
	916	mesh['region_tags'].append(tags[i].tostring().strip())
	917	except KeyError:
	918	for ob in mesh['regions']:
	919	mesh['region_tags'].append('')
	920	try:
	921	mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
	922	except KeyError:
	923	mesh['region_max_areas'] = array([])
	924	#mesh[''] = fid.variables[''][:]
	925
	926	try:
	927	geo_reference = Geo_reference(NetCDFObject=fid)
	928	mesh['geo_reference'] = geo_reference
	929	except AttributeError, e:
	930	#geo_ref not compulsory
	931	mesh['geo_reference'] = None
	932
	933	fid.close()
	934
	935	return mesh
	936
	937
	938	## used by alpha shapes
	939
	940	def export_boundary_file( file_name, points, title, delimiter = ','):
	941	"""
	942	export a file, ofile, with the format
	943
	944	First line: Title variable
	945	Following lines: [point index][delimiter][point index]
	946
	947	file_name - the name of the new file
	948	points - List of point index pairs [[p1, p2],[p3, p4]..]
	949	title - info to write in the first line
	950	"""
	951
	952	fd = open(file_name,'w')
	953
	954	fd.write(title+"\n")
	955	#[point index][delimiter][point index]
	956	for point in points:
	957	fd.write( str(point[0]) + delimiter
	958	+ str(point[1]) + "\n")
	959	fd.close()
	960
	961
	962	###
	963	# IMPORT/EXPORT POINTS FILES
	964	###
	965
[2480]	966	def export_points_file(ofile, point_dict):
[2253]	967	"""
[2558]	968	write a points file, ofile, as a text (.xya) or binary (.pts) file
[2257]	969
	970	ofile is the file name, including the extension
	971
	972	The point_dict is defined at the top of this file.
[2253]	973	"""
	974	#this was done for all keys in the mesh file.
	975	#if not mesh_dict.has_key('points'):
	976	# mesh_dict['points'] = []
	977	if (ofile[-4:] == ".xya"):
	978	_write_xya_file(ofile, point_dict)
	979	elif (ofile[-4:] == ".pts"):
	980	_write_pts_file(ofile, point_dict)
	981	else:
	982	msg = 'Unknown file type %s ' %ofile
	983	raise IOError, msg
	984
	985	def import_points_file(ofile, delimiter = None, verbose = False):
	986	"""
	987	load an .xya or .pts file
	988
	989	Note: will throw an IOError if it can't load the file.
	990	Catch these!
	991	"""
	992
[2459]	993
	994	#FIXME (Ole): This function should really return a Geospatial_data object. #FIXME (DSG): Do you know it does, in the points dic?
[2941]	995	#No. Where this function is used, geospatial objects should be used.
[2459]	996
[2253]	997	if ofile[-4:]== ".xya":
	998	try:
	999	if delimiter == None:
	1000	try:
	1001	fd = open(ofile)
	1002	points_dict = _read_xya_file(fd, ',')
[2941]	1003	except TitleError: # this is catching the error thrown by geo_ref
[2253]	1004	fd.close()
	1005	fd = open(ofile)
	1006	points_dict = _read_xya_file(fd, ' ')
	1007	else:
	1008	fd = open(ofile)
	1009	points_dict = _read_xya_file(fd, delimiter)
	1010	fd.close()
	1011	except (IndexError,ValueError,SyntaxError):
[2941]	1012	fd.close()
[2253]	1013	msg = 'Could not open file %s ' %ofile
	1014	raise IOError, msg
[2941]	1015	except TitleError:
	1016	print "reclassifying title error"
	1017	fd.close()
	1018	msg = 'Could not open file %s ' %ofile
	1019	raise IOError, msg
[2253]	1020	except IOError:
	1021	# Catch this to add an error message
	1022	msg = 'Could not open file %s ' %ofile
	1023	raise IOError, msg
	1024
	1025	elif ofile[-4:]== ".pts":
	1026	try:
	1027	points_dict = _read_pts_file(ofile, verbose)
	1028	except IOError, e:
	1029	msg = 'Could not open file %s ' %ofile
	1030	raise IOError, msg
	1031	else:
	1032	msg = 'Extension %s is unknown' %ofile[-4:]
	1033	raise IOError, msg
	1034	return points_dict
	1035
	1036	def extent_point_atts(point_atts):
	1037	"""
	1038	Returns 4 points representing the extent
	1039	This losses attribute info.
	1040	"""
	1041	point_atts['pointlist'] = extent(point_atts['pointlist'])
	1042	point_atts['attributelist'] = {}
	1043	return point_atts
	1044
	1045	def extent(points):
	1046	points = array(points).astype(Float)
	1047	max_x = min_x = points[0][0]
	1048	max_y = min_y = points[0][1]
	1049	for point in points[1:]:
	1050	x = point[0]
	1051	if x > max_x: max_x = x
	1052	if x < min_x: min_x = x
	1053	y = point[1]
	1054	if y > max_y: max_y = y
	1055	if y < min_y: min_y = y
	1056	extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
	1057	#print "extent",extent
	1058	return extent
	1059
	1060	def reduce_pts(infile, outfile, max_points, verbose = False):
	1061	"""
	1062	reduces a points file by removong every second point until the # of points
	1063	is less than max_points.
	1064	"""
[2888]	1065	# check out pts2rectangular in least squares, and the use of reduction.
	1066	# Maybe it does the same sort of thing?
[2253]	1067	point_atts = _read_pts_file(infile)
	1068	while point_atts['pointlist'].shape[0] > max_points:
	1069	if verbose: print "point_atts['pointlist'].shape[0]"
	1070	point_atts = half_pts(point_atts)
	1071	export_points_file(outfile, point_atts)
	1072
	1073	def produce_half_point_files(infile, max_points, delimiter, verbose = False):
	1074	point_atts = _read_pts_file(infile)
	1075	root, ext = splitext(infile)
	1076	outfiles = []
	1077	if verbose: print "# of points", point_atts['pointlist'].shape[0]
	1078	while point_atts['pointlist'].shape[0] > max_points:
	1079	point_atts = half_pts(point_atts)
	1080	if verbose: print "# of points", point_atts['pointlist'].shape[0]
	1081	outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
	1082	outfiles.append(outfile)
	1083	export_points_file(outfile,
	1084	point_atts)
	1085	return outfiles
	1086
	1087	def point_atts2array(point_atts):
	1088	point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
	1089
	1090	for key in point_atts['attributelist'].keys():
	1091	point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
	1092	return point_atts
	1093
	1094	def half_pts(point_atts):
	1095	point_atts2array(point_atts)
	1096	point_atts['pointlist'] = point_atts['pointlist'][::2]
	1097
	1098	for key in point_atts['attributelist'].keys():
	1099	point_atts['attributelist'][key]= point_atts['attributelist'][key] [::2]
	1100	return point_atts
	1101
	1102	def concatinate_attributelist(dic):
	1103	"""
	1104	giving a dic[attribute title] = attribute
	1105	return list of attribute titles, array of attributes
	1106	"""
	1107	point_attributes = array([]).astype(Float)
	1108	keys = dic.keys()
	1109	key = keys.pop(0)
	1110	point_attributes = reshape(dic[key],(dic[key].shape[0],1))
	1111	for key in keys:
	1112	#point_attributes = concatenate([point_attributes, dic[key]], axis=1)
	1113	reshaped = reshape(dic[key],(dic[key].shape[0],1))
	1114	point_attributes = concatenate([point_attributes, reshaped], axis=1)
	1115	return dic.keys(), point_attributes
	1116
	1117	def _read_pts_file(file_name, verbose = False):
	1118	"""Read .pts NetCDF file
	1119
	1120	Return a dic of array of points, and dic of array of attribute
	1121
	1122	eg
	1123	dic['points'] = [[1.0,2.0],[3.0,5.0]]
	1124	dic['attributelist']['elevation'] = [[7.0,5.0]
	1125	"""
	1126	#FIXME: (DSG)This format has issues.
	1127	# There can't be an attribute called points
	1128	# consider format change
	1129
	1130
	1131	if verbose: print 'Reading ', file_name
	1132
	1133
	1134	# see if the file is there. Throw a QUIET IO error if it isn't
	1135	fd = open(file_name,'r')
	1136	fd.close()
	1137
	1138	#throws prints to screen if file not present
	1139	fid = NetCDFFile(file_name, 'r')
	1140
	1141	point_atts = {}
	1142	# Get the variables
	1143	point_atts['pointlist'] = array(fid.variables['points'])
	1144	keys = fid.variables.keys()
	1145	if verbose: print 'Got %d variables: %s' %(len(keys), keys)
	1146	try:
	1147	keys.remove('points')
	1148	except IOError, e:
	1149	fid.close()
	1150	msg = 'Expected keyword "points" but could not find it'
	1151	raise IOError, msg
	1152
	1153	attributes = {}
	1154	for key in keys:
	1155	if verbose: print "reading attribute '%s'" %key
	1156
	1157	attributes[key] = array(fid.variables[key])
	1158
	1159	point_atts['attributelist'] = attributes
	1160
	1161	try:
	1162	geo_reference = Geo_reference(NetCDFObject=fid)
	1163	point_atts['geo_reference'] = geo_reference
	1164	except AttributeError, e:
	1165	#geo_ref not compulsory
	1166	point_atts['geo_reference'] = None
	1167
	1168
	1169	fid.close()
	1170
	1171	#print "point_atts",point_atts
	1172	return point_atts
	1173
	1174	def _write_pts_file(file_name, point_atts):
	1175	"""Write .pts NetCDF file
	1176
	1177	WARNING: This function mangles the point_atts data structure
	1178	"""
	1179	#FIXME: (DSG)This format has issues.
	1180	# There can't be an attribute called points
	1181	# consider format change
[2262]	1182
	1183
	1184	legal_keys = ['pointlist', 'attributelist', 'geo_reference']
	1185	for key in point_atts.keys():
	1186	msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
	1187	assert key in legal_keys, msg
[2253]	1188
	1189	point_atts2array(point_atts)
	1190	# NetCDF file definition
	1191	outfile = NetCDFFile(file_name, 'w')
	1192
	1193	#Create new file
	1194	outfile.institution = 'Geoscience Australia'
	1195	outfile.description = 'NetCDF format for compact and portable storage ' +\
	1196	'of spatial point data'
	1197
	1198	# dimension definitions
	1199	shape = point_atts['pointlist'].shape[0]
	1200	outfile.createDimension('number_of_points', shape)
	1201	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
	1202
	1203	# variable definition
	1204	outfile.createVariable('points', Float, ('number_of_points',
	1205	'number_of_dimensions'))
	1206
	1207	#create variables
	1208	outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
	1209	for key in point_atts['attributelist'].keys():
	1210	outfile.createVariable(key, Float, ('number_of_points',))
	1211	outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
	1212
	1213	if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
	1214	point_atts['geo_reference'].write_NetCDF(outfile)
	1215
	1216	outfile.close()
	1217
	1218	def _read_xya_file(fd, delimiter):
	1219	#lines = fd.readlines()
	1220	points = []
	1221	pointattributes = []
	1222	#if len(lines) <= 1:
	1223	# raise SyntaxError
	1224	title = fd.readline()
	1225	#title = lines.pop(0) # the first (title) line
	1226	att_names = clean_line(title,delimiter)
	1227
	1228	att_dict = {}
	1229	line = fd.readline()
	1230	numbers = clean_line(line,delimiter)
	1231	#for line in lines:
	1232	while len(numbers) > 1:
	1233	#print "line >%s" %line
	1234	#numbers = clean_line(line,delimiter)
	1235	#print "numbers >%s<" %numbers
	1236	#if len(numbers) < 2 and numbers != []:
	1237
	1238	# A line without two numbers
	1239	# or a bad delimiter
	1240	#FIXME dsg-dsg change error that is raised.
	1241	#raise SyntaxError
	1242	if numbers != []:
	1243	try:
	1244	x = float(numbers[0])
	1245	y = float(numbers[1])
	1246	points.append([x,y])
	1247	numbers.pop(0)
	1248	numbers.pop(0)
	1249	#attributes = []
	1250	#print "att_names",att_names
	1251	#print "numbers",numbers
	1252	if len(att_names) != len(numbers):
	1253	fd.close()
	1254	# It might not be a problem with the title
	1255	#raise TitleAmountError
	1256	raise IOError
	1257	for i,num in enumerate(numbers):
	1258	num.strip()
	1259	if num != '\n' and num != '':
	1260	#attributes.append(float(num))
	1261	att_dict.setdefault(att_names[i],[]).append(float(num))
	1262
	1263	except ValueError:
	1264	raise SyntaxError
	1265	line = fd.readline()
	1266	numbers = clean_line(line,delimiter)
	1267
	1268	if line == '':
	1269	# end of file
	1270	geo_reference = None
	1271	else:
	1272	geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
	1273
	1274	xya_dict = {}
	1275	xya_dict['pointlist'] = array(points).astype(Float)
	1276
	1277	for key in att_dict.keys():
	1278	att_dict[key] = array(att_dict[key]).astype(Float)
	1279	xya_dict['attributelist'] = att_dict
	1280	xya_dict['geo_reference'] = geo_reference
	1281	#print "xya_dict",xya_dict
	1282	return xya_dict
	1283
	1284	#FIXME(dsg), turn this dict plus methods into a class?
	1285	def take_points(dict,indices_to_keep):
	1286	dict = point_atts2array(dict)
	1287	#FIXME maybe the points data structure should become a class?
	1288	dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
	1289
	1290	for key in dict['attributelist'].keys():
	1291	dict['attributelist'][key]= take(dict['attributelist'][key],
	1292	indices_to_keep)
	1293	return dict
	1294
	1295	def add_point_dictionaries (dict1, dict2):
	1296	"""
	1297	"""
	1298	dict1 = point_atts2array(dict1)
	1299	dict2 = point_atts2array(dict2)
	1300
	1301	combined = {}
	1302	combined['pointlist'] = concatenate((dict2['pointlist'],
	1303	dict1['pointlist']),axis=0)
	1304	atts = {}
	1305	for key in dict2['attributelist'].keys():
	1306	atts[key]= concatenate((dict2['attributelist'][key],
	1307	dict1['attributelist'][key]), axis=0)
	1308	combined['attributelist']=atts
	1309	combined['geo_reference'] = dict1['geo_reference']
	1310	return combined
	1311
	1312	def _write_xya_file( file_name, xya_dict, delimiter = ','):
	1313	"""
	1314	export a file, ofile, with the xya format
	1315
	1316	"""
	1317	points = xya_dict['pointlist']
	1318	pointattributes = xya_dict['attributelist']
	1319
	1320	fd = open(file_name,'w')
	1321
	1322	titlelist = ""
	1323	for title in pointattributes.keys():
	1324	titlelist = titlelist + title + delimiter
	1325	titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
	1326	fd.write(titlelist+"\n")
	1327	#<vertex #> <x> <y> [attributes]
	1328	for i,vert in enumerate( points):
	1329
	1330	attlist = ","
	1331	for att in pointattributes.keys():
	1332	attlist = attlist + str(pointattributes[att][i])+ delimiter
	1333	attlist = attlist[0:-len(delimiter)] # remove the last delimiter
	1334	attlist.strip()
	1335	fd.write( str(vert[0]) + delimiter
	1336	+ str(vert[1])
	1337	+ attlist + "\n")
	1338
	1339	# geo_reference info
	1340	if xya_dict.has_key('geo_reference') and \
	1341	not xya_dict['geo_reference'] is None:
	1342	xya_dict['geo_reference'].write_ASCII(fd)
	1343	fd.close()
	1344
	1345	if __name__ == "__main__":
	1346	pass

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

Context Navigation

source: inundation/load_mesh/loadASCII.py @ 3381

Download in other formats: