Issue
I am trying to plot two arrays with coordinates (around 4.8 millions elements) using the projection aitoff
like:
fig = plt.figure(figsize=(12,10))
ax = fig.add_subplot(111, projection="aitoff")
ax.grid(alpha=0.3)
ax.hist2d(ra, dec, bins = 100)
But I get the following error:
TypeError: Changing axes limits of a geographic projection is not supported. Please consider using Cartopy.
The coordinates (ra
and dec
) are already in radians.
I also tried with the following hexbin:
ax.hexbin(ra, dec, bins = 100)
and that works. Can you please help?
Solution
I found an answer. Here is a function to manually project all the points using degrees instead of radians:
degrad = np.pi/180.
def project(li,bi,lz):
sa = li-lz
if len(sa) == 1:
sa = np.zeros(1)+sa
x180 = np.where(sa >= 180.0)
sa = sa
sa[x180] = sa[x180]-360.0*abs(np.cos(lz*degrad/2.))#uncomment b=0
alpha2 = sa*degrad/2.
delta = bi*degrad
r2 = np.sqrt(2.)
f = 2.*r2/np.pi
cdec = np.cos(delta)
denom = np.sqrt(1.+cdec*np.cos(alpha2))
xx = cdec*np.sin(alpha2)*2.*r2/denom
yy = np.sin(delta)*r2/denom
xx = xx*(1./degrad)/f
yy = yy*(1./degrad)/f
return xx,yy
And then
ra, dec = project(np.array(ra), np.array(dec), 180)
binner1 = np.linspace(-180, 180, 1000)
binner2 = np.linspace(-90, 90, 1000)
img, xbins,ybins = np.histogram2d(ra, dec, bins=(binner1,binner2))
fig, ax = plt.subplots(figsize=(10.,6.))
plot = ax.imshow(img.T, origin='lower', extent=[xbins[0],xbins[-1],ybins[0],ybins[-1]],
cmap=plt.cm.viridis, interpolation='gaussian', aspect='auto',
vmax=40,vmin=0)
Answered By - Alessandro Peca
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.