Saturday, 31 August 2013

Circular mean for directional data

When working with angles or time (minutes, days, months or years) the mean (average) of a list of values can't be determined in the same fashion as we normally do.  This example show one how to calculate the values of angles (azimuths specifically).  The function can be used in any script once the CircularMean.py script is imported.  I will get around to other descriptive measures for circular data.

'''
CircularMean

Author
Dan Patterson
 
obidangis@gmail.com
  Dan_Patterson@carleton.ca
 Date  Aug 2013
Purpose:
  To calculate the mean (average) for points along a circle using the
  angular measure to these points from the center of the circle

Notes:
  If the angles are given in degrees then there is no need to convert
  If the units are time, such as hours (24 hour clock) or months (12) then
  these need to be converted by
    angle in degrees = time quantity / number of equal intervals
    ie 0.5 = (12 hours / 24 hours per day)
    or 180 degrees
required modules
arcpy, math

'''
def circularMean(angleList):
  '''
  angles to be given in degrees
  math module needs to be imported in the main module
  '''
  cosList = []
  cosSum = 0.0
  sinSum = 0.0
  for i in angleList:
    theCos = math.cos(math.radians(float(i)))
    theSin = math.sin(math.radians(float(i)))
    cosSum += theCos
    sinSum += theSin
  N = len(angleList)
  C = cosSum/N
  S = sinSum/N
  theMean = math.atan2(S,C)
  if theMean < 0.0:
    theMean += math.radians(360.0)
  return theMean

if __name__ == "__main__":
  import math
  sampleData = [ [350,10],
                 [350, 5],
                 [355, 10],
                 [90,270],
                 [180,0],
                 [180,360],
                 [0,360],
                 [90,180,270]
               ]
  print "\nCircular Mean Demo"
  print  "%20s %20s %20s" % ("angle List", "Mean radians", "Mean degrees")

  for angleList in sampleData:
    theMean = circularMean(angleList)
    print  "%20s %20.8f %20.8f" % (angleList, theMean, math.degrees(theMean))

No comments:

Post a Comment