NDData算法

介绍

NDDataRef 实现以下算术运算：

• 附加： add()

• 减法： subtract()

• 乘法： multiply()

• 师： divide()

使用基本的算术方法

使用标准算术方法要求第一个操作数是 NDDataRef 实例：

>>> from astropy.nddata import NDDataRef
>>> from astropy.wcs import WCS
>>> import numpy as np
>>> ndd1 = NDDataRef([1, 2, 3, 4])

而对第二个操作数的要求是它必须可以转换为第一个操作数。它可以是一个数字：

>>> ndd1.add(3)
NDDataRef([4, 5, 6, 7])

或者 list ：：

>>> ndd1.subtract([1,1,1,1])
NDDataRef([0, 1, 2, 3])

或者 numpy.ndarray ：：

>>> ndd1.multiply(np.arange(4, 8))
NDDataRef([ 4, 10, 18, 28])
>>> ndd1.divide(np.arange(1,13).reshape(3,4))  # a 3 x 4 numpy array  
NDDataRef([[1.        , 1.        , 1.        , 1.        ],
           [0.2       , 0.33333333, 0.42857143, 0.5       ],
           [0.11111111, 0.2       , 0.27272727, 0.33333333]])

在这里，广播照顾到了不同的维度。还可以接受其他几种类型的操作数。

使用算术类方法

这里两个操作数不必都是 NDDataRef -比如：

>>> NDDataRef.add(1, 3)
NDDataRef(4)

在两个量之间包装算术运算的结果：

>>> import astropy.units as u
>>> ndd = NDDataRef.multiply([1,2] * u.m, [10, 20] * u.cm)
>>> ndd  
NDDataRef([10., 40.], unit='cm m')
>>> ndd.unit
Unit("cm m")

或者取一个 NDDataRef 对象：

>>> NDDataRef.divide(1, ndd1)  
NDDataRef([1.        , 0.5       , 0.33333333, 0.25      ])

可能的操作数

操作数的可能输入类型有：

• 任何类型的标量

• 包含数字的列表（或嵌套列表）

• numpy 数组

• numpy 屏蔽数组

• astropy 量

• 其他 nddata 类或子类

高级选项

普通的Python操作符 + ， - ，等没有实现，因为这些方法提供了几个关于如何处理附加属性的选项。

数据和单位

为了 data 和 unit 没有参数。每个算术运算都可以 astropy.units.Quantity -框架评估结果或失败并中止操作。

加二 NDData 单位工程相同的对象：

>>> ndd1 = NDDataRef([1,2,3,4,5], unit='m')
>>> ndd2 = NDDataRef([100,150,200,50,500], unit='m')

>>> ndd = ndd1.add(ndd2)
>>> ndd.data  
array([101., 152., 203., 54., 505.])
>>> ndd.unit
Unit("m")

加二 NDData 具有兼容单位的对象也可以工作：

>>> ndd1 = NDDataRef(ndd1, unit='pc')
INFO: overwriting NDData's current unit with specified unit. [astropy.nddata.nddata]
>>> ndd2 = NDDataRef(ndd2, unit='lyr')
INFO: overwriting NDData's current unit with specified unit. [astropy.nddata.nddata]

>>> ndd = ndd1.subtract(ndd2)
>>> ndd.data  
array([ -29.66013938,  -43.99020907,  -58.32027876,  -11.33006969,
       -148.30069689])
>>> ndd.unit
Unit("pc")

这将默认保留第一个操作数的单位。但是，在划分过程中不会分解单元：

>>> ndd = ndd2.divide(ndd1)
>>> ndd.data  
array([100. , 75. , 66.66666667, 12.5 , 100. ])
>>> ndd.unit
Unit("lyr / pc")

面具

这个 handle_mask 用于算术运算的参数实现结果掩码的内容。有几种选择。

• None ，结果将没有 mask ：：

  >>> ndd1 = NDDataRef(1, mask=True)
  >>> ndd2 = NDDataRef(1, mask=False)
  >>> ndd1.add(ndd2, handle_mask=None).mask is None
  True
  

• "first_found" 或 "ff" ，结果将具有 mask 第一个操作数的 None , the mask 第二个操作数：

  >>> ndd1 = NDDataRef(1, mask=True)
  >>> ndd2 = NDDataRef(1, mask=False)
  >>> ndd1.add(ndd2, handle_mask="first_found").mask
  True
  >>> ndd3 = NDDataRef(1)
  >>> ndd3.add(ndd2, handle_mask="first_found").mask
  False
  

• 至少有两个参数的函数（或任意可调用函数）。例如， numpy.logical_or 是默认值：

  >>> ndd1 = NDDataRef(1, mask=np.array([True, False, True, False]))
  >>> ndd2 = NDDataRef(1, mask=np.array([True, False, False, True]))
  >>> ndd1.add(ndd2).mask
  array([ True, False,  True,  True]...)
  

  默认为 "first_found" 万一只有一个 mask 不是没有：

  >>> ndd1 = NDDataRef(1)
  >>> ndd2 = NDDataRef(1, mask=np.array([True, False, False, True]))
  >>> ndd1.add(ndd2).mask
  array([ True, False, False,  True]...)
  

  也可以自定义函数：

  >>> def take_alternating_values(mask1, mask2, start=0):
  ...     result = np.zeros(mask1.shape, dtype=np.bool_)
  ...     result[start::2] = mask1[start::2]
  ...     result[start+1::2] = mask2[start+1::2]
  ...     return result
  

  这个函数是无稽之谈，但我们仍然可以看到它是如何执行的：

  >>> ndd1 = NDDataRef(1, mask=np.array([True, False, True, False]))
  >>> ndd2 = NDDataRef(1, mask=np.array([True, False, False, True]))
  >>> ndd1.add(ndd2, handle_mask=take_alternating_values).mask
  array([ True, False,  True,  True]...)
  

  附加参数可以通过在它们前面加上前缀来给出 mask_ （在传递给函数之前将被剥离）：

  >>> ndd1.add(ndd2, handle_mask=take_alternating_values, mask_start=1).mask
  array([False, False, False, False]...)
  >>> ndd1.add(ndd2, handle_mask=take_alternating_values, mask_start=2).mask
  array([False, False,  True,  True]...)
  

元

这个 handle_meta 用于算术运算的参数实现结果 meta 将。选项与 mask ：

• 如果 None 结果 meta 会是空的 collections.OrderedDict .

  >>> ndd1 = NDDataRef(1, meta={'object': 'sun'})
  >>> ndd2 = NDDataRef(1, meta={'object': 'moon'})
  >>> ndd1.add(ndd2, handle_meta=None).meta
  OrderedDict()
  

  为了 meta 这是默认值，因此在这种情况下不需要传递它：：

  >>> ndd1.add(ndd2).meta
  OrderedDict()
  

• 如果 "first_found" 或 "ff" ，结果 meta 将是 meta 如果第一个操作数不包含键，则 meta 取第二个操作数的。

  >>> ndd1 = NDDataRef(1, meta={'object': 'sun'})
  >>> ndd2 = NDDataRef(1, meta={'object': 'moon'})
  >>> ndd1.add(ndd2, handle_meta='ff').meta
  {'object': 'sun'}
  

• 如果是 callable 它至少需要两个参数。两者兼而有之 meta 属性将传递给此函数（即使其中一个或两个都为空），并且可调用函数计算结果的 meta . 例如，合并这两个函数的函数：

  >>> # It's expected with arithmetic that the result is not a reference,
  >>> # so we need to copy
  >>> from copy import deepcopy
  
  >>> def combine_meta(meta1, meta2):
  ...     if not meta1:
  ...         return deepcopy(meta2)
  ...     elif not meta2:
  ...         return deepcopy(meta1)
  ...     else:
  ...         meta_final = deepcopy(meta1)
  ...         meta_final.update(meta2)
  ...         return meta_final
  
  >>> ndd1 = NDDataRef(1, meta={'time': 'today'})
  >>> ndd2 = NDDataRef(1, meta={'object': 'moon'})
  >>> ndd1.subtract(ndd2, handle_meta=combine_meta).meta 
  {'object': 'moon', 'time': 'today'}
  

  这里同样可以使用前缀传递函数的其他参数 meta_ （在传递给此函数之前将被剥离）。有关详细信息，请参见mask属性的说明。

世界坐标系（WCS）

这个 compare_wcs 参数将决定结果是什么 wcs 否则将禁止操作。可能的值与 mask 和 meta ：

• 如果 None 结果 wcs 会是空的 None .

  >>> ndd1 = NDDataRef(1, wcs=None)
  >>> ndd2 = NDDataRef(1, wcs=WCS())
  >>> ndd1.add(ndd2, compare_wcs=None).wcs is None
  True
  

• 如果 "first_found" 或 "ff" 结果 wcs 将是 wcs 第一个操作数的 None , the meta 取第二个操作数的。

  >>> wcs = WCS()
  >>> ndd1 = NDDataRef(1, wcs=wcs)
  >>> ndd2 = NDDataRef(1, wcs=None)
  >>> str(ndd1.add(ndd2, compare_wcs='ff').wcs) == str(wcs)
  True
  

• 如果是 callable 它至少需要两个参数。两者兼而有之 wcs 属性将传递给此函数（即使其中一个或两个 None )回电话的人应该回来了 True 如果这些 wcs 相同（足够）以允许算术运算或 False 如果算术运算应该用 ValueError .如果 True 这个 wcs 相同，第一个用于结果：

  >>> def compare_wcs_scalar(wcs1, wcs2, allowed_deviation=0.1):
  ...     if wcs1 is None and wcs2 is None:
  ...         return True  # both have no WCS so they are identical
  ...     if wcs1 is None or wcs2 is None:
  ...         return False  # one has WCS, the other doesn't not possible
  ...     else:
  ...         # Consider wcs close if centers are close enough
  ...         return all(abs(wcs1.wcs.crpix - wcs2.wcs.crpix) < allowed_deviation)
  
  >>> ndd1 = NDDataRef(1, wcs=None)
  >>> ndd2 = NDDataRef(1, wcs=None)
  >>> ndd1.subtract(ndd2, compare_wcs=compare_wcs_scalar).wcs
  

  附加参数可以在前缀中传递 wcs_ （此前缀在传递给函数之前将被删除）：：

  >>> ndd1 = NDDataRef(1, wcs=WCS())
  >>> ndd1.wcs.wcs.crpix = [1, 1]
  >>> ndd2 = NDDataRef(1, wcs=WCS())
  >>> ndd1.subtract(ndd2, compare_wcs=compare_wcs_scalar, wcs_allowed_deviation=2).wcs.wcs.crpix
  array([1., 1.])
  

  如果您正在使用 WCS 对象，一个非常方便的函数可能是：

  >>> def wcs_compare(wcs1, wcs2, *args, **kwargs):
  ...     return wcs1.wcs.compare(wcs2.wcs, *args, **kwargs)
  

  见 astropy.wcs.Wcsprm.compare() 对于这个比较允许的参数。

不确定性

这个 propagate_uncertainties 参数可用于打开或关闭不确定性的传播。

• 如果 None 结果没有不确定性：

  >>> from astropy.nddata import StdDevUncertainty
  >>> ndd1 = NDDataRef(1, uncertainty=StdDevUncertainty(0))
  >>> ndd2 = NDDataRef(1, uncertainty=StdDevUncertainty(1))
  >>> ndd1.add(ndd2, propagate_uncertainties=None).uncertainty is None
  True
  

• 如果 False 结果将有第一次发现的不确定性。

  备注

  设置 propagate_uncertainties=False 一般不推荐。

• 如果 True 两种不确定性都必须 NDUncertainty 实现传播的子类。这是可能的 StdDevUncertainty ：：

  >>> ndd1 = NDDataRef(1, uncertainty=StdDevUncertainty([10]))
  >>> ndd2 = NDDataRef(1, uncertainty=StdDevUncertainty([10]))
  >>> ndd1.add(ndd2, propagate_uncertainties=True).uncertainty  
  StdDevUncertainty([14.14213562])
  

相关不确定度

如果 propagate_uncertainties 是 True 你也可以为 uncertainty_correlation . StdDevUncertainty 它本身无法跟踪其相关性，但它可以评估正确的结果不确定性，如果正确的话 correlation 给出。

默认值 (0 )表示不相关的while 1 均值相关和 -1 反相关。如果给予 numpy.ndarray 它应该代表元素的相关系数。

实例

无相关性，减去 NDDataRef 实例本身导致非零不确定性：

>>> ndd1 = NDDataRef(1, uncertainty=StdDevUncertainty([10]))
>>> ndd1.subtract(ndd1, propagate_uncertainties=True).uncertainty  
StdDevUncertainty([14.14213562])

给出了 1 （因为它们明显相关）给出了正确的不确定性 0 ：：

>>> ndd1 = NDDataRef(1, uncertainty=StdDevUncertainty([10]))
>>> ndd1.subtract(ndd1, propagate_uncertainties=True,
...               uncertainty_correlation=1).uncertainty  
StdDevUncertainty([0.])

与等效操作一致 ndd1 * 0 ：：

>>> ndd1.multiply(0, propagate_uncertainties=True).uncertainty 
StdDevUncertainty([0.])

警告

用户需要手动计算或知道适当的值或数组，并将其传递给 uncertainty_correlation . 实现遵循一般的一阶误差传播公式。例如，参见： Wikipedia .

你也可以给出元素相关性：

>>> ndd1 = NDDataRef([1,1,1,1], uncertainty=StdDevUncertainty([1,1,1,1]))
>>> ndd2 = NDDataRef([2,2,2,2], uncertainty=StdDevUncertainty([2,2,2,2]))
>>> ndd1.add(ndd2,uncertainty_correlation=np.array([1,0.5,0,-1])).uncertainty  
StdDevUncertainty([3.        , 2.64575131, 2.23606798, 1.        ])

相关性 np.array([1, 0.5, 0, -1]) 表明第一个元素是完全相关的，第二个元素是部分相关的，而第三个元素是不相关的，第四个元素是反相关的。

单位不确定性

StdDevUncertainty 即使数据单位与不确定度单位不同，也能实现正确的误差传播：

>>> ndd1 = NDDataRef([10], unit='m', uncertainty=StdDevUncertainty([10], unit='cm'))
>>> ndd2 = NDDataRef([20], unit='m', uncertainty=StdDevUncertainty([10]))
>>> ndd1.subtract(ndd2, propagate_uncertainties=True).uncertainty  
StdDevUncertainty([10.00049999])

但它需要转换成数据的单位。