在NumpyArray中解释特定位

Question

我试图解释HDF文件，或者更精确地解释文件中的光栅带。有一个用于质量评估的频带，它将位信息表示为相应的整数（例如，01000000表示为64）

根据具体的位，我想得到一个计数，例如，位位置5上有多少像素是1。如果以前的计数考虑了它们，则不应再次获得计数。现在，我正在根据自己的优先级列表更改每个单元格中的条目。这需要很长时间，我确信有更快的方法，因为我以前从未使用过bits

这是我的密码：

from osgeo import gdal
import numpy as np
QA_Band = gdal.Open(hdf.GetSubDatasets()[B][0],gdal.GA_ReadOnly)
QA = QA_Band.ReadAsArray()

# Calculate Bit-Representation of QA-Band
bin_vec = np.vectorize(np.binary_repr)
QAB = bin_vec(QA, width = 8)

# Filter by Bit-Values
QAB = np.where(QAB == '11111111', "OutOfSwath", QAB)

for i in range(QAB.shape[0]):
    for j in range(QAB.shape[1]):
        if QAB[i,j][6] == '1':
            QAB[i,j] = "Cloud"
            Cloud = (QAB == "Cloud").sum()
        elif QAB[i,j][4] == '1':
            QAB[i,j] = "Cloud Shadow"      
            Shadow = (QAB == "Cloud Shadow").sum()
        elif QAB[i,j][5] == '1':
            QAB[i,j] = "Adjacent Cloud"
            AC = (QAB == "Adjacent Cloud").sum()

        elif QAB[i,j][7] == '1':
            QAB[i,j] = "Cirrus"
            Cirrus = (QAB == "Cirrus").sum()

        elif QAB[i,j][3] == '1':
            QAB[i,j] = "Snow/Ice"
            SnowC = (QAB == "Snow/Ice").sum()
        elif QAB[i,j][2] == '1':
            QAB[i,j] = "Water"
            WaterC = (QAB == "Water").sum()

        elif QAB[i,j][0:1] == '11':
            QAB[i,j] = "High Aerosol"        
        elif QAB[i,j][0:1] == '10':
            QAB[i,j] = "Avrg. Aerosol"        
        elif QAB[i,j][0:1] == '01':
            QAB[i,j] = "Low Aerosol"
        elif QAB[i,j][0:1] == '00':
            QAB[i,j] = "Aerosol Climatology"

这将产生一个表示不同事物的字符串数组，但正如前面所说，这需要时间

有关如何访问位表示的任何帮助都会很有帮助：）

Answer 1

程序的逻辑实际上精确地映射到numpy.select上，元素根据条件（布尔）数组列表从数组列表中选择，第一个匹配获胜。所以你可以简洁地写下

conditions = QAB&(128>>np.arange(8))[:,None,None]
values = ["OutOfSwath","Cloud","Cloud Shadow","Adjacent Cloud","Cirrus","Snow/Ice",
          "Water","High Aerosol","Avrg. Aerosol","Low Aerosol","Aerosol Climatology"]

np.select([QAB==255,*conditions[[6,4,5,7,3,2]],*(QAB==np.arange(0,256,64)[::-1,None,None])],values)

Answer 2

您可以使用numpy函数unpackbits来处理位。对于numpy，我们更喜欢使用numpy方法和函数，而不是python for loop—通常更快。因此，您可以将每个数字解压到第三个轴，然后像QAB[i, j][5] == '1'这样的条件变为result[bits[5]]。我颠倒了elif子句的顺序，就好像你QAB[i, j][6] == '1'，然后它将其设置为"Cloud"，并且从不重写子句，所以如果我们运行每个条件，这个例子应该是最后一个重写的例子。另外，您的最后一个案例（如QAB[i,j][0:1] == '11'）从未触发，因为左侧的长度始终为1。所以我用你说的QAB[i,j][0:2] == ...来重写

"""
>>> print(bits)
[[[False False False]
  [False False  True]]
<BLANKLINE>
 [[False False False]
  [False False  True]]
<BLANKLINE>
 [[False False False]
  [False False  True]]
<BLANKLINE>
 [[False False False]
  [False False  True]]
<BLANKLINE>
 [[False False False]
  [False  True  True]]
<BLANKLINE>
 [[False False False]
  [False False  True]]
<BLANKLINE>
 [[False  True  True]
  [ True False  True]]
<BLANKLINE>
 [[ True False  True]
  [ True  True  True]]]
>>> print(result)
[['Cirrus' 'Cloud' 'Cloud']
 ['Cloud' 'Cloud Shadow' 'OutOfSwath']]
>>> Cloud
3
"""
import numpy as np
QA = [[1, 2, 3], [3, 9, 255]]

bits = np.unpackbits(np.array(QA, dtype=np.uint8, ndmin=3), axis=0).astype(bool)
result = np.array([[""] * bits.shape[2] for _ in range(bits.shape[1])], dtype=object)
result[~bits[0] & ~bits[1]] = "Aerosol Climatology"
result[~bits[0] & bits[1]] = "Low Aerosol"
result[bits[0] & ~bits[1]] = "Avrg. Aerosol"
result[bits[0] & bits[1]] = "High Aerosol"
result[bits[2]] = "Water"
result[bits[3]] = "Snow/Ice"
result[bits[7]] = "Cirrus"
result[bits[5]] = "Adjacent Cloud"
result[bits[4]] = "Cloud Shadow"
result[bits[6]] = "Cloud"
result[bits.all(axis=0)] = "OutOfSwath"
Cloud = (result == "Cloud").sum()
Shadow = (result == "Cloud Shadow").sum()
AC = (result == "Adjacent Cloud").sum()
Cirrus = (result == "Cirrus").sum()
SnowC = (result == "Snow/Ice").sum()
WaterC = (result == "Water").sum()