We previously analyzed the equirectangular projection of great circles using sampling and simulations. This was done mostly to visualize the shape of the curve and understand its behavior. This time, let’s try to find an analytical solution to describe the curve of a great circle in an equirectangular map.
A great circle is the intersection between a sphere and a plane, with the plane passing by the center of the sphere. For simplicity, we will assume the sphere is centered on the origin. We don’t lose generality with this assumption, since the equirectangular coordinates are centered on the sphere center.
The coordinates in the equirectangular projection are the longitude \(\theta\) and the latitude \(\phi\).
Our goal is to find the function representing the intersection of a plane and a sphere in those coordinates. The function of a plane centered on the origin in Cartesian coordinates is
$$ax + by + cz = 0$$
where \(\vec{n} = (a,b,c)\) is the normal of the plane. The function of a sphere centered on the origin in spherical coordinates is
$$ x = r \cos(\theta)\cos( \phi) $$
$$ y = r \sin(\theta)\cos(\phi) $$
$$ z = r \sin(\phi) $$
where \(r\) is the radius of the sphere. By replacing the variables in the plane equation with the ones from the sphere equations, we get the function of the intersection
$$ a r \cos(\theta)\cos(\phi) + b r \sin(\theta)\cos(\phi) + c r \sin(\phi) = 0 .$$
Dividing by \(r\), we can remove the radius from the equation and we obtain
$$ a \cos(\theta)\cos(\phi) + b\sin(\theta)\cos(\phi) + c \sin(\phi) = 0 $$
meaning that the sphere radius as no effect on the shape of the great circle in the equirectangular projection. With some algebraic manipulations, we get
$$ \frac{a \cos(\theta) + b\sin(\theta)}{-c} = \frac{\sin(\phi)}{\cos(\phi)}$$
which can be simplified again with trigonometric identities to get
$$ \frac{a \cos(\theta) + b\sin(\theta)}{-c} = \tan(\phi).$$
Finally, isolating \(\phi\) we get
$$ \phi = \arctan\left(\frac{a \cos(\theta) + b\sin(\theta)}{-c}\right) = gc(\theta)$$
where \(gc(\theta)\) is the function mapping longitude \(\theta\) to latitude \(\phi\) along a great circle with parameters \((a,b,c)\) in an equirectangular projection.
Remember that \(\vec{n} = (a,b,c)\) is the normal of the plane. This unit vector can be re-parameterized using two variables representing the orientation of the plane in spherical coordinates : azimuth \(\alpha\) and inclination \(\beta\). The equations mapping the new parameters to the original parameters are
$$ a = \sin(\beta) \sin(\alpha) $$
$$b = -\sin(\beta) \cos(\alpha) $$
$$c = \cos(\beta) $$
We could directly replace \(a\) , \(b\) and \(c\) in \(gc(\theta)\) using these three equations, but the function would be a bit unwieldy. There is an additional simplification we can apply to obtain a cleaner function. If we analyze the function when the azimuth \(\alpha\) of the great circle is 0, i.e. when the great circle is crossing the equator at longitude 0°, we get
$$ a =0 $$
$$b = -\sin(\beta) $$
$$c = \cos(\beta) $$
and if we insert those in \(gc(\theta)\), we obtain
$$ gc_{\alpha=0}(\theta) = \arctan\left(\tan(\beta)\sin(\theta)\right).$$
While \(gc_{\alpha=0}(\theta)\) is only valid for \(\alpha=0\), it is possible to apply a rotation of \(\alpha\) to the sphere to get an equivalent great circle. After applying this rotation, we get the final form of \(gc(\theta)\)
$$ gc(\theta) = \arctan\left(\tan(\beta)\sin(\theta – \alpha)\right) = \phi.$$
Using this function, we can plot any great circle in an equirectangular projection using two parameters : azimuth \(\alpha\) and inclination \(\beta\). First, we plot a series of great circles with various inclinations, but with a fixed azimuth of 0°. The result is presented in figure 1.
As expected, this figure is very similar to figure 3 in the previous article. The only difference is a symmetry along the longitude axis. This is due to a change in the convention used for the inclination angle.
Then, we repeat the previous experiment, but with a fixed inclination of 85° and a series of different azimuths. The results are displayed in figure 2.
Finally, let’s plot the trajectory of the International Space Station in an equirectangular map using \(gc(\theta)\). The ISS has an inclination of 51.6°. Its azimuth is constantly changing due to the rotation of the Earth, but let’s use an azimuth of -130° because that makes the trajectory pass over Montreal. After filling these inclination and azimuth values in \(gc(\theta)\), we get the trajectory displayed in figure 3. The Earth image used for the background was obtained from NASA’s Blue Marble.
If you compare this synthetic trajectory to the real one presented in the previous article, we see that they are very similar, but the real trajectory doesn’t return to the same point over the planet after completing an orbit. This is due to the planet rotation.
The Python code used to generate the figures is available.