forked from jasonfleming/pputils
-
Notifications
You must be signed in to change notification settings - Fork 0
/
adcirc2asc.py
executable file
·104 lines (89 loc) · 3.51 KB
/
adcirc2asc.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
#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
# #
# adcirc2asc.py #
# #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
#
# Date: June 29, 2015
#
# Purpose: Script takes in a tin in ADCIRC format, and generates an ESRI *.asc
# file for easy visualization by a GIS.
#
# Uses: Python 2 or 3, Matplotlib, Numpy
#
# Example:
#
# python adcirc2asc.py -i tin.14 -s 10 -o tin.asc
# where:
# -i input adcirc mesh file
# -s spacing (in m) of the *.asc grid file
# -o generated *.asc 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
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# MAIN
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
curdir = os.getcwd()
#
# I/O
if len(sys.argv) != 7 :
print('Wrong number of Arguments, stopping now...')
print('Usage:')
print('python adcirc2asc.py -i tin.14 -s 10 -o tin.asc')
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 *.asc grid
# to create the output file
fout = open(output_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
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
# write the header string
header_str = "NCOLS " + str(interp_zz.shape[1]) + "\n"
header_str = header_str + "NROWS " + str(interp_zz.shape[0]) + "\n"
header_str = header_str + "XLLCORNER " + str(x_regs[0]) + "\n"
header_str = header_str + "YLLCORNER " + str(y_regs[0]) + "\n"
header_str = header_str + "CELLSIZE " + str(spacing) + "\n"
header_str = header_str + "NODATA_VALUE " + str(-999.00)
#np.savetxt("temp.out", z, fmt='%.2f', delimiter='') # this has no max spaces
np.savetxt(output_file, np.flipud(interp_zz), fmt='%10.3f', header = header_str,
comments = '', delimiter='') # this has 10 char spaces, 2 after decimal
print("All done!")