diff --git a/cave/com.raytheon.viz.gfe/localization/gfe/userPython/procedures/TCMWindTool.py b/cave/com.raytheon.viz.gfe/localization/gfe/userPython/procedures/TCMWindTool.py index d268b70202..3bc8a64d85 100644 --- a/cave/com.raytheon.viz.gfe/localization/gfe/userPython/procedures/TCMWindTool.py +++ b/cave/com.raytheon.viz.gfe/localization/gfe/userPython/procedures/TCMWindTool.py @@ -1,22 +1,3 @@ -## -# This software was developed and / or modified by Raytheon Company, -# pursuant to Contract DG133W-05-CQ-1067 with the US Government. -# -# U.S. EXPORT CONTROLLED TECHNICAL DATA -# This software product contains export-restricted data whose -# export/transfer/disclosure is restricted by U.S. law. Dissemination -# to non-U.S. persons whether in the United States or abroad requires -# an export license or other authorization. -# -# Contractor Name: Raytheon Company -# Contractor Address: 6825 Pine Street, Suite 340 -# Mail Stop B8 -# Omaha, NE 68106 -# 402.291.0100 -# -# See the AWIPS II Master Rights File ("Master Rights File.pdf") for -# further licensing information. -## # ---------------------------------------------------------------------------- # 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 @@ -24,19 +5,29 @@ # # TCMWindTool # -# Version 2.4 July 18, 2006 # Version 2.7.1 2 Sept 2010 Modified to Fix RCL error # Version 2.7.2 01 Feb 2011 Fixed Pie Slice Algorithm/Added Backgroun Options +# Version Last 14 Apr 2014 Added User-editable max winds +# Modified On 22 May 2014 Introduced option for handling asymetry +# in inner core (RMW), option to use 85th wind radii reduction based on +# 2009 paper by DeMaria or the NCST Bias Correction scheme outside radius +# of MaxWind, corrected problem with ring of lower wind value introduced +# at times in the transition between 34 knots wind radii and background +# field, and introduced option to use preliminary TCM message being +# pushed to all offices text databases beginning with 2014 season. +# +# Modified: 1 Jun 2014 to fix bugs with Lat grids and add option +# to use WindReductionFactor grids for Mid Atlantic offices. +# Modified: 6 June 2014 to fix bugs with large reduction factors over land. +# Modified: 9 June 2014 to fix GUI option to run or not over Selected Time Range. +# +# Last Modified: 2 July to fix decoding of PRE TCM files +# Whatever options are needed should be carefully coordinated among +# offices. # # Author: Tom LeFebvre # Contributor: Pablo Santos # ---------------------------------------------------------------------------- -# SOFTWARE HISTORY -# Date Ticket# Engineer Description -# ------------ ---------- ----------- -------------------------- -# 3/6/2013 15658 ryu Merge in change from AWIPS I DR 21414, which fixed -# makeMaxWindGrid() for when center is outside domain. -# Mar 13, 2013 1793 bsteffen Performance improvements for TCMWindTool # The MenuItems list defines the GFE menu item(s) under which the # Procedure is to appear. @@ -47,29 +38,40 @@ VariableList = [("Product to\ndecode:", [], "check", ["preTCM","WRKTCM","TCMAT1", "TCMAT2", "TCMAT3", "TCMAT4", "TCMAT5", "TCMEP1", "TCMEP2", "TCMEP3", "TCMEP4", "TCMEP5"]), ("Product to\n decode:", [], "check", - ["TCMCP1", "TCMCP2", "TCMCP3", "TCMCP4", "TCMCP5", - "TCPWP1", "TCPWP2", "TCPWP3", "TCPWP4", "TCPWP5"]), - ("Background\nModel:", "Fcst", "radio", ["GFS40", "NAM12", "ECMWFHiRes", "Fcst"]), - ("Number of Pie Slices?", "4", "radio", ["4", "8", "12", "16"]), + ["PREAT1", "PREAT2", "PREAT3", "PREAT4", "PREAT5", + "PREEP1", "PREEP2", "PREEP3", "PREEP4", "PREEP5"]), + ("Background\nModel:", "Fcst", "radio", ["GFS0p5degGbl", "UKMET", "ECMWFHiRes", "Fcst"]), + ("Number of Pie Slices?", "16", "radio", ["4", "8", "12", "16", "24"]), ("Eye Diameter:", 0, "scale", [0, 100], 1), ("34 knot radius at 3 days (NM):", 100, "scale", [0, 1000], 10), ("34 knot radius at 4 days (NM):", 100, "scale", [0, 1000], 10), ("34 knot radius at 5 days (NM):", 100, "scale", [0, 1000], 10), - ("Decrease Wind over Land by (%):", 10, "scale", [-20, 50], 1), - ("Make Grids over Selected Time Only:", "No", "radio", ["Yes", "No"]), - ("MaxWind Swath for TCWindThreat?", "No", "radio", ["Yes", "No"]), + ("Decrease Wind over Land by (%):", 15, "scale", [-20, 50], 1), + ("Make Grids over \nSelected Time Only:", "No", "radio", ["Yes", "No"]), + ("MaxWind Swath for \nTCWindThreat?", "No", "radio", ["Yes", "No"]), + ("Define Asymmetrical \nMax Winds?", "No", "radio", ["Yes", "No"]), + ("Reduce Radii by 15% or \n NC State Bias Correction", "Reduce by 15%", + "radio", ["Reduce by 15%", "NC State Bias Correction"]), + ("Constant Land\nReduction (Slider Bar)\nor Wind Reduction\nFactor Grid?", + "Constant", "radio", ["Constant", "Grid"]), ] -import TimeRange -import AbsTime +try: # See if this is the AWIPS I environment + from Numeric import * + import AFPS + AWIPS_ENVIRON = "AWIPS1" +except: # Must be the AWIPS II environment + from numpy import * + import AbsTime + import TimeRange + AWIPS_ENVIRON = "AWIPS2" import SmartScript +import DefineMaxWindGUI +import MetLib -import string, time -from numpy import * - - -## For available commands, see SmartScript +import popen2, string, time, os, cPickle +import Exceptions, types, copy class TCMDecoder: def __init__(self): @@ -95,26 +97,26 @@ class TCMDecoder: "RADIUS OF 050 KT WINDS" : self.decodeJTWCRadii, "RADIUS OF 064 KT WINDS" : self.decodeJTWCRadii, "RADIUS OF 100 KT WINDS" : self.decodeJTWCRadii, - " ---" : self.endJTWCWindForecast, + " ---" : self.endJTWCWindForecast, "REMARKS:" : self.stopDecodingJTWC, } self.fcstList = [] # a place to store all of the forecasts self.text = [] # the text product - + self.currentFcst = {} # the current forecast we are docoding self.baseProductTime = 0 self.foundEyeDiameter = 0 - - self.AltFileName = "" + self.AltFileName = "" + def calcEyeDiameter(self, center, maxWind): lat = center[0] # latitude in degrees - maxWindC = maxWind / 1.944 # convert to meters per second - rmw = 46.29 * exp(-0.0153 * maxWindC + 0.0166 * lat) + maxWind = maxWind / 1.944 # convert to meters per second + rmw = 46.29 * exp(-0.0153 * maxWind + 0.0166 * lat) # convert to diameter and convert from km to nm ed = rmw * 2.0 / 1.852 @@ -122,7 +124,7 @@ class TCMDecoder: def stripText(self): endStr = chr(13) + chr(13) + chr(10) - for i in xrange(len(self.text)): + for i in range(len(self.text)): self.text[i] = string.replace(self.text[i], endStr, "") return @@ -131,15 +133,15 @@ class TCMDecoder: def getBaseProductTime(self): return self.baseProductTime - + def getAltInfoFileName(self): return self.AltFileName - + def currentLine(self): return self.text[self.pos] def nextLine(self): - self.pos += 1 + self.pos = self.pos + 1 if self.pos < len(self.text): return self.text[self.pos] else: @@ -152,7 +154,7 @@ class TCMDecoder: try: return monthList.index(monthStr) + 1 except ValueError: - return 0 + return 0 def convertBaseTime(self, timeStr): # timeStr format: "HHMM UTC DAY MON DD YYYY" @@ -197,10 +199,10 @@ class TCMDecoder: month = tupleTime[1] # see if we crossed over to a new month if tupleTime[2] > day: - month += 1 + month = month + 1 if month > 12: month = 1 - year += 1 + year = year + 1 newTuple = (year, month, day, hour, minute, tupleTime[5], tupleTime[6], tupleTime[7], tupleTime[8]) @@ -220,17 +222,17 @@ class TCMDecoder: for c in hhmm: if not c in string.digits: return - + baseTime = self.convertBaseTime(timeStr) self.baseProductTime = baseTime return - + def decodeAltFileName(self): nameStr = self.currentLine() parts = string.split(nameStr) self.AltFileName = parts[-1] # grab the last string token - #print "I AM HERE AND AltFileName is: ", self.AltFileName + return def decodeCenterLocation(self): @@ -238,7 +240,7 @@ class TCMDecoder: # check for the repeat center....don't want this one if string.find(locStr, "REPEAT") >= 0: return - + keyWord = "NEAR" pos = string.find(locStr, keyWord) if pos > -1: # found it @@ -265,9 +267,9 @@ class TCMDecoder: pos = string.find(presStr, keyWord) if pos > -1: # found it presStr = presStr[pos + len(keyWord):] - + return - + def decodeMaxSustainedWinds(self): keyWord = "MAX SUSTAINED WINDS" windStr = self.currentLine() @@ -275,7 +277,7 @@ class TCMDecoder: if pos > -1: # found it windList = [] tokenList = string.split(windStr) - for i in xrange(len(tokenList)): + for i in range(len(tokenList)): if string.find(tokenList[i], "KT") >= 0: windList.append(float(tokenList[i - 1])) @@ -300,7 +302,7 @@ class TCMDecoder: self.currentFcst['maxWind']) else: # otherwise use what's been defined or read from the text self.currentFcst['eyeDiameter'] = self.defaultEyeDiameter - + return def decodeMaxWind(self): @@ -356,17 +358,17 @@ class TCMDecoder: # store the radii info self.currentFcst['radii'][radiiWindValue] = radiusList - return + return def decodeWindForecast(self): # if we're decoding a new forecast, save the old one first if self.currentFcst != {}: self.fcstList.append(self.currentFcst) self.currentFcst = {} # reset - + str = self.currentLine() str = string.replace(str, '...', ' ') # remove ... - + tokenList = string.split(str) # decode the validTime validTime = self.convert_ddhhmm(tokenList[2], self.baseProductTime) @@ -397,7 +399,7 @@ class TCMDecoder: self.currentFcst['eyeDiameter'] = diameter - # Since we found it in the product, set the default diameter + # Since we found it in the procuct, set the default diameter self.defaultEyeDiameter = diameter self.foundEyeDiameter = 1 # mark that we found it return @@ -409,7 +411,7 @@ class TCMDecoder: self.defaultEyeDiameter = eyeDiameter self.stripText() - + try: while self.pos < len(TCMProduct): line = self.currentLine() @@ -418,16 +420,16 @@ class TCMDecoder: self.keyWordDict[k]() break self.pos = self.pos + 1 - + # store the last forecast in the list of forecasts if self.currentFcst != {}: self.fcstList.append(self.currentFcst) self.currentFcst = {} # reset except: - # Some problem occurred during the decoding process so return an empty fcst + # Some problem occured during the decoding process so return an empty fcst self.baseProductTime = 0 self.fcstList = {} # reset - + return def decodeLatLonToken(self, latLonStr): @@ -456,7 +458,7 @@ class TCMDecoder: self.baseProductTime = int(self.baseProductTime / 3600) * 3600 return None - + def decodeJTWCTimeCenter(self): line = self.nextLine() tokenList = string.split(line) @@ -468,7 +470,7 @@ class TCMDecoder: lat = self.decodeLatLonToken(latStr) lon = self.decodeLatLonToken(lonStr) if lon > 0: - lon -= 360.0 + lon = lon - 360.0 productTime = self.convert_ddhhmm(dateTimeStr, self.baseProductTime) # make a new fcst object to store the analysis @@ -477,7 +479,7 @@ class TCMDecoder: self.currentFcst['centerLocation'] = (lat, lon) self.currentFcst['radii'] = {} self.currentFcst['eyeDiameter'] = self.defaultEyeDiameter - + def decodeJTWCWindForecast(self): line = self.nextLine() @@ -497,7 +499,7 @@ class TCMDecoder: lat = self.decodeLatLonToken(latStr) lon = self.decodeLatLonToken(lonStr) if lon > 0: - lon -= 360.0 + lon = lon - 360.0 # make a new currentFcst and store the info self.currentFcst = {} @@ -505,7 +507,7 @@ class TCMDecoder: self.currentFcst['centerLocation'] = (lat, lon) self.currentFcst['radii'] = {} self.currentFcst['eyeDiameter'] = self.defaultEyeDiameter - return + return def decodeJTWCRadii(self): line = self.currentLine() @@ -528,10 +530,10 @@ class TCMDecoder: radius = float(tokenList[6]) radList = [radius] else: # no RADIUS found so maybe a QUADRANT line - if string.find(line, "QUADRANT") >= 0: + if string.find(line, "QUADRANT") >= 0: radius = float(tokenList[0]) radList.append(radius) - + line = self.nextLine() # save the last radii info @@ -542,22 +544,22 @@ class TCMDecoder: self.fcstList.append(self.currentFcst) self.currentFcst = {} - return + return def endJTWCWindForecast(self): - + if self.currentFcst != {}: self.fcstList.append(self.currentFcst) self.currentFcst = {} - return - + return + def stopDecodingJTWC(self): line = "ZZZZZ" while line != "": line = self.nextLine() return - + # end class TCMDecoder # begin class CircleEA @@ -570,18 +572,18 @@ class CircleEA(SmartScript.SmartScript): self.xDist = (lonGrid - center[1]) * 111.1 * cosLatGrid self.yDist = (latGrid - center[0]) * 111.1 self.distGrid = sqrt(pow(self.xDist, 2)+ pow(self.yDist, 2)) - + self.tanGrid = arctan2(-self.xDist, -self.yDist) # mask off all but the specified quadrant. self.quadList = [] - for quad in xrange(1, slices + 1): + for quad in range(1, slices + 1): minValue = -pi + (quad - 1) * 2 * pi / slices maxValue = -pi + quad * 2 * pi / slices quadrant = logical_and(greater_equal(self.tanGrid, minValue), less(self.tanGrid, maxValue)) self.quadList.append(quadrant) - + return # Return an edit area for just one quadrant. @@ -599,13 +601,62 @@ class CircleEA(SmartScript.SmartScript): def getXYDistGrids(self): return self.xDist, self.yDist - + # end class CircleEA ------------------------------------------------------- class Procedure (SmartScript.SmartScript): def __init__(self, dbss): SmartScript.SmartScript.__init__(self, dbss) + self._dbss = dbss + + # Make a timeRange based on the start and end int times + def makeTimeRange(self, start=0, end=0): + + if AWIPS_ENVIRON == "AWIPS1": + + if start == 0 and end == 0: + return AFPS.TimeRange.allTimes() + + startTime = AFPS.AbsTime(start) + endTime = AFPS.AbsTime(end) + + tr = AFPS.TimeRange(startTime, endTime) + + elif AWIPS_ENVIRON == "AWIPS2": + if start == 0 and end == 0: + startTime = AbsTime.AbsTime(start) + endTime = AbsTime.maxFutureTime() + else: + startTime = AbsTime.AbsTime(start) + endTime = AbsTime.AbsTime(end) + + tr = TimeRange.TimeRange(startTime, endTime) + else: + self.statusBarMsg("Unknown AWIPS version", "U") + tr = None + + return tr + + def getParmTimeConstraints(self, weName, dbName): + + parm = self.getParm(dbName, weName, "SFC") + + if AWIPS_ENVIRON == "AWIPS1": + parmStart = parm.timeConstraints().startTime() + parmDuration = parm.timeConstraints().duration() + parmRepeat = parm.timeConstraints().repeatInterval() + + elif AWIPS_ENVIRON == "AWIPS2": + parmStart = parm.getGridInfo().getTimeConstraints().getStartTime() + parmDuration = parm.getGridInfo().getTimeConstraints().getDuration() + parmRepeat = parm.getGridInfo().getTimeConstraints().getRepeatInterval() + else: + self.statusBarMsg("Unknown AWIPS version", "U") + return None, None, None + + return parmStart, parmDuration, parmRepeat + # Use this method if you have no luck getting products # directly from the text database @@ -620,10 +671,30 @@ class Procedure (SmartScript.SmartScript): f.close() return textList - def printFcst(self, f, baseTime): + # Retrieves a text product from the text database + def getTextProductFromDB(self, productID): + + cmd = "textdb -r " + productID + + # if your path does not include FXA_HOME/bin, + # this line may work instead of the above line. +# cmd = "/awips/fxa/bin/textdb -r " + productID + + (stdout, stdin, stderr) = popen2.popen3(cmd) + + textList = [] + line = stdout.readline() + textList.append(line) + while line != "": + line = stdout.readline() + textList.append(line) + return textList + + def printFcst(self, f, baseTime=None): print "==============================================================" print "Time:", time.asctime(time.gmtime(f['validTime'])), - print "LeadTime:", (f['validTime'] - baseTime) / 3600 + 3 + if baseTime is not None: + print "LeadTime:", (f['validTime'] - baseTime) / 3600 + 3 print "Center:", f['centerLocation'] print "Eye:", f['eyeDiameter'] if f.has_key('maxWind'): @@ -635,15 +706,15 @@ class Procedure (SmartScript.SmartScript): print r, "kts:", f['radii'][r] def getWEInventory(self, modelName, WEName, level): - yesterday = self._gmtime() - (2 * 24 * 3600) # two days ago - later = self._gmtime() + 10 * 24 * 3600 # 10 days from now - allTimes = TimeRange.TimeRange(yesterday, later) - parm = self.getParm(modelName, WEName, level); - inv = parm.getGridInventory(allTimes.toJavaObj()) + allTimes = self.makeTimeRange(0, 0) + gridInfo = self.getGridInfo(modelName, WEName, level, allTimes) trList = [] - for gd in inv: - tr = TimeRange.TimeRange(gd.getGridTime()) + for g in gridInfo: + start = g.gridTime().startTime().unixTime() + end = g.gridTime().endTime().unixTime() + tr = self.makeTimeRange(start, end) trList.append(tr) + return trList def timeRangeSort(self, a, b): @@ -651,15 +722,15 @@ class Procedure (SmartScript.SmartScript): return -1 else: return 1 - + # returns a wind grid from the specified model most closely matched in # time def getClosestWindGrid(self, modelName, bgDict, timeTarget): topo = self.getTopo() calmGrid = self.makeWindGrid(0.0, 0.0, topo.shape) - + if len(bgDict.keys()) == 0: - print "No background grids available...Using calm grid." +# print "No background grids available...Using calm grid." return calmGrid minDiff = 3600 * 24 * 365 # just a large number @@ -669,12 +740,12 @@ class Procedure (SmartScript.SmartScript): # sort the keys by time so we get consistent behavior bgKeys = bgDict.keys() bgKeys.sort(self.timeRangeSort) - + targetTR = self.makeTimeRange(timeTarget, timeTarget + 3600) # figure out which grid is closest in time for invTR in bgKeys: # if we have an exact match, we're done - if invTR.contains(AbsTime.AbsTime(timeTarget)): + if invTR.overlaps(targetTR): tr = invTR # set the tr minDiff = 0 break @@ -711,7 +782,7 @@ class Procedure (SmartScript.SmartScript): convV = cycU u = cycU * cycWt + convU * convWt v = cycV * cycWt + convV * convWt - mag, dir = self.UVToMagDir(u, v) + mag, dir = self.UVToMagDir(u, v) return dir @@ -739,7 +810,7 @@ class Procedure (SmartScript.SmartScript): while len(f2Radii[r]) < 4: f2Radii[r].append(0) - for i in xrange(4): + for i in range(4): r1 = f1Radii[r][i] r2 = f2Radii[r][i] radius = r1 + (r2 - r1) * (newTime - t1) / (t2 - t1) @@ -748,11 +819,12 @@ class Procedure (SmartScript.SmartScript): return newRadii - # interpolates the wind forecasts inbetween the two specified forecasts. + # interpolate the wind forecasts inbetween the two specified forecasts. # interval is assumed to be specified in hours. # returns a new list of forecasts with f1 at the front of the list # and f2 not present at all in the list. def interpolateWindFcst(self, f1, f2, interval): + intSecs = 3600 * interval t1 = f1['validTime'] t2 = f2['validTime'] @@ -776,13 +848,22 @@ class Procedure (SmartScript.SmartScript): dMaxWind = (f2MaxWind - f1MaxWind) / timeSlots f1Radii = f1['radii'] f2Radii = f2['radii'] + + if f1.has_key("editedMaxWinds"): + emw1 = array(f1["editedMaxWinds"]) + emw2 = array(f2["editedMaxWinds"]) + demw = (emw2 - emw1) / timeSlots + fcstList = [f1] # include the first fcst in the list - for i in xrange(1, timeSlots): + for i in range(1, timeSlots): newTime = t1 + (i * intSecs) newLat = f1Lat + (i * dLat) newLon = f1Lon + (i * dLon) newEye = f1Eye + (i * dEye) newMaxWind = f1MaxWind + (i * dMaxWind) + if f1.has_key("editedMaxWinds"): + newEMW = emw1 + (i * demw) + newRadii = self.interpRadii(t1, t2, newTime, f1Radii, f2Radii) f = {} f['centerLocation'] = (newLat, newLon) @@ -790,28 +871,33 @@ class Procedure (SmartScript.SmartScript): f['validTime'] = newTime f['maxWind'] = newMaxWind f['radii'] = newRadii + if f1.has_key("editedMaxWinds"): + f['editedMaxWinds'] = list(newEMW) fcstList.append(f) return fcstList def calcRadiusList(self, maxWind, rmw, rad34, newRadii): - for i in xrange(len(newRadii)): + for i in range(len(newRadii)): # linearly interpolate newRadii[i] = rmw + ((rmw - rad34) / (maxWind - 34.0)) / (64.0 - maxWind) if newRadii[i] < 0: newRadii[i] = 0 return newRadii + # This method fills in radii/wind values one way for the 36-72 hour period # and another way for the 72-120 hour period. The concept is to add more # data values to the wind field so that the wind grids look more realistic. def extrapolateRadii(self, fcstList, baseDecodedTime, radiiFactor): - for i in xrange(1, len(fcstList)): + for i in range(1, len(fcstList)): fcst = fcstList[i] prevFcst = fcstList[i-1] # calc the lead time in hours leadTime = (fcst['validTime'] - baseDecodedTime) / 3600 + 3 + + extRadius = self.getOutlookRadius(leadTime) zeroRadius = extRadius * radiiFactor @@ -836,23 +922,24 @@ class Procedure (SmartScript.SmartScript): prev64 = prevFcst['radii'][64] fcst50 = fcst['radii'][50] newRadii = [0, 0, 0, 0] - for i in xrange(len(prev50)): + for i in range(len(prev50)): if prev50[i] == 0: continue newRadii[i] = fcst50[i] / prev50[i] * prev64[i] - fcst['radii'][64] = newRadii + if not fcst['radii'].has_key(64): + fcst['radii'][64] = newRadii # add in a 5 knot radius for better blending fcst['radii'][5.0] = [zeroRadius, zeroRadius, zeroRadius, zeroRadius] - + elif leadTime > 72: # different algorithm for beyond 72 hours # if there are radii already defined, don't extrapolate new radii if fcst.has_key("radii"): if len(fcst["radii"]) > 0: continue - + # Stuff radii into the rDict to make a cyclone maxWind = 0 if fcst.has_key("maxWind"): @@ -860,24 +947,23 @@ class Procedure (SmartScript.SmartScript): rDict = {} - # add the radii for maxWind at the rmw + # add the radii for maxWind at the rmw if maxWind > 0: # calculate an rmw lat = fcst["centerLocation"][0] # latitude in degrees rmw = 46.29 * exp(-0.0153 * (maxWind / 1.944) + 0.0166 * lat) - rmw /= 1.852 # convert to nautical miles + rmw = rmw / 1.852 # convert to nautical miles for ws in [64.0, 50.0]: newRadii = [0, 0, 0, 0] if ws < maxWind: newRadii = self.calcRadiusList(maxWind, rmw, extRadius, newRadii) rDict[ws] = newRadii - + rDict[34.0] = [extRadius, extRadius, extRadius, extRadius] rDict[5.0] = [zeroRadius, zeroRadius, zeroRadius, zeroRadius] fcst['radii'] = rDict -# print "From extrapolateRadii added rDict:", rDict - + return fcstList # Smooths the specified grid by the specified factor @@ -888,13 +974,19 @@ class Procedure (SmartScript.SmartScript): # factors of less than 3 are useless or dangerous if factor < 3: return grid + + # Specifying the grid type depends on the environment + + typecode = float64 + st = time.time() half = int(factor)/ 2 - sg = zeros(grid.shape,float64) - count = zeros(grid.shape,float64) - gridOfOnes = ones(grid.shape,float64) - for y in xrange(-half, half + 1): - for x in xrange(-half, half + 1): + sg = zeros(grid.shape, typecode) + count = zeros(grid.shape, typecode) + gridOfOnes = ones(grid.shape, typecode) + + for y in range(-half, half + 1): + for x in range(-half, half + 1): if y < 0: yTargetSlice = slice(-y, None, None) ySrcSlice = slice(0, y, None) @@ -916,46 +1008,90 @@ class Procedure (SmartScript.SmartScript): target = [yTargetSlice, xTargetSlice] src = [ySrcSlice, xSrcSlice] - sg[target] += grid[src] - count[target] += gridOfOnes[src] + sg[target] = sg[target] + grid[src] + count[target] = count[target] + gridOfOnes[src] return sg / count # Smooths the direction grid without regard to the magnitude - def smoothDirectionGrid(self, dir, factor): - mag = ones(dir.shape, float) # 1.0 everywhere - u, v = self.MagDirToUV(mag, dir) + def smoothDirectionGrid(self, dirGrid, factor): + mag = ones(dirGrid.shape, float32) # 1.0 everywhere + u, v = self.MagDirToUV(mag, dirGrid) u = self.smoothGrid(u, factor) v = self.smoothGrid(v, factor) - mag, dir = self.UVToMagDir(u, v) - return dir - + mag, dirGrid = self.UVToMagDir(u, v) + return dirGrid + def makeWindGrid(self, mag, direct, gridShape): - mag = ones(gridShape, float) * mag - direct = ones(gridShape, float) * direct + mag = ones(gridShape, float32) * mag + direct = ones(gridShape, float32) * direct return mag, direct - def decreaseWindOverLand(self, grid, fraction, Topo): - mask = greater(Topo, 0.0) - + def decreaseWindOverLand(self, grid, fraction, Topo, timeRange): + + if self.lessOverLandGrid == "Grid": + + windFactorGrid = self.getWindReductionFactorGrid("Fcst", timeRange) + if windFactorGrid is not None: + # Restrict reduction to the cyclone winds defined by the TCM + grid = where(self._cycloneMask, grid * (1 - windFactorGrid), grid) + return grid + else: + # If no grid was found just return the standard reduction + self.statusBarMsg("Wind Reduction Factor grid not found. Using standard reduction." , "S") + # If area over which you desire to apply land correction you prefer be # based on Edit Are instead of areas with Topo greater than zero then # uncomment the next two lines and specify Edit Area to use. - #editArea = self.getEditArea("LAND_EDIT_AREA_NAME_HERE") + #editArea = self.getEditArea("LAND_EDIT_ARE_NAME_HERE") #mask = self.encodeEditArea(editArea) - gridC = where(mask, grid * fraction, grid) - return gridC + # Restrict reduction to the cyclone winds defined by the TCM + mask = logical_and(greater(Topo, 0.0), self._cycloneMask) + + grid = where(mask, grid * fraction, grid) + + return grid + + # fetches and returns all of the wind reduction factor grids in Fcst DB. + def getWindReductionFactorGrid(self, modelName, timeRange): + try: + inv = self.getWEInventory(modelName, "WindReductionFactor", "SFC") + for tr in inv: + if tr.overlaps(timeRange): + WindRedGrid = self.getGrids(modelName, "WindReductionFactor", "SFC", + timeRange, mode="First") + return WindRedGrid + # If no overlapping grids, return None + return None + except: + return None + def getTimeConstraintDuration(self, element): - return self.getParm("Fcst", element, "SFC").getGridInfo()\ - .getTimeConstraints().getDuration() + + parmStart, parmDuration, parmRepeat = self.getParmTimeConstraints(element, "Fcst") + return parmDuration + def getParmMinMaxLimits(self, modelName, weName): + + parm = self.getParm(modelName, weName, "SFC") + + if AWIPS_ENVIRON == "AWIPS1": + return parm.minLimit(), parm.maxLimit() + elif AWIPS_ENVIRON == "AWIPS2": + return parm.getGridInfo().getMinValue(), parm.getGridInfo().getMaxValue() + else: + self.statusBarMsg("Unknown AWIPS version", "U") + return None, None + + return + # returns the maximum allowable wind speed based on NWS directives def getMaxAllowableWind(self, maxWind): - parm = self.getParm("Fcst", "Wind", "SFC") - maxAllowable = parm.getGridInfo().getMaxValue() - return min(maxWind, maxAllowable) + minAllowable, maxAllowable = self.getParmMinMaxLimits("Fcst", "Wind") + return min(maxWind, maxAllowable) + # returns an interpolated radius based on input radii def getOutlookRadius(self, leadTime): leadTimeList = [72, 96, 120] @@ -963,8 +1099,8 @@ class Procedure (SmartScript.SmartScript): if leadTime < leadTimeList[0]: return radiusList[0] - - for i in xrange(1, len(leadTimeList)): + + for i in range(1, len(leadTimeList)): if leadTime < leadTimeList[i]: dt = leadTimeList[i] - leadTimeList[i - 1] dr = radiusList[i] - radiusList[i - 1] @@ -975,98 +1111,54 @@ class Procedure (SmartScript.SmartScript): # Blends the specified grid together def blendGrids(self, windGrid, bgGrid): - # make a mask around the edge - windMag = windGrid[0] - backMag = bgGrid[0] - mag = windMag.copy() + # Combine the two grids using the windGrid for the cyclone and the + # background grid everywhere else. - # if we have a calm cyclone, return the background grid - windMask = greater(mag, 1.0) - if sum(sum(windMask)) == 0: - return bgGrid + windMag, windDir = windGrid + bgMag, bgDir = bgGrid - # now check the background grid - bgMask = greater(backMag, 1.0) - if sum(sum(bgMask)) == 0: - return windGrid + mask = greater_equal(windMag, 34.0) - # make a weightingGrid - upper = 33.9 # start blending at this mag - # stop blending at the average background mag - lower = sum(sum(backMag)) / sum(sum(ones(backMag.shape))) - if lower >= upper: - print "Problem calculating lower and upper ring thresholds." - print "lower = ", lower, "upper = ", upper - return bgGrid + # No background winds inside any defined wind radii + # Add in the point inside the defined wind radii + mask = logical_or(mask, self._cycloneMask) - # calculate the average value over the area where blending will occur - ringMask = logical_and(less_equal(mag, upper), greater_equal(mag, lower)) + magGrid = where(mask, windMag, bgMag) + dirGrid = where(mask, windDir, bgDir) - ringMaskSum = sum(sum(ringMask)) - if ringMaskSum == 0: - print "Problem calculating ringMask. No blending for this grid." - windMag = self.smoothGrid(windMag, 9) - ringMask = logical_and(less_equal(mag, upper), greater_equal(mag, lower)) - windMask = greater(windMag, 5.0) - magGrid = where(windMask, windMag, backMag) - dirGrid = where(windMask, windGrid[1], bgGrid[1]) - return (magGrid, dirGrid) - - avgGrid = where(ringMask, backMag, 0.0) - lower = sum(sum(avgGrid)) / sum(sum(ringMask)) - - # a nearly calm grid means no blending required so return the cyclone - if lower < 1.0: - return windGrid - - wtGrid = where(greater(mag, upper), 1.0, 0.0) - ringMask = logical_and(less(mag, upper), greater(mag, lower)) - wtGrid = where(ringMask, (mag - lower) / (upper - lower), wtGrid) - wtGrid[less(mag, lower)] = 0.0 - wtGrid = self.smoothGrid(wtGrid, 5) - - # calculate the new mag grid - mag *= wtGrid - mag += backMag * (1 - wtGrid) - - # calculate direction grid - onesGrid = ones(mag.shape) - gridU, gridV = self.MagDirToUV(onesGrid, windGrid[1]) - bgU, bgV = self.MagDirToUV(onesGrid, bgGrid[1]) - gridU *= wtGrid - gridU += bgU * (1 - wtGrid) - gridV *= wtGrid - gridV += bgV * (1 - wtGrid) - - # get the dirGrid and toss out the magnitude - magGrid, dirGrid = self.UVToMagDir(gridU, gridV) - - return mag, dirGrid + return magGrid, dirGrid def getLatLonGrids(self): # Try to get them from the fcst database to save time - startTime = AbsTime.current() - 86400 - endTime = AbsTime.current() + 86400 # 1 day - timeRange = TimeRange.TimeRange(startTime, endTime) - latGrid = self.getGrids("Fcst", "latGrid", "SFC", timeRange, - mode = "First", noDataError = 0) - lonGrid = self.getGrids("Fcst", "lonGrid", "SFC", timeRange, - mode = "First", noDataError = 0) - if latGrid != None and lonGrid != None: - return latGrid, lonGrid + try: + trList = self.getWEInventory("Fcst", "latGrid", "SFC") + except: + trList= [] - # make the latGird and lonGrid - latGrid, lonGrid = SmartScript.SmartScript.getLatLonGrids(self) - - # Temporarliy save them in the forecast database - startTime = AbsTime.current() - endTime = AbsTime.current() + 86400 * 7 # 7 days - timeRange = TimeRange.TimeRange(startTime, endTime) + if len(trList) > 0: + timeRange = trList[0] + latGrid = self.getGrids("Fcst", "latGrid", "SFC", timeRange, + mode = "First", noDataError = 0) + lonGrid = self.getGrids("Fcst", "lonGrid", "SFC", timeRange, + mode = "First", noDataError = 0) + if latGrid != None and lonGrid != None: + return latGrid, lonGrid + + # make the lat and lon grids + gridLoc = self.getGridLoc() + + latGrid, lonGrid = MetLib.getLatLonGrids(gridLoc) + + start = int(time.time() / (24 * 3600)) * 24 * 3600 + end = start + (24 * 3600) + timeRange = self.makeTimeRange(start, end) + + # Temporarily save them in the forecast database self.createGrid("Fcst", "latGrid", "SCALAR", latGrid, timeRange, descriptiveName=None, timeConstraints=None, - precision=1, minAllowedValue=0.0, + precision=1, minAllowedValue=-90.0, maxAllowedValue=90.0) - + self.createGrid("Fcst", "lonGrid", "SCALAR", lonGrid, timeRange, descriptiveName=None, timeConstraints=None, precision=1, minAllowedValue=-360.0, @@ -1088,11 +1180,11 @@ class Procedure (SmartScript.SmartScript): interpFactor = pieSlices / len(rList) newList = [] - for i in xrange(-1, len(rList) -1): + for i in range(-1, len(rList) -1): minVal = rList[i] maxVal = rList[i + 1] dVal = (maxVal - minVal) / interpFactor - for f in xrange(interpFactor): + for f in range(interpFactor): radius = minVal + dVal * f # make sure we never exceed the forecast radius ## if radius > minVal: @@ -1103,13 +1195,16 @@ class Procedure (SmartScript.SmartScript): # the list so that it starts at North to conform to convention shift = int(pieSlices / 4) shiftedList = newList[shift:] - shiftedList += newList[:shift] + shiftedList = shiftedList + newList[:shift] newDict[k] = shiftedList return newDict + # fetches and returns all of the wind grids specified by the model # name. Should be called before any new wind grids are created def getBackgroundGrids(self, modelName): + bgDict = {} + siteID = self.getSiteID() if modelName == "Fcst": level = "SFC" @@ -1123,64 +1218,149 @@ class Procedure (SmartScript.SmartScript): level = "FHAG10" elementName = "wind" - inv = self.getWEInventory(modelName, elementName, level) - bgDict = self.getGrids(modelName, elementName, level, inv, - mode="First") + for tr in inv: + bgDict[tr] = self.getGrids(modelName, elementName, level, + tr, mode="First") return bgDict - - def secondsToYYYYMMDDHH(self, baseTime): - gTime = time.gmtime(baseTime) - return time.strftime("%Y%m%d%H", gTime) + def secondsToYYYYMMDDHH(self, baseTime): + # convert the base time to a string + gTime = time.gmtime(baseTime) + yearStr = str(gTime.tm_year) + monthStr = str(gTime.tm_mon) + dayStr = str(gTime.tm_mday) + hourStr = str(gTime.tm_hour) + while len(monthStr) < 2: + monthStr = "0" + monthStr + while len(dayStr) < 2: + dayStr = "0" + dayStr + while len(hourStr) < 2: + hourStr = "0" + hourStr + + baseTimeStr = yearStr + monthStr + dayStr + hourStr + + return baseTimeStr + + # returns the index corresponding to the specified timeStr and fcstHour def findFcst(self, fcstList, fcstHour): - for i in xrange(len(fcstList)): + for i in range(len(fcstList)): validTime = fcstList[i]["validTime"] leadTime = (validTime - self.baseDecodedTime) / 3600 if fcstHour == leadTime: return i - return None + return None + + # Accepts the number of slices to interpolate and a list of defined + # wind values. Returns a new list of length slices with the specified + # windList interpolated to the new resolution. + def interpWindMax(self, slices, windList): + + maxWindList = [0.0] * slices + + quads = len(windList) + ratio = slices / quads + intOffset = int(ratio / 2) + floatOffset = float(intOffset) / ratio + sliceMap = [] + windPos = [0] * len(windList) + + # Figure out the left and right positions for each new slice + for i in range(slices): + left = int((i - int(ratio/2)) / ratio) + if i % ratio == int(ratio/2): + right = left + windPos[left] = i + else: + right = left + 1 + + if right >= quads: + right = right - quads + + sliceMap.append((left, right)) + + # Do the actual interpolation based on the above positions + interpWindList = [] + for i in range(slices): + left, right = sliceMap[i] + + if left == right: + val = windList[left] + absDist = 1.1111 + elif windPos[left] > windPos[right]: + absDist = slices - abs(windPos[right] - windPos[left]) + else: + absDist = abs(windPos[right] - windPos[left]) + + diff = i - windPos[left] + if diff < 0: + diff = slices + diff + val = windList[left] + diff * ((windList[right] - windList[left]) / absDist) + interpWindList.append(val) + + return interpWindList + + + # Calculate the radius of the maxWind based on teh specified eyeDiameter + def maxWindRadius(self, eyeDiameter=None): + + if eyeDiameter is None: + return 12.5 + + rmw = (eyeDiameter / 2.0) + 8.0 + + return rmw + + def adjustMaxWind(self, outSpeed, inSpeed, outRadius, inRadius, + globalMaxWind, maxWindList, maxWindRadius, quad): + + maxWind = maxWindList[quad] + + # check which speed/radius should be modified + if outSpeed == globalMaxWind: + outSpd = maxWind + outRad = maxWindRadius + inSpd = inSpeed + inRad = inRadius + elif inSpeed == globalMaxWind: + inSpd = maxWind + inRad = maxWindRadius + outSpd = outSpeed + outRad = outRadius + else: + print "ERROR!!! Neither inSpeed or outSpeed is max!!!" + + return outSpd, inSpd, outRad, inRad, maxWind # Makes a Rankine Vortex wind speed grid that decreases exponentially # from the known values at known radii. Inside the Radius of maximum # wind the wind decreases linearly toward the center - def makeRankine(self, f, latGrid, lonGrid, pieSlices, radiiFactor): + def makeRankine(self, f, latGrid, lonGrid, pieSlices, radiiFactor, timeRange): st = time.time() grid = zeros(latGrid.shape) rDict = f['radii'] - -## print "rDict before interpolating Quads:" -## for r in rDict: -## print rDict[r] rDict = self.interpolateQuadrants(rDict, pieSlices) - -## print "rDict after interpolating Quads:" -## for r in rDict: -## print rDict[r] -## -## return (None, None) validTime = f['validTime'] center = f['centerLocation'] maxWind = f['maxWind'] + circleEA = CircleEA(latGrid, lonGrid, center, pieSlices) # make a list that contains the highest non-zero radius speed -# centerWindList = [0, 0, 0, 0] centerWindList = [0] * pieSlices for k in rDict.keys(): - for i in xrange(len(rDict[k])): + for i in range(len(rDict[k])): if rDict[k][i] > 0 and k > centerWindList[i]: centerWindList[i] = k - for k in rDict.keys(): -# if rDict[k] == [0, 0, 0, 0]: - if rDict[k] == [0] * pieSlices: - del rDict[k] + for k in rDict.keys(): + if rDict[k] == [0] * pieSlices: + del rDict[k] # make a list of lowest wind speed found with zero radius # and save the next lowest wind speed for later if rDict.has_key(100.0): @@ -1188,82 +1368,112 @@ class Procedure (SmartScript.SmartScript): else: speedList = [None, 64.0, 50.0, 34.0, 5.0] -# zeroRadList = [999, 999, 999, 999] -# validRadList = [999, 999, 999, 999] zeroRadList = [999] * pieSlices validRadList = [999] * pieSlices - for s in xrange(len(speedList) - 1): + for s in range(len(speedList) - 1): speed = speedList[s] nextSpeed = speedList[s + 1] if not rDict.has_key(speed): -# zeroRadList = [speed, speed, speed, speed] -# validRadList = [nextSpeed, nextSpeed, nextSpeed, nextSpeed] zeroRadList = [speed] * pieSlices validRadList = [nextSpeed] * pieSlices else: - for i in xrange(len(rDict[speed])): + for i in range(len(rDict[speed])): if rDict[speed][i] == 0: zeroRadList[i] = speed validRadList[i] = nextSpeed # get the distance grid and make sure it's never zero anywhere distanceGrid = circleEA.getDistanceGrid() / 1.852 # dist in NM - distanceGrid[equal(distanceGrid, 0)] = 0.01 + distanceGrid[distanceGrid == 0] = 0.01 # make a grid into which we will define the wind speeds - grid = zeros(latGrid.shape, float) + grid = zeros(latGrid.shape, float32) # The cyclone algorithm depends on the forecast lead time fcstLeadTime = (validTime - self.baseDecodedTime) / 3600 + # add the radius for maxWind for interpolation if f.has_key('eyeDiameter'): - maxRad = f['eyeDiameter'] / 2.0 + 8.0 + eyeDiameter = f['eyeDiameter'] else: print "Error --- no eye diameter found." - maxRad = 12.5 # half of default 25 nm eye diameter + eyeDiameter = None -# rDict[maxWind] = [maxRad, maxRad, maxRad, maxRad] + maxWindRadius = self.maxWindRadius(eyeDiameter) + + maxWindList = [] + # add the edited maxWind values, if any + if f.has_key("editedMaxWinds"): + # First interpolate based on pie slices + maxWindList = self.interpWindMax(pieSlices, f["editedMaxWinds"]) + + # Add in the maxWind and radius as a point if not rDict.has_key('maxWind'): - rDict[maxWind] = [maxRad] * pieSlices - - # make sure no wind exceeds the maximum allowable windspeed - # remove any rDict entries that are higher than this - + rDict[maxWind] = [maxWindRadius] * pieSlices # extract the list and sort it wsList = rDict.keys() wsList.sort() - + # insert a dummy wind near the center and append so it's done last -# rDict[1] = [1, 1, 1, 1] rDict[1] = [1] * pieSlices - wsList.append(1.0) + wsList.append(1.0) # insert an artificial 5 knot radius at a distance proportional # to the 34 knot radius for that quadrant -# tenKnotRadiusList = [108.0, 108.0, 108.0, 108.0] # 200 km tenKnotRadiusList = [108.0] * pieSlices + if rDict.has_key(34.0): tenKnotRadList = [] radList34 = rDict[34.0] for r in radList34: tenKnotRadList.append(r * radiiFactor) - - # insert the 10 knot radius at the beginning so is made first + + # insert the 5 knot radius at the beginning so is made first rDict[5.0] = tenKnotRadList wsList.insert(0, 5.0) + insideRMWMask = zeros(latGrid.shape) + self._cycloneMask = zeros(latGrid.shape, bool) # for each rDict record and quadrant, make the grid one piece at a time - for i in xrange(len(wsList) - 1): + for i in range(len(wsList) - 1): + self.lastRadius = [None] * pieSlices if not rDict.has_key(wsList[i]): continue radiusList = rDict[wsList[i]] nextRadiusList = rDict[wsList[i + 1]] - for quad in xrange(len(radiusList)): + + maxRadius = maxWindRadius # temp copy + for quad in range(len(radiusList)): + + maxRadius = maxWindRadius # temp copy + maxWind = f['maxWind'] # reset maxWind as we may fiddle with it + + # fetch the speeds and radii we'll need outSpeed = float(wsList[i]) inSpeed = float(wsList[i + 1]) outRadius = float(radiusList[quad]) inRadius = float(nextRadiusList[quad]) - # reset the speeds if they exceed the max non-zero wind + + # Here's where the speeds and radii are adjusted based + # on the edited values but only if they have been edited + # and only if we're working on the maxWind. + if f.has_key("editedMaxWinds"): + if maxWind == wsList[i] or maxWind == wsList[i+1]: + outSpeed, inSpeed, outRadius, inRadius, maxWind = \ + self.adjustMaxWind(outSpeed, inSpeed, outRadius, + inRadius, maxWind, maxWindList, + maxWindRadius, quad) + + # Some cases require we adjust the maxWindRadius + if outSpeed in [64.0, 50.0, 34.0] and outRadius <= maxWindRadius: + inRadius = outRadius * 0.9 + self.lastRadius[quad] = outRadius + elif inSpeed == 1.0 and self.lastRadius[quad] is not None: + outRadius = self.lastRadius[quad] * 0.9 + #print "Adjusting MaxWindRadius at:", inSpeed, "kts" + self.lastRadius[quad] = None + + # reset the speeds if they exceed the maxWind if fcstLeadTime <= 72 and zeroRadList[quad] is not None: if inSpeed >= zeroRadList[quad]: inSpeed = validRadList[quad] @@ -1286,7 +1496,6 @@ class Procedure (SmartScript.SmartScript): inRadius = 0.1 if outRadius == 0.0: outRadius = 0.1 - # no wind speed can never exceed the maximum allowable wind speed if inSpeed > maxWind: inSpeed = maxWind @@ -1299,32 +1508,48 @@ class Procedure (SmartScript.SmartScript): if inRadius > outRadius: continue + if inSpeed == 0.0 or outSpeed == 0.0: + continue # calculate the exponent so that we exactly fit the next radius denom = log10(inRadius / outRadius) if denom == 0: exponent = 1.0 - else: + else: exponent = (log10(outSpeed) - log10(inSpeed)) / denom - + # make sure the exponent behaves itself if exponent > 10.0: exponent = 10.0 # inside RMW gets a linear slope to largest of max wind forecasts if inRadius <= 1.0: dSdR = (outSpeed - inSpeed) / (outRadius - inRadius) - grid[mask] = inSpeed + (dSdR * distanceGrid[mask]) + grid = where(mask, inSpeed + (dSdR * distanceGrid), grid) + insideRMWMask = logical_or(insideRMWMask, mask) else: # outside RMW - grid[mask] = inSpeed * power((inRadius / distanceGrid[mask]), exponent) -# grid = clip(grid, 0.0, 200.0) + grid = where(mask, inSpeed * power((inRadius / distanceGrid), exponent), + grid) + if outSpeed >= 34.0 and inSpeed >= 34.0: + self._cycloneMask = logical_or(self._cycloneMask, mask) + + # Apply the NC State correction outside the RMW + if self._applyNCSCorrection: + corrGrid = self.makeCorrectionGrid(latGrid, lonGrid, center) +## self.createGrid("Fcst", "NCSCorr", "SCALAR", corrGrid, self._timeRange, +## precision=3, minAllowedValue=-1.0, maxAllowedValue=1.0) + + grid = where(logical_not(insideRMWMask), grid * (1 - corrGrid), grid) + + maxWind = f['maxWind'] # reset again before clipping + dirGrid = self.makeDirectionGrid(latGrid, lonGrid, center[0], center[1]) + # clip values between zero and the maximum allowable wind speed maxWind = self.getMaxAllowableWind(maxWind) grid = clip(grid, 0.0, maxWind) # apply the wind reduction over land fraction = 1.0 - (self.lessOverLand / 100.0) - grid = self.decreaseWindOverLand(grid, fraction, self.elevation) - + grid = self.decreaseWindOverLand(grid, fraction, self.elevation, timeRange) return (grid, dirGrid) def makeMaxWindGrid(self, interpFcstList, interval, latGrid, lonGrid, pieSlices, @@ -1333,22 +1558,19 @@ class Procedure (SmartScript.SmartScript): ## if len(interpFcstList) == 0: ## return - maxWindGrid = zeros_like(self.getTopo()) + maxWindGrid = zeros(self.getTopo().shape) startTime = interpFcstList[0]["validTime"] endTime = startTime + (123 * 3600) # 123 hours later - - start = AbsTime.AbsTime(startTime) - end = AbsTime.AbsTime(endTime) - timeRange = TimeRange.TimeRange(start, end) - + timeRange = self.makeTimeRange(startTime, endTime) + # Used getGrids to calculate the maximum wind grid. # # Fetch the max of the wind grids just generated as this is very fast. maxWindGrid, maxDirGrid = self.getGrids("Fcst", "Wind", "SFC", timeRange, mode="Max") - maxWindGrid = self.smoothGrid(maxWindGrid, 3) + maxWindGrid = self.smoothGrid(maxWindGrid,3) self.createGrid("Fcst", "TCMMaxWindComposite", "SCALAR", maxWindGrid, timeRange, precision=1, minAllowedValue=0.0, maxAllowedValue=200.0) @@ -1358,24 +1580,7 @@ class Procedure (SmartScript.SmartScript): return - def validateCycloneForecast(self, fcstList, baseTime): - print "Validating Forecasts:" -## -## This section is overly restrictive, so it is commented out pending -## better logic that compares the TCM to the RCL forecasts -## -## leadTimeList = [12, 24, 36, 48, 72, 96, 120] -## fcstHourList = [] -## for f in fcstList: -## # calc the leadTime in hours -## leadTime = (f['validTime'] - baseTime) / 3600 + 3 -## fcstHourList.append(int(leadTime)) -## -## # Make sure all of the forecasts are there -## for leadTime in leadTimeList: -## if leadTime not in fcstHourList: -## return False # Now check each forecast to make sure that we have a radius for any # standard wind values less than the maxWind @@ -1394,15 +1599,93 @@ class Procedure (SmartScript.SmartScript): return True + # Returns a dictionary that lists the min and max allowed wind for each hour + def makeWindDict(self, fcstList): + + windDict = {} + + + for f in fcstList: + windValues = f["radii"].keys() + hour = (f["validTime"] - self.baseDecodedTime) / 3600 + maxWind = f["maxWind"] + minWind = 999999.0 + if len(f["radii"].keys()) == 0: + minWind = 0.0 + + # Grab the first (highest) forecast wind speed value + if len(windValues) > 0: + minWind = windValues[0] + else: + minWind = 0.0 + + windDict[hour] = (minWind, maxWind) + + return windDict + + # Pop up a GUI that will maxWind values for each quadrant and time + def launchMaxWindGUI(self, fcstList): + + windDict = self.makeWindDict(fcstList) + if AWIPS_ENVIRON == "AWIPS1": + eaMgr = self.eaMgr() + else: + eaMgr = None + + self._maxWindGUI = DefineMaxWindGUI.DefineMaxWindGUI(self._dbss, eaMgr) + + newMaxWinds = self._maxWindGUI.displayGUI(windDict) + + if newMaxWinds is not None: + + hourList = newMaxWinds.keys() + hourList.sort() + + self._maxWindGUI.cancelCommand() + + return newMaxWinds + + # Make the NCState bais correction grid based on the forecast. + def makeCorrectionGrid(self, latGrid, lonGrid, center): + + + # structure to hold the polynomial coefficients + coeff = [[1.282e-011, -3.067e-008, 2.16e-005, -5.258e-003, 3.794e-001], + [3.768e-011, -4.729e-008, 2.097e-005, -3.904e-003, 2.722e-001], + [4.692e-011, -5.832e-008, 2.565e-005, -4.673e-003, 2.952e-001], + [3.869e-011, -4.486e-008, 1.84e-005, -3.331e-003, 2.738e-001]] + + # make the circle edit area and distance grid + pieSlices = 4 + circleEA = CircleEA(latGrid, lonGrid, center, pieSlices) + + dist = circleEA.getDistanceGrid() # dist in km + + corrGrid = zeros(dist.shape) + + for quad in range(pieSlices): + + ea = circleEA.getQuadrant(quad + 1, 500.0) + grid = coeff[quad][0] * pow(dist, 4) + coeff[quad][1] * pow(dist, 3) + \ + coeff[quad][2] * pow(dist, 2) + coeff[quad][3] * dist + \ + coeff[quad][4] + + corrGrid = where(ea, grid, corrGrid) + + return corrGrid + def execute(self, varDict, timeRange): + + RADII_FACTOR = 4.5 + self.setToolType("numeric") self.toolTimeRange = timeRange # define the default eye diameter for bulletins where they are missing eyeStr = varDict["Eye Diameter:"] self.dialogEyeDiameter = float(eyeStr) - maxwindswath = varDict["MaxWind Swath for TCWindThreat?"] - + maxwindswath = varDict["MaxWind Swath for \nTCWindThreat?"] + Topo = self.getTopo() tcDuration = self.getTimeConstraintDuration("Wind") @@ -1413,16 +1696,16 @@ class Procedure (SmartScript.SmartScript): # get the product ID productList1 = varDict["Product to\ndecode:"] productList2 = varDict["Product to\n decode:"] - productList1C = productList1 + productList2 # concatenate - if len(productList1C) != 1: + productList1 = productList1 + productList2 # concatenate + if len(productList1) != 1: self.statusBarMsg("Please select one TCM bulletin only.", "S") return None - productID = productList1C[0] + productID = productList1[0] # get the ID for this site siteID = self.getSiteID() - + bgModelName = varDict["Background\nModel:"] self.day3Radius = varDict["34 knot radius at 3 days (NM):"] self.day4Radius = varDict["34 knot radius at 4 days (NM):"] @@ -1435,29 +1718,29 @@ class Procedure (SmartScript.SmartScript): # forecast into more radial pieces. Recommended alternative values: # 12, 20, 36, 72. pieSlices = int(varDict["Number of Pie Slices?"]) - print "Number of Pie Slices: ", pieSlices - # pieSlices = 8 # define radii factor - may make this configurable - radiiFactor = 1.5 + # Multiply 3-5 day radius by this factor to get the zero radius. + # Smaller values ramp the cyclone down to zero more quickly. + self.lessOverLand = int(varDict["Decrease Wind over Land by (%):"]) + self.lessOverLandGrid = varDict["Constant Land\nReduction (Slider Bar)\nor Wind Reduction\nFactor Grid?"] self.elevation = Topo - rclDecoder = TCMDecoder() tcmDecoder = TCMDecoder() msg = "" -# # Use this method to fetch a product from the text database + # Fetch the text product if productID == "preTCM": - textProduct = self.getTextProductFromFile("/tmp/Isaac.txt") + textProduct = self.getTextProductFromFile("/tmp/Wilma.txt") decoder = TCMDecoder() decoder.decodeTCMProduct(textProduct, self.dialogEyeDiameter) fcstList = decoder.getFcstList() baseTime = decoder.getBaseProductTime() #elif productID == "WRKTCM": - # textProduct = self.getTextProductFromFile("/data/local/research/TPCWindProb/WRKTCM") + # textProduct = self.getTextProductFromFile("/data/local/research/TPCWindProb/WRKTCM") else: # try fetching the RCL first. rclProductID = "MIARCL" + productID[3:] @@ -1473,8 +1756,12 @@ class Procedure (SmartScript.SmartScript): rclFcstList = rclDecoder.getFcstList() rclBaseTime = rclDecoder.getBaseProductTime() completeFcst = self.validateCycloneForecast(rclFcstList, rclBaseTime) - + + if productID[:3] == "PRE": + productID = "MIA" + productID + tcmTextProduct = self.getTextProductFromDB(productID) + if len(tcmTextProduct) < 5: msg = productID + " could not be retrieved from the text database." self.statusBarMsg(msg, "S") @@ -1484,10 +1771,10 @@ class Procedure (SmartScript.SmartScript): tcmFcstList = tcmDecoder.getFcstList() tcmBaseTime = tcmDecoder.getBaseProductTime() - print "TCM and RCL Base Times are: ", tcmBaseTime, rclBaseTime + #print "TCM and RCL Base Times are: ", tcmBaseTime, rclBaseTime if not completeFcst or rclBaseTime != tcmBaseTime: msg = "Problem decoding " + rclProductID + " Used TCM to make cyclone.\n" - msg += " Used GUI sliders for 3, 4, 5 day forecast." + msg = msg + " Used GUI sliders for 3, 4, 5 day forecast." #self.statusBarMsg(msg, "S") fcstList = tcmFcstList baseTime = tcmBaseTime @@ -1496,51 +1783,58 @@ class Procedure (SmartScript.SmartScript): fcstList = rclFcstList baseTime = rclBaseTime productID = rclProductID - -## # Use this method to fetch a text product from a file -## textProduct = self.getTextProductFromFile("/tmp/Irene.txt") -## -## decoder = TCMDecoder() -## decoder.decodeTCMProduct(textProduct, self.dialogEyeDiameter) -## fcstList = decoder.getFcstList() -## baseTime = decoder.getBaseProductTime() - + print "Decoded:", len(fcstList), " forecasts." # Set the baseDecodedTime - validTime of first entry - 3 hours if len(fcstList) > 0: self.baseDecodedTime = fcstList[0]['validTime'] - 3 * 3600 - - fcstList = self.extrapolateRadii(fcstList, baseTime, radiiFactor) - # See if the decoded fcst is close to the current time. This is needed - # so the tool will work on archived data sets (testMode) - testMode = 0 - if abs(time.time() - self.baseDecodedTime) > 2 * 24 * 3600: # older than 2 days - testMode = 1 + if varDict["Define Asymmetrical \nMax Winds?"] == "Yes": + + newMaxWinds = self.launchMaxWindGUI(fcstList) + for i in range(len(fcstList)): + fcstHour = (fcstList[i]['validTime'] - baseTime) / 3600 + 3 + maxList = newMaxWinds[fcstHour] + fcstList[i]["editedMaxWinds"] = maxList + + fcstList = self.extrapolateRadii(fcstList, baseTime, RADII_FACTOR) + +## # See if the decoded fcst is close to the current time. This is needed +## # so the tool will work on archived data sets (testMode) +## testMode = False +## if abs(time.time() - self.baseDecodedTime) > 2 * 24 * 3600: # older than 2 days +## testMode = True + # restrict grids to the selected time period if option is selected. - restrictAnswer = varDict["Make Grids over Selected Time Only:"] + testMode = False + restrictAnswer = varDict["Make Grids over \nSelected Time Only:"] if restrictAnswer == "Yes": - testMode = 1 + testMode = True # Turn off testMode if the selected timeRange is less than an hour in duration if self.toolTimeRange.duration() < 3600: - testMode = 0 + testMode = False # interpolate the wind forecasts we got from the decoder - print "Interpolating wind forecasts:" selectedStartTime = self.toolTimeRange.startTime().unixTime() selectedEndTime = self.toolTimeRange.endTime().unixTime() - interpFcstList = [] - for i in xrange(len(fcstList) - 1): + for i in range(len(fcstList) - 1): newFcstList = self.interpolateWindFcst(fcstList[i], fcstList[i+1], interval) + # Make sure the fcst is within the selected time range or we're in testMode for f in newFcstList: if (testMode and (f['validTime'] >= selectedStartTime and \ - f['validTime'] < selectedEndTime)) or testMode == 0: + f['validTime'] < selectedEndTime)) or (not testMode): interpFcstList.append(f) + + # append the very last forecast on to the end of the interpolated list + if len(fcstList) > 0: + if (testMode and (f['validTime'] >= selectedStartTime and \ + f['validTime'] < selectedEndTime)) or (not testMode): + interpFcstList.append(fcstList[-1]) if len(fcstList) == 1: interpFcstList = fcstList @@ -1549,47 +1843,74 @@ class Procedure (SmartScript.SmartScript): self.statusBarMsg("No cyclone forecasts found within the Selected TimeRange", "S") else: + # If the wind grids are more than 3 hours long, the first grid ends up being double + # duration. So, add an extra duplicate forecast at the beginning and reset + # the validTime + print "tcHours:", tcHours + if tcHours > 3: + interpFcstList.insert(0, copy.deepcopy(interpFcstList[0])) + interpFcstList[0]["validTime"] = (int(interpFcstList[0]["validTime"] / tcDuration) \ + * tcDuration) + interpFcstList[1]["validTime"] = (int(interpFcstList[0]["validTime"] / tcDuration) \ + * tcDuration) + tcDuration + print "Adjusted time for first forecast" print "Generating", len(interpFcstList), "wind grids" - + # get the lat, lon grids latGrid, lonGrid = self.getLatLonGrids() - + + self._applyNCSCorrection = False + if varDict["Reduce Radii by 15% or \n NC State Bias Correction"] == "Reduce by 15%": + # Reduce the extent of the wind radii per Mark De Maria's research + # Loop through each wind radius and modify in place + for f in interpFcstList: + for windValue in f["radii"]: + for i in range(len(f["radii"][windValue])): + f["radii"][windValue][i] = f["radii"][windValue][i] * 0.85 + elif varDict["Reduce Radii by 15% or \n NC State Bias Correction"] == "NC State Bias Correction": + self._applyNCSCorrection = True + # make a grid for each interpolate forecast gridCount = 0 for f in interpFcstList: - t1 = time.time() - windGrid = self.makeRankine(f, latGrid, lonGrid, pieSlices, radiiFactor) - print "time to makeRankine:", time.time() - t1 + + self._timeRange = timeRange validTime = int(f['validTime'] / 3600) * 3600 bgGrid = self.getClosestWindGrid(bgModelName, bgDict, validTime) + startTime = validTime + endTime = validTime + (interval * 3600) + timeRange = self.makeTimeRange(startTime, endTime) + self._cycloneTimeRange = timeRange - start = AbsTime.AbsTime(int(validTime)) - timeRange = TimeRange.TimeRange(start, start + interval * 3600) + t1 = time.time() + windGrid = self.makeRankine(f, latGrid, lonGrid, pieSlices, RADII_FACTOR, timeRange) + print "Time to makeRankine:", time.time() - t1 - grid = self.blendGrids(windGrid, bgGrid) + magGrid, dirGrid = self.blendGrids(windGrid, bgGrid) + magGrid = self.smoothGrid(magGrid, 5) + dirGrid = self.smoothDirectionGrid(dirGrid, 5) name = "Wind" - self.createGrid("Fcst", name, "VECTOR", grid, timeRange, + self.createGrid("Fcst", name, "VECTOR", (magGrid, dirGrid), timeRange, descriptiveName=None, timeConstraints=None, precision=1, minAllowedValue=0.0, maxAllowedValue=200.0) - - gridCount += 1 - print "TCMWindTool:", productID, "- Generated",gridCount, \ - "out of", len(interpFcstList), "grids", time.asctime(time.gmtime(timeRange.startTime().unixTime())) + + gridCount = gridCount + 1 + print "TCMWindTool:", productID, "- Generated", gridCount, \ + "out of", len(interpFcstList), "grids", \ + time.asctime(time.gmtime(timeRange.startTime().unixTime())) # interpolate through forecast period to very high resolution and make # a composite maxWind grid from those wind grids if maxwindswath is "Yes": t1 = time.time() self.makeMaxWindGrid(interpFcstList, interval, latGrid, lonGrid, pieSlices, - radiiFactor) + RADII_FACTOR) print time.time() - t1, "seconds to generate Max wind composite." if msg != "": self.statusBarMsg(msg, "S") - + return None - - diff --git a/cave/com.raytheon.viz.gfe/localization/gfe/userPython/utilities/DefineMaxWindGUI.py b/cave/com.raytheon.viz.gfe/localization/gfe/userPython/utilities/DefineMaxWindGUI.py new file mode 100644 index 0000000000..df96b09244 --- /dev/null +++ b/cave/com.raytheon.viz.gfe/localization/gfe/userPython/utilities/DefineMaxWindGUI.py @@ -0,0 +1,225 @@ +# ---------------------------------------------------------------------------- +# 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. +# +# DefineMaxWindGUI +# +# Author: +# ---------------------------------------------------------------------------- + +try: # See if this is the AWIPS I environment + from Numeric import * + import AFPS + AWIPS_ENVIRON = "AWIPS1" +except: # Must be the AWIPS II environment + from numpy import * + import AbsTime + import TimeRange + AWIPS_ENVIRON = "AWIPS2" + +import SmartScript +import Tkinter +import tkFont + + +class DefineMaxWindGUI(SmartScript.SmartScript): + + def __init__(self, dbss, eaMgr=None): + SmartScript.SmartScript.__init__(self, dbss) + + if AWIPS_ENVIRON == "AWIPS1": + self.setUp(eaMgr) + + def reportQuadError(self): + self.statusBarMsg("Only three quadants at a time may be reduced.\n" + \ + "...Please toggle another quadrant off first.","S") + return + + def toggleButton(self, buttonLabel): + b = self._buttonLabels.index(buttonLabel) + if self._buttonState[b]: # button was on + self._buttonState[b] = False + self._buttonList[b].configure(background="gray", activebackground="gray") + else: # button was off + if sum(self._buttonState) > 2: # only three at a time allowed + self.reportQuadError() + return + self._buttonState[b] = True + self._buttonList[b].configure(background="green", activebackground="green") + + + def NEClick(self): + self.toggleButton("NE") + return + + def SEClick(self): + self.toggleButton("SE") + return + + def SWClick(self): + self.toggleButton("SW") + return + + def NWClick(self): + self.toggleButton("NW") + return + + def makeLabel(self, frame): + + label = Tkinter.Label(frame, fg="red", font=self._boldFont, + text="Max winds will be reduced by\n 20% in selected quadrants") + label.grid(row=0) + + return + + def makeBottomButtons(self, frame): + # Create the Execute/Cancel buttons + self._doneButton = Tkinter.Button(frame, text="Done", + command=self.doneCommand) + self._doneButton.grid(row=3, column=0, padx=20, pady=5, sticky=Tkinter.W) + + self._cancelButton = Tkinter.Button(frame, text="Cancel", + command=self.cancelCommand) + self._cancelButton.grid(row=3, column=2, padx=20, pady=5, sticky=Tkinter.E) + + frame.grid(columnspan=3, sticky=Tkinter.EW) + + return + + def makeQuadButtons(self, frame): + # Create the quadrant buttons + commandList = [self.NWClick, self.SWClick, self.SEClick, self.NEClick] + self._buttonLabels = ["NW", "SW", "SE", "NE"] + + # define the button position in geometric order + buttonPos = [(0, 0), (1, 0), (1, 1), (0, 1)] + for b in range(len(self._buttonLabels)): + label = self._buttonLabels[b] + + self._buttonList[b] = Tkinter.Button(frame, text=label, + command=commandList[b], + font=self._bigFont, width=3) + rowPos, colPos = buttonPos[b] + self._buttonList[b].grid(row=rowPos, column=colPos, padx=30, pady=10) + + return + + def setUpUI(self): + + self._labelFrame = Tkinter.Frame(self._master) + + self._labelFrame.grid(row=0) + + + self._buttonFrame = Tkinter.Frame(self._master, borderwidth=3, + relief=Tkinter.RIDGE, bd=2, pady=5) + self._buttonFrame.grid(row=1, column=0,padx=25, + sticky=Tkinter.E+Tkinter.W, pady=5) + + + self._bottomFrame = Tkinter.Frame(self._master, borderwidth=3, + relief=Tkinter.RIDGE, bd=2) + self._bottomFrame.grid(row=2, column=0, columnspan=2, padx=25) + + self._master.title("Reduce Max Wind by Quadrant") + + self.makeLabel(self._labelFrame) + + self.makeQuadButtons(self._buttonFrame) + + self.makeBottomButtons(self._bottomFrame) + + return + + def doneCommand(self): + self._master.quit() + + quadCount = 4 + reducePct = 0.80 + + # Gather up the maxWind info to return to the main tool + self._maxWindDict = {} + for h in range(len(self._hourList)): + windList = [] + for quad in range(quadCount): + + # Reduce the value if that quadrant was selected + if self._buttonState[quad]: + windValue = self._maxWind[quad][h] * self._allTimeMaxWind * reducePct + else: + windValue = self._maxWind[quad][h] * self._allTimeMaxWind + + windList.append(windValue) + + windList.reverse() + self._maxWindDict[self._hourList[h]] = windList + + return + + def cancelCommand(self): + self._master.destroy() + + return None + + def displayGUI(self, windDict): + + self._windDict = windDict + self._maxWindDict = None + self._quadCount = 4 + + hourKeys = self._windDict.keys() + hourKeys.sort() + self._hourList = hourKeys + self._initialMinWind = [] + self._initialMaxWind = [] + for hour in hourKeys: + minWind, maxWind = windDict[hour] + self._initialMinWind.append(minWind) + self._initialMaxWind.append(maxWind) + + self._hourCount = len(hourKeys) + + if AWIPS_ENVIRON == "AWIPS1": + self._allTimeMaxWind = max(self._initialMaxWind) + else: # numpy uses "amax" + self._allTimeMaxWind = amax(self._initialMaxWind) + + + + # Make the topLevel window - different for A1 and A2 + if AWIPS_ENVIRON == 'AWIPS1': + self._master = Tkinter.Toplevel(self.eaMgr().root()) + self._master.transient(self.eaMgr().root()) # always draw on top of GFE + else: + self._tkmaster = Tkinter.Tk() + self._master = Tkinter.Toplevel(self._tkmaster) + self._tkmaster.withdraw() + + self._buttonLabels = ["NW", "SW", "SE", "NE"] + + self._buttonState = [False, False, False, False] + self._buttonList = [None, None, None, None] + + self._boldFont = tkFont.Font(family="Helvetica", size=12, weight="bold") + self._bigFont = tkFont.Font(family="Helvetica", size=16) + + self.setUpUI() + + self._maxWind = zeros((self._quadCount, self._hourCount)) * 1.0 + + for hour in range(len(hourKeys)): + for quad in range(self._quadCount): + minWind, maxWind = self._windDict[hourKeys[hour]] + self._maxWind[quad][hour] = maxWind / self._allTimeMaxWind + + + #self.updateDisplay() # Draws everything + + + self._master.mainloop() + + self._master.withdraw() + self._master.destroy() + + return self._maxWindDict