Numpy值错误：使用序列设置数组元素

Question

我正在使用numpy和astropy用Python编写代码。在代码中，我希望创建与我的数据集大体相似的随机numpy数组。之后，我想把这些球坐标的随机数组转换成笛卡尔坐标。不幸的是，我不断得到一个值错误，我完全困惑于为什么会发生这种情况，我曾试图做一些虚拟测试，如如果他们是相同的形状，所有数组的值是合理的，相同的类型，等等，但我卡住了。这是我的密码：

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

R   = 445 + np.random.randn(262615) 
print(np.shape(R))
dec = 2 + np.random.randn(262615)
print(np.shape(dec))
ra  = 150 + np.random.randn(262615)
print(np.shape(ra))
c   = np.zeros(262615) 
print(np.shape(c))


for i in range(262615):
    c[i] = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc)
    print(c[i])      

这是我的错误信息：

PS C:\Users\sirep\Documents\C++ scripts> cd 'c:\Users\sirep\Documents\C++ scripts'; ${env:PYTHONIOENCODING}='UTF-8'; ${env:PYTHONUNBUFFERED}='1'; & 'C:\Users\sirep\Anaconda3\python.exe' 'c:\Users\sirep\.vscode\extensions\ms-python.python-2018.5.0\pythonFiles\PythonTools\visualstudio_py_launcher.py' 'c:\Users\sirep\Documents\C++ scripts' '57764' '34806ad9-833a-4524-8cd6-18ca4aa74f14' 'RedirectOutput,RedirectOutput' 'c:\Users\sirep\Documents\Python Scripts\sph2cart.py'
(262615,)
(262615,)
(262615,)
(262615,)
Traceback (most recent call last):
  File "c:\Users\sirep\Documents\Python Scripts\sph2cart.py", line 16, in <module>
    c[i] = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc)
ValueError: setting an array element with a sequence.

谢谢大家抽出时间！你知道吗

Answer 1

我发现了我的错误。SkyCoord返回3个值：x坐标、y坐标和z坐标。我试着给数组的一个元素赋三个值。为了解决这个问题，我必须首先为每个坐标创建3个独立的数组，然后确保每个值在输入到各自的数组时是无量纲的：

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


R   = 445 + np.random.randn(262615) 
print(np.shape(R))
dec = 2 + np.random.randn(262615)
print(np.shape(dec))
ra  = 150 + np.random.randn(262615)
print(np.shape(ra))
cx   = np.zeros(262615) 
cy   = np.zeros(262615)
cz   = np.zeros(262615)
print(np.shape(cx))


for i in range(262):
    c = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc)
    cx[i] = c.cartesian.x.value
    cy[i] = c.cartesian.y.value
    cz[i] = c.cartesian.z.value
    print(cx[i],cy[i],cz[i])

Answer 2

我想我应该把我的comment扩展成一个较长的答案，因为这里有一些事情值得解释和澄清，以供将来的读者阅读。你知道吗

你在your answer中写道：

I've figured out my error. SkyCoord returns 3 values an x coordinate, y coordinate and z coordinate. I was trying to assign three values to a single element of an array.

这当然是正确的，但并不完全正确。在原始代码中，您有如下内容：

c = np.zeros(262615)

这已经有点麻烦了，因为您没有指定数据类型，但是默认情况下，数据类型是float64，这可能是许多应用程序所需要的（对于这个应用程序来说确实如此）。在任何情况下，Numpy数组类型化意味着，如果您像在原始代码中那样为数组的单个元素赋值：

c[i] = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc)

指定的值最好是一个浮点数，或者至少是其他一些可以明确转换为float的数值类型（比如int）。对SkyCoord来说不是这样的，因为正如你所说的，它是一个三维的多重集。我的观点是，一般来说，如果您使用Numpy数组，您需要注意它的dtype是什么，以及您试图分配给它的元素什么。对于更随意的对象，您可能会得到更清晰的错误，如：

>>> c[0] = object()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: float() argument must be a string or a number

仍然不是很好，但至少它证明了它正试图调用float()将参数转换为float。但是SkyCoord会得到不同的结果，因为SkyCoord可以是一个包含许多坐标的数组的容器，Numpy会看到这一点，并尝试将其视为将一系列值赋给一个标量，这就是所得到的错误。你知道吗

顺便说一句，在Numpy中也可以使用structured arrays创建更复杂的数组类型。这允许您创建一个(x, y, z)坐标数组，例如：

>>> c = np.zeros(262615, dtype=[('x', 'f8'), ('y', 'f8'), ('z', 'f8')])
>>> c
array([(0.0, 0.0, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0), ...,
       (0.0, 0.0, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0)],
      dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])
>>> c[0]
(0.0, 0.0, 0.0)

尽管您不能将SkyCoord直接赋给这些值中的一个（我认为从技术上讲SkyCoord被视为无坐标，不管您使用什么坐标系来实例化它，但我可能错了），但是您可以指定例如：

>>> c[0] = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc).cartesian.xyz

但是，这仍然是没有必要的，因为正如我在评论中提到并进一步解释的in the docsa SkyCoord可以表示一个坐标数组，如：

>>> coords = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, distance=R*u.mpc)

你可以一次把这些转换成笛卡尔坐标，然后检索x，y，z坐标的独立数组，比如：

>>> x, y, z = coords.cartesian.xyz

这还有一个额外的优点，即使用最合适的长度维度将坐标返回为Quantity（在本例中是Mpc，因为这是您给出的距离）。然而，coords.cartesian本身实际上已经是一个(x, y, z)坐标的数组，非常像我上面的结构化数组示例（从技术上讲，它不是Numpy数组，但它有许多相同的方法，可以转换为类似的方法）：

>>> coords.cartesian._values
array([(0.19718680211339326, 0.002173755110841713, 0.0021735811221131776),
       (0.6853033697941637, 0.005924402286034272, 0.004262079913938389),
       ...
      dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])

但是这是一个不应该使用的未记录的内部属性（尽管我不知道为什么这个接口没有公开，因为它可能有用…）

最后，我要补充一点，使用这个接口要快得多，因为所有的循环都是向量化的数组操作，大部分都是用C语言编写的。任何时候在Python级别执行诸如分配给数组（c[i] = ...）或属性访问（c.cartesian.x.value）之类的操作时，都会导致严重的性能损失，因为需要将值从C转换为Python，然后再转换回C。使用矢量化操作可以避免所有这些。所以当我制作一个SkyCoord数组时，我得到：

In [7]: %%timeit
   ...: c = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, distance=R*u.mpc)
   ...: c.cartesian.xyz
   ...:
10 loops, best of 3: 111 ms per loop

或111ms表示262615坐标，如原始示例中所示。然而，用“天真”的方式去做会让我：

In [11]: %%timeit
    ...: for i in range(262615):
    ...:     c = SkyCoord(ra=ra[i]*u.degree,dec=dec[i]*u.degree,distance=R[i]*u.mpc)
    ...:     cx[i] = c.cartesian.x.value
    ...:     cy[i] = c.cartesian.y.value
    ...:     cz[i] = c.cartesian.z.value
    ...:

1 loop, best of 3: 18min 26s per loop