«IndexError: слишком много индексов для массива» при объединении VIPERS и PRIMUS - PullRequest
1 голос
/ 06 марта 2020

Привет. Я пытаюсь извлечь RA, De c и информацию о красном смещении из двух опросов (PRIMUS и VIPERS) и собрать их в один массив nd. Код выглядит следующим образом:

from astropy.io import fits
import numpy as np

hdulist_PRIMUS = fits.open('data/PRIMUS_2013_zcat_v1.fits')

data_PRIMUS = hdulist_PRIMUS[1].data
data_PRIMUS = np.column_stack((data_PRIMUS['RA'], data_PRIMUS['DEC'],
    data_PRIMUS['Z'], data_PRIMUS['FIELD']))
data_PRIMUS = np.array(filter(lambda x: x[3].strip() == 'xmm', data_PRIMUS))[:, :3]
data_PRIMUS = np.array(map(lambda x: [float(x[0]), float(x[1]), float(x[2])], data_PRIMUS))

hdulist_VIPERS = fits.open('data/VIPERS_W1_SPECTRO_PDR2.fits')
data_VIPERS = hdulist_VIPERS[1].data
data_VIPERS = np.column_stack((data_VIPERS['alpha'], data_VIPERS['delta'], data_VIPERS['zspec']))

from astropy import units as u
from astropy.coordinates import SkyCoord

PRIMUS_catalog = SkyCoord(ra=data_PRIMUS[:, 0]*u.degree, dec =data_PRIMUS[:, 1]*u.degree)
VIPERS_catalog = SkyCoord(ra=data_VIPERS[:, 0]*u.degree, dec=data_VIPERS [:, 1]*u.degree)

idx, d2d, d3d = PRIMUS_catalog.match_to_catalog_sky(VIPERS_catalog)
feasible_indices = np.array(map(
    lambda x: x[0],
    filter(lambda x: x[1].value > 1e-3, zip(idx, d2d))))
data_VIPERS = data_VIPERS[feasible_indices]
data_HZ = np.vstack((data_PRIMUS, data_VIPERS))

Когда я запускаю это, я получаю «IndexError: слишком много индексов для массива»

Наборы данных: Каталог PRIMUS Redshift - https://primus.ucsd.edu/version1.html Каталог VIPS Redshift - https://projects.ift.uam-csic.es/skies-universes/VIPERS/photometry/

1 Ответ

1 голос
/ 08 марта 2020

Я думаю, что есть несколько способов сделать это, когда вы усложняете себе задачу, не эффективно используя существующие доступные инструменты. Например, поскольку вы работаете с табличными данными из файла 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)

Пожалуйста, дайте мне знать, если есть какие-то части этого вопроса, у вас есть вопросы.

...