forked from jasonfleming/pputils
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathadcirc2flt.py
executable file
·136 lines (116 loc) · 4.51 KB
/
adcirc2flt.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
#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
# #
# adcirc2flt.py #
# #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
#
# Date: Feb 16, 2016
#
# Purpose: Script takes in a tin in ADCIRC format, and generates an ESRI *.flt
# file for easy visualization by a GIS. It works exactly as my adcirc2asc.py
# script, except that it produces binary files instead of ascii files.
#
# Uses: Python 2 or 3, Matplotlib, Numpy
#
# Example:
#
# python adcirc2flt.py -i tin.grd -s 10 -o tin.flt
# where:
# -i input adcirc mesh file
# -s spacing (in m) of the *.flt grid file
# -o generated *.flt grid file
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Global Imports
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import os,sys # system parameters
import matplotlib.tri as mtri # matplotlib triangulations
import numpy as np # numpy
from ppmodules.readMesh import * # to get all readMesh functions
import struct # to write binary data to file
# uses Rick van Hattem's progressbar
# https://github.com/WoLpH/python-progressbar
from progressbar import ProgressBar, Bar, Percentage, ETA
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# MAIN
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
curdir = os.getcwd()
#
# I/O
if len(sys.argv) != 7 :
print('Wrong number of Arguments, stopping now...')
print('Usage:')
print('python adcirc2flt.py -i tin.grd -s 10 -o tin.flt')
sys.exit()
dummy1 = sys.argv[1]
adcirc_file = sys.argv[2]
dummy2 = sys.argv[3]
spacing = sys.argv[4]
spacing = float(spacing)
dummy3 = sys.argv[5]
output_file = sys.argv[6] # output *.flt grid
# to create the output *.flt file
fout = open(output_file,"wb")
# the name of the header file (this is asci format)
# use reverse split!
header_file = output_file.rsplit('.',1)[0] + '.hdr'
fhdr = open(header_file,"w")
# read the adcirc file
n,e,x,y,z,ikle = readAdcirc(adcirc_file)
# create tin triangulation object using matplotlib
tin = mtri.Triangulation(x, y, ikle)
# determine the spacing of the regular grid
range_in_x = x.max() - x.min()
range_in_y = y.max() - y.min()
max_range = max(range_in_x, range_in_y)
# first index is integer divider, second is remainder
num_x_pts = divmod(range_in_x, spacing)
num_y_pts = divmod(range_in_y, spacing)
print("Size of output matrix is : " + str(int(num_x_pts[0])) + " x " + str(int(num_y_pts[0])))
print("Grid resolution is : " + str(spacing) + " m")
# creates the regular grid
xreg, yreg = np.meshgrid(np.linspace(x.min(), x.max(), int(num_x_pts[0])),
np.linspace(y.min(), y.max(), int(num_y_pts[0])))
x_regs = xreg[1,:]
y_regs = yreg[:,1]
# perform the triangulation
print('Interpolating ...')
interpolator = mtri.LinearTriInterpolator(tin, z)
interp_zz = interpolator(xreg, yreg)
#print "Shape of array z: " + str(interp_zz.shape[0])
#print "Shape of arrays xreg and yreg: " + str(x_regs.shape) + " " + str(y_regs.shape)
where_are_NaNs = np.isnan(interp_zz)
interp_zz[where_are_NaNs] = -999.0
#np.savetxt("temp.out", z, fmt='%.2f', delimiter='') # this has no max spaces
#np.savetxt('temp.out', np.flipud(interp_zz), fmt='%10.3f', delimiter='') # this has 10 char spaces, 2 after decimal
# open the output *.hdr file, and write the header info
fhdr.write("NCOLS " + str(interp_zz.shape[1]) + "\n")
fhdr.write("NROWS " + str(interp_zz.shape[0]) + "\n")
fhdr.write("XLLCORNER " + str(x_regs[0]) + "\n")
fhdr.write("YLLCORNER " + str(y_regs[0]) + "\n")
fhdr.write("CELLSIZE " + str(spacing) + "\n")
fhdr.write("NODATA_VALUE " + str(-999.00) + "\n")
fhdr.write("BYTEORDER LSBFIRST " + "\n")
print('Writing binary data file ...')
'''
# this also works too, but can't use progress bar
num_pts = int(num_x_pts[0]) * int(num_y_pts[0])
interp_zz_reshaped = np.reshape(np.flipud(interp_zz),num_pts)
s = struct.pack('f'*num_pts,*interp_zz_reshaped)
fout.write(s)
fout.close()
'''
w = [Percentage(), Bar(), ETA()]
pbar = ProgressBar(widgets=w, maxval=interp_zz.shape[0]).start()
for i in range(interp_zz.shape[0]):
s = struct.pack('f'*interp_zz.shape[1], *np.flipud(interp_zz)[i,:])
fout.write(s)
pbar.update(i+1)
fout.close()
pbar.finish()
print('All done!')