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 ------------------------------------------------------------------------------ All the data continuously generated in your IT infrastructure contains a definitive record of customers, application performance, security threats, fraudulent activity, and more. Splunk takes this data and makes sense of it. IT sense. And common sense. http://p.sf.net/sfu/splunk-novd2d _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
