' Name: GPS.CalcCEP
'
' Author: Timothy N. Loesch
'         Minnesota Department of Natural Resources
'         GIS Applications Coordinator
'         500 Lafayette Road, Box 11
'         St. Paul, MN 55155
'         tim.loesch@dnr.state.mn.us
'         (651) 296-0654
'
' Date: Thu Jan 17 10:18:13 2002
' Revised by: 
' Revision Date: 
' Revisions: 
' -------------------------------------------------------------------
' Description: This script calculates Circular Error of Probability rings for 
'              a specific dataset.
' Requires: 
' Runs: 
' Run by: 
' Self: 
' Returns: 
'-------------------------------------------------------------------

theview = av.getactivedoc
theprj = theView.GetProjection
MapUnits = TheView.GetDisplay.GetUnits
DistUnits = theView.GetDisplay.GetDistanceUnits

if ((theprj.IsNull) and (mapUnits = #UNITS_LINEAR_DEGREES)) then
  msgbox.error("Unable to calculate CEP for unprojected decimal degree data. Set a projection for your view and try again.",Script.The.GetName)
  return nil
end
thethemes = theView.GetActiveThemes
if ( (thethemes.count = 1) and (thethemes.get(0).getftab.GetShapeClass.GetClassName = "Point") ) then
  thetheme = thethemes.get(0)
 else
  thethemes = theView.GetThemes
  if (thethemes.count < 1) then
    msgbox.error("No themes. Unable to comply",Script.The.GetName)
    return nil
  end
  pthms = {}
  for each athm in thethemes
    if (athm.getftab.GetShapeClass.GetClassName = "Point") then
      pthms.add(athm)
    end
  end
  if (pthms.count < 1) then
    msgbox.error("No point themes. Unable to comply",Script.The.GetName)
    return nil
  elseif (pthms.count = 1) then
    thetheme = pthms.get(0)
  else
    thetheme = msgbox.ListAsString("Select Point Theme",Script.The.GetName)
    if (thetheme = nil) then
      return nil
    end
  end
end

theftab = thetheme.getftab
shpfld = theftab.FindField("SHAPE")
if (shpfld = nil) then
  msgbox.Error("unable to find shape field. Aborting....",Script.The.GetName)
  return nil
end
basepoint = Msgbox.MultiInput("Enter Known Coordinates - CANCEL for NONE","Base Coordinate",{"X-Coord","Y-Coord"},{"",""})
if (basepoint.count <1) then
  DOAVERAGE = msgbox.yesno("Do you want to use the Average X,Y as the base point?","Base Coordinate",TRUE)
  if (DOAVERAGE.Not) then
    return nil
  end
 else
  DOAVERAGE = FALSE
  x = Basepoint.get(0)
  y = BasePoint.Get(1)
  if ((x.IsNumber.Not) or (y.IsNumber.Not) ) then
    msgbox.Error("Invalid Coordinate. Unable to continue.",Script.The.GetName)
    return nil
  end
  basepoint = point.make(x.AsNumber,y.AsNumber)
end

thePrecision = "d.dddddddddd"
xSum = 0
ysum = 0

theCount = 0
xmin = nil
xmax = nil
ymin = nil
ymax = nil

pointDict = Dictionary.Make(theFtab.GetNumRecords)

for each rec in theftab
  theshp = theftab.returnvalue(shpfld,rec)
  if (theprj.IsNull.not) then
    theshp = theshp.ReturnProjected(theprj)
  end
  pointDict.Add(rec.clone,theshp)
  xVal = theshp.getx
  yval = theshp.gety
  if ( xmin = nil ) then
    xmin = xVal
    xmax = xVal
  else
    xmin = xMin min xval
    xmax = xMax max xval
  end
  if ( ymin = nil ) then
    ymin = yVal
    ymax = yVal
  else
    ymin = yMin min yval
    ymax = yMax max yval
  end
  xSum = xval + xSum
  ySum = yval + ySum
  theCount = theCount + 1
end

xMean = xSum / thecount
yMean = YSum / thecount

xSumSqDev = 0
ySumSqDev = 0
for each rec in theftab
  theshp = theftab.returnvalue(shpfld,rec)
  xValue = theshp.getx
  yvalue = theshp.gety
  if ( not ( xValue.IsNull ) ) then
    xSqDev = ( xValue - xMean ) * ( xValue - xMean )
    xSumSqDev = xSqDev + xSumSqDev
  end
  if ( not ( yValue.IsNull ) ) then
    ySqDev = ( yValue - yMean ) * ( yValue - yMean )
    ySumSqDev = ySqDev + ySumSqDev
  end
end

if (theCount > 1) then
  xVariance = xSumsqdev / (theCount - 1)
  xStdDev = xVariance.Sqrt
  yVariance = ySumsqdev / (theCount - 1)
  yStdDev = yVariance.Sqrt
else
  theVariance = 0
  theStdDev = 0
end

avpnt = point.make(xmean,ymean)

if (DOAVERAGE) then
  thepoint = avpnt
 else
  thepoint = basepoint
end
thepntgr = GraphicShape.Make(thepoint)
theView.GetGraphics.Add(thepntgr)
thetheme.GetGraphics.add(thepntgr)

DistList = {}
for each apntId in pointDict.ReturnKeys
  theDist = thePoint.Distance(PointDict.Get(apntID)).clone
  DistList.Add(theDist)
end



' now sort the list of distances
DistList.Sort(TRUE)
cep50 = DistList.Get( (DistList.count * 0.50).Truncate)
cep90 = DistList.Get( (DistList.count * 0.90).Truncate)
cep95 = DistList.Get( (DistList.count * 0.95).Truncate)
cep98 = DistList.Get( (DistList.Count * 0.98).Truncate)

thecirc = Circle.Make(thepoint,cep50)
CEP50gr = graphicshape.make(Circle.Make(thepoint,cep50))
theView.GetGraphics.Add(CEP50gr)
thetheme.getgraphics.add(CEP50gr)

CEP90gr = graphicshape.make(Circle.Make(thePoint,cep90))
theView.GetGraphics.Add(CEP90gr)
thetheme.getgraphics.add(CEP90gr)

CEP95gr = GraphicShape.make(Circle.Make(thePoint,cep95))
theView.GetGraphics.Add(CEP95gr)
theTheme.GetGraphics.Add(CEP95gr)

CEP98gr = GraphicShape.Make(Circle.Make(thePoint,cep98))
theView.GetGraphics.Add(CEP98gr)
thetheme.getgraphics.add(CEP98gr)

thetheme.Clearselection

theView.Invalidate

If (MapUnits = #UNITS_LINEAR_UNKNOWN) then
  convFact = 1
else
  if (DistUnits <> #UNITS_LINEAR_UNKNOWN) then
    convFact = Units.Convert(1,MapUnits,DistUnits)
   else
    convfact = 1
  end
end

msg = "Statistics"+NL
msg = msg + "Average = "+xmean.AsString++ymean.AsString + nl
msg = msg + "SD      = "+xStdDev.AsString ++ yStdDev.AsString + nl
msg = msg + "Circular Error Probabilities (CEP) - "+Units.GetFullUnitString(DistUnits)+NL
msg = msg + "    50% = "+(CEP50 * convFact).AsString +NL
msg = msg + "    90% = "+(CEP90 * convFact).AsString +NL
msg = msg + "    95% = "+(CEP95 * convFact).AsString +NL
msg = msg + "    98% = "+(CEP98 * convFact).AsString +NL
msgbox.report(msg,"")
""theView.Invalidate