3

I am trying to interpret HDF-files or more precisely the Rasterbands inside the files. There is a Band that is used for Quality Assessment which represents Bit-Information as corresponding integers (e.g. 01000000 as 64).

Depending on specific bits I want to get a count, e.g. how many of the Pixels are 1 on bit position 5. If they are taken into account by a previous count, the should not get count again. Right now I am changing the entries in each cell depending on a priority-list by myself. This takes ages and I am really sure that there are faster ways since I have never worked with bits before.

Here is my code:

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"

This will lead to an array of strings representing the different things, but takes ages as said before.

Any help in how to access the Bit-Representation would be helpful :)

Improve this question

2 Answers 2

Reset to default
2

You can use numpy function unpackbits to handle bits. With numpy we prefer to use numpy methods and functions instead of python for loop - generally it is faster. So you can unpack each number to the third axis and then conditions like QAB[i, j][5] == '1' turn to result[bits[5]]. I reverse the order of you elif clauses as if you QAB[i, j][6] == '1' then it set it to "Cloud" and never run over clauses so this case should be last to rewrite overs if we run each condition. Also your last cases like QAB[i,j][0:1] == '11' never triggered as left hand side always has the length 1. So I rewrite with the consumption that you mean 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()
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

1

The logic of your program actually maps precisely onto numpy.select which element-wise selects from a list of arrays according to a list of condition (boolean) arrays, first match wins. So you could compactly write something like

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)
Improve this answer

6 Comments

I tried you solution but get an error: TypeError: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
@Felix What type is your (original) QAB? It should be an intenger type. Could you please share the output of QAB.dtype?
QA is "Array of uint8", QAB was the result of bin_vec = np.vectorize(np.binary_repr) QAB = bin_vec(QA, width = 8), it is an "array of str256". I used ur code-snippet on QAB, but maybe should have used it on QA instead if I interpret ur answer right.
@Felix or rather I shouöld have spotted the difference between QA and QAB. Anyways, does it work with QA?
We had holiday yesterday and i was not in office today to test it. I will check on monday and update you :)
|

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.