# How to create a QQ plot of azimuths to test rotational symmetry of a spherical point dataset?

I am trying to do some Q-Q plots and the Kuiper test for rotational symmetry about the mean direction. These are unit vector, spherical data. What I am struggling to understand is the rotation from my own directions, where φ is the azimuth direction (0 to 360) and θ is the dip up/down direction (0 to 180), to the rotated/modified φ' and θ'. I'll limit this posting to the Q-Q plot question.

Here is a small example set of 10 azimuth/φ angles: 146.6 198.0 182.5 148.5 183.0 130.5 344.9 355.5 82.3 166.2

So, first the graphical method: Q-Q plot. Images, including #1 below, are from Fisher, Lewis and Embleton 1987 (reference below). I have the φ values, but the plot requires the rotated/modified φ' values. Now below in Image #2 is the method for doing the rotation, equation 3.8. I'm thinking I can set the ψ0 to zero, so I could use the "simpler" equation 3.9. But equations 3.8 or 3.9 require equations 3.5 and 3.6, which are below in Images #3 and #4 respectively. I'm able to calculate the Resultant length (R) and x-,y- and z-hat, so I have those, and which allows me to get φ-hat and φ-hat. So, I'm not sure what to do with equation 3.9 or 3.9, or how to rotate my numbers. Any suggestions? Apologies for the long post of images.  Fisher, N.I., T. Lewis and B.J.J. Embleton. 1987. Statistical Analysis of Spherical Data. Cambridge, Cambridge University Press.

## 1 Answer

You just want to study the azimuths of a set of spherical points $P_i$ relative to their spherical mean $\bar P$. The most straightforward solution solves the spherical triangles $(N,\bar P, P_i)$ where $N$ is the North Pole.

Let the co-latitudes of the points $P_i$ and $\bar P$ (angles from North) be $a$ and $b$ respectively. Let $\gamma$ be the angle between them: it's just the difference between the longitudes of the same two points. The azimuth, with due East being zero and orienting angles counterclockwise, is determined by

$$\arctan_2(\sin(b)\cos(a) - \cos(b)\sin(a)\cos(\gamma),\ \sin(a)\sin(\gamma))$$

where $\arctan_2(y,x)$ is the angle of a point $(x,y)$ in the plane. (This is supposed to be a numerically stable version of the formula, but I haven't tested it extensively.)

In this example, $100$ points were generated according to a (symmetric) Fisher-von Mises distribution distributed throughout the southern and western hemispheres, along with another $50$ points focused in the south and east. The resulting distribution is not symmetric. The mean point is shown as a red triangle.

Relative to the mean point, there is a cluster of points to its right (East) and upward (North), creating a swath of azimuths in the QQ plot between $0$ and $1$ (expressed in radians). The diffuse cluster to its west creates a broader swath of azimuths between $3$ and $5$. The QQ plot is clearly not uniform (for otherwise it would lie close to the dashed diagonal line), reflecting the bimodality of the spherical point distribution.

The R code that produced this example can be used to generate azimuthal QQ plots for any data. It assumes the spherical coordinates are provided as rows in an array; the relevant rows are indexed by "phi" and "theta".

#
# Spherical triangle, two sides and included angle given.
# Returns the angle alpha opposite a, in radians between 0 and 2*pi.
#
SAS <- function(a, gamma, b) {
atan2(sin(b)*cos(a) - cos(b)*sin(a)*cos(gamma), sin(a)*sin(gamma)) %% (2*pi)
}
#
# Cartesian coordinate conversion (for generating points).
#
xyz.to.spherical <- function(xyz) {
xyz <- matrix(xyz, nrow=3)
x <- xyz[1,]; y <- xyz[2,]; z <- xyz[3,]
r2 <- x^2 + y^2
rho <- sqrt(r2 + z^2)
theta <- pi/2 - atan2(z, sqrt(r2))
phi <- atan2(y, x)
theta[x==0 && y==0] <- sign(z) * pi/2
return (rbind(rho, theta, phi))
}
#
# Generate random points on the sphere.
#
library(MASS)
set.seed(17)
n.1 <- 100
n.2 <- 50
mu.1 <- c(0,-1,-1/4) * 2        # Center of first distribution
mu.2 <- c(1,1,-1/2) * 5         # Center of the second distribution
Sigma <- outer(1:3, 1:3, "==")  # Identity covariance matrix
xyz.1 <- t(mvrnorm(n.1, mu.1, Sigma)) # Each column is a point
xyz.2 <- t(mvrnorm(n.2, mu.2, Sigma))
xyz <- cbind(xyz.1, xyz.2)      # The Cartesian coordinates
rtf <- xyz.to.spherical(xyz)    # The spherical coordinates (also in columns)
#
# Compute the spherical mean and the azimuths relative to that mean.
#
mean.rtf <- xyz.to.spherical(rowMeans(xyz))
a <- SAS(rtf["theta",], rtf["phi",]-mean.rtf["phi",], mean.rtf["theta",])
#
# Plot the data and a QQ plot of the azimuths.
#
par(mfrow=c(1,2))
plot(c(-pi, pi), c(-1,1), type="n",
xlab="Phi", ylab="Cos(theta)", main="Sample points and their mean")
abline(h=0, col="Gray") # The Equator
abline(v=0, col="Gray") # The Prime Meridian
points(rtf["phi",], cos(rtf["theta", ]), col="#00000080")
points(mean.rtf["phi",], cos(mean.rtf["theta",]), bg="Red", pch=24, cex=1.25)

plot(c(0,1), c(0,2*pi), type="n",
xlab="Quantile", ylab="Azimuth", main="Azimuthal QQ Plot")
abline(c(0, 2*pi), lty=3, lwd=2, col="Gray")
points(seq(0, 1, along.with=a), sort(a))

• Thank you. I'm still confused though. What is the subscript 2 on the equation beginning with arctan signifying? What do I do with that comma in the middle of the same equation? Or, what do I do with the two separate parts within the arctan calculation? My azimuths have 0 at the top, and then proceed clockwise around the equator. The Q-Q plot I'm interested in only uses the azimuths, but for my altitude/dip/co-latitude angles I use -90 to 0 (S hemisphere) and 0 to 90 (N hemisphere) rather than 0-180, for esoteric archaeological reasons. I've been using R since Dec. & the code is over my head – mdhunstiger Mar 23 '15 at 5:18
• Regarding the R code, so I would put my angles in two rows, say in a csv file made in excel? Say, call cell A1 phi and A2 theta, with the actual data following horizontally in each row? I can never (initially) understand how to simply set up data in R so whatever function I'm using will use it. It's trial and error. Where is the name of the file (array name?) that I need to call my data so that R knows where to find theta and phi in them? In other instances I've done it where a dataset is called something, say x, and then the column with the relevant numbers is given as x\$example. – mdhunstiger Mar 23 '15 at 5:28
• Arctangent functions come in two flavours, with one and two arguments. A two argument version is needed to give results covering the entire circle. @Whuber's notation with subscript 2 echoes one often used in software where some function name like ATAN2 is used. For more see e.g. mathworld.wolfram.com/InverseTangent.html – Nick Cox Mar 23 '15 at 10:52
• Your specific questions in comments about basic difficulties in setting this up in R are pointing off-topic. There is a fine line here: @whuber's answer explained principles and provided code as back-up, but once the focus is entirely on code Cross Validated is no longer the best forum. See the Help Center for guidance on software-related questions. – Nick Cox Mar 23 '15 at 10:56
• Indeed the R name atan2() encapsulates the same distinction. – Nick Cox Mar 23 '15 at 10:59