3

Measurements of angles (e.g. compass bearings, aspect) are circular in nature (e.g. 0 and 360 degrees coincide) and therefore cannot be meaningfully described by ordinary arithmetic means. I have a shapefile of 22 locations containing aspect value for each point. I want to calculate mean aspect of the locations using circular statistics in R. "circstats" is one of the packages that are suggested to me. Question: Do I need to transform aspect data before feeding into R (aspect data is calculated in GRASS)? How I can calculate the mean direction of aspect in R?

7

In R, the package CircStats is old and of rather limited scope and has been replaced by the more complete Circular package. There are tutorials and a book, Circular Statistics with R (2013, A. Pewsey, M. Neuhäuser, and G. D. Ruxton, Oxford University Press, 208 pp.) which explains how to use it (The R scripts can be downloaded from the resources site of the book)

1) If the data are directional (= circular with values in the entire 0–360◦ range and 0-360° equivalent), the solution is to use the angular values as vectors an apply the basic vector operations. The mean direction is then found by computing the vector sum of the vectors that represent the various directions in the data.

enter image description here

Since the coordinates of these vectors (normalized) are Xi =cosθi and Yi =sinθi, the mean direction is given by atan2(sum(sinθi),sum(cosθi)) ( -> the solution of RadouxJu)

2) If the data are axial (two-headed vectors with a 180◦ ambiguity, as a line) there are others solutions.

Therefore for circular data, you have many solutions with GRASS GIS:

  • run R in the GRASS shell and open the GRASS vector with the spgrass6 (GRASS 6.4.x) or rgrass7 (GRASS 7.x) packages (look at spgrass6 and GIS_LOCK, for vectors, the command is readVECT6) and use the Circular package

    library(circular)
    angles = c(341.0, 359.0, 334.0, 15.0, 330.0, 301.0, 299.0, 9.0, 7.0, 353.0, 28.0, 25.0, 23.0, 30.0, 350.0, 25.0, 22.0, 8.0, 356.0, 27.0)
    anglecir =  circular(angles, type="angles", units="degrees",modulo="2pi", template='geographics')
    summary(anglecir)
         n     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      Rho 
    20.0000  30.0000  23.5000   7.5000   1.0610 -12.2500 -61.0000   0.8962 
    mean(anglecir)
    Circular Data: 
    Type = angles 
    Units = degrees 
    Template = geographics 
    Modulo = 2pi
    Zero = 1.570796 
    Rotation = clock 
    [1] 1.060902
    
  • run Python from the GRASS shell (with grass.script or pygrass modules) or from the Python shell and use the solution of radouxju, very easy to implement in Python.

    import numpy as np
    angles = np.array([341.0, 359.0, 334.0, 15.0, 330.0, 301.0, 299.0, 9.0, 7.0, 353.0, 28.0, 25.0, 23.0, 30.0, 350.0, 25.0, 22.0, 8.0, 356.0, 27.0])
    def getCircularMean(angles):
        n = len(angles)
        sineMean = np.divide(np.sum(np.sin(np.radians(angles))), n)
        cosineMean = np.divide(np.sum(np.cos(np.radians(angles))), n)
        vectorMean = np.arctan2(sineMean, cosineMean)
        return np.degrees(vectorMean)
    print "{:.6f}".format(getCircularMean(angles))
    1.060902
    
|improve this answer|||||
  • Good references. However, you ought to update your arctan solution, which is invalid (as well as being syntactically incorrect in R), with a correct one based on atan2. – whuber May 13 '15 at 17:34
  • done, thanks, but it was pure formalism here because the real trigonometric function is Arctanwith atan2 instead of atan for computations and correct results. – gene May 13 '15 at 18:03
  • gene... and of course this works in ArcMap's python shell as well, for completeness – user681 May 14 '15 at 18:23
  • ArcGIS has also the command Linear Directional Mean (Spatial Statistics), but for lines only (not points with attributes values) – gene May 14 '15 at 18:55
2

you can convert your aspects into the sine and cosine, compute the mean of the sine's and the mean of the cosine's, then turn it back to aspect using atan2(sine,cosine). For more details, see Wikipedia

|improve this answer|||||
  • Do I need also to transform aspect degrees to radius before changing to sin/cos? – Mah May 13 '15 at 14:20
  • by using atan(sine/cosine), I have a negative number above 1. Is it suppose to be between -1 and 1? – Mah May 13 '15 at 15:32
  • 1
    An arctangent is an angle--it can be anything. However, you shouldn't be using atan at all--use the atan2 version (called atan(x,y) in the GRASS manual; its value is in degrees). Otherwise you will only get a 180 degree range of angles rather than the full 360 degrees. – whuber May 13 '15 at 15:42
  • @whuber Thank you, with atan2(x,y) you don't need to test the sign of sine and cos. I've updated my answer. – radouxju May 13 '15 at 17:30
  • @Mah yes, R works with angles in radian, so you need to transform from degree to radian before sin and cos, then from radian to degree after atan2 – radouxju May 13 '15 at 17:31
1

Most textbooks suggest using atan2(Sigma(sin(x)), Sigma(cos(x))), however this is not always the right thing to do. For example, the average of 0, 0 and 90 degrees is atan( (sin(0)+sin(0)+sin(90)) / (cos(0)+cos(0)+cos(90)) ) = atan(1/2)= 26.56 deg, and not 30 deg as one may expect.

Take a look at my article on CodeProject "Circular Values Math and Statistics with C++11", especially section 26 - Averaging n circular values, and section 29 - Circular parameter estimation based on noisy measurements, were I suggest a method that is consistent with the general definition on an average, and also more accurate when the data has a wrapped normal distribution.

An equivalent method is suggested in Edwin Olson's paper "On computing the average orientation of vectors and lines" (Algorithm 1).

|improve this answer|||||
  • This is an interesting and valuable observation (+1). Ultimately, the choice of how to average directional values depends on the purpose of the analysis as well as the nature of the random errors (the "noise"). Aspects are computed from DEMs, typically, and as such tend to be functions of least squares estimates of locally fitted surfaces. Although the distributions of those parameter estimates tend to be approximately Normal, the aspects (as angles) will not be. In effect, the atan2(sin,cos) procedure appears to be a more appropriate mean than Olson's method. – whuber May 14 '15 at 19:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.