I would like to suggest the following fixes for the
mpl_toolkits.basemap module to improve its treatment of celestial
(rather than geographic) coordinates.
The first one, posted at
https://github.com/mollyswanson/basemap/commit/23db4bbebf4d7fe6ca202b5dad50b6a2054dd685
changes the call function in basemap's init.py to correctly transform
lat/lon values into xy map coordinates in the case of a cyclic or
polycyclic projection with lon_0 not equal to 0.
The second, posted at
https://github.com/mollyswanson/basemap/commit/35470b51523e9429d26cefc911dca843264581b9
changes one line in the drawparallels. This line is to avoid drawing
lines between points on the parallel that span the whole map. However,
the old version uses a fixed value for the distance between the points
rather than scaling it to the radius of the sphere used in the
projection, so if you use a non-default radius (such as 180/pi, so
your x-y values are in degrees on the sky instead of meters on the
earth) it won't work. This fix scales the cutoff value to the radius
of the projection sphere.
The following example illustrates the issues that are addressed here:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
figure(1)
#make a basemap centered on longitude of 90
m=Basemap(celestial=False,lon_0=90,projection='hammer')
#draw map boundary and grid
m.drawmapboundary()
m.drawparallels(np.arange(-90.,91.,30.),labels=[1,0,0,0])
m.drawmeridians(np.arange(-90.,271.,30.),labels=[0,0,0,0])
#define a test polygon - a triangle with corners at
[lon,lat]=[90,30],[120,60],[120,30]
polygon=array([[90,30],[120,60],[120,30],[90,30]])
#convert to map coordinates
polyxy=m(polygon[:,0],polygon[:,1])
plt.plot(polyxy[0],polyxy[1])
plt.savefig('basemap1.png')
figure(2)
#make a celestial basemap centered on longitude of 90
m=Basemap(celestial=True,lon_0=90,projection='hammer',rsphere=180./pi)
#draw map boundary and grid
m.drawmapboundary()
m.drawparallels(np.arange(-90.,91.,30.),labels=[1,0,0,0])
m.drawmeridians(np.arange(-90.,271.,30.),labels=[0,0,0,0])
#define a test polygon - a triangle with corners at
[lon,lat]=[90,30],[120,60],[120,30]
polygon=array([[90,30],[120,60],[120,30],[90,30]])
#convert to map coordinates
polyxy=m(polygon[:,0],polygon[:,1])
plt.plot(polyxy[0],polyxy[1])
plt.savefig('celestial_basemap1.png')
Thank you!
Molly Swanson