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

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

• coordinates in 3D or in some projected 2D space? – nbubis Mar 11 '13 at 18:58
• This is for a 2D GIS application but I would like to use 3D space to avoid GIS projection issues – htulipe Mar 11 '13 at 20:18
• Is the radius of the small circle measured around the surface of the earth, straight through the earth from the given center point, or perpendicularly from a point within the earth that is the true center of the circle? For small circles it won't matter. then are the output in $x,y,z$ or lat/lon? – Ross Millikan Mar 11 '13 at 21:19
• The idea is to stick to sphere surface (hence the Great Circle distance) and the output would be in lat/lon for direct use in GIS application – htulipe Mar 11 '13 at 21:24

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.

• Can't I say that the pole is my circle center? That would save a matrix transform. Anyway, what would be $n$ and $N$ in your $\phi$ definition? – htulipe Mar 11 '13 at 21:32
• $N$ - number of points, $n$ the index. – nbubis Mar 11 '13 at 23:49
• What are R2 and r2? – Charles Holbrow Apr 13 '15 at 19:07
• @CharlesHolbrow - $R$ is the sphere radius, $r$ the circle radius. – nbubis Apr 14 '15 at 7:27

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.