Я пытался использовать пакет pv-extractor для извлечения скорости по определенному пути.затем преобразуйте положение в координате пикселя в мировую координату, а частоту в скорость.вот код, который я написал:
###############################
## import packages
###############################
import matplotlib
matplotlib.use('Qt4Agg')
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy import wcs
from astropy import coordinates
from astropy import units as u
import os
import aplpy
import numpy as np
from pvextractor import extract_pv_slice
from pvextractor.geometry import Path
from astropy import constants as const
import array as arr
###############################
fitsfile = fits.open('Image.fits', cache=True)
hdu = fitsfile[0]
data = hdu.data[0,:,:,:]
header = hdu.header
w = wcs.WCS(hdu.header)
###############################
## pv extract
###############################
from pvextractor import Path
list = [(1120.40,780.97), (881.37,741.97)]
path1 = Path(list, width=20 )
pv1 = extract_pv_slice(data, path1)
fig = plt.figure( )
F1 = aplpy.FITSFigure(pv1, figure=fig)
F1.show_colorscale(cmap='YlGnBu')
plt.show()
##################################
## define functions
##################################
restf = 8.6754290e10
f0 = header['CRVAL3']
f0 = np.float(f0)
df = header['CDELT3']
df = np.float(df)
ref = header['CRPIX3']
ref = np.float(ref)
chan_max = header['NAXIS3']
chan_max = np.float(chan_max)
channels = np.arange(0, chan_max) # generate array with number of channels
##################################
## velocity
##################################
f0=np.float(f0)
df=np.float(df)
freq = []
freq = f0 + ( ( channels - ref ) * df )
v = []
for i in xrange(0 , len(freq)):
v.append(const.c.value * (1 - (freq[i] / restf)))
##################################
## distance
##################################
pixcrd = np.array(list, np.float)
world = w.wcs_pix2world(pixcrd,1)
ListDist =[]
for i in renge (0, len(list)):
Pixel_x_1 =world[i][0]
Pixel_y_1 =world[i][1]
print Pixel_x_1
print Pixel_y_1
try:
Pixel_x_2 =world[i+1][0]
Pixel_y_2 =world[i+1][1]
except IndexError:
break
dist_dc = abs(Pixel_y_2 - Pixel_y_1)
dist_ra = abs(Pixel_x_2 - Pixel_x_1) * np.cos(())
ListDist.append(dist_connect)
print ListDist
dist_total = sum(ListDist) #it should be in arcsec
print dist_total
###############################
Он отлично извлекает путь и показывает карту, также я мог бы также извлечь скорость, но на последнем шаге я получаю следующую ошибку, которую я получаюне знаю как это решить:
>>Traceback (most recent call last):
>> File "distance.py", line 43, in <module>
>> world = w.wcs_pix2world(pixcrd,1)
>> File "/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.py", line 1355, in wcs_pix2world
>> 'output', *args, **kwargs)
>> File "/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.py", line 1257, in _array_converter
>> return _return_single_array(xy, origin)
>> File "/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.py", line 1238, in _return_single_array
>> "of shape (N, {0})".format(self.naxis))
>>ValueError: When providing two arguments, the array must be of shape (N, 4)