Circular mean for directional data

4891
0
01-10-2016 11:32 AM
Labels (1)
DanPatterson_Retired
MVP Emeritus
0 0 4,891

Original:    Aug 2013

Modified:  Sept 2015  to include Numpy version

Reposted:  Jan 2016

General references

-   http://webspace.ship.edu/pgmarr/Geo441/Lectures/Lec%2016%20-%20Directional%20Statistics.pdf

-   http://ncss.wpengine.netdna-cdn.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Circular_Data_Analysi...

-   SciPy implementation  scipy/morestats.py at master · scipy/scipy ·GitHublines 2601 - 2642  using complex numbers

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 circular_mean.py script is imported.

I will get around to other descriptive measures for circular data but they are covered in the references

"""
CircularMean
Author: Dan.Patterson@carleton.ca
Date:     Aug 2013
Modified: Sept 2015
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
numpy, math
"""
import math
import numpy as np
def circ_mean(angles):
    """angles in degrees"
    """
    cosList = []
    cosSum = 0.0
    sinSum = 0.0
    for i in angles:
        theCos = math.cos(math.radians(float(i)))
        theSin = math.sin(math.radians(float(i)))
        cosSum += theCos
        sinSum += theSin
    N = len(angles)
    C = cosSum/N
    S = sinSum/N
    theMean = math.atan2(S,C)
    if theMean < 0.0:
        theMean += math.radians(360.0)
    return theMean
def circ_mean_np(angles,azimuth=True):
    """ numpy version of above"""
    rads = np.deg2rad(angles)
    av_sin = np.mean(np.sin(rads))
    av_cos = np.mean(np.cos(rads))
    ang_rad = np.arctan2(av_sin,av_cos)
    ang_deg = np.rad2deg(ang_rad)
    if azimuth:
        ang_deg = np.mod(ang_deg,360.)
    return ang_rad, ang_deg
if __name__ == "__main__":

    sampleData = [ [350,10],
                   [350, 5],
                   [355, 10],
                   [90,270],
                   [180,0],
                   [180,360],
                   [0,360],
                   [90,180,270]
                 ]
    print("\nCircular Mean Demo...python version\n")
    frmt ="{:15s}{:15s}{:15s}"
    print(frmt.format("angles", "Mean radians", "Mean degrees"))
    for angles in sampleData:
        theMean = circ_mean(angles)
        args = (angles, theMean, math.degrees(theMean))
        print("{!s:>15}{:>12.4}{:>15.4}".format(*args))
    print("\nCircular Mean Demo...numpy version\n")
    frmt ="{:15s}{:15s}{:15s}"
    print(frmt.format("angles", "Mean radians", "Mean degrees"))
    #
    for angles in sampleData:
        result = circ_mean_np(angles)
        print("{!s:>15}{:>12.4f}{:>15.1f}".format(angles, result[0],result[1]))

Results from both

angles         Mean radians   Mean degrees   
      [350, 10]       6.283          360.0
       [350, 5]        6.24          357.5
      [355, 10]     0.04363            2.5
      [90, 270]       3.142          180.0
       [180, 0]       1.571           90.0
     [180, 360]       4.712          270.0
       [0, 360]       6.283          360.0
 [90, 180, 270]       3.142          180.0
About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...
Labels