Я думаю, что есть несколько способов сделать это, когда вы усложняете себе задачу, не эффективно используя существующие доступные инструменты. Например, поскольку вы работаете с табличными данными из файла FITS, вы можете воспользоваться интерфейсом Astropy Table :
>>> from astropy.table import Table
>>> primus = Table.read('PRIMUS_2013_zcat_v1.fits')
(для этого конкретного файла я получил несколько предупреждений о некоторых заголовки в таблице нестандартны, но это можно игнорировать).
Если вы хотите выполнить некоторые операции над несколькими столбцами таблицы, вы можете легко это сделать. Например, вместо того, чтобы делать то, что вы сделали, выделив несколько столбцов вместе, а затем сложив их в новый массив
np.column_stack((data_PRIMUS['RA'], data_PRIMUS['DEC'],
data_PRIMUS['Z'], data_PRIMUS['FIELD']))
, вы можете выбрать подмножество столбцов из таблицы следующим образом:
>>> primus[['RA', 'DEC', 'Z', 'FIELD']]
<Table length=213696>
RA DEC Z FIELD
degree degree
float64 float64 float32 bytes13
------------------ ------------------- ---------- -------------
52.892275339281994 -27.833172368069615 0.3420992 calib
52.88448889270391 -27.85252305560996 0.4824943 calib
52.880363885710295 -27.86221750021335 0.33976158 calib
52.88334306466262 -27.86937808271639 0.6134631 calib
52.8866138857103 -27.871773055662942 0.58744365 calib
52.885607068267845 -27.889578785511922 0.26873255 calib
... ... ... ...
34.54856 -4.5544 0.8544105 xmm
34.56942 -4.57564 0.6331108 xmm
34.567412432719756 -4.572718190305209 1.1456184 xmm
34.57134 -4.56414 0.6346616 xmm
34.58088 -4.56804 1.081143 xmm
34.58686 -4.57449 0.7471819 xmm
Затем кажется, что вы выбираете столбцы RA
, DEC
и Z
, где поле равно xmm
, с помощью функции filter
, но, поскольку это массивы Numpy, вы можете использовать возможности фильтрации, встроенные в индексирование массива Numpy, а также индексирование таблиц. Единственная сложность заключается в том, что, поскольку это поля символов фиксированной ширины, вам все равно необходимо правильно выполнять сравнения. Для этого вы можете использовать строковые функции Numpy, например np.char.startswith
:
>>> primus = primus[np.char.startswith(primus['FIELD'], b'xmm')]
В процессе сравнения производительности я понял, что в этой строке вы вероятно, получая ошибку IndexError: too many indices for array
:
>>> np.array(filter(lambda x: x[3].strip() == 'xmm', primus))
array(<filter object at 0x7f5170981940>, dtype=object)
В Python 3 функция filter
возвращает итерацию, поэтому ее перенос в np.array()
просто создает массив 0-D, содержащий этот Python объект; это, вероятно, не то, что вы намеревались, поэтому он терпит неудачу здесь (именно здесь можно было бы посмотреть на трассировку). Даже если вы поместите вызов filter()
в list()
, он не будет работать, потому что np.array()
обычно принимает только однородные массивы. Так что подход, подобный тому, который я дал, вполне достаточен (хотя могут быть и несколько более эффективные способы). Это также делает следующую строку:
np.array(map(lambda x: [float(x[0]), float(x[1]), float(x[2])], data_PRIMUS))
ненужной. В частности, первые три столбца уже представлены в формате с плавающей запятой, поэтому в любом случае в этом нет необходимости.
Некоторые аналогичные рекомендации применимы и к другим частям вашего кода. Я бы написал так примерно так:
import numpy as np
from astropy.table import Table, vstack
from astropy import units as u
from astropy.coordinates import SkyCoord
primus = Table.read('PRIMUS_2013_zcat_v1.fits')
primus_field = primus['FIELD']
primus = primus[['RA', 'DEC', 'Z']]
primus = primus[np.char.startswith(primus_field, b'xmm')]
vipers = Table.read('VIPERS_W1_SPECTRO_PDR2.fits')[['alpha', 'delta', 'zspec']]
primus_catalog = SkyCoord(ra=primus['RA']*u.degree, dec=primus['DEC']*u.degree)
vipers_catalog = SkyCoord(ra=vipers['alpha']*u.degree, dec=vipers['delta']*u.degree)
idx, d2d, d3d = primus_catalog.match_to_catalog_sky(vipers_catalog)
feasible_indices = idx[d2d > 1e-3]
vipers = vipers[feasible_indices]
vipers.rename_columns(['alpha', 'delta', 'zspec'], ['RA', 'DEC', 'Z'])
hz = vstack(primus, vipers)
Пожалуйста, дайте мне знать, если есть какие-то части этого вопроса, у вас есть вопросы.