Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

fixes to basemap using celestial coordinates #6

Closed
wants to merge 2 commits into from

2 participants

@mollyswanson

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
Collaborator

applied to master - thanks!

@jswhit jswhit closed this
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Nov 29, 2011
  1. @mollyswanson

    Fixes issue with using lon_0 not equal to 0 for celestial coordinates…

    mollyswanson authored
    … in a cyclic or polycyclic projection.
  2. @mollyswanson

    fixes bug in drawparallels when the radius of the sphere of the proje…

    mollyswanson authored
    …ction is set to something other than the default.
This page is out of date. Refresh to see the latest.
Showing with 9 additions and 6 deletions.
  1. +9 −6 lib/mpl_toolkits/basemap/__init__.py
View
15 lib/mpl_toolkits/basemap/__init__.py
@@ -835,7 +835,7 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None,
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
dist = np.sqrt(xd+yd)
- split = dist > 5000000.
+ split = dist > 5000000./6370997.0*self.projparams['R']
if np.sum(split) and self.projection not in _cylproj:
ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
iprev = 0
@@ -887,20 +887,23 @@ def __call__(self,x,y,inverse=False):
Input arguments lon, lat can be either scalar floats,
sequences, or numpy arrays.
"""
+ if (self.projection in _pseudocyl) or (self.projection in _cylproj):
+ lon_0=self.projparams['lon_0']
+ else:
+ lon_0=0.
if self.celestial and not inverse:
try:
- x = -x
+ x = 2.*lon_0-x
except TypeError:
- x = [-xx for xx in x]
+ x = [2.*lon_0-xx for xx in x]
xout,yout = self.projtran(x,y,inverse=inverse)
if self.celestial and inverse:
try:
- xout = -xout
+ xout = 2.*lon_0-xout
except:
- xout = [-xx for xx in xout]
+ xout = [2.*lon_0-xx for xx in xout]
return xout,yout
-
def makegrid(self,nx,ny,returnxy=False):
"""
return arrays of shape (ny,nx) containing lon,lat coordinates of
Something went wrong with that request. Please try again.