使用投影时出现问题“类型错误：不支持更改地理投影的轴限制。”

发布于 2025-01-16 14:49:14 字数 490 浏览 1 评论 0原文

我正在尝试使用投影 aitoff 绘制两个带有坐标的数组（大约 480 万个元素），例如：

fig = plt.figure(figsize=(12,10))
ax = fig.add_subplot(111, projection="aitoff")
ax.grid(alpha=0.3)
ax.hist2d(ra, dec, bins = 100)

但出现以下错误：

TypeError: Changing axes limits of a geographic projection is not supported.  Please consider using Cartopy.

坐标（ra 和 dec) 已经是弧度了。我还尝试了以下 hexbin：

ax.hexbin(ra, dec, bins = 100)

，效果很好。你能帮忙吗？

原文

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?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

那支青花 2025-01-23 14:49:14

我找到了答案。这是一个使用度数而不是弧度手动投影所有点的函数：

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

然后

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)

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)

回复收藏 0 原文

~没有更多了~