drawgreatcircle, directions, and equator

[derek-dalle]

It appears that drawgreatcircle() has some sign issues, especially when an equator crossing is included.

I've attached a picture showing the resulting paths created by drawgreatcircle() using a few major cities. I called great circle in both directions for each pair of city. It appears that it only works if the direction is west-to-east and no equator crossing is involved. The Tokyo-to-LA path shows the expected behavior, so crossing 180 longitude is ok.

Ignore the red to in the middle of the Atlantic Ocean, which I put there manually for testing.

My basemap version is whatever is in the Ubuntu 12.10 standard repositories.

[jswhit]

Could you post a script here that triggers the problem? I can't reproduce it here.

[derek-dalle]

Indeed I had some trouble reproducing this problem with a simple script. I took this example from http://matplotlib.org/basemap/users/examples.html and did it in the opposite direction expecting to get the error shown above. But it turns out that it only has a problem if the latitudes and longitudes are not regular floats. The example below works normally if the np.float64 calls are removed.

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# create new figure, axes instances.
fig=plt.figure()
ax=fig.add_axes([0.1,0.1,0.8,0.8])
# setup mercator map projection.
m = Basemap(llcrnrlon=-180.,llcrnrlat=-60.,\
            urcrnrlon=180.,urcrnrlat=60.,\
            projection='merc',\
            lat_0=40.,lon_0=-20.,lat_ts=20.)
m.drawcoastlines()
m.fillcontinents()
# nylat, nylon are lat/lon of New York
nylat = np.float64(40.78)
nylon = np.float64(-73.98)
# lonlat, lonlon are lat/lon of London.
lonlat = np.float64(51.53)
lonlon = np.float64(0.08)
# draw great circle route between NY and London
m.drawgreatcircle(nylon,nylat,lonlon,lonlat,linewidth=2,color='b')
# now go the other direction
m.drawgreatcircle(lonlon,lonlat,nylon,nylat,linewidth=2,color='b')
# draw parallels
m.drawparallels(np.arange(-60,60,30))
# draw meridians
m.drawmeridians(np.arange(-180,180,30))
plt.show()

[jswhit]

Your example works with the most recent basemap. The great circle computation code in pyproj has been reworked, and in the process the bug you mention has been fixed.