Recently, I've been working on orbital planes intersecting stationary locations on the ground. It occurred to me that it might be helpful to plot those intersections. In doing so, I believe I encountered an inconsistency in how mpl_toolkits.basemap handles coordinates passed to Basemap.plot(). When using the ortho or nsper projection, it does not plot the same path as, for example, moll or mill , which both produce the expected output.
Here's a minimal working example:
import math
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
# Utilities:
RAD2DEG = 180 / math.pi
DEG2RAD = math.pi / 180
# Angles and coordinates:
i = 51.6390 * DEG2RAD # Orbital inclination of approx. 51 degrees matching ISS in radians
observer_lat = 28.5728722 * DEG2RAD # Latitude matching KSC in radians
observer_lon = -80.6489808 * DEG2RAD # Longitude matching KSC in radians
plt.figure(figsize=(14,14))
offset = -math.acos(-1/(math.tan(i) * math.tan(math.pi/2 - observer_lat))) # Offset for intersecting KSC
map = Basemap(projection='ortho', lat_0 = observer_lat * RAD2DEG, lon_0 = observer_lon * RAD2DEG)
map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color = 'whitesmoke')
map.drawmapboundary()
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
def path(lon): # Latitude as a function of longitude for the plane
return math.atan(-1/(math.tan(i)*(math.cos(lon + offset)))) + math.pi/2
lons = np.linspace(0, 2 * math.pi, 1000)
map.plot(lons * RAD2DEG, (np.vectorize(path)(lons) * RAD2DEG + 90) % 180 - 90, color='r', latlon=True)
map.plot(observer_lon * RAD2DEG, observer_lat * RAD2DEG, color='blue', marker='.', latlon=True)
plt.show()
The above produces this unexpected output. Changing only the projection to moll produces this output. The input points have not changed - nor has anything else - yet the outputs differ. Note how the Mollweide projection has the red line intersecting the blue dot, whereas the orthogonal projection doesn't come close to intersecting.
I have scoured the documentation, but can only find information on the latlon flag. I also tried removing that flag altogether by using an approach similar to the second example here. Nevertheless, I consistently produce non-matching outputs between projections.
I'd love to correctly plot using either the ortho or nsper projection. Any help is much appreciated.