gis2gmsh_kd.py

#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#                                                                       #
#                                 gis2gmsh_kd.py                        # 
#                                                                       #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
# 
# Date: August 16, 2015
#
# Modified: Feb 20, 2016
# Made it work for python 2 and 3
#
# Purpose: Script takes in a text file of the geometry generated in qgis
# (or any other gis or cad package) and produces geometry files used by
# gmsh mesh generator program.
#
# Same as gis2gmsh.py, except that it uses kd tree searching algorithm for
# node searching. This is far superior than my node searching algorithm!
#
# Revised: Nov 21, 2016
# Changed KDTree to cKDTree to improve performance.
#
# Revised: Nov 23, 2016
# For some test cases cKDTree crashed, while KDTree went to completion.
# Therefore, revert back to using KDTree.
#
# Uses: Python2.7.9, Numpy v1.8.2
#
# Example:
#
# python gis2gmsh.py -n nodes.csv -b boundary.csv -l lines.csv -h holes.csv -o out.geo
#
# where:
#       --> -n is the file listing of all nodes (incl. embedded nodes
#                        if any). The nodes file consist of x,y,z or x,y,z,size;
#                        The size parameter is an optional input, and is used 
#                        by gmsh as an extra parameter that forces element 
#                        size around a particular node. The nodes file must
#                        be comma separated, and have no header lines. 
#
#       --> -b is the node listing of the outer boundary for the mesh.
#                        The boundary file is generated by snapping lines
#                        to the nodes from the nodes.csv file. The boundary file 
#                        consists of shapeid,x,y of all the lines in the file.
#                        Boundary has to be a closed shape, where first and last 
#                        nodes are identical. Shapeid is a integer, where the
#                        boundary is defined with a distict id (i.e., shapeid 
#                        of 0). 
#
#       --> -l is the node listing of the constraint lines for the mesh.
#                        The lines file can include open or closed polylines.
#                        The file listing has shapeid,x,y, where x,y have to 
#                        reasonable match that of the nodes.csv file. Each distinct
#                        line has to have an individual (integer) shapeid. If no 
#                        constraint lines in the mesh, enter 'none' without the quotes.
#
#       --> -h is the node listing of the holes in the mesh.
#                        The holes file must include closed polylines. The
#                        file listing has shapeid,x,y, where x,y have to reasonably
#                        match that of the nodes.csv file. Each distinct hole has to 
#                        have an individual (integer) shapeid. If no holes 
#                        (islands) in the mesh, enter 'none' without the quotes.
#
#      --> -o is the output gmsh geometry format. To generate the mesh, launch 
#                        gmsh and go to: Modules --> Mesh --> 2D
#
#      --> -d is an optional flag to ignore removal of duplicate nodes in the 
#                        nodes file. By default, duplicate nodes are removed 
#                        from the nodes.csv file, and the user need not set
#                        this flag. 
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Global Imports
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import os,sys                              # system parameters
import numpy             as np             # numpy
from collections import OrderedDict        # for removal of duplicate nodes
from scipy import spatial                  # kd tree for searching coords
curdir = os.getcwd()
#
# I/O
if len(sys.argv) == 11 :
	dummy1 =  sys.argv[1]
	nodes_file = sys.argv[2]
	dummy2 =  sys.argv[3]
	boundary_file = sys.argv[4]
	dummy3 =  sys.argv[5]
	lines_file = sys.argv[6]
	dummy4 =  sys.argv[7]
	holes_file = sys.argv[8]
	dummy5 =  sys.argv[9]
	output_file = sys.argv[10]
	dummy6 = ' '
	duplicates_flag = 1 # removal of duplicate nodes is on by default
elif (len(sys.argv) == 13):
	dummy1 =  sys.argv[1]
	nodes_file = sys.argv[2]
	dummy2 =  sys.argv[3]
	boundary_file = sys.argv[4]
	dummy3 =  sys.argv[5]
	lines_file = sys.argv[6]
	dummy4 =  sys.argv[7]
	holes_file = sys.argv[8]
	dummy5 =  sys.argv[9]
	output_file = sys.argv[10]
	dummy6 = sys.argv[11]
	duplicates_flag = sys.argv[12]
else:
	print('Wrong number of Arguments, stopping now...')
	print('Usage:')
	print('python gis2gmsh.py -n nodes.csv -b boundary.csv -l lines.csv -h holes.csv -o out.geo')
	#print 'or, if wanting to turn off duplicate removal algorithm'
	#print 'python gis2gmsh.py -n nodes.csv -b boundary.csv -l lines.csv -h holes.csv -o out.geo -d 0'
	sys.exit()

# find out if the nodes file is x,y,z or x,y,x,size
with open(nodes_file, 'r') as f:
    line = next(f) # read 1 line
    n_attr = len(line.split(','))

# to create the output file
fout = open(output_file,"w")

# use numpy to read the file
# each column in the file is a row in data read by no.loadtxt method
nodes_data = np.loadtxt(nodes_file, delimiter=',',skiprows=0,unpack=True)
boundary_data = np.loadtxt(boundary_file, delimiter=',',skiprows=0,unpack=True)

if (lines_file != 'none'):
	lines_data = np.loadtxt(lines_file, delimiter=',',skiprows=0,unpack=True)

if (holes_file != 'none'):
	holes_data = np.loadtxt(holes_file, delimiter=',',skiprows=0,unpack=True)

# master nodes in the file (from the nodes file)
x = nodes_data[0,:]
y = nodes_data[1,:]
z = nodes_data[2,:]
if (n_attr == 4):
	size = nodes_data[3,:]
else:
	size = np.zeros(len(x))

# n is the number of nodes
n = len(x)

# creates node numbers from the nodes file
node = np.zeros(n,dtype=np.int32)

# to check for duplicate nodes
# crop all the points to three decimals only
x = np.around(x,decimals=3)
y = np.around(y,decimals=3)
z = np.around(z,decimals=3)
size = np.around(size,decimals=3)

# this piece of code uses OrderedDict to remove duplicate nodes
# source "http://stackoverflow.com/questions/12698987"
# ###################################################################
tmp = OrderedDict()
for point in zip(x, y, z, size):
  tmp.setdefault(point[:2], point)

# in python 3 tmp.values() is a view object that needs to be 
# converted to a list
mypoints = list(tmp.values()) 
# ###################################################################
n_rev = len(mypoints)

# replace x,y,z,size and n with their unique equivalents
if (duplicates_flag == 1):
	for i in range(n_rev):
		x[i] = mypoints[i][0]
		y[i] = mypoints[i][1]
		z[i] = mypoints[i][2]
		size[i] = mypoints[i][3]
	n = n_rev

# when I made the change to python 3, had to use np.column_stack
# http://stackoverflow.com/questions/28551279/error-running-scipy-kdtree-example

# to create the tuples of the master points
points = np.column_stack((x,y))
tree = spatial.KDTree(points)

# if node is part of boundary or lines, then it is not embedded
is_node_emb = np.zeros(n,dtype=np.int32)
for i in range(0,n):
	node[i] = i+1
	is_node_emb[i] = 1

# boundary data
shapeid_bnd = boundary_data[0,:]
x_bnd = boundary_data[1,:]
y_bnd = boundary_data[2,:]

# round boundary nodes to three decimals
x_bnd = np.around(x_bnd,decimals=3)
y_bnd = np.around(y_bnd,decimals=3)

# number of nodes in the boundary file
n_bnd = len(x_bnd)

# count lines from boundary lines
count_bnd = 0

# lines data
if (lines_file != 'none'):
	shapeid_lns = lines_data[0,:]
	x_lns = lines_data[1,:]
	y_lns = lines_data[2,:]
	
	# round lines nodes to three decimals
	x_lns = np.around(x_lns,decimals=3)
	y_lns = np.around(y_lns,decimals=3)
	
	# number of nodes in the lines file
	n_lns = len(x_lns)

# holes data
if (holes_file != 'none'):
	shapeid_hls = holes_data[0,:]
	x_hls = holes_data[1,:]
	y_hls = holes_data[2,:]
	
	# round lines nodes to three decimals
	x_hls = np.around(x_hls,decimals=3)
	y_hls = np.around(y_hls,decimals=3)
	
	# number of nodes in the holes file
	n_hls = len(x_hls)

count_lns = 0

# writes the nodes in gmsh format
for i in range(0,n):
	fout.write("Point(" + str(i+1) + ") = {" + str("{:.3f}".format(x[i])) +
		str(", ") + str("{:.3f}".format(y[i])) + str(", ") + 
		str("{:.3f}".format(z[i])) + str(", ") + str("{:.3f}".format(size[i])) + 
		str("};") + "\n")
	
# BOUNDARY LINES
# index of the minimum, for each boundary node
minidx = np.zeros(n_bnd,dtype=np.int32) -1

pt_bnd = list()

for i in range(0,n_bnd):
	pt_bnd.append(x_bnd[i])
	pt_bnd.append(y_bnd[i])
	
	# find the index of the boundary point from the nodes file
	minidx_temp = tree.query_ball_point(pt_bnd, 0.01)
	
	if (len(minidx_temp) > 0):
		minidx[i] = minidx_temp[0]
	else:
		print('Boundary node ' + str(x_bnd[i]) + ' ' + str(y_bnd[i]) + ' not found')
	
	# fill in the is_node_emb array 
	is_node_emb[minidx[i]] = 0
	
	# remove the node to search for
	pt_bnd.remove(x_bnd[i])
	pt_bnd.remove(y_bnd[i])
	
# write the boundary in gmsh format
for i in range(0,n_bnd-1):
	if (i == 0) :
		fout.write("Line(" + str(i+1) + str(") = {") + str(node[minidx[0]])
			+ str(", ") + str(node[minidx[1]]) + str("};") + "\n")
		count_bnd =count_bnd +1
	else:
		fout.write("Line(" + str(i+1) + str(") = {") + str(node[minidx[i]])
			+ str(", ") + str(node[minidx[i+1]]) + str("};") + "\n")
		count_bnd =count_bnd +1

# the lines numbering continues from the boundary numbering
count_lns = count_bnd + 1

# CONSTRAINT LINES
if (lines_file != 'none'):
	####################################
	# index for the minimum, for each lines node
	minidx_lns = np.zeros(n_lns,dtype=np.int32) -1
	
	pt_lns = list()
	
	for i in range(0,n_lns):
		pt_lns.append(x_lns[i])
		pt_lns.append(y_lns[i])
		
		# find the index of the lines from the nodes file
		minidx_lns_temp = tree.query_ball_point(pt_lns, 0.01)
		
		if (len(minidx_lns_temp) > 0):
			minidx_lns[i] = minidx_lns_temp[0]
		else:
			print('Lines node ' + str(x_lns[i]) + ' ' + str(y_lns[i]) + ' not found')
			print('Exiting ...')
			sys.exit()

		#fout.write(str(i) + " " + str(minidx_lns[i]) + "\n")
		
		# fill in the is_node_emb array 
		is_node_emb[minidx_lns[i]] = 0
		
		# to remove the node to search for
		pt_lns.remove(x_lns[i])
		pt_lns.remove(y_lns[i])	
		
	cur_lns_shapeid = shapeid_lns[0]
	prev_lns_shapeid = shapeid_lns[0]	
	
	# write the constraint lines
	for i in range(0,n_lns):
		if (i>0):
			cur_lns_shapeid = shapeid_lns[i]
			prev_lns_shapeid = shapeid_lns[i-1]
			if (cur_lns_shapeid - prev_lns_shapeid < 0.001):
				#fout.write(str(cur_lns_shapeid) + " " + str(prev_lns_shapeid) + " ")
				fout.write("Line(" + str(count_lns) + str(") = {") + 
					str(node[minidx_lns[i-1]]) + str(", ") + str(node[minidx_lns[i]]) + str("};") + "\n")
				count_lns = count_lns + 1
	
####################################
# holes
count_hls = count_lns +1

hole_nodes = list()
if (holes_file != 'none'):
	# index for the minimum, for each lines node
	minidx_hls = np.zeros(n_hls,dtype=np.int32) - 1
	
	pt_hls = list()
	
	for i in range(0,n_hls):
		pt_hls.append(x_hls[i])
		pt_hls.append(y_hls[i])
		
		# find the index of each holes point
		minidx_hls_temp = tree.query_ball_point(pt_hls, 0.01)
		
		if (len(minidx_hls_temp) > 0):
			minidx_hls[i] = minidx_hls_temp[0]
		else:
			print('Holes node ' + str(x_hls[i]) + ' ' + str(y_hls[i]) + ' not found')
			print('Exiting ...')
			sys.exit()		
		
		#fout.write(str(i) + " " + str(minidx_hls[i]) + "\n")
		
		# fill in the is_node_emb array 
		is_node_emb[minidx_hls[i]] = 0
		
		# to remove the node to search for
		pt_hls.remove(x_hls[i])
		pt_hls.remove(y_hls[i])	
		
	cur_hls_shapeid = shapeid_hls[0]
	prev_hls_shapeid = shapeid_hls[0]	
	
	# write the constraint lines
	for i in range(0,n_hls):
		if (i>0):
			cur_hls_shapeid = shapeid_hls[i]
			prev_hls_shapeid = shapeid_hls[i-1]
			if (cur_hls_shapeid - prev_hls_shapeid < 0.001):
				#fout.write(str(cur_hls_shapeid) + " " + str(prev_hls_shapeid) + " ")
				hole_nodes.append(count_hls)
				fout.write("Line(" + str(count_hls) + str(") = {") + 
					str(node[minidx_hls[i-1]]) + str(", ") + str(node[minidx_hls[i]]) + str("};") + "\n")
				count_hls = count_hls + 1
	####################################

#fout.write("HOLES" + str(hole_nodes) + '\n')

n_holes = len(hole_nodes)

# writes the line loop and the plane surface
fout.write("Line Loop(1) = {1:"+ str(count_bnd))
if (holes_file != 'none'):
	fout.write(', ')
	for i in range(n_holes-1):
		fout.write(str(hole_nodes[i]*-1) + ", ")
	fout.write(str(-1*hole_nodes[n_holes-1]))	
fout.write(str("};") + "\n")

fout.write("Physical Line(1) = {1:"+ str(count_bnd))
if (holes_file != 'none'):
	fout.write(', ')
	for i in range(n_holes-1):
		fout.write(str(hole_nodes[i]*-1) + ", ")
	fout.write(str(-1*hole_nodes[n_holes-1]))	
fout.write(str("};") + "\n")


fout.write("Plane Surface(1) = {1};" + "\n")
fout.write("Physical Surface(1) = {1};" + "\n")

if (lines_file != 'none'):
	# write the embedded lines
	# re-set the count_lns back to what it was before
	count_lns = count_bnd + 1
	for i in range(0,n_lns):
		if (i>0):
			cur_lns_shapeid = shapeid_lns[i]
			prev_lns_shapeid = shapeid_lns[i-1]
			if (cur_lns_shapeid - prev_lns_shapeid < 0.001):
				#fout.write(str(cur_lns_shapeid) + " " + str(prev_lns_shapeid) + " ")
				fout.write(str("Line {") + str(count_lns) + "} In Surface {1};" + "\n")
				count_lns = count_lns + 1
	
# if there are embedded nodes, write them to the file
# embedded nodes should be used in tin applications, not in mesh generation
for i in range(0,n):
	if (is_node_emb[i] == 1):
		fout.write(str("Point {") + str(node[i]) + "} In Surface {1};" + "\n")

# gmsh option to make sure the elements size extend from boundary
# write zero when doing a TIN; write one when doing a mesh!
#fout.write(str("Mesh.CharacteristicLengthExtendFromBoundary = 0;") + "\n")