NDData算法

介绍

NDDataRef 实现以下算术运算:

使用基本的算术方法

使用标准算术方法要求第一个操作数是 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.])
>>> ndd.unit
Unit("cm m")

或者取一个 NDDataRef 对象:

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

可能的操作数

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

  • 任何类型的标量

  • 包含数字的列表(或嵌套列表)

  • numpy 数组

  • numpy 屏蔽数组

  • astropy

  • 其他 nddata 类或子类

高级选项

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

数据和单位

为了 dataunit 没有参数。每个算术运算都可以 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 arithmetics 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 否则将禁止操作。可能的值与 maskmeta

  • 如果 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_uncertaintiesTrue 你也可以为 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])

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