-
Notifications
You must be signed in to change notification settings - Fork 13
/
gis2tin.py
executable file
·192 lines (178 loc) · 7.17 KB
/
gis2tin.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
# #
# gis2tin.py #
# #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
#
# Date: February 2, 2016
#
# Purpose: Same inputs as gis2triangle.py; calls gis2triangle.py, then
# uses python's subprocess module to call triangle to create the tin, then calls
# triangle2adcirc.py to produce an adcirc tin file.
#
# Modified: Feb 20, 2016
# Made it work for python 2 and 3
#
# Revised: Apr 29, 2017
# Changed how different system architectures are called; made it work
# for the raspberry pi system.
#
# Revised: May 6, 2017
# Placed a call to processor type inside the posix if statement.
#
# Uses: Python 2 or 3, Numpy
#
# Example:
#
# python gis2tin.py -n nodes.csv -b boundary.csv -l lines.csv -h none -o out.grd
#
# 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. For triangle, it has
# no meaning. 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 listing of the holes in the mesh. The holes file is
# generated by placing a single node marker inside a
# closed line constraint. The holes file must include a
# x,y within the hole boundary. If no holes (islands)
# in the mesh, enter 'none' without the quotes. Note
# that for triangle, the format of the holes file is
# different than for gmsh!!!
#
# --> -o is the output adcirc file that is the TIN.
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Global Imports
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import os,sys # system parameters
import numpy as np # numpy
from collections import OrderedDict # for removal of duplicate nodes
import struct # to determine sys architecture
import subprocess # to execute binaries
from ppmodules.readMesh import * # to get all readMesh functions
import matplotlib.tri as mtri # matplotlib triangulations
#
curdir = os.getcwd()
#
try:
# this only works when the paths are sourced!
pputils_path = os.environ['PPUTILS']
except:
pputils_path = curdir
# this is to maintain legacy support
if (sys.version_info > (3, 0)):
version = 3
pystr = 'python3'
elif (sys.version_info > (2, 7)):
version = 2
pystr = 'python'
#
# 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]
else:
print('Wrong number of Arguments, stopping now...')
print('Usage:')
print('python gis2tin.py -n nodes.csv -b boundary.csv -l lines.csv -h holes.csv -o out.grd')
sys.exit()
# to determine if the system is 32 or 64 bit
archtype = struct.calcsize("P") * 8
# Tell the user what is going on
print('Constructing Triangle poly file ...')
# call gis2triangle.py
# have to deal with cases when paths are sources, and when they are not
try:
# this only works when the paths are sourced!
pputils_path = os.environ['PPUTILS']
subprocess.call(['gis2triangle_kd.py', '-n', nodes_file,
'-b', boundary_file, '-l', lines_file, '-h', holes_file,
'-o', 'tin.poly'])
except:
subprocess.call([pystr, 'gis2triangle_kd.py', '-n', nodes_file,
'-b', boundary_file, '-l', lines_file, '-h', holes_file,
'-o', 'tin.poly'])
print('Generating TIN using Triangle ...')
if (os.name == 'posix'):
# determines processor type
proctype = os.uname()[4][:]
# for linux32 its i686
# for linux64 its x86_64
# for raspberry pi 32 its armv7l
# this assumes chmod +x has already been applied to the binaries
if (proctype == 'i686'):
subprocess.call( [pputils_path + '/triangle/bin/triangle_32', 'tin.poly' ] )
elif (proctype == 'x86_64'):
subprocess.call( [pputils_path + '/triangle/bin/triangle_64', 'tin.poly' ] )
elif (proctype == 'armv7l'):
subprocess.call( [pputils_path + '/triangle/bin/triangle_pi32', 'tin.poly' ] )
elif (os.name == 'nt'):
subprocess.call( ['.\\triangle\\bin\\triangle_32.exe', 'tin.poly' ] )
else:
print('OS not supported!')
print('Exiting!')
sys.exit()
# call triangle2adcirc.py
try:
# this only works when the paths are sourced!
subprocess.call(['gis2triangle_kd.py', '-n', nodes_file,
'-b', boundary_file, '-l', lines_file, '-h', holes_file,
'-o', 'tin.poly'])
except:
subprocess.call([pystr, 'gis2triangle_kd.py', '-n', nodes_file,
'-b', boundary_file, '-l', lines_file, '-h', holes_file,
'-o', 'tin.poly'])
print('Converting TIN from Triangle to ADCIRC format ...')
try:
subprocess.call(['triangle2adcirc.py', '-n', 'tin.1.node',
'-e', 'tin.1.ele', '-o', output_file])
except:
subprocess.call([pystr, 'triangle2adcirc.py', '-n', 'tin.1.node',
'-e', 'tin.1.ele', '-o', output_file])
# construct the output wkt file
wkt_file = output_file.rsplit('.',1)[0] + '_WKT.csv'
# now convert the *.grd file to a *.wkt file by calling adcirc2wkt.py
print('Converting TIN from ADCIRC to WKT format ...')
try:
subprocess.call(['adcirc2wkt.py', '-i', output_file, '-o', wkt_file])
except:
subprocess.call([pystr, 'adcirc2wkt.py', '-i', output_file, '-o', wkt_file])
# to remove the temporary files
os.remove('tin.poly')
os.remove('tin.1.poly')
os.remove('tin.1.node')
os.remove('tin.1.ele')
print('All done!')