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))