Discretize a circle on a sphere with a given center and radius

Question

I would like to draw a discretized circle on the surface of a sphere (the Earth in that case).

The input of the algorithm would be the center of the circle (expressed as longitude and latitude), the radius (expressed in meters), and the number of points to calculate.

The output would be the coordinates of each calculated point.

I am guessing this involve the Great Circle formula but I can't see how. Any help would be appreciated.

Regards

Answer 1

The easiest way to do this would probably be to first find the coordinates of a circle at one of the poles, and then use rotation matrix transforms to move the circle to the desired location. At a pole, the points have coordinates: $$(r\cos\phi,r\sin\phi,\sqrt{R^2-r^2})$$ Which you can evenly space by using $\phi=2n\pi/N$. Once you have an array of points, you want to transform each point by rotating it around the $X$ axis $\pi-\theta_{lat}$ degrees, and then around the $Z$ axis $\phi_{long}$ degrees.

Answer 2

Let $p(\theta, \phi) = r(\cos \phi \cos \theta, -\cos \phi \sin \theta,\sin \phi)$ (latitude is $\phi$ measured north, longitude is $\theta$ measured west, the point $0°$, $0 °$ corresponds to $(r,0,0)$, the north pole is $(0,0,r)$). This gives the $(x,y,z)$ coordinates of the point, where $(0,0,0)$ is the center of the earth.

You are given a center $x_c = p(\theta_c, \phi_c)$, and a distance $d$. The desired circle consists of all points that are a distance $d$ from $x_c$, lying on the surface, ie, $C = \{x \in \mathbb{R}^3 | \|x\|=r,\ \|x-x_c\| = d \}$.

To parameterize points in $C$, we will express points as $x = \lambda x_c + p$, where $p \bot x_c$. Since $p \bot x_c$, we have $\lambda = \frac{\langle x , x_c \rangle}{\|x_c\|^2} = \frac{\langle x , x_c \rangle}{r^2}$, and since $d^2 = \|x-x_c\|^2 = \|x\|^2 + \|x_c\|^2 - 2 \langle x , x_c \rangle $, we have $\langle x , x_c \rangle = \frac{1}{2}(2 r^2 -d^2)$, hence $\lambda = 1 - \frac{1}{2}(\frac{d}{r})^2$. Since $p \bot x_c$, we also have $\|x\|^2 = |\lambda|^2 \|x_c\|^2 + \|p\|^2$, or $\|p\| = d_\bot = r\sqrt{1-\lambda^2} = d \sqrt{1-(\frac{d}{2r})^2}$.

Now suppose the subspace perpendicular to $x_c$ is given by $\{x_c\}^\bot = \operatorname{sp} \{ p_1, p_2 \}$, where the $p_i$ are orthonormal. Then $C$ can be parameterized by $C = \{\lambda x_c + d_\bot((\cos \alpha) p_1 + (\sin \alpha) p_2) \}_{\alpha \in [0,2 \pi)}$.

To actually do the computation, compute $x_c = p(\theta_c, \phi_c)$, compute one orthonormal vector $p_1 = \frac{p(\theta_c, \phi_c+\frac{\pi}{2})}{\| p(\theta_c, \phi_c+\frac{\pi}{2}) \|}$, and compute $p_2 = \frac{x_c}{r} \times p_1$ ($\times$ is the cross product). To compute $n$ points on the circle, compute $\lambda x_c + d_\bot((\cos (2 \pi \frac{k}{n})) p_1 + (\sin (2 \pi \frac{k}{n})) p_2)$ for $k = 0,...,n-1$.

(To compute the latitude & longitude of a point (x,y,z) on the surface, use $\phi = \arcsin \frac{z}{r}$, $\theta = \operatorname{atan2}(y,x)$.)

Note: The circle radius ($d$ above) is the Euclidean distance. If you want the great circle distance instead, say $\gamma$, then use $d = 2 r \sin \frac{\gamma}{2r}$ in the above computation.