1234 lines
51 KiB
Python
1234 lines
51 KiB
Python
##
|
|
# 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
|
|
# any purpose.
|
|
#
|
|
# ISC_Utility - Version 3.05 (Tim Barker - SOO Boise, ID)
|
|
#
|
|
# Supports new routines that mimic NDFD algorithms that were changed
|
|
# in late 2005 or early 2006. Algorithms now have potentially different
|
|
# thresholds for every 'border pair', based on topography, and the values
|
|
# of the grids themselves (i.e. large values of waveheight have more
|
|
# leniant' thresholds). Algorithms now consider a border as 'discrepant'
|
|
# if the average absolute difference along the border is larger than the
|
|
# average threshold along that border. Some tools/algorithms will also show
|
|
# which individual pairs violate their particular threshold.
|
|
#
|
|
# Author: barker
|
|
# 2017-12-13 - dgilling - fix direction checking in _checkGridBorders. Fixes
|
|
# Show_ISC_Info tool.
|
|
# 2008-11-19 - Barker - Version 3.05. Added code to check for 'office type'
|
|
# of editAreas, so that it checks only ISC_xxxx areas for wfos -
|
|
# not rfcs. Also added a check for ISC_xxxx editAreas without
|
|
# corresponding xxxx editArea.
|
|
# 2006-01-23 - Barker - Version 3.04. Add thresholds for PoP12hr, QPF12hr
|
|
# SnowAmt12hr (standard in ER), WindChill, HeatIndex (instead
|
|
# of AppT), and PoP6, PoP12 (common in SR).
|
|
# 2006-01-19 - Barker - Version 3.03. Another typo for non-square grids.
|
|
# 2006-01-17 - Barker - Version 3.02. Fix problem for non-square grids.
|
|
# 2006-01-13 - Barker - Version 3.01. Changed for new NDFD algorithm
|
|
# Thresholds now vary at each gridpoint - overall average
|
|
# difference along border must be less than average threshold
|
|
# along that border (a much better algorithm!).
|
|
#
|
|
# ----------------------------------------------------------------------------
|
|
|
|
##
|
|
# This is a base file that is not intended to be overridden.
|
|
##
|
|
|
|
|
|
|
|
import numpy
|
|
import SmartScript
|
|
import time
|
|
import TimeRange
|
|
from com.raytheon.uf.common.dataplugin.gfe.db.objects import GridParmInfo
|
|
GridType = GridParmInfo.GridType
|
|
|
|
class ISC_Utility(SmartScript.SmartScript):
|
|
def __init__(self, dbss, eaMgr, mdMode=None, toolType="numeric"):
|
|
SmartScript.SmartScript.__init__(self, dbss)
|
|
self.setToolType(toolType)
|
|
self._dbss = dbss
|
|
# self.setUp(eaMgr, mdMode, toolType)
|
|
self.configuration()
|
|
#
|
|
# Always check to see if BorderPairs are current
|
|
#
|
|
refresh=7200 # seconds between refresh of borders
|
|
self._debug=0 # set to 1 or 5 or 10 for increasing info
|
|
self.list=""
|
|
#
|
|
# get Border Pair info - either from cache, or by calculating
|
|
#
|
|
self.pairInfo=self._getCachedPairs("BorderPairs","ISCPairs",refresh)
|
|
if self.pairInfo is None:
|
|
if self._debug>=1:
|
|
self.statusBarMsg("Calculating Pairs","R")
|
|
self.pairInfo=self._getPairInfo(self.getTopo())
|
|
if self._debug>=1:
|
|
self.statusBarMsg("Calculating Pairs Done","R")
|
|
self._cachePairs("BorderPairs",self.pairInfo,"ISCPairs")
|
|
else:
|
|
if self._debug>=1:
|
|
self.statusBarMsg("Pair info obtained from cache","R")
|
|
|
|
|
|
#-------------------------------------------------------------------------
|
|
#
|
|
# C O N F I G U R A T I O N S E C T I O N F O R A L L I S C TOOLS
|
|
#
|
|
# this function is intended to be overridden and sets up the default
|
|
# configuration for the set of ISC tools. Copy and place into the
|
|
# ISC_Utility_Local_New and modify as needed.
|
|
def configuration(self):
|
|
#
|
|
# points which have an elevation difference greater than this will NOT
|
|
# be considered in ISC statistics (in feet). NDFD sets this to 1000ft.
|
|
#
|
|
self.MAXTOPODIFF=1000.0
|
|
#
|
|
# NDFD checks are not performed when one side of a border is a land
|
|
# point and the other side is an ocean point. To do this, an EditArea
|
|
# with land/sea points needs to be calculated. With LANDEDITAREA set
|
|
# to None - the code will calculate the land area by a 'union' of all
|
|
# points found in the CWA editareas named XXX, where the XXX values
|
|
# are taken from all the editareas name ISC_XXX. If you have not
|
|
# overridden the ISC_XXX editarea or XXX edit areas, then this will
|
|
# work fine. If you HAVE overridden these edit area - use the
|
|
# LANDEDITAREA to specify the name of an editarea that contains just
|
|
# land points (all others are assumed to be ocean points).
|
|
#
|
|
self.LANDEDITAREA=None # or string with name of EditArea containing land
|
|
#
|
|
#--------------------------------------------------------------------
|
|
# These configuration items for Show_ISC_Info and Show_ISC_Highlights.
|
|
#
|
|
# If you want the check for a particular parm to ACTUALLY check other
|
|
# parms, then list them here. Vector parms need not be listed - but
|
|
# the threshold for Vector parms in GFE is assumed to be the threshold
|
|
# for the magnitude part - and the threshold for the direction part is
|
|
# hard-coded below
|
|
#
|
|
self.MultiParms={"MaxRH":("MinT","TdMrn","MaxRH"),
|
|
"MinRH":("MaxT","TdAft","MinRH"),
|
|
"RH": ("T","Td","RH"),
|
|
}
|
|
|
|
# Minimum number of points along a border before it considers a
|
|
# failing average threshold "significant" (to get rid of short borders)
|
|
#
|
|
self.MINPOINTS=10
|
|
#
|
|
#------------------------------------------------------------------
|
|
#
|
|
# NDFD thresholds - should not need to be modified.
|
|
#
|
|
# Each entry in THRESHOLDS contains a tuple ( parmnames, thresholdinfo)
|
|
# parmnames can be a tuple with many parms listed that use the
|
|
# same threshold
|
|
# thresholdinfo contains (thresholdtype,thresholdvalues,
|
|
# conditions,dirflag) where:
|
|
# thresholdtype="contant","topo" or "graduated"
|
|
# thresholdvalues=
|
|
# for "constant" type: value
|
|
#
|
|
# differences greater than value are considered
|
|
# discrepant
|
|
#
|
|
# for "topo" type: (elev,lowvalue,highvalue)
|
|
#
|
|
# if the elevation difference between points is
|
|
# less than elev, then the lowvalue is used as
|
|
# the threshold value. Otherwise the highvalue
|
|
# is used for the threshold value
|
|
#
|
|
# for "graduated" type: (bigvalue,(lessthan,value),(lessthan,value),...)
|
|
#
|
|
# bigvalue is the default threshold value. However
|
|
# if the lowest of the two pair points is less than the
|
|
# 'lessthan', then that 'value' is used for the
|
|
# threshold instead. All 'lessthan' values are checked,
|
|
# so they should be listed in decreasing order.
|
|
#
|
|
self.DEFAULT_THRESHOLD=("constant",5,("none",0,0),0)
|
|
self.THRESHOLDS=[
|
|
(("T","Td","MaxT","MinT","TdAft","TdMrn"),
|
|
("topo",(500,5,7),("none",0,0),0)),
|
|
(("HeatIndex","WindChill"),
|
|
("topo",(500,7,9),("none",0,0),0)),
|
|
(("PoP","PoP12","PoP6","PoP12hr","PoP6hr"),
|
|
("constant",20,("none",0,0),0)),
|
|
(("WindSpd","TransWindSpd","WindGust"),
|
|
("graduated",(15,(20,10)),("greater_equal",12,0),0)),
|
|
(("WindDirec","TransWindDirec"),
|
|
("topo",(500,45,90),("greater_equal",12,1),1)),
|
|
(("Sky"),
|
|
("topo",(500,25,35),("none",0,0),0)),
|
|
(("QPF","QPF6hr"),
|
|
("graduated",(1.0,(3.0,0.5),(1.5,0.25)),("greater",0.25,0),0)),
|
|
(("SnowAmt","SnowAmt6hr"),
|
|
("graduated",(6,(12,4),(6,2)),("greater",2,0),0)),
|
|
(("SnowLevel","FzLevel","MixHgt"),
|
|
("constant",1000,("none",0,0),0)),
|
|
(("RH","MaxRH","MinRH"),
|
|
("graduated",(25,(75,20),(50,15),(25,10)),("none",0,0),0)),
|
|
(("WaveHeight"),
|
|
("graduated",(10,(36,9),(32,8),(28,7),(24,6),(20,5),(16,4),(12,3),(6,2)),("greater",0,1),0)),
|
|
(("CWR"),
|
|
("constant",10,("none",0,0),0)),
|
|
(("Haines"),
|
|
("constant",1,("none",0,0),0)),
|
|
]
|
|
# --------- E N D C O N F I G U R A T I O N S E C T I O N ----------
|
|
|
|
|
|
|
|
#=================================================================
|
|
# _getPairInfo - calculate pairInfo data structure from edit areas
|
|
#
|
|
# each entry in pairInfo is a tuple with (label,pairlist)
|
|
# where label=string with CWA neighbor name
|
|
# pairlist=list of tuples
|
|
# where each entry in pairlist is a tuple with: (insidex,insidey,outsidex,outsidey,TopoDiff)
|
|
# where insidex,insidey = pair coordinates inside CWA
|
|
# outsidex,outsidey = pair coordinates outside CWA
|
|
# absTopoDiff = topography difference (always positive)
|
|
#
|
|
def _getPairInfo(self,Topo):
|
|
pairInfo=[]
|
|
#
|
|
# Find coordinates of neighboring pairs
|
|
# first setup eah=home edit area
|
|
#
|
|
homeSite=self.getSiteID()
|
|
name="ISC_"+homeSite
|
|
eah=self.encodeEditArea(self.getEditArea(name))
|
|
homeType=self.myOfficeType()
|
|
#
|
|
# Get grids with home edit area shifted 1 pixel
|
|
# each direction...and topo difference when shifted
|
|
# 1 pixel in each direction
|
|
#
|
|
xshift=( 1,-1, 0, 0)
|
|
yshift=( 0, 0, 1,-1)
|
|
homeshifts=[]
|
|
topodiffs=[]
|
|
for i in range(4):
|
|
homeshifts.append(self.offset(eah,xshift[i],yshift[i]))
|
|
topodiffs.append(Topo-self.offset(Topo,xshift[i],yshift[i]))
|
|
landGrid=self._getLandEditArea()
|
|
#
|
|
# Loop through other ISC_xxx edit areas - except the one
|
|
# for the home edit area
|
|
#
|
|
eanames=self.editAreaList()
|
|
|
|
for eaname in eanames:
|
|
if (len(eaname)==7)and(eaname[0:4]=="ISC_")and(eaname[4:]!=homeSite):
|
|
siteName=eaname[4:]
|
|
sType=self.officeType(siteName)
|
|
if sType is None:
|
|
continue
|
|
if sType==homeType:
|
|
ean=self.encodeEditArea(self.getEditArea(eaname))
|
|
#
|
|
# Compare edit area to the shifted home edit areas
|
|
# looking for common points. Also check that topo
|
|
# difference is less than MAXTOPODIFF
|
|
#
|
|
pairs=[]
|
|
for k in range(4):
|
|
pair=numpy.logical_and(homeshifts[k],ean)
|
|
topodiff=topodiffs[k]
|
|
|
|
for x in range(eah.shape[1]):
|
|
for y in range(eah.shape[0]):
|
|
if pair[y,x]==1: # common point
|
|
tdiff=numpy.abs(topodiff[y,x])
|
|
homepointx=x-xshift[k]
|
|
homepointy=y-yshift[k]
|
|
land1=landGrid[homepointy,homepointx]
|
|
land2=landGrid[y,x]
|
|
coast=numpy.not_equal(land1,land2) # 0 if both land or both water
|
|
pairs.append((homepointx,homepointy,x,y,tdiff,coast))
|
|
|
|
if len(pairs)>0:
|
|
label=siteName
|
|
pairInfo.append((label,pairs))
|
|
|
|
|
|
return pairInfo
|
|
#=================================================================
|
|
# _getCachedPairs - See if cached pairInfo structure is still
|
|
# valid. If so, return it, otherwise return
|
|
# None.
|
|
#
|
|
def _getCachedPairs(self, name, category, timeLimit):
|
|
try:
|
|
cachedObject = self.getObject(name+"_"+self._dbss.getSiteID(), category)
|
|
pairInfo, timeWritten, geoInfo = cachedObject
|
|
if timeLimit != 0 and time.time() - timeWritten > timeLimit:
|
|
return None #too old
|
|
|
|
# validate geoinfo
|
|
geo = self._dbss.getSiteID()
|
|
if not geoInfo == geo:
|
|
return None #different geoinfo
|
|
return pairInfo
|
|
except:
|
|
return None
|
|
#=================================================================
|
|
# _cachePairs - save pairInfo structure data, along with time and
|
|
# grid location info to IFPS database - so it can be
|
|
# grabbed quickly later without re-calculating
|
|
#
|
|
def _cachePairs(self, name, pairInfo, category):
|
|
objectToSave = pairInfo, time.time(), self._dbss.getSiteID()
|
|
self.saveObject(name+"_"+self._dbss.getSiteID(), objectToSave, category)
|
|
#=================================================================
|
|
# _getLandEditArea - get Land editArea grid, calculating from
|
|
# cwa edit areas if not specified via
|
|
# configuration
|
|
#
|
|
# Returns grid of 0/1 for land. Everything not land is assumed
|
|
# to be water.
|
|
#
|
|
def _getLandEditArea(self):
|
|
#
|
|
# Get points that are land - either by the specified edit
|
|
# area in configuration, or by looking for ISC_xxx edit
|
|
# areas, and adding up all the corresponding xxx areas.
|
|
#
|
|
# This does NOT work if the default CWA edit areas named 'xxx'
|
|
# have been overridden with other names or different areas, or
|
|
# if new ISC_xxx edit areas have been added to the system and
|
|
# have a corresponding xxx edit area.
|
|
#
|
|
landGrid=None
|
|
if (self.LANDEDITAREA is not None):
|
|
landea=self.getEditArea(self.LANDEDITAREA)
|
|
if (landea is not None):
|
|
landGrid=self.encodeEditArea(landea)
|
|
if landGrid is None:
|
|
landGrid = self.empty(numpy.bool)
|
|
eanames=self.editAreaList()
|
|
for eaname in eanames:
|
|
if ((len(eaname)==7)and(eaname[0:4]=="ISC_")):
|
|
name=eaname[4:]
|
|
if name in eanames:
|
|
ea=self.getEditArea(name)
|
|
if ea is not None:
|
|
grid=self.encodeEditArea(ea)
|
|
landGrid |= grid
|
|
return landGrid
|
|
#========================================================================
|
|
# _getThresholdInfo - return thresholdInfo structure for the
|
|
# specified parm.
|
|
#
|
|
def _getThresholdInfo(self,parmName):
|
|
thresholdInfo=self.DEFAULT_THRESHOLD
|
|
for (names,threshold) in self.THRESHOLDS:
|
|
if parmName in names:
|
|
thresholdInfo=threshold
|
|
break
|
|
return thresholdInfo
|
|
#=================================================================
|
|
# _getListing - get text detailing most recent checks
|
|
#
|
|
def _getListing(self):
|
|
return self.list
|
|
#=================================================================
|
|
# _checkParmBorders - check ISC borders for this parm - which may
|
|
# mean checking more than one grid (see the
|
|
# MultiParms configuration section)
|
|
#
|
|
def _checkParmBorders(self,WEname,GridTimeRange,listing=0):
|
|
timetext=self.makeTimeMsg(GridTimeRange)
|
|
if listing==1:
|
|
self.list="ISC Discrepancies Check for %s %s:\n\n"%(WEname,timetext)
|
|
#
|
|
# Get list of parms to check - which might be more than one
|
|
#
|
|
if WEname in self.MultiParms:
|
|
parmlist=self.MultiParms[WEname]
|
|
else:
|
|
parmlist=(WEname,)
|
|
#
|
|
# Loop over each parm - put status messages for reading problems
|
|
#
|
|
totalviolate=0
|
|
totalwarning=0
|
|
totalchecked=0
|
|
for parmname in parmlist:
|
|
results=self._checkGridBorders(parmname, GridTimeRange, listing=listing)
|
|
(status,numchecked,violate,warning)=results
|
|
if status==1:
|
|
msg="No %s GridInfo for %s"%(WEname,timetext)
|
|
self.statusBarMsg(msg,"A")
|
|
continue
|
|
if status==2:
|
|
msg="No %s ISC data for %s"%(WEname,timetext)
|
|
self.statusBarMsg(msg,"A")
|
|
continue
|
|
if status==3:
|
|
msg="No border checks for weather or discrete elements: %s"%WEname
|
|
self.statusBarMsg(msg,"A")
|
|
continue
|
|
totalviolate=totalviolate+violate
|
|
totalwarning=totalwarning+warning
|
|
totalchecked=totalchecked+numchecked
|
|
return(totalchecked,totalviolate,totalwarning)
|
|
#=================================================================
|
|
# _checkGridBorders - check Borders for a single grid (which
|
|
# might be a vector - so it might actually
|
|
# check two grids)
|
|
#
|
|
# returns (status,numchecked,violate,warning)
|
|
#
|
|
# if status is non-zero, then there was a problem reading the
|
|
# grid and the values are meaningless). numchecked is the total
|
|
# number of pairs checked - if zero - then no checks were done
|
|
# because no points met the conditions. That is different than
|
|
# no violations/warning when lots of points were checked.
|
|
#
|
|
def _checkGridBorders(self,parmName,GridTimeRange,listing=0):
|
|
#
|
|
# Get grid info - return status=1 if no gridInfo
|
|
#
|
|
try:
|
|
gridInfoList=self.getGridInfo(self.mutableID(), parmName, "SFC",
|
|
GridTimeRange)
|
|
except:
|
|
return(1,0,0,0)
|
|
#
|
|
# Return status=1 if gridInfo is empty
|
|
#
|
|
if (len(gridInfoList)<1):
|
|
return (1,0,0,0)
|
|
gridInfo=gridInfoList[0]
|
|
#
|
|
# get the ISC data - return status=2 if could not read ISC data
|
|
#
|
|
bits,isc=self._getBitsAndISC(parmName, gridInfo, GridTimeRange)
|
|
if ((bits is None) or (isc is None)):
|
|
return (2,0,0,0)
|
|
#
|
|
# Only know how to check for Scalar or Vector grids
|
|
#
|
|
WEtype=gridInfo.type()
|
|
if (not GridType.SCALAR.equals(WEtype)) and (not GridType.VECTOR.equals(WEtype)):
|
|
return (3,0,0,0)
|
|
|
|
#
|
|
# If a vector - then check each part seperately - but the
|
|
# condition always depends on speed. For scalar grids the
|
|
# condition is always the same as the grid itself.
|
|
#
|
|
gridList=[]
|
|
if GridType.SCALAR.equals(WEtype):
|
|
mag=isc
|
|
gridList.append((parmName,mag,mag))
|
|
else:
|
|
name1=parmName+"Spd"
|
|
(mag,direc)=isc
|
|
gridList.append((name1,mag,mag))
|
|
name2=parmName+"Direc"
|
|
gridList.append((name2,direc,mag))
|
|
#
|
|
# Loop over grids (usually one, but possibly two for vectors)
|
|
#
|
|
totalchecked=0
|
|
totalviolate=0
|
|
totalwarning=0
|
|
for (pname,grid,condgrid) in gridList:
|
|
if self._debug>=10:
|
|
print("checking {!s} grid".format(pname))
|
|
if listing==1:
|
|
self.list=self.list+"For %s:\n"%pname
|
|
(numchecked,violate,warning)=self._checkAllBorders(pname,grid,condgrid,bits,listing=listing)
|
|
totalchecked=totalchecked+numchecked
|
|
totalviolate=totalviolate+violate
|
|
totalwarning=totalwarning+warning
|
|
if listing==1:
|
|
self.list=self.list+"\n"
|
|
return(0,totalchecked,totalviolate,totalwarning)
|
|
#=================================================================
|
|
# _checkAllBorders - check all borders for a single scalar grid
|
|
# that has a condition grid of congrid and bits
|
|
# indicates where ISC data was available.
|
|
#
|
|
# if listing=1, then adds text to self.list
|
|
# that shows some stats for each border
|
|
#
|
|
# returns:
|
|
# numchecked - total number of pairs checked
|
|
# over all borders
|
|
# numviolations - total number of borders that
|
|
# were in violation
|
|
# numwarnings - total number of borders that
|
|
# were in warning (violating but
|
|
# having less than MINPOINTS
|
|
# pairs)
|
|
#
|
|
def _checkAllBorders(self,parmName,grid,condgrid,bits,listing=0):
|
|
#
|
|
# get Threshold info for this parm
|
|
#
|
|
thresholdInfo=self._getThresholdInfo(parmName)
|
|
if self._debug>=10:
|
|
print("thresholdInfo=", thresholdInfo)
|
|
(thresholdType,thresholdValues,conditions,dirflag)=thresholdInfo
|
|
#
|
|
# Loop over each neighbors border
|
|
#
|
|
violate=0
|
|
warning=0
|
|
totalchecks=0
|
|
for (label,pairList) in self.pairInfo:
|
|
if self._debug>=5:
|
|
print("Checking borders with {!s}".format(label))
|
|
results=self._checkOneBorder(grid,condgrid,bits,pairList,thresholdInfo)
|
|
(returnvalue,totalnum,numchecked,numviolate,
|
|
avgbias,avgdiff,avgthresh)=results
|
|
if self._debug>=5:
|
|
print(" totalpoints :", totalnum)
|
|
print(" numberchecked:", numchecked)
|
|
print(" numviolate :", numviolate)
|
|
if (numchecked>0):
|
|
print(" bias :", avgbias)
|
|
print(" diff :", avgdiff)
|
|
print(" threshold :", avgthresh)
|
|
if listing==1:
|
|
self.list=self.list+" Avg Diff for %s is %7.2f (limit %7.2f) [%4d pairs - %4d failed] - "%(label,avgdiff,avgthresh,numchecked,numviolate)
|
|
if avgdiff>avgthresh:
|
|
if numchecked<self.MINPOINTS:
|
|
self.list=self.list+"IGNORED"
|
|
else:
|
|
self.list=self.list+"FAILED"
|
|
else:
|
|
self.list=self.list+"OK"
|
|
self.list=self.list+"\n"
|
|
totalchecks=totalchecks+numchecked
|
|
if returnvalue==1:
|
|
warning=warning+1
|
|
if returnvalue==2:
|
|
violate=violate+1
|
|
return (numchecked,violate,warning)
|
|
|
|
#=================================================================
|
|
# _meetCondiditions - returns 1 if conditions for checking against
|
|
# the threshold are met - or 0 if the conditions
|
|
# are not met.
|
|
#
|
|
def _meetConditions(self, conditions, x1, y1, x2, y2, condgrid):
|
|
(conditionType, conditionValue, bothFlag) = conditions
|
|
if "none" == conditionType:
|
|
return 1
|
|
value1 = condgrid[y1,x1]
|
|
value2 = condgrid[y2,x2]
|
|
if "greater" == conditionType:
|
|
cond1 = (value1>conditionValue)
|
|
cond2 = (value2>conditionValue)
|
|
elif "greater_equal" == conditionType:
|
|
cond1 = (value1>=conditionValue)
|
|
cond2 = (value2>=conditionValue)
|
|
elif "less" == conditionType:
|
|
cond1 = (value1<conditionValue)
|
|
cond2 = (value2<conditionValue)
|
|
elif "less_equal" == conditionType:
|
|
cond1 = (value1<=conditionValue)
|
|
cond2 = (value2<=conditionValue)
|
|
if bothFlag:
|
|
return (cond1 and cond2)
|
|
else:
|
|
return (cond1 or cond2)
|
|
|
|
#=================================================================
|
|
# _getDiffGrid - check Borders for a single grid - and return
|
|
# a grid with the max difference between the
|
|
# point and any neighboring points. If vectordir
|
|
# =1, and grid is a vector - then make grid for
|
|
# direction - otherwise get grid for magnitude part.
|
|
# If MaskNonViolators=1 then do not set the grid
|
|
# value unless it violates the threshold for these
|
|
# points.
|
|
#
|
|
# returns grid of non-zero values, or None if there were problems
|
|
# reading the grids.
|
|
#
|
|
#
|
|
def _getDiffGrid(self,parmName,GridTimeRange,vectordir=0,maskNonViolators=0):
|
|
#
|
|
# Get grid info - return None if no gridInfo
|
|
#
|
|
diffGrid=self.empty()
|
|
|
|
try:
|
|
gridInfoList=self.getGridInfo(self.mutableID(), parmName, "SFC",
|
|
GridTimeRange)
|
|
except:
|
|
return None
|
|
#
|
|
# Return None if an empty list of gridInfos
|
|
#
|
|
if (len(gridInfoList)<1):
|
|
return None
|
|
gridInfo=gridInfoList[0]
|
|
#
|
|
# get the ISC data - return None if it could not read data
|
|
#
|
|
bits,isc=self._getBitsAndISC(parmName, gridInfo, GridTimeRange)
|
|
if ((bits is None) or (isc is None)):
|
|
return None
|
|
#
|
|
# Only know how to check for Scalar or Vector grids
|
|
#
|
|
WEtype=gridInfo.type()
|
|
if (not GridType.SCALAR.equals(WEtype)) and (not GridType.VECTOR.equals(WEtype)):
|
|
return None
|
|
#
|
|
# If a vector - then check only magnitude part
|
|
#
|
|
if GridType.SCALAR.equals(WEtype):
|
|
grid=isc
|
|
condgrid=isc
|
|
checkName=parmName
|
|
else:
|
|
(mag,direc)=isc
|
|
if vectordir==0:
|
|
checkName=parmName+"Spd"
|
|
grid=mag
|
|
condgrid=mag
|
|
else:
|
|
checkName=parmName+"Direc"
|
|
grid=direc
|
|
condgrid=mag
|
|
#
|
|
# get Threshold info for this parm
|
|
#
|
|
thresholdInfo=self._getThresholdInfo(checkName)
|
|
if self._debug>=10:
|
|
print("thresholdInfo=", thresholdInfo)
|
|
(thresholdType,thresholdValues,conditions,dirflag)=thresholdInfo
|
|
#
|
|
# Loop over each neighbors border
|
|
#
|
|
for (label,pairList) in self.pairInfo:
|
|
if self._debug>=5:
|
|
print("Checking borders with {!s}", label)
|
|
for pair in pairList:
|
|
(x1,y1,x2,y2,topodiff,coast)=pair
|
|
if self._debug>=10:
|
|
print(" point {:3d},{:3d}-->{:3d},{:3d}".format(x1, y1, x2, y2))
|
|
#
|
|
# get values across the border
|
|
#
|
|
value1=grid[y1,x1]
|
|
value2=grid[y2,x2]
|
|
#
|
|
# Get the difference across the border
|
|
#
|
|
diff=value1-value2
|
|
if dirflag==1:
|
|
if diff>180.0:
|
|
diff=diff-360.0
|
|
if diff<-180.0:
|
|
diff=diff+360.0
|
|
absdiff=abs(diff)
|
|
#
|
|
# Get the threshold (which might depend on the values)
|
|
#
|
|
if thresholdType=="constant":
|
|
thresh=thresholdValues
|
|
elif thresholdType=="topo":
|
|
(elevation,low,high)=thresholdValues
|
|
if (topodiff<=elevation):
|
|
thresh=low
|
|
else:
|
|
thresh=high
|
|
elif thresholdType=="graduated":
|
|
minvalue=min(value1,value2)
|
|
thresh=thresholdValues[0]
|
|
for i in range(1,len(thresholdValues)):
|
|
(lessthan,newthresh)=thresholdValues[i]
|
|
if minvalue<lessthan:
|
|
thresh=newthresh
|
|
#
|
|
# Calculate whether this particular point violates the threshold
|
|
#
|
|
if absdiff<=thresh:
|
|
itPasses=1
|
|
else:
|
|
itPasses=0
|
|
if self._debug>=10:
|
|
print(" {:f} {:f} diff:{:f} threshold:{:f} itPasses:{:1d}".format(value1, value2, diff, thresh, itPasses))
|
|
#
|
|
# If too big a topodiff - or a coast - then it passes anyway
|
|
#
|
|
if ((topodiff>self.MAXTOPODIFF) or (coast==1)):
|
|
itPasses=1
|
|
#
|
|
# If no ISC data - then it passes anyway
|
|
# available in ISC
|
|
#
|
|
if (bits[y2,x2]<0.5):
|
|
itPasses=1
|
|
#
|
|
# Make sure conditions for checking are met
|
|
#
|
|
if (not self._meetConditions(conditions,x1,y1,x2,y2,condgrid)):
|
|
itPasses=1
|
|
#
|
|
# Do nothing if this point passes and configuration
|
|
# indicates that nonViolators will not be displayed
|
|
#
|
|
if (maskNonViolators and itPasses):
|
|
continue
|
|
#
|
|
if absdiff>abs(diffGrid[y1,x1]):
|
|
diffGrid[y1,x1]=diff
|
|
return diffGrid
|
|
#=================================================================
|
|
#
|
|
# checkOneBorder - given a list of points in pairList, the threshold
|
|
# info, the grid and condition grid, and the bits
|
|
# that indicates where ISC data is available:
|
|
# then for each point:
|
|
# get the grid values,
|
|
# check that conditions are met,
|
|
# get difference,
|
|
# get threshold,
|
|
# check difference comapred to threshold and
|
|
# add to total differences, biases, etc.
|
|
#
|
|
# return:
|
|
# code = 0 if border passed (even if no points)
|
|
# 1 if border violated but had MINPOINTS
|
|
# or fewer pairs
|
|
# 2 if border violated and had more than
|
|
# MINPOINTS pairs
|
|
# totalnum = total number of pairs along the
|
|
# border
|
|
# numchecked = total pairs checked along the
|
|
# border after skipping ones without
|
|
# ISC data (bits=0), topo difference
|
|
# higher than limit, or condgrid
|
|
# values not meeting condition
|
|
# numviolate = number of individual pairs that
|
|
# violated their threshold
|
|
# avgbias = average difference among pairs
|
|
# avgdiff = average absolute difference among pairs
|
|
# avgthresh = average threshold along border
|
|
#
|
|
#
|
|
def _checkOneBorder(self,grid,condgrid,bits,pairList,thresholdInfo):
|
|
(thresholdType,thresholdValues,conditions,dirflag)=thresholdInfo
|
|
totalnum=len(pairList)
|
|
numchecked=0
|
|
numviolate=0
|
|
biastotal=0.0
|
|
difftotal=0.0
|
|
threshtotal=0.0
|
|
|
|
for pair in pairList:
|
|
(x1,y1,x2,y2,topodiff,coast)=pair
|
|
if self._debug>=10:
|
|
print(" point {:3d},{:3d}-->{:3d},{:3d}".format(x1, y1, x2, y2))
|
|
#
|
|
# no tests if too big of topodiff, or if a coastline pair
|
|
#
|
|
if ((topodiff>self.MAXTOPODIFF) or (coast==1)):
|
|
if self._debug>=10:
|
|
print(" skipped because of elevation or coastline")
|
|
continue
|
|
#
|
|
# no tests if bits of outside point indicate it is NOT
|
|
# available in ISC
|
|
#
|
|
if (bits[y2,x2]<0.5):
|
|
if self._debug>=10:
|
|
print(" skipped because ISC data not available")
|
|
continue
|
|
#
|
|
# get values across the border
|
|
#
|
|
value1=grid[y1,x1]
|
|
value2=grid[y2,x2]
|
|
#
|
|
# Make sure conditions for checking are met
|
|
#
|
|
if (not self._meetConditions(conditions,x1,y1,x2,y2,condgrid)):
|
|
if self._debug>=10:
|
|
print(" skipped because {:f} and {:f} did not meet the conditions:{!s}".format(value1, value2, conditions))
|
|
continue
|
|
numchecked=numchecked+1
|
|
#
|
|
# Get the difference across the border - add it to the
|
|
# total difference (biastotal) and total of absolute
|
|
# value of differences (difftotal)
|
|
#
|
|
diff=value1-value2
|
|
if dirflag==1:
|
|
if diff>180.0:
|
|
diff=diff-360.0
|
|
if diff<-180.0:
|
|
diff=diff+360.0
|
|
absdiff=abs(diff)
|
|
biastotal=biastotal+diff
|
|
difftotal=difftotal+absdiff
|
|
#
|
|
# Get the threshold (which might depend on the values)
|
|
# Add it to the total of thresholds
|
|
#
|
|
if thresholdType=="constant":
|
|
thresh=thresholdValues
|
|
elif thresholdType=="topo":
|
|
(elevation,low,high)=thresholdValues
|
|
if (topodiff<=elevation):
|
|
thresh=low
|
|
else:
|
|
thresh=high
|
|
elif thresholdType=="graduated":
|
|
minvalue=min(value1,value2)
|
|
thresh=thresholdValues[0]
|
|
for i in range(1,len(thresholdValues)):
|
|
(lessthan,newthresh)=thresholdValues[i]
|
|
if minvalue<lessthan:
|
|
thresh=newthresh
|
|
threshtotal=threshtotal+float(thresh)
|
|
#
|
|
# Calculate whether this particular point violates the threshold
|
|
#
|
|
if absdiff<=thresh:
|
|
itPasses=1
|
|
else:
|
|
itPasses=0
|
|
numviolate=numviolate+1
|
|
if self._debug>=10:
|
|
print(" {:f} {:f} diff:{:f} threshold:{:f} itPasses:{:1d}".format(value1, value2, diff, thresh, itPasses))
|
|
#
|
|
# Calculate average bias, average abs difference, average threshold
|
|
#
|
|
if numchecked>0:
|
|
avgbias=biastotal/numchecked
|
|
avgdiff=difftotal/numchecked
|
|
avgthresh=threshtotal/numchecked
|
|
else:
|
|
avgbias=0.0
|
|
avgdiff=0.0
|
|
avgthresh=0.0
|
|
#
|
|
# Setup return value (0 if OK, 1 if failed but only a few points, 2 if failed)
|
|
#
|
|
returnvalue=0
|
|
if ((avgdiff>avgthresh) and (avgthresh>0.0)):
|
|
if numchecked>self.MINPOINTS:
|
|
returnvalue=2
|
|
else:
|
|
returnvalue=1
|
|
return(returnvalue,totalnum,numchecked,numviolate,avgbias,avgdiff,avgthresh)
|
|
#=========================================================================
|
|
#
|
|
# _getBitsAndISC - a routine to get the ISC composite - broken out
|
|
# from the neighboringPoints routine so that it need not be called
|
|
# several times when looping over different edit areas (for each
|
|
# of the neighboring CWAs).
|
|
#
|
|
def _getBitsAndISC(self,WEname,GridInfo,GridTimeRange):
|
|
isc = self._getBetterComposite(WEname, GridTimeRange)
|
|
if isc is None:
|
|
return None, None
|
|
#
|
|
# See if we are working with a Scalar or Vector element
|
|
#
|
|
wxType = GridInfo.type()
|
|
if GridType.SCALAR.equals(wxType):
|
|
bits, isc = isc
|
|
return bits,isc
|
|
elif GridType.VECTOR.equals(wxType):
|
|
bits, isc, direc = isc
|
|
return bits,(isc,direc)
|
|
else:
|
|
return None, None
|
|
#========================================================================
|
|
#
|
|
# Essentially the same as the SmartScript getComposite routine
|
|
# but correctly handles multiple ISC grids within the timeRange
|
|
# of the grid you want. Can return None if no ISC or specified
|
|
# grids lie within the TimeRange specified.
|
|
#
|
|
# 2005-01-24 - Changed again because accumulative parms return
|
|
# different values from getComposite after IFPS16.
|
|
#
|
|
def _getBetterComposite(self,parmName, timeRange):
|
|
#
|
|
# Get the type, rateParm flag, and limits
|
|
# for the parameter name passed in.
|
|
#
|
|
mutableID=self.mutableID()
|
|
baseGrid=self.getGrids(mutableID,parmName,"SFC",timeRange,noDataError=0)
|
|
if baseGrid is None:
|
|
return None
|
|
gridInfoList=self.getGridInfo(mutableID,parmName,"SFC",timeRange)
|
|
if (len(gridInfoList)<1):
|
|
return None
|
|
for gridInfo in gridInfoList:
|
|
wxType=gridInfo.type()
|
|
minlimit=gridInfo.minLimit()
|
|
#
|
|
# Make sure ISC grids exist for this parm
|
|
#
|
|
parm=self.getParm("ISC",parmName,"SFC")
|
|
if parm is None:
|
|
return None
|
|
#
|
|
# Get list of all ISC time-blocks that fit in the
|
|
# timerange of the specified GridTimeRange grid
|
|
#
|
|
iscInfos=self.getGridInfo("ISC",parmName,"SFC",timeRange)
|
|
if (len(iscInfos)<1):
|
|
return None
|
|
alltrs=[]
|
|
for info in iscInfos:
|
|
tr=info.gridTime()
|
|
alltrs.append(tr)
|
|
#
|
|
# setup sumGrid/counter for average
|
|
#
|
|
if ((parmName=="MaxT")or(parmName=="PoP")):
|
|
sumGrid=self.newGrid(-150.0)
|
|
elif (parmName=="MinT"):
|
|
sumGrid=self.newGrid(150.0)
|
|
else:
|
|
sumGrid=self.empty()
|
|
if GridType.VECTOR.equals(wxType):
|
|
sumv=self.empty()
|
|
cnt = self.empty()
|
|
#
|
|
# foreach time range...get the ISC composite for
|
|
# that hour
|
|
#
|
|
for tr in alltrs:
|
|
comp=self.getComposite(parmName,tr)
|
|
if comp[0].shape!=sumGrid.shape:
|
|
continue
|
|
#
|
|
# Add to sums, or min/max
|
|
#
|
|
if GridType.SCALAR.equals(wxType):
|
|
bits,isc=comp
|
|
#isc=self.getGrids("ISC",parmName,"SFC",tr)
|
|
if parmName in ["MaxT", "PoP"]:
|
|
sumGrid[bits] = numpy.maximum(isc,sumGrid)[bits]
|
|
cnt[bits] = 1
|
|
elif parmName=="MinT":
|
|
sumGrid[bits] = numpy.minimum(isc,sumGrid)[bits]
|
|
cnt[bits] = 1
|
|
else:
|
|
sumGrid[bits] += isc[bits]
|
|
cnt[bits] += 1
|
|
if GridType.VECTOR.equals(wxType):
|
|
bits,mag,direc = comp
|
|
#(mag,direc)=self.getGrids("ISC",parmName,"SFC",tr)
|
|
(u,v)=self.MagDirToUV(mag,direc)
|
|
sumGrid[bits] += u[bits]
|
|
sumv[bits] += v[bits]
|
|
cnt[bits] += 1
|
|
if GridType.WEATHER.equals(wxType):
|
|
bits = comp
|
|
bits,keys,strings=comp
|
|
#(keys,strings)=self.getGrids("ISC",parmName,"SFC",tr)
|
|
#
|
|
# now calculate average/max/min, etc.
|
|
# (count is always 1 for max/min)
|
|
#
|
|
noISC=numpy.less(cnt,0.5)
|
|
bits=numpy.greater(cnt,0.5)
|
|
if GridType.SCALAR.equals(wxType) or GridType.VECTOR.equals(wxType):
|
|
cnt[numpy.less(cnt, 1)] = 1
|
|
if GridType.VECTOR.equals(wxType):
|
|
sumGrid /= cnt
|
|
sumGrid[noISC]= minlimit
|
|
|
|
sumv /= cnt
|
|
sumv[noISC] = minlimit
|
|
|
|
(mag,direc)=self.UVToMagDir(sumGrid,sumv)
|
|
(baseMag,baseDir)=baseGrid
|
|
mag[noISC] = baseMag[noISC]
|
|
direc[noISC] = baseDir[noISC]
|
|
return bits,mag,direc
|
|
else:
|
|
sumGrid /= cnt
|
|
sumGrid[noISC] = baseGrid[noISC]
|
|
return bits,sumGrid
|
|
else:
|
|
return bits,keys,strings
|
|
#=================================================================
|
|
# makeTimeMsg - Make short string with time of this grid, usually
|
|
# something like "Mon (6/5) 12-18Z:", but can get
|
|
# complicated like "Mon (6/5) 18Z - Tue (6/6) 03Z:"
|
|
#
|
|
def makeTimeMsg(self,tr):
|
|
DAYS=("Monday","Tuesday","Wednesday","Thursday",
|
|
"Friday","Saturday","Sunday")
|
|
swdy=DAYS[time.gmtime(tr.startTime().unixTime())[6]]
|
|
sday=tr.startTime().day
|
|
smon=tr.startTime().month
|
|
sdate="%s (%d/%d)"%(swdy,smon,sday)
|
|
shou=tr.startTime().hour
|
|
ewdy=DAYS[time.gmtime(tr.endTime().unixTime())[6]]
|
|
eday=tr.endTime().day
|
|
emon=tr.endTime().month
|
|
edate="%s (%d/%d)"%(ewdy,emon,eday)
|
|
ehou=tr.endTime().hour
|
|
if (sdate==edate):
|
|
msg="%s %2.2d-%2.2dZ:"%(sdate,shou,ehou)
|
|
else:
|
|
msg="%s %2.2dZ - %s %2.2dZ:"%(sdate,shou,edate,ehou)
|
|
return msg
|
|
#================================================================
|
|
# _getElementList - get sorted list of currently displayed mutable
|
|
# model elements (by default it excludes nonSCALAR and non-
|
|
# VECTOR elements - but you can include them with the flag)
|
|
#
|
|
def _getElementList(self,excludeWxDiscrete=1):
|
|
mutableModel=self.mutableID().modelName()
|
|
parmList = self._dbss.getParmManager().getDisplayedParms()
|
|
elementList = []
|
|
for parm in parmList:
|
|
model = parm.getParmID().getDbId().getModelName()
|
|
if model == mutableModel:
|
|
wxType = parm.getGridInfo().getGridType()
|
|
if ((excludeWxDiscrete == 1) and (GridType.SCALAR.equals(wxType) or \
|
|
GridType.VECTOR.equals(wxType))):
|
|
elementList.append(parm.expressionName())
|
|
elementList.sort()
|
|
return elementList
|
|
#================================================================
|
|
# _getTimeRangeList - get list of all PublishTimes,
|
|
# plus "All Grids" and "Selected Time" added at the top of
|
|
# the list
|
|
#
|
|
def _getTimeRangeList(self):
|
|
trList = ["All Grids", "Selected Time"]
|
|
publishTimes = self.getConfigItem("PublishTimes",[])
|
|
#inv = self._dbss.dataManager().selectTRMgr().inventory()
|
|
inv = self._dbss.getSelectTimeRangeManager().inventory()
|
|
for t in publishTimes:
|
|
if t in inv:
|
|
trList.append(t)
|
|
return trList
|
|
#================================================================
|
|
# _convertTimeRange - given a timeRangeName - return the
|
|
# timeRange (including the bogus names "All Grids" and
|
|
# "Selected Time").
|
|
#
|
|
def _convertTimeRange(self, trName):
|
|
if trName == "All Grids":
|
|
curTime = self.gmtime()
|
|
startHour = curTime[3]
|
|
timeRange = self.createTimeRange(startHour, 204, mode="Zulu")
|
|
#timeRange = TimeRange.allTimes()
|
|
# timeRange = AFPS.TimeRange(AFPS.AbsTime(0),
|
|
# AFPS.AbsTime_maxFutureTime())
|
|
elif trName == "Selected Time":
|
|
selectedTime = self._dbss.getParmOp().getSelectionTimeRange()
|
|
if selectedTime is None:
|
|
return None
|
|
else:
|
|
tr = self._dbss.getParmOp().getSelectionTimeRange()
|
|
return TimeRange.TimeRange(tr.getStart(), tr.getEnd())
|
|
else:
|
|
timeRange = self.getTimeRange(trName)
|
|
return timeRange
|
|
|
|
def _getCachedGrid(self, name, category, timeLimit):
|
|
try:
|
|
cachedObject = self.getObject(name, category)
|
|
grid, timeWritten, geoInfo = cachedObject
|
|
if timeLimit != 0 and time.time() - timeWritten > timeLimit:
|
|
return None #too old
|
|
|
|
# validate geoinfo
|
|
geo = self._dbss.getParmManager().compositeGridLocation()
|
|
if geoInfo != repr(geo):
|
|
return None #different geoinfo
|
|
return grid
|
|
except:
|
|
return None
|
|
|
|
def _cacheGrid(self, name, grid, category):
|
|
geo = self._dbss.getParmManager().compositeGridLocation()
|
|
objectToSave = grid, time.time(), repr(geo)
|
|
self.saveObject(name, objectToSave, category)
|
|
|
|
#========================================================================
|
|
#
|
|
# _checkViolate makes a mask of points INSIDE the SITE_EDITAREA that
|
|
# have a difference with a neighbors gridpoints (defined by areamask)
|
|
# that have a magnitude greater than threshold). Uses bits and criteria
|
|
# just like other calcuations
|
|
#
|
|
def _checkViolate(self, bits, criteria, areamask, discGrid, threshold):
|
|
violate = self.empty(bool)
|
|
for i in range(4): # range(8) to consider diagonal neighbors
|
|
#
|
|
# make sure data exists for both points
|
|
#
|
|
bitshift=self.offset(bits,self._ishift[i],self._jshift[i])
|
|
exist=logical_and(bitshift,bits)
|
|
#
|
|
# Make sure at least one of the points meets the criteria
|
|
#
|
|
critshift=self.offset(criteria,self._ishift[i],self._jshift[i])
|
|
meetcrit=logical_or(critshift,criteria)
|
|
#
|
|
# Make sure it borders the specified area
|
|
#
|
|
areashift=self.offset(areamask,self._ishift[i],self._jshift[i])
|
|
onborder=logical_and(areashift,self._siteAreaMask)
|
|
#
|
|
# Make sure it meets all criteria: exist, meetcrit, onborder
|
|
# and checkNDFD (meets topo thresholds)
|
|
#
|
|
mask=logical_and(logical_and(logical_and(exist,meetcrit),onborder),less(self._topodiff[i],self.MAXTOPODIFF))
|
|
#
|
|
violate[logical_and(mask, greater(abs(discGrid), threshold))] = True
|
|
return violate
|
|
|
|
|
|
#========================================================================
|
|
#
|
|
# _nbmask
|
|
#
|
|
# Replicate the NDFD neighboring points algorithm. It checks 'border
|
|
# pairs' where one point is 'inside' the CWA and the immediate
|
|
# neighbor (east-west-north-south) is 'outside' the CWA. A point
|
|
# inside the CWA might be used in several 'border pairs' - since it
|
|
# might have a neighbor outside the CWA to both the north and the west,
|
|
# etc.
|
|
#
|
|
# A 'border pair' is compared if:
|
|
# . The elevation difference between the two points is less
|
|
# than 1000 feet.
|
|
# . Both points have ISC data available.
|
|
# . The point 'outside' must be inside the area specified by
|
|
# areamask (usually the mask of the CWA you are checking against
|
|
#
|
|
# After finding the differences for each of the remaining border pairs,
|
|
# It averages them, and returns the average border pair difference, the
|
|
# average elevation difference for the pairs, the overall average of the
|
|
# values, the minimum of the average of the values 'inside' or the values
|
|
# 'outside', and the number of pairs used in this average.
|
|
#
|
|
def _nbmask(self, bits, isc, criteria, areamask, dirtype=0):
|
|
totalpts=0
|
|
totaldiff=0
|
|
totaltopo=0
|
|
totalvalue=0
|
|
totalvaluein=0
|
|
totalvalueout=0
|
|
for i in range(4): # range(8) to consider diagonal neighbors
|
|
#
|
|
# make sure data exists for both points
|
|
#
|
|
bitshift=self.offset(bits,self._ishift[i],self._jshift[i])
|
|
exist=logical_and(bitshift,bits)
|
|
#
|
|
# Make sure at least one of the points meets the criteria
|
|
#
|
|
critshift=self.offset(criteria,self._ishift[i],self._jshift[i])
|
|
meetcrit=logical_or(critshift,criteria)
|
|
#
|
|
# Make sure it borders the specified area
|
|
#
|
|
areashift=self.offset(areamask,self._ishift[i],self._jshift[i])
|
|
onborder=logical_and(areashift,self._siteAreaMask)
|
|
#
|
|
# Make sure it meets all criteria: exist, meetcrit, onborder
|
|
# and elevation difference less than MAXTOPODIFF (meets topo thresholds)
|
|
#
|
|
mask=logical_and(logical_and(logical_and(exist,meetcrit),onborder),less(self._topodiff[i],self.MAXTOPODIFF))
|
|
#
|
|
# Shift the data the directions
|
|
#
|
|
shift=self.offset(isc,self._ishift[i],self._jshift[i])
|
|
#
|
|
# Get the difference (different for directions)
|
|
#
|
|
if dirtype==1:
|
|
d1=self._dirdiff(isc,shift)
|
|
else:
|
|
d1=abs(isc-shift)
|
|
#
|
|
# Get number of points
|
|
#
|
|
pts=add.reduce(add.reduce(mask))
|
|
#
|
|
# get sum of differences for each shift direction
|
|
#
|
|
diff=sum(compress(mask.flat,d1.flat))
|
|
#
|
|
# get sum of topo differences
|
|
#
|
|
topo=sum(compress(mask.flat,self._topodiff[i].flat))
|
|
#
|
|
# get sums for points inside/outside
|
|
#
|
|
totalout=sum(compress(mask.flat,shift.flat))
|
|
totalin=sum(compress(mask.flat,isc.flat))
|
|
totalpts=totalpts+pts
|
|
totaldiff=totaldiff+diff
|
|
totaltopo=totaltopo+topo
|
|
totalvalue=totalvalue+totalout+totalin
|
|
totalvaluein=totalvaluein+totalin
|
|
totalvalueout=totalvalueout+totalout
|
|
#
|
|
# calculate total average
|
|
#
|
|
if (totalpts>0):
|
|
avg=totaldiff/totalpts
|
|
topoavg=totaltopo/totalpts
|
|
avgvalue=totalvalue/(totalpts*2)
|
|
avgvaluein=totalvaluein/totalpts
|
|
avgvalueout=totalvalueout/totalpts
|
|
minavgvalue=min(avgvaluein,avgvalueout)
|
|
else:
|
|
avg=0
|
|
topoavg=0
|
|
avgvalue=0
|
|
minavgvalue=0
|
|
return(avg,topoavg,avgvalue,minavgvalue,totalpts)
|
|
|
|
#========================================================================
|
|
#
|
|
# Get the smallest direction difference between two directions
|
|
#
|
|
def _dirdiff(self,dir1,dir2):
|
|
diff=abs(dir1-dir2)
|
|
less1=less(dir1,dir2)
|
|
diff=where(greater(diff,180.0),where(less1,dir1+360.0-dir2,dir2+360.0-dir1),diff)
|
|
return diff
|
|
|