Context Navigation

source: anuga_core/source_numpy_conversion/anuga/load_mesh/loadASCII.py @ 5946

Last change on this file since 5946 was 5946, checked in by rwilson, 15 years ago
More NumPy? changes.
File size: 41.2 KB

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

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

Download in other formats: