fixes to basemap using celestial coordinates #6
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
I would like to suggest the attached fixes for the basemap module to improve its treatment of celestial (rather than geographic) coordinates.
The first 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 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')