Reputation: 101
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 :)
Upvotes: 3
Views: 649
Reputation: 53079
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)
Upvotes: 1
Reputation: 2719
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()
Upvotes: 2