Does mpl_toolkits.basemap handle coordinates consistently throughout projections?

41 views Asked by At

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.

0

There are 0 answers