【问题标题】:How could I convert a contour plot (matplotlib) to FITS format with header?如何将等高线图(matplotlib)转换为带有标题的 FITS 格式?
【发布时间】:2014-02-07 19:39:29
【问题描述】:

我需要制作等高线图并覆盖图像上的等高线。我使用 aplpy 库来覆盖天文图像上的轮廓。 我已经在 APlpy 网站(https://github.com/aplpy/aplpy-examples/tree/master/data) 中下载了 2MASS 数据并编写了以下代码:

import numpy as np
import aplpy
import atpy
from pyavm import AVM
import scipy
import matplotlib.pyplot as plt
import scipy.ndimage
from matplotlib import cm
import montage_wrapper
import matplotlib.colors as colors
from scipy import stats

def get_contour_verts(cn):
    contours = []
    # for each contour line
    for cc in cn.collections:
        paths = []
        # for each separate section of the contour line
        for pp in cc.get_paths():
            xy = []
            # for each segment of that section
            for vv in pp.iter_segments():
                xy.append(vv[0])
            paths.append(np.vstack(xy))
        contours.append(paths)
    return contours

# Convert all images to common projection
aplpy.make_rgb_cube(['2MASS_h.fits', '2MASS_j.fits', '2MASS_k.fits'], '2MASS_rgb.fits')

# Customize it
aplpy.make_rgb_image('2MASS_rgb.fits','2MASS_rgb_contour.png',embed_avm_tags=True)

# Launch APLpy figure of 2D cube
img = aplpy.FITSFigure('2MASS_rgb_contour.png')
img.show_rgb()

# Modify the tick labels for precision and format
img.tick_labels.set_xformat('hhmm')
img.tick_labels.set_yformat('ddmm')
# Move the tick labels
img.tick_labels.set_xposition('top')
img.tick_labels.set_yposition('right')
img.tick_labels.set_font(size='small', family='sans-serif', style='normal')

data = scipy.loadtxt('sources.txt')
m1=data[:,0]
m2=data[:,1]
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()
#Gridding the data 

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values)
Z = np.reshape(kernel(positions).T, X.shape)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r, extent=[xmin, xmax, ymin, ymax])
ax.plot(m1, m2, 'k.', markersize=2)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])

CS = plt.contour(X, Y, Z)
Contour_arrays=get_contour_verts(CS)

#Adding contour plots
for contours_at_level in Contour_arrays:
    radec = [img.pixel2world(*verts.T) for verts in contours_at_level]
    new_radec=[]
    for coor in radec:
        #new_radec.append(np.column_stack((coor[0], coor[1])))
        new_radec.append(np.vstack((coor[0], coor[1])))
    print new_radec
    img.show_lines(new_radec,color='red', zorder=100)

img.save('tutorial.png')

它似乎仍然无法正常工作。

【问题讨论】:

  • 关于使用 matplotlib 绘制等高线图,您有很多关于 StackOverflow 的示例。现在,如果要将绘图保存为 FITS 文件,您首先需要将绘图转换为 numpy 数组(例如,按照here 描述的过程),然后使用 PyFITS 创建包含数组的 FITS 文件。
  • @RicardoCárdenes 当我将 matplotlib 等高线图转换为数组时,如何从赤经和赤纬坐标传输信息?
  • @Dalek - 更改您的循环以匹配我编辑的循环。我犯了一个错误;您创建的轮廓在 ra/dec 坐标中已经,因此您不需要进行像素->世界转换。请注意,链接示例(github.com/aplpy/aplpy-examples/pull/1)的版本正确

标签: python matplotlib pyfits astropy aplpy


【解决方案1】:

您不能直接“将轮廓另存为 FITS 文件”,但您可以尝试其他方法。

您可以使用matplotlib._cntr 工具,如下所述:Python: find contour lines from matplotlib.pyplot.contour() 获取图形坐标中的端点,然后使用 WCS 在像素坐标和世界坐标之间进行转换。 aplpy.FITSFigures 有方便的函数,F.world2pixelF.pixel2world,每个接受 2 个数组:

F.pixel2world(arange(5),arange(5))

因此,如果您使用的网格与FITSFigure 窗口中显示的网格相同,您可以将点转换为世界坐标并使用show_lines 绘制它们:

ra,dec = F.pixel2world(xpix,ypix)
F.show_lines([[ra,dec]])

或者对于更真实的轮廓案例,现在复制链接文章中的代码:

import numpy as np
import aplpy

def get_contour_verts(cn):
    contours = []
    # for each contour line
    for cc in cn.collections:
        paths = []
        # for each separate section of the contour line
        for pp in cc.get_paths():
            xy = []
            # for each segment of that section
            for vv in pp.iter_segments():
                xy.append(vv[0])
            paths.append(np.vstack(xy))
        contours.append(paths)

    return contours

# This line is copied from the question's code, and assumes that has been run previously
CS = plt.contour(Xgrid, Ygrid, Hsmooth,levels=loglvl,extent=extent,norm=LogNorm())

contours = get_contour_verts(CS) # use the linked code


gc = aplpy.FITSFigure('2MASS.fits',figsize=(10,9))

# each level comes with a different set of vertices, so you have to loop over them
for contours_at_level in Contour_arrays:
    clines = [cl.T for cl in contours_at_level]
    img.show_lines(clines,color='red', zorder=100)

仍然是最简单的方法是,如果您有一个 FITS 文件,您可以从中生成所述轮廓,则直接使用该文件。或者,如果您没有 FITS 文件,您可以使用 pyfits 通过创建自己的标头来制作。

【讨论】:

  • 您能否详细说明我如何使用 WCS(它来自 astropy)以便将等高线图中的网格数据转换为 WCS?
  • Dalek:添加了详细答案
  • 对我来说仍然很模糊,我应该如何实现它。如果我使用 c = cntr.Cntr(Xgrid, Ygrid, Hsmooth) 然后 c.get_cdata() 给我一个 (nx,ny) 数组,但它是每个网格点的数量,另一方面 c.trace(z)轮廓中位于 ra 和 dec 坐标中的点列表。哪一个应该转换为适合文件,你能解释一下哪一个可以被 aplpy.FITSFigure() 中的 show_contour() 读取吗?谢谢
  • 您不能为此使用show_contour,您需要使用show_lines。我已经更新了示例,以明确您在问题中概述的情况
  • 我在运行您的答案中写的最后一个循环时收到此错误消息:错误:TypeError:pixel2world() 恰好需要 3 个参数(给定 2 个)[IPython.core.interactiveshell] 不是 pixel2world () 应该只得到 x 和 y 坐标?
猜你喜欢
  • 1970-01-01
  • 2015-04-30
  • 2020-03-29
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多