Context Navigation

source: inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py @ 1376

Last change on this file since 1376 was 1376, checked in by duncan, 20 years ago
refactoring / changing method names
File size: 47.9 KB

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

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

Download in other formats: