Как преобразовать контурный график (matplotlib) в формат FITS с заголовком?

Мне нужно сделать контурный график и наложить контуры на изображение. Я использовал библиотеку aplpy, чтобы наложить контуры на астрономическое изображение. Я скачал данные 2MASS на веб-сайте APlpy(https://github.com/aplpy/aplpy-examples/tree/master/data) и написал следующий фрагмент кода:

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

Кажется, он до сих пор вообще не работает.


person Dalek    schedule 07.02.2014    source источник
comment
О создании контурных графиков с помощью matplotlib у вас есть множество примеров в StackOverflow. Теперь, если вы хотите сохранить график в виде файла FITS, вам нужно сначала преобразовать график в массив numpy (например, следуя описанному процессу здесь), а затем используйте PyFITS для создания файла FITS с массивом.   -  person Ricardo Cárdenes    schedule 07.02.2014
comment
@RicardoCárdenes Как я могу передать информацию из координат прямого восхождения и склонения, когда я конвертирую контурный график matplotlib в массив?   -  person Dalek    schedule 08.02.2014
comment
@Dalek - измени свой цикл, чтобы он соответствовал моему отредактированному. я сделал ошибку; контуры, которые вы создали, уже находятся в координатах ra/dec, поэтому вам не нужно выполнять преобразование пикселей в мир. Обратите внимание, что связанный пример (github.com/aplpy/aplpy-examples/pull/1) имеет правильную версию.   -  person keflavich    schedule 12.02.2014


Ответы (1)


Вы не можете «сохранить контуры как файл FITS» напрямую, но есть и другие подходы, которые вы можете попробовать.

Вы можете использовать инструмент matplotlib._cntr, как описано здесь: Python: найти контурные линии из matplotlib.pyplot.contour(), чтобы получить конечные точки в координатах фигуры, затем используйте WCS для преобразования между пикселями и мировыми координатами. aplpy.FITSFigure имеют вспомогательные функции, F.world2pixel и F.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, создав собственный заголовок.

person keflavich    schedule 07.02.2014
comment
не могли бы вы подробнее рассказать, как я могу использовать WCS (это из астропии), чтобы преобразовать данные с координатной сеткой на контурном графике в WCS? - person Dalek; 08.02.2014
comment
Далек: добавлен подробный ответ - person keflavich; 08.02.2014
comment
для меня все еще очень расплывчато, как я должен это реализовать. если я использую c = cntr.Cntr(Xgrid, Ygrid, Hsmooth), то c.get_cdata() дает мне массив (nx,ny), но это количество каждой точки сетки, а с другой стороны c.trace(z) список точек контура, которые находятся в координатах ra и dec. Какой из них следует преобразовать в файл fits, и не могли бы вы немного объяснить, какой из них можно прочитать с помощью show_contour() в aplpy.FITSFigure()? Спасибо - person Dalek; 08.02.2014
comment
Вы не можете использовать show_contour для этой цели, вместо этого вам нужно использовать show_lines. Я обновил пример, чтобы он был явным для случая, который вы указали в вопросе. - person keflavich; 09.02.2014
comment
Я получил это сообщение об ошибке, когда запустил последний цикл, написанный в вашем ответе: ОШИБКА: TypeError: pixel2world() принимает ровно 3 аргумента (2 задано) [IPython.core.interactiveshell] Разве pixel2world() не должен получать только x и твои координаты? - person Dalek; 09.02.2014
comment
Да, но это также занимает self. Должно быть gc.pixel2world(*verts.T) - person keflavich; 10.02.2014
comment
Спасибо! Я сталкиваюсь с новой проблемой, работающей построчно. Извините, что публикую слишком много. Я получил это сообщение об ошибке для последней строки кода: ----> 3 img.show_lines(radec) 1256 for line in line_list: -> 1257 xp, yp = wcs_util.world2pix(self._wcs, line[0, :], line[1, :]) 1258 lines.append(np.column_stack((xp, yp))) 1259 TypeError: list indices must be integers, not tuple - person Dalek; 10.02.2014
comment
Это очень расстраивает, потому что я не получаю никаких сообщений об ошибках, но в конце он не рисует контуры! - person Dalek; 10.02.2014
comment
Быстрая проверка работоспособности: есть ли все эти координаты ra,dec на вашем изображении? Если это так, попробуйте поиграть с ключевым словом zorder: например, установите zorder=100. - person keflavich; 11.02.2014
comment
Он все еще не показывает контуры! Есть ли какая-либо другая проверка работоспособности, чтобы убедиться, что функция получила правильные входные данные и может построить график? Или вы можете опубликовать пример здесь или на веб-странице приложения? - person Dalek; 11.02.2014
comment
Далек: Похоже, это ошибка! Добавьте color='r', и они должны появиться. См. этот новый пример: github.com/aplpy/aplpy-examples/pull/1 - person keflavich; 11.02.2014
comment
Я изменил свой вопрос, но собрал кусочки из ваших ответов и снова опубликовал как исходный вопрос! Не могли бы вы взглянуть и уточнить для меня, где я сделал ошибки, чтобы сделать контурные графики? - person Dalek; 12.02.2014
comment
Я пытаюсь использовать свой старый код, немного изменить его и добавить что-то похожее на clabel(clines[3:], inline=1, fontsize=20., fmt=ticker.LogFormatterMathtext()) в matplotlib для aplpy, но это не работает, вы не представляете, как я могу заставить его работать! Заранее спасибо~ - person Dalek; 10.10.2014