-
Notifications
You must be signed in to change notification settings - Fork 12
/
Copy pathextract_line_t.py
executable file
·160 lines (136 loc) · 4.84 KB
/
extract_line_t.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
#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
# #
# extract_line_t.py #
# #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
#
# Date: Feb 22, 2017
#
# Purpose: Script designed to open 2D telemac binary file and write the
# values of all variables (for a particular time step) for the nodes
# along the given line.
#
# Revised: Mar 9, 2017
# Deleted the trailing comma in the text output.
#
# Uses: Python 2 or 3, Numpy
#
# Usage:
# python extract_line_t.py -i res.slf -t 8 -l line.csv -o line_t.csv
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Global Imports
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import sys
import numpy as np
import matplotlib.tri as mtri
from ppmodules.selafin_io_pp import *
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# MAIN
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# I/O
if len(sys.argv) == 9:
input_file = sys.argv[2] # input *.slf file
t = int(sys.argv[4]) # index of time record
lines_file = sys.argv[6] # input lines file
output_file = sys.argv[8] # output lines file
else:
print('Wrong number of arguments ... stopping now ...')
print('Usage:')
print('python extract_line_t.py -i res.slf -t 8 -l line.csv -o line_t.csv')
sys.exit()
# use Numpy to read the lines file
# all data is read as floats
lines_data = np.loadtxt(lines_file, delimiter=',',skiprows=0,unpack=True)
shapeid = lines_data[0,:]
shapeid = shapeid.astype(np.int32)
lnx = lines_data[1,:]
lny = lines_data[2,:]
# Read the header of the selafin result file and get geometry and
# variable names and units
print('Reading result file headers ...')
# use selafin_io_pp class ppSELAFIN
slf = ppSELAFIN(input_file)
slf.readHeader()
slf.readTimes()
# get the mesh properties from the resultfile
NELEM, NPOIN, NDP, IKLE, IPOBO, x, y = slf.getMesh()
# the IKLE array starts at element 1, but matplotlib needs it to start at zero
# note that doing this invalidates the IPOBO array; but we don't need IPOBO here
IKLE[:,:] = IKLE[:,:] - 1
# use get methods from ppSELAFIN class
times = slf.getTimes()
vnames = slf.getVarNames()
vunits = slf.getVarUnits()
NPLAN = slf.getNPLAN()
# to warn the user that this applies for 2d SELAFIN files only
if (NPLAN > 1):
print('It does not make sense to run this script for 3d SELAFIN files.')
print('Use POSTEL-3D or Paraview instead!')
sys.exit(0)
# number of variables
numvars = len(vnames)
# create the new output variables for the output file
ln_interp = np.zeros((numvars, len(lnx)))
sta = np.zeros(len(lnx))
tempid = np.zeros(len(lnx))
dist = np.zeros(len(lnx))
# to create the sta array
sta[0] = 0.0
tempid = shapeid
dist[0] = 0.0
for i in range(1,len(lnx)):
if (tempid[i] - shapeid[i-1] < 0.001):
xdist = lnx[i] - lnx[i-1]
ydist = lny[i] - lny[i-1]
dist[i] = np.sqrt(np.power(xdist,2.0) + np.power(ydist,2.0))
sta[i] = sta[i-1] + dist[i]
# create a result triangulation object
triang = mtri.Triangulation(x, y, IKLE)
# read the variables for the specified time step t
slf.readVariables(t)
results = slf.getVarValues()
print('Interpolating ...')
# now to interpolate the result file to the nodes of the lines file, for
# each variable in the file (for the specified time step only).
for j in range(numvars):
interpolator = mtri.LinearTriInterpolator(triang, results[j,:])
ln_interp[j,:] = interpolator(lnx,lny)
print('Transposing ...')
# tranpose the interpolated result, so that they will print nicely
ln_interp_tr = np.transpose(ln_interp)
where_are_NaNs = np.isnan(ln_interp_tr)
if (np.sum(where_are_NaNs) > 0):
print('#####################################################')
print('')
print('WARNING: Some line nodes are outside of the mesh boundary!!!')
print('')
print('#####################################################')
print('Writing output to file ...')
# to write the output file
fout = open(output_file, 'w')
# to write variable names header (this is not in true style of pputils
# lines data, but it will be easier on the eyes ...)
fout.write(str('id, x, y, sta, '))
for j in range(numvars):
if (j < len(vnames)-1):
fout.write(str(vnames[j]) + ', ')
else:
fout.write(str(vnames[j]))
fout.write('\n')
# to write the data
for i in range(len(lnx)):
fout.write(str(shapeid[i]) + ', ' + str(lnx[i]) + ', ' + str(lny[i]) + ', ')
fout.write(str(sta[i]) + ', ')
for j in range(numvars):
if (j < len(vnames)-1):
fout.write(str(ln_interp_tr[i,j]) + ', ')
else:
fout.write(str(ln_interp_tr[i,j]))
fout.write('\n')