separate the Fmask from HLS into the different masks that compose it

Fmask_bit_band_separation.py

# separate the Fmask from HLS into the different masks that compose it
# meaning of each output band:
# cirrus, cloud, adjacent cloud, cloud shadow, snow or ice, water, aerosol
# by Luisa V. Lucchese 
# June 2023

# MIT License

#Copyright (c) 2023 Luisa V. Lucchese

#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.

from qgis.core import *
import processing
import numpy as np
from osgeo import gdal

def decodefmaskvalue(value):
    #simple function to know the bits that are 1 in a byte
    #in a byte, each bit means:
    bitcomp=[0,1,2,4,8,16,32,64,128]
    #initialization for the remainders vector
    bitcomp_mod=[0,1,2,4,8,16,32,64,128]
    #initialization for the subtraction vector
    subtraction=[0,1,2,4,8,16,32,64]
    bit_val=[0,0,0,0,0,0,0,0]
    for i in range(0,8):
        bitcomp_mod[i]=value%bitcomp[i+1]
    for i in range(0,8):
        subtraction[i]=bitcomp_mod[i]-bitcomp_mod[i+1]
        if (subtraction[i]<0):
            bit_val[i]=1
    return bit_val


def Fmask_bit_band_separation(rasterloc,saveloc):
    #generates a 7-band raster saved in saveloc
    #each band is a different mask from Fmask
    
    #open the fmask file and read the array
    openedraster = gdal.Open(rasterloc, gdal.GA_ReadOnly)
    
    #get geotransformation
    GT_input = openedraster.GetGeoTransform()
    
    #read the one and only band from Fmask
    raster_input = openedraster.GetRasterBand(1)
    fmask_array = raster_input.ReadAsArray()
    #save the shape of the array to use it in loops
    size1,size2=fmask_array.shape 

    bands=np.zeros((size1,size2,9)) #initialization

    #running the function for every point on Fmask
    for i in range(0,(size1)): #
        for j in range(0,(size2)): #
            bands[i,j,0:8]=decodefmaskvalue(fmask_array[i,j])

    #last two bands mean actually just one thing, aerosol
    for i in range(0,(size1)): #
        for j in range(0,(size2)): #
            if (bands[i,j,7]==0 and bands[i,j,8]==0):
                bands[i,j,7]=0 #climatology
            elif (bands[i,j,7]==0 and bands[i,j,8]==1):
                bands[i,j,7]=1 #low
            elif (bands[i,j,7]==1 and bands[i,j,8]==0):
                bands[i,j,7]=2 #average
            elif (bands[i,j,7]==1 and bands[i,j,8]==1):
                bands[i,j,7]=3 #high
            else:
                bands[i,j,7]=4 #somethings off

    #save this as a raster in Geotiff format
    driver = gdal.GetDriverByName( 'GTiff' )
    dst_ds=driver.Create(saveloc,size2,size1,7,gdal.GDT_Float64)
    dst_ds.SetGeoTransform(GT_input)

    for iband in range(0,7):
        dst_ds.GetRasterBand(iband+1).WriteArray(bands[0:size1,0:size2,iband])
    
    dst_ds=None    

    #open file automatically
    fmask_layer=iface.addRasterLayer(saveloc, "Fmask")
    return fmask_layer