awips2/cave/com.raytheon.viz.gfe/localization/gfe/userPython/procedures/TCWindThreat.py
2018-06-20 17:39:08 -06:00

614 lines
23 KiB
Python

# ----------------------------------------------------------------------------
# This software is in the public domain, furnished "as is", without technical
# support, and with no warranty, express or implied, as to its usefulness for
# any purpose.
#
# A2TCWindThreatFinal_3
# This version is for use with MLB's Graphical HLS.
# This version modifies the gHLS WindThreat element.
#
# Author: lefebvre 9/11/06
# Modified by: Volkmer/MLB and Santos/MFL 5/23/07
# More modifications: by LeFebvre/Santos/Sharp 05/7/09
# LeFebvre/Santos 07/20/10
# Modified: by Santos 04/26/2011 to accomodate alternate logic for NE Coast and hide bars
# Modified: by Santos 04/20/2012 to get rid off Very Low and tweak GUI.
# Modified: by Shannon/Pablo 06/19/2012 to make A2 and DRT compatible
# Modified: by Pablo/Shannon on 05/21/2014 to fix bug introduced in 2012 when Very Low was eliminated
# Modified: By Santos/Lefebvre 09/17/2014 to make changes for official implementation in 2015.
# Modified: By Belk 07/15/2016 to make efficiency improvements, and
# refactor to make use of a utility containing common methods with other tools
# CHECKED IN FOR 17.1.1: By LeFebvre 09/23/2016 finish converstion to numpy conventions.
# Modified: By LeFebvre 09/26/16 - Removed commented out code to pass code review.
# Modified: By LeFebvre 10/31/16 - Added more code to ensure only one cyclone center point is calculated
# Modified: By LeFebvre 07/18/17 - Added option to populate based on Preliminary or Official prob guidance.
#----------------------------------------------------------------------------
##
# This is an absolute override file, indicating that a higher priority version
# of the file will completely replace a lower priority version of the file.
##
#
# THIS CODE SHALL NOT BE CHANGED WITHOUT CONSULTING ORIGINAL DEVELOPERS.
#
# The MenuItems list defines the GFE menu item(s) under which the
# Procedure is to appear.
# Possible items are: Populate, Edit, Consistency, Verify, Hazards
MenuItems = ["Populate","Edit"]
# The ToolList is optional, but recommended, if you are calling
# Smart Tools from your Script.
# If present, it can be used to show which grids will be
# modified by the Script.
VariableList = [("Probabilistic Wind Source?", "Official", "radio", ["Official", "Preliminary"]),
("Forecast Confidence?", "Typical (Combined; For ill defined or for most systems anytime within 48 hours of landfall or closest approach)",
"radio", ["Typical (Combined; For ill defined or for most systems anytime within 48 hours of landfall or closest approach)",
"High (Combined; for well-behaved systems within 12 hours of landfall or closest approach)",
"Higher (Combined; for very well-behaved systems within 6 hours of landfall or closest approach)",
"Highest (Deterministic-only; assumes a perfect forecast)"]),
("WindMax Source?", "WFO Mesoscale (default)", "radio", ["WFO Mesoscale (default)","NHC/TCM Synoptic"]),
]
import TropicalUtility
import string
import TimeRange
import AbsTime
import time
import MetLib
import sys
import math
import numpy as np
class Procedure (TropicalUtility.TropicalUtility):
def __init__(self, dbss):
TropicalUtility.TropicalUtility.__init__(self, dbss)
# Finds the TPC prob database, extracts all of the grids and returns
# the last one of each type
def getProbGrids(self):
ap = self.availableParms()
searchStr = self.getSiteID() + "_GRID_D2D_" + self._probWindModelSource
probModel = ""
for elem, level, model in ap:
if searchStr in model.modelIdentifier():
probModel = model.modelIdentifier()
break
if probModel == "":
self.statusBarMsg("TPC Wind probability grids not found.", "S")
return None, None, None
# make a big timeRange
timeRange = TimeRange.allTimes()
# get the TPC prob grids, and keep the last one
prob34List = self.getGrids(probModel, "prob34", "FHAG10", timeRange,
mode = "List")
prob34Grid = prob34List[-1]
prob50List = self.getGrids(probModel, "prob50", "FHAG10", timeRange,
mode = "List")
prob50Grid = prob50List[-1]
prob64List = self.getGrids(probModel, "prob64", "FHAG10", timeRange,
mode = "List")
prob64Grid = prob64List[-1]
return prob34Grid, prob50Grid, prob64Grid
# Make a timeRange spanning from the current hour to 6 hours hence
def make6hrTimeRange(self):
cTime = int(self._gmtime().unixTime()/ 3600) * 3600
startTime = AbsTime.AbsTime(cTime)
endTime = AbsTime.AbsTime(cTime + (6 * 3600)) # 6 hours
tr = TimeRange.TimeRange(startTime, endTime)
return tr
# Make a timeRange spanning the usual -24 to 10-day db
def makeDBTimeRange(self):
cTime = int(self._gmtime().unixTime()/ 3600) * 3600
startTime = AbsTime.AbsTime(cTime - (24*3600))
endTime = AbsTime.AbsTime(cTime + (240 * 3600)) # 10 days
tr = TimeRange.TimeRange(startTime, endTime)
return tr
# Removes all grids from a selected set of temporary parms
def removeTempGrids(self):
elementList = ["Prob34", "Prob50", "Prob64", "WindMax"]
ap = self.availableParms()
tr = TimeRange.allTimes()
for elem, level, model in ap:
modelName = model.modelIdentifier()
if elem in elementList and level == "SFC":
self.deleteCmd([elem], tr)
return
# Returns the distance between two points
def pointDist(self, a, b):
return math.sqrt((a[0] - a[1]) * (a[0] - a[1])) + ((b[0] - b[1]) * (b[0] - b[1]))
def distSort(self, a, b):
return cmp(a, b)
# Calculates the average location of the specified points
def avgLocation(self, pointList):
# calc avg
xSum = 0.0
ySum = 0.0
for x, y in pointList:
xSum = xSum + x
ySum = ySum + y
avgX = xSum / len(pointList)
avgY = ySum / len(pointList)
avgLoc = (avgX, avgY)
return avgLoc
# Given a list of (x, y) pairs, determine the true center of "the cluster" of
# points.
def selectSingleCenter(self, pointList):
# If there's only one point, we're done.
if len(pointList) == 1:
return pointList[0]
# calc the distance from each point to the average location of all the points
avgLoc = self.avgLocation(pointList)
#print "avgLoc:", avgLoc
distList = []
for x, y in pointList:
dist = self.pointDist(avgLoc, (x, y))
distList.append((dist, x, y))
# sort smallest to largest distance
distList.sort(self.distSort)
#print "DistanceList after sort:", distList
# Now collect the closest n points
closestPoint = (distList[0][1], distList[0][2])
#print "closestPoint:", closestPoint
# See if the closest point
# if self.pointDist(avgLoc, closestPoint) > 20:
# self.statusBarMsg("Warning! No cluster of points found. Center not reliable", "U")
cluster = [closestPoint]
for d, x, y in distList:
dist = self.pointDist(closestPoint, (x, y))
if dist < 5:
cluster.append((x, y))
avgX, avgY = self.avgLocation(cluster)
avgX = int(avgX + 0.5)
avgY = int(avgY + 0.5)
return avgX, avgY
# Calculates the grid position of the storm center
def getStormCenter(self, windGrid, tr):
u, v = self.MagDirToUV(windGrid[0], windGrid[1])
vortGrid = MetLib.vorticity((u, v))
vortGrid = self.GM_smoothGrid(vortGrid, 3)
vortGrid = np.clip(vortGrid, 0.0, 1000000)
# Find the max vorticity magnitude
#maxVort = np.amax(vortGrid)
maxVort = np.max(vortGrid.flat)
if maxVort < 10:
return None, None
# Find the location of the max (i.e. the center)
mask = (vortGrid == maxVort)
centerList = np.nonzero(mask.flat)[0]
# Convert the list of nonzero points to x, y coordinates and append to a list
pointList = []
for c in centerList:
row = c / mask.shape[1]
col = c % mask.shape[1]
pointList.append((col, row))
# Call this method to select a single center
centerCol, centerRow = self.selectSingleCenter(pointList)
return centerCol, centerRow
# Using the center, a distance grid from that center and a mask
# indicating the area that includes the maxWind field. This
# method returns the radius of that area in gridpoint units.
# Note this method adds 2 gridpoints since the area tends to
# fall a little short of the actual max wind value.
def calcRadius(self, center, distGrid, mask):
maskPoints = np.nonzero(mask)
if len(maskPoints) == 0:
return 0
histDict = {}
yCoords, xCoords = maskPoints
for i in range(len(yCoords)):
dist = int(distGrid[yCoords[i], xCoords[i]])
if dist not in histDict:
histDict[dist] = 1
else:
histDict[dist] = histDict[dist] + 1
histKeys = histDict.keys()
histKeys.sort()
lastKey = 0
for key in histKeys:
if key - lastKey >=2:
return lastKey + 2
lastKey = key
return lastKey + 2
# Makes a grid of distance from the specified center
# which must be expressed as a tuple (x, y)
def makeDistGrid(self, center):
iGrid = np.indices(self.getGridShape())
yDist = center[1] - iGrid[0]
xDist = center[0] - iGrid[1]
distGrid = np.sqrt(xDist**2 + yDist**2)
return distGrid
# Given a specified wind speed and direction grid along with
# a distance grid, this method computes a mask that covers
# the area defined by the ring of max winds.
def makeCoreMask(self, ws, wd, distGrid):
# Find the area inside the max wind ring
u, v = self.MagDirToUV(ws, wd)
windGrad = MetLib.gradient(ws)
distGrad = MetLib.gradient(distGrid)
dotGrid = MetLib.dot(windGrad, distGrad)
# smooth this grid to remove noise
dotGrid = self.GM_smoothGrid(dotGrid, 9)
mask = dotGrid > 0.0
return mask
# This method returns the max wind speed value found in the
# specified wind speed grid. This method restricts the search
# to an area just inside and just outside the specified radius.
def getMaxWindValue(self, distGrid, radius, ws, tr):
maxWS = 0.0
lessMask = distGrid < (radius + 2)
moreMask = distGrid > (radius - 7)
mask = lessMask & moreMask
## self.createGrid("Fcst", "MaxMask", "SCALAR", mask, tr)
ringGrid = ws * mask
thisMax = np.amax(ringGrid)
if thisMax > maxWS:
maxWS = thisMax
return maxWS
# This method adjusts and returns the specified maxWindGrid by
# temporally interpolating each wind grid found in the current
# inventory. A vorticity field is used to determine the storm
# center. Then the wind speed and direction grids are used to
# calculate the area surrounded by the ring of max winds. This
# area is used to calculate the radius of max wind and the all
# of the attributes are used to calculate the max wind field
# over the entire event by interpolating temporally. This max
# wind field is later used to calculate the WindThreat.
def adjustWindMax(self, maxWindGrid):
trList = self.GM_getWEInventory("Wind")
if len(trList) == 0:
return
adjustedWindMax = maxWindGrid
stormMaxWindValue = -1
modifiedMask = np.zeros(self.getGridShape(), np.bool)
lastXPos = None
lastYPos = None
lastRadius = None
lastWindMax = None
for tr in trList:
ws, wd = self.getGrids("Fcst", "Wind", "SFC", tr, mode="First")
xPos, yPos = self.getStormCenter((ws, wd), tr)
# See if the storm is outside the GFE domain. If it is, just update
# the adjustedWindMax and the stormMaxWindValue
if xPos is None or yPos is None:
mask = ws > adjustedWindMax
adjustedWindMax[mask] = ws[mask]
# Update the overall max
gridWindMax = np.amax(ws) # Max value over the whole grid
stormMaxWindValue = max(gridWindMax, stormMaxWindValue)
continue
# calc change in wind speed as a function of radius
distGrid = self.makeDistGrid((xPos, yPos))
coreMask = self.makeCoreMask(ws, wd, distGrid)
# first time through just store the position
if lastXPos == None or lastYPos == None:
lastXPos = xPos
lastYPos = yPos
lastRadius = self.calcRadius((xPos, yPos), distGrid, coreMask)
lastWindMax = self.getMaxWindValue(distGrid, lastRadius, ws, tr)
continue
# get the maxWindRadius in grid cell coordinates
radius = self.calcRadius((xPos, yPos), distGrid, coreMask)
if radius is None:
continue
# Get the max wind value, but restrict it to near the radius
maxWindValue = self.getMaxWindValue(distGrid, radius, ws, tr)
#print "unrestricted max Wind Value: ", maxWindValue
dr = radius - lastRadius
# Calculate distance from last point
dx = xPos - lastXPos
dy = yPos - lastYPos
dWind = maxWindValue - lastWindMax
# One iteration per grid point
steps = abs(max(dx, dy))
for i in range(steps):
x = lastXPos + (float(dx) / steps) * i
y = lastYPos + (float(dy) / steps) * i
r = lastRadius + (float(dr) / steps) * i
windValue = lastWindMax + (float(dWind) / steps) * i
distGrid = self.makeDistGrid((x, y))
# Change grids at the intersection of a circle defined by
# the center and radius and points that whose value is
# less than the current maxWind value
distMask = distGrid <= r
lessMask = adjustedWindMax < windValue
mask = distMask & lessMask # union of dist and less masks
adjustedWindMax[mask] = windValue
modifiedMask = mask | modifiedMask
lastXPos = xPos
lastYPos = yPos
lastRadius = radius
lastWindMax = maxWindValue
# Calculate the maximum wind value over the modified area
stormMaxWindValue = np.amax(adjustedWindMax * modifiedMask)
#print "stormMaxWindValue: ", stormMaxWindValue
# smooth out moderate values to remove "quadrant effect"
adjustedWindMax = self.GM_smoothGrid(adjustedWindMax, 5)
return adjustedWindMax, stormMaxWindValue
# fetch a grid that represents the maxWind everywhere for the next 5 days
def getWindMax(self, timeRange):
cTime = int(self._gmtime().unixTime()/ 3600) * 3600
startTime = AbsTime.AbsTime(cTime)
endTime = startTime + (5 * 24 * 3600) # 5 days from the startTime
tr = TimeRange.TimeRange(startTime, endTime)
try:
windMax, dir = self.getGrids("Fcst", "Wind", "SFC", tr, mode = "Max")
except:
windMax = None
self.statusBarMsg("No Wind grids found. Please define Wind grids first.",
"S")
return windMax
def getMaxWindGrid(self):
try:
grid = self.getObject("TCMMaxWindGrid", "WindGrid")
maxWindValue = np.amax(grid)
return (grid, maxWindValue)
except:
self.statusBarMsg("No TCMMaxWindGrid found.", "S")
return (None, None)
def execute(self, varDict):
# Fetch the model source and define the model name
sourceDB = varDict["Probabilistic Wind Source?"]
if sourceDB == "Official":
self._probWindModelSource = "TPCProb"
elif sourceDB == "Preliminary":
self._probWindModelSource = "TPCProbPrelim"
else:
self.statusBarMsg("Unknown model source selected. Aborting.", "U")
return
# Get confidence value from the dialog
confidenceStr = varDict["Forecast Confidence?"]
tcmwindmax = varDict["WindMax Source?"]
# define a couple of boolean flags that will come in handy later
# allProb = confidenceStr == "Low (Probability-only; 10% Default-05% NE Coast)"
# allWind = confidenceStr == "Highest (Deterministic-only; MaxWind Composite)"
#allProb = confidenceStr == "Low (Probability-only; for ill-behaved systems)"
allWind = confidenceStr == "Highest (Deterministic-only; assumes a perfect forecast)"
typical = confidenceStr == "Typical (Combined; For ill defined or for most systems anytime within 48 hours of landfall or closest approach)"
high = confidenceStr == "High (Combined; for well-behaved systems within 12 hours of landfall or closest approach)"
higher = confidenceStr == "Higher (Combined; for very well-behaved systems within 6 hours of landfall or closest approach)"
#print "allProb and allWind are: ", allProb, allWind
# extract the percent value from this string
# pctPos = confidenceStr.find("% Default")
# pctStr = confidenceStr[pctPos - 2:pctPos]
pctStr="10"
if high:
pctStr="20"
elif higher:
pctStr="30"
# Percent thresholds for each confidence category
threatDict = {"10" : [10.0, 10.0, 10.0, 20.0, 30.0],
"20" : [20.0, 20.0, 20.0, 30.0, 40.0],
"30" : [30.0, 30.0, 30.0, 40.0, 50.0],
}
# wind thresholds for each threat category
windDict = {'None' : 0.0,
# 'Very Low' : 34.0,
#'Elevated' : 50.0,
'Elevated': 34.0,
'Mod' : 50.0,
'High1' : 64.0,
'High2' : 89.0,
'Extreme' : 96.0,
}
# Make sure the string is valid. If not then assign any value since the user
# has indicted Highest confidence in the wind field and is not using the
# probabilistic grids at all.
#if not threatDict.has_key(pctStr):
# pctStr = "10"
#print "pctStr is: ", pctStr
# Extract the proper list and assign thresholds
thresholdList = threatDict[pctStr]
#
t34TS1 = thresholdList[0]
t50TS2 = thresholdList[1]
t64Cat1 = thresholdList[2]
t64Cat2 = thresholdList[3]
t64Cat3 = thresholdList[4]
timeRange = self.make6hrTimeRange()
# Remove previous version of grids.
self.removeTempGrids()
# set up the indices for the discrete keys
keys = self.getDiscreteKeys("WindThreat")
## print "WindThreat keys are:", keys
noneIndex = self.getIndex("None", keys)
lowIndex = self.getIndex("Elevated", keys)
modIndex = self.getIndex("Mod", keys)
highIndex = self.getIndex("High", keys)
extremeIndex = self.getIndex("Extreme", keys)
# Initialize the threat grid
threatGrid = self.empty(np.int8) # a grid of zeros
# Attempt to get the grid from the server.
windMax, maxWindValue = self.getMaxWindGrid()
#print "MAXWIND ACROSS DOMAIN IS: ", maxWindValue
# Use the old-fashioned method
if windMax is None or tcmwindmax == "WFO Mesoscale (default)":
# Get and adjust a grid of maximum wind over the entire storm
windMax = self.getWindMax(timeRange)
windMax, maxWindValue = self.adjustWindMax(windMax)
# Create this grid, if it is valid
if windMax is not None:
self.createGrid("Fcst", "WindMax", "SCALAR", windMax, timeRange,
minAllowedValue=0, maxAllowedValue=200)
# Assign values to the grid based on the probability grids
if allWind:
threatGrid[windMax >= windDict["Elevated"]] = lowIndex
threatGrid[windMax >= windDict["Mod"]] = modIndex
threatGrid[windMax >= windDict["High1"]] = highIndex
threatGrid[windMax >= windDict["Extreme"]] = extremeIndex
else:
# high and extreme threats require maxWind to meet particular windMax criteria
# Fetch the probabilistic grids
prob34Grid, prob50Grid, prob64Grid = self.getProbGrids()
#self.createGrid("Fcst", "Prob34", "SCALAR", prob34Grid, timeRange)
#self.createGrid("Fcst", "Prob50", "SCALAR", prob50Grid, timeRange)
#self.createGrid("Fcst", "Prob64", "SCALAR", prob64Grid, timeRange)
#print "MAXWIND IS: ", maxWindValue
threatGrid[prob34Grid >= t34TS1] = lowIndex
threatGrid[prob50Grid >= t50TS2] = modIndex
threatGrid[prob64Grid >= t64Cat1] = highIndex
if maxWindValue >= windDict['High2']:
threatGrid[prob64Grid >= t64Cat3] = extremeIndex
if maxWindValue >= windDict['Extreme']:
threatGrid[prob64Grid >= t64Cat2] = extremeIndex
#===================================================================
# Upgrade windThreat based on windMax grid
##
# Upgrade None to Elevated
windMask = ((windMax >= windDict['Elevated']) &
(windMax < windDict['Mod']))
threatMask = threatGrid < lowIndex
threatGrid[windMask & threatMask] = lowIndex
# Upgrade Elevated to Med
windMask = ((windMax >= windDict['Mod']) &
(windMax < windDict['High1']))
threatMask = threatGrid < modIndex
threatGrid[windMask & threatMask] = modIndex
# Upgrade Med to High
windMask = ((windMax >= windDict['High1']) &
(windMax < windDict['Extreme']))
threatMask = threatGrid < highIndex
threatGrid[windMask & threatMask] =highIndex
# Upgrade High to Extreme
windMask = windMax >= windDict['Extreme']
threatMask = threatGrid < extremeIndex
threatGrid[windMask & threatMask] = extremeIndex
# Remove previous version of grid.
dbTimes = self.makeDBTimeRange()
self.deleteCmd(['WindThreat'], dbTimes)
# create the threat grid
self.createGrid("Fcst", "WindThreat", "DISCRETE", (threatGrid, keys),
timeRange, discreteKeys=keys, discreteAuxDataLength=5,
discreteOverlap=0)