Skip to content

fixes to basemap using celestial coordinates #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 2 commits into from

Conversation

mollyswanson
Copy link

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')

@jswhit
Copy link

jswhit commented Nov 29, 2011

applied to master - thanks!

@jswhit jswhit closed this Nov 29, 2011
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants