forked from wterwey/SPOUT
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain.py
270 lines (235 loc) · 11.5 KB
/
main.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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
import netCDF4
import math, os, pickle
import numpy as np
import matplotlib.pyplot as plt
import updraftClasses as up
import modelReadIn as modelread
import globalVariables as GV
import found2Dupdraft as f2D
from matplotlib import cm
global modelParams, TrackerParams, modelData
def is_netcdf(filename):
try:
netCDF4.Dataset(filename, 'r')
return True
except:
return False
def count_frames_in_netcdf(filename):
return len(netCDF4.Dataset(filename, 'r').variables["XTIME"])
def list_WRFout_in_path(dir_path):
files = os.listdir(dir_path)
wrfs = []
for f in files:
fpath = os.path.join(dir_path, f)
if not f.startswith('wrfinput'): # ignore wrfinput files
if is_netcdf(fpath):
wrfs.append(fpath)
return sorted(wrfs)
def count_frames_in_wrffiles(fileList):
"""Count number of time frames. Input absolute paths."""
frames = 0
for f in fileList:
frames += count_frames_in_netcdf(f)
return frames
updrafts = []
fileList = list_WRFout_in_path(GV.modelParams['INPUTROOT'])
frames = count_frames_in_wrffiles(fileList)
if frames != len(fileList):
print('File with multiple frames found (WRF: `frames_per_outfile` > 1). '
'Assuming all frames are within first file.')
assume_allframes_firstfile = True
os.makedirs(GV.modelParams['OUTPUTROOT'], exist_ok=True)
outputfile = open(GV.modelParams['OUTPUTROOT'] + '/finalUpdrafts.dat', 'w')
idNum = 0
skipframes = GV.TrackerParams.get('SKIPFRAMES', 0) # ignore first n frames
for t in range(skipframes, frames):
print('time index', t, 'of', frames)
if assume_allframes_firstfile:
currFile = fileList[0]
modelread.modelData_ReadInFromFile_WRF(currFile, time_index=t)
else: # each frame in separate file
currFile = fileList[t]
modelread.modelData_ReadInFromFile_WRF(currFile, time_index=0)
print('Establishing current time updrafts...')
curr3D = []
#-----------------------------------------------------
# Begin searching through the W array for 2D updrafts,
# then look for connections to make the collection of
# 3D updrafts.
#-----------------------------------------------------
# TODO: why does it need to run 20 times with the same level index?
for m in GV.TrackerParams['checkLevels'][::-1]:
print('level', m)
testArray = GV.W[:, :, m]
mVal = m - GV.TrackerParams['MINUPHEIGHT'] + 1
checkVal = GV.TrackerParams['checkValues'][mVal]
wTh = GV.TrackerParams['WTHRES']
for i in range(2, GV.modelParams['NX']-2):
for j in range(2, GV.modelParams['NY']-2):
if (testArray[i, j] >= checkVal):
tempTotal = 0
for x in range(-1, 2):
for y in range(-1, 2):
if (testArray[i+x, j+y] > checkVal):
tempTotal = tempTotal + 1
if (tempTotal >= 5):
print('inner')
if ((testArray[i, j] >= testArray[i-2, j] + wTh) and \
(testArray[i, j] >= testArray[i-1, j] + wTh) and \
(testArray[i, j] >= testArray[i+1, j] + wTh) and \
(testArray[i, j] >= testArray[i+2, j] + wTh) and \
(testArray[i, j] >= testArray[i-1, j-1] + wTh) and \
(testArray[i, j] >= testArray[i, j-1] + wTh) and \
(testArray[i, j] >= testArray[i+1, j-1] + wTh) and \
(testArray[i, j] >= testArray[i, j-2] + wTh) and \
(testArray[i, j] >= testArray[i-1, j+1] + wTh) and \
(testArray[i, j] >= testArray[i, j+1] + wTh) and \
(testArray[i, j] >= testArray[i+1, j+1] + wTh) and \
(testArray[i, j] >= testArray[i, j+2] + wTh)):
tempData = f2D.found2Dupdraft(i, j, m)
tempDraft = up.updraft2D(i, j, m, t, tempData)
if (m == max(GV.TrackerParams['checkLevels'])):
curr3D.append(up.updraft3D([tempDraft], m))
# If at the top level, start new 3D updrafts; if
# not, look for connections above.
else:
connection = 0
tempRank = 1000.0
print('vertrange')
for upTemp3 in curr3D:
if (((upTemp3.vertRange())[0] == m + 1)):
for upTemp2 in upTemp3.linked2DUpdrafts:
if (upTemp2.zPos == m+1):
tempRangeX = abs(i - upTemp2.xPos)
tempRangeY = abs(j - upTemp2.yPos)
if ((tempRangeX <= GV.TrackerParams['SPACETHRES']) and \
(tempRangeY <= GV.TrackerParams['SPACETHRES'])):
currRank = (tempRangeX ** 2.0 + tempRangeY ** 2.0) ** 0.5
if (currRank < tempRank):
tempRank = currRank
temp3Updraft = upTemp3
connection = 1
# if within SPACETHRES, add this point to the 3D updraft, if
# none already matched up.
if (connection == 0):
curr3D.append(up.updraft3D([tempDraft], m))
# If not connected, start a new 3D updraft
else:
temp3Updraft.linked2DUpdrafts.append(tempDraft)
temp3Updraft.bottomHeight = m
print('secondpart')
# Complete the 3D updraft finder for the current time.
newCurr3D = []
for upTemp3 in curr3D:
tempRange = upTemp3.vertRange()
if (tempRange[1]-tempRange[0]+1 >= GV.TrackerParams['MINUPHEIGHT']):
newCurr3D.append(upTemp3)
if (t == 0):
for temp3D in newCurr3D:
updrafts.append(up.updraft4D([temp3D], t, idNum))
idNum = idNum + 1
else:
for temp3D in newCurr3D:
minRank1 = 1000.0
minRank2 = 1001.0
minRank3 = 1002.0
minUp1 = []
minUp2 = []
minUp3 = []
minPts1 = 0
minPts2 = 0
minPts3 = 0
for temp4D in updrafts:
up3Dprior = temp4D.findTime(t-1)
up3Dcurr = temp4D.findTime(t)
if ((up3Dcurr == -1) and (up3Dprior != -1)):
priorRange = up3Dprior.vertRange()
totalUadv = 0.0
totalVadv = 0.0
tempCount = 0
for temp2D in up3Dprior.linked2DUpdrafts:
totalUadv = totalUadv + temp2D.savedData['nearUMean']
totalVadv = totalVadv + temp2D.savedData['nearVMean']
tempCount = tempCount + 1
avgUadv = totalUadv / tempCount
avgVadv = totalVadv / tempCount
totalXadv = 0.0
totalYadv = 0.0
tempCount = 0
gridCorrX = GV.modelParams['DT'] * 3600.0 / GV.modelParams['DX']
gridCorrY = GV.modelParams['DT'] * 3600.0 / GV.modelParams['DY']
for temp2D in up3Dprior.linked2DUpdrafts:
totalXadv = (temp2D.xPos + GV.TrackerParams['velCorrect'] * avgUadv * gridCorrX) + \
totalXadv
totalYadv = (temp2D.yPos + GV.TrackerParams['velCorrect'] * avgVadv * gridCorrY) + \
totalYadv
tempCount = tempCount + 1
advectedX = totalXadv / tempCount
advectedY = totalYadv / tempCount
totalXcur = 0.0
totalYcur = 0.0
tempCount = 0
for temp2D in temp3D.linked2DUpdrafts:
totalXcur = temp2D.xPos + totalXcur
totalYcur = temp2D.yPos + totalYcur
tempCount = tempCount + 1
currentX = totalXcur / tempCount
currentY = totalYcur / tempCount
checkRank = ((currentX - advectedX) ** 2 + (currentY - advectedY) ** 2) ** 0.5
currentRange = temp3D.vertRange()
numRankPts = min([currentRange[1], priorRange[1]]) - \
max([currentRange[0], priorRange[0]]) + 1
if (numRankPts <= 0):
numRankPts = 0
if (checkRank < minRank1):
minRank3 = minRank2
minRank2 = minRank1
minRank1 = checkRank
minUp3 = minUp2
minUp2 = minUp1
minUp1 = temp4D
minPts3 = minPts2
minPts2 = minPts1
minPts1 = numRankPts
elif (checkRank < minRank2):
minRank3 = minRank2
minRank2 = checkRank
minUp3 = minUp2
minUp2 = temp4D
minPts3 = minPts2
minPts2 = numRankPts
elif (checkRank < minRank3):
minRank3 = checkRank
minUp3 = temp4D
minPts3 = numRankPts
#print (minRank1, minPts1), (minRank2, minPts2), (minRank3, minPts3)
if ((minRank1 <= GV.TrackerParams['RANKTHRES']) and \
(minPts1 >= GV.TrackerParams['POINTTHRES'])):
minUp1.linked3DUpdrafts.append(temp3D)
minUp1.endTime = t
elif ((minRank2 <= GV.TrackerParams['RANKTHRES']) and \
(minPts2 >= GV.TrackerParams['POINTTHRES'])):
minUp2.linked3DUpdrafts.append(temp3D)
minUp2.endTime = t
elif ((minRank3 <= GV.TrackerParams['RANKTHRES']) and \
(minPts3 >= GV.TrackerParams['POINTTHRES'])):
minUp3.linked3DUpdrafts.append(temp3D)
minUp3.endTime = t
else:
updrafts.append(up.updraft4D([temp3D], t, idNum))
idNum = idNum + 1
currUpdrafts = []
for temp4D in updrafts:
if (temp4D.endTime != t):
if (temp4D.endTime - temp4D.beginTime + 1 >= GV.TrackerParams['OUTPUTTIMETHRES']):
pickle.dump(temp4D, outputfile)
else:
currUpdrafts.append(temp4D)
updrafts = currUpdrafts
#-----------------------------------------------
# End time/file loop...
#-----------------------------------------------
for temp4D in updrafts:
if(temp4D.endTime - temp4D.beginTime + 1 >= GV.TrackerParams['OUTPUTTIMETHRES']):
pickle.dump(temp4D, outputfile)
outputfile.close()