Skip to content
Snippets Groups Projects
main.py 23.5 KiB
Newer Older
import ruamel_yaml as yml
import numpy as np
import read_data as rd
# import ocean_day_tests as odt
# import restoral
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0]
_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1454.001.*.uwssec.nc')[0]
_threshold_file = '/home/pveglio/mvcm/thresholds.mvcm.snpp.v0.0.1.yaml'
_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4'
_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4'
_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1430.V01.nc4'
_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.V01.nc4'
_geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
_ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
_sst_file = 'oisst.20220622'
# #################################################################### #
Paolo Veglio's avatar
Paolo Veglio committed


def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
         img02=_fname_img02, img03=_fname_img03, threshold_file=_threshold_file,
         geos_atm_1=_geos_atm_1, geos_atm_2=_geos_atm_2, geos_land=_geos_land,
         geos_ocean=_geos_ocean, geos_constants=_geos_constants, ndvi_file=_ndvi_file, sst_file=_sst_file):

    # datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
    # fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
    # fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
    # thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'

    file_names = {'MOD02': f'{mod02}',
                  'MOD03': f'{mod03}',
                  'IMG02': f'{img02}',
                  'IMG03': f'{img03}',
                  'GEOS_atm_1': f'{geos_atm_1}',
                  'GEOS_atm_2': f'{geos_atm_2}',
                  'GEOS_land': f'{geos_land}',
                  'GEOS_ocean': f'{geos_ocean}',
                  'GEOS_constants': f'{geos_constants}',
                  'NDVI': f'{ndvi_file}',
                  'SST': f'{sst_file}',
                  'ANC_DIR': f'{data_path}/ancillary'
                  }
        text = f.read()
    thresholds = yml.safe_load(text)

    sunglint_angle = thresholds['Sun_Glint']['bounds'][3]

    viirs_data = rd.get_data(file_names, sunglint_angle)

    # scene_xr = xr.Dataset()
    # for s in scn._scene_list:
    #    scene_xr[s] = (('number_of_lines', 'number_of_pixels'), scn.scene_id[s])
    # scene_xr['latitude'] = viirs_xr.latitude
    # scene_xr['longitude'] = viirs_xr.longitude
    #
    # viirs_data = xr.Dataset(viirs_xr, coords=scene_xr)
    # viirs_data.drop_vars(['latitude', 'longitude'])

    cmin_G1 = np.ones(viirs_data.M01.shape)
    cmin_G2 = np.ones(viirs_data.M01.shape)
    cmin_G3 = np.ones(viirs_data.M01.shape)
    cmin_G4 = np.ones(viirs_data.M01.shape)
    cmin_G5 = np.ones(viirs_data.M01.shape)

#    cmin_test = {'Ocean_Day': np.ones(viirs_data.M01.shape),
#                 'Polar_Ocean_Day': np.ones(viirs_data.M01.shape),
#                 'Polar_Ocean_Night': np.ones(viirs_data.M01.shape)
#                 }
#    cmin2 = np.ones(viirs_data.M01.shape)
#    cmin3 = np.ones(viirs_data.M01.shape)
#    cmin4 = np.ones(viirs_data.M01.shape)

    # ------------------- #
    # ### Testing Setup ###
    # ------------------- #
    perform = {'11um BT Test': False,
               'CO2 High Clouds Test': False,
               'Water Vapor High Clouds Test': False,
               'Surface Temperature Test': False,
Paolo Veglio's avatar
Paolo Veglio committed
               'SST Test': False,
               '8.6-11um BT Difference Test': False,
               '11-12um BTD Thin Cirrus Test': False,
               '11-4um BT Difference Test': False,
               '7.3-11um BT Difference Mid-level Clouds': False,
               'Water Vapor Cloud Test': False,
               '11um BT Variability Test': False,
               '11-4um BTD Oceanic Stratus': False,
Paolo Veglio's avatar
Paolo Veglio committed
               'NIR Reflectance Test': True,
               'Vis/NIR Ratio Test': False,
               '1.6um or 2.1um NIR Reflectance Test': False,
               'Visible Reflectance Test': False,
               'GEMI Test': False,
               '1.38um High Cloud Test': False,
               '4-12um BTD Thin Cirrus Test': False
               }

    # --------------------- #
    # ### Group 1 Tests ### #
    # --------------------- #

    if perform['11um BT Test'] is True:
        # 11um BT Test
        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G1 = SceneType.single_threshold_test('11um_Test', 'M15', cmin_G1)

    if perform['CO2 High Clouds Test'] is True:
        # CO2 High Clouds Test
        for scene_name in ['Land_Day', 'Land_Night', 'Land_Day_Coast', 'Land_Day_Desert',
                           'Land_Day_Desert_Coast', 'Ocean_Day', 'Ocean_Night', 'Day_Snow', 'Night_Snow']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G1 = SceneType.single_threshold_test('CO2_High_Clouds_tests', 'bad_data', cmin_G1)

    if perform['Water Vapor High Clouds Test'] is True:
        # Water Vapor High Clouds Test
        for scene_name in ['Land_Day', 'Land_Night', 'Land_Day_Coast', 'Land_Day_Desert',
                           'Land_Day_Desert_Coast', 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land',
                           'Polar_Night_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
                           'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow', 'Polar_Ocean_Day',
                           'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow', 'Antarctic_Day']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G1 = SceneType.single_threshold_test('Water_Vapor_High_Clouds_tests', 'bad_data', cmin_G1)

    if perform['Surface Temperature Test'] is True:
        # Surface Temperature Test
        # NOTE: this test requires the thresholds to be preprocessed
        for scene_name in ['Land_Night', 'Polar_Night_Land']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G1 = SceneType.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)

    if perform['SST Test'] is True:
        # SST Test
        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G1 = SceneType.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)

    # --------------------- #
    # ### Group 2 tests ### #
    # --------------------- #

    if perform['8.6-11um BT Difference Test'] is True:
        # 8.6-11um BT Difference Test
        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G2 = SceneType.single_threshold_test('8.6-11um_Test', 'M15-M14', cmin_G2)

    if perform['11-12um BTD Thin Cirrus Test'] is True:
        # 11-12um BT BTD Transmissive Cirrus Test
        # NOTE: some of the tests have some differences in how the thresholds are derived
        # The commented list scene_name is the complete list, the one currently in use is to test that
        # the template code works at least with a subset of scenes that have the same way of deriving the
        # thresholds
        # for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
        #                    'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land', 'Polar_Day_Coast',
        #                    'Polar_Day_Desert', 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow',
        #                    'Polar_Ocean_Day', 'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow']:
        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
                           'Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G2 = SceneType.single_threshold_test('11-12um_Cirrus_Test', 'M15-M16', cmin_G2)

    if perform['11-4um BT Difference Test'] is True:
        # 11-4um BT Difference Test
        for scene_name in ['Ocean_Night', 'Polar_Ocean_Night']:
            pass

    if perform['7.3-11um BT Difference Mid-level Clouds'] is True:
        for scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
            pass

    if perform['Water Vapor Cloud Test'] is True:
        for scene_name in ['Polar_Night_Ocean']:
            pass

    if perform['11um BT Variability Test'] is True:
        for scene_name in ['Polar_Night_Ocean']:
            pass

    if perform['11-4um BTD Oceanic Stratus'] is True:
        # 11-4um BT Difference for oceanic stratus (low emissivity water cloud) Test
        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
                           'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land', 'Polar_Day_Coast',
                           'Polar_Day_Desert', 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow',
                           'Polar_Ocean_Day', 'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow']:
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G2 = SceneType.single_threshold_test('11-4um_BTD_Oceanic_Stratus_Test', 'M15-M13', cmin_G2)

    # --------------------- #
    # ### Group 3 tests ### #
    # --------------------- #

    if perform['NIR Reflectance Test'] is True:
        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
Paolo Veglio's avatar
Paolo Veglio committed
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G3 = SceneType.single_thredhold_test('NIR_Reflectance_Test', 'M07', cmin_G3)

    if perform['Vis/NIR Ratio Test'] is True:
        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
            pass

    if perform['1.6um or 2.1um NIR Reflectance Test'] is True:
        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
            pass

    if perform['Visible Reflectance Test'] is True:
        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
                           'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', 'Polar_Day_Desert_Coast']:
            pass

    if perform['GEMI Test'] is True:
        for scene_name in ['Land_Day_Desert', 'Polar_Day_Desert']:
            pass

    # --------------------- #
    # ### Group 4 tests ### #
    # --------------------- #

    if perform['1.38um High Cloud Test'] is True:
        # for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
        #                    'Ocean_Day', 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
        #                    'Polar_Day_Desert_Coast', 'Polar_Ocean_Day', 'Day_Snow']:
        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
                           'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
                           'Polar_Day_Desert_Coast', 'Day_Snow']:
            # The current loop is missing Ocean_Day and Polar_Ocean_Day because they need some
            # preprocessing to compute the thresholds. Once things are implemented I can use the commented
            # loop
            SceneType = CloudTests(viirs_data, scene_name, thresholds)
            cmin_G4 = SceneType.single_threshold_test('1.38um_High_Cloud_Test', 'M09', cmin_G4)

    # --------------------- #
    # ### Group 5 tests ### #
    # --------------------- #

    if perform['4-12um BTD Thin Cirrus Test'] is True:
        for scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
            pass

    cmin_total = cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5
    return cmin_total
    Land_Day = CloudTests(viirs_data, 'Land_Day', thresholds)
    Land_Night = CloudTests(viirs_data, 'Land_Night', thresholds)
    Land_Day_Coast = CloudTests(viirs_data, 'Land_Day_Coast', thresholds)
    Land_Day_Desert = CloudTests(viirs_data, 'Land_Day_Desert', thresholds)
    Land_Day_Desert_Coast = CloudTests(viirs_data, 'Land_Day_Desert_Coast', thresholds)
    Ocean_Day = CloudTests(viirs_data, 'Ocean_Day', thresholds)
    Ocean_Night = CloudTests(viirs_data, 'Ocean_Night', thresholds)
    Polar_Day_Land = CloudTests(viirs_data, 'Polar_Day_Land', thresholds)
    Polar_Night_Land = CloudTests(viirs_data, 'Polar_Night_Land', thresholds)
    Polar_Day_Coast = CloudTests(viirs_data, 'Polar_Day_Coast', thresholds)
    Polar_Day_Desert = CloudTests(viirs_data, 'Polar_Day_Desert', thresholds)
    Polar_Day_Desert_Coast = CloudTests(viirs_data, 'Polar_Day_Desert_Coast', thresholds)
    Polar_Day_Snow = CloudTests(viirs_data, 'Polar_Day_Snow', thresholds)
    Polar_Night_Snow = CloudTests(viirs_data, 'Polar_Night_Snow', thresholds)
    Polar_Ocean_Day = CloudTests(viirs_data, 'Polar_Ocean_Day', thresholds)
    Polar_Ocean_Night = CloudTests(viirs_data, 'Polar_Ocean_Night', thresholds)
    Day_Snow = CloudTests(viirs_data, 'Day_Snow', thresholds)
    Night_Snow = CloudTests(viirs_data, 'Night_Snow', thresholds)
    Antarctic_Day = CloudTests(viirs_data, 'Antarctic_Day', thresholds)
    cmin_G1 = Ocean_Day.single_threshold_test('11um_Test', 'M15', cmin_G1)
    cmin_G1 = Polar_Ocean_Day.single_threshold_test('11um_Test', 'M15', cmin_G1)
    cmin_G1 = Polar_Ocean_Night.single_threshold_test('11um_Test', 'M15', cmin_G1)

    # CO2 High Clouds Test
    # NOTE: VIIRS doesn't have the MODIS equivalent of B35 so this test is not performed
    cmin_G1 = Land_Day.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Night.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Coast.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Desert.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Desert_Coast.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Ocean_Day.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Ocean_Night.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Day_Snow.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Night_Snow.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)

    # Water Vapor High Clouds Test
    # NOTE: VIIRS doesn't have the MODIS equivalent of B27 so this test is not performed
    cmin_G1 = Land_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Desert.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Land_Day_Desert_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Ocean_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Ocean_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Day_Land.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Night_Land.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Day_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Day_Desert.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Day_Desert_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Day_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Night_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Ocean_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Polar_Ocean_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Day_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Night_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
    cmin_G1 = Antarctic_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)

    # Surface Temperature Test
    # ## NOTE: This requires some calculations for the thresholds.
    # Moreover this test is using the 11um - 12um difference instead of a single channel
    # Also, look at the test carefully for these two cases. Polar_Night_Land uses some hardcoded coeffs
    # (i.e. *_df1 and *_df2) that might be worth moving in the thresholds file for consistency
    cmin_G1 = Land_Night.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)
    cmin_G1 = Polar_Night_Land.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)

    # SST Test
    # ## NOTE: This requires some calculations for the thresholds.
    # Moreover this test is using the 11um - 12um difference instead of a single channel
    cmin_G1 = Ocean_Day.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
    cmin_G1 = Polar_Ocean_Day.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
    cmin_G1 = Polar_Ocean_Night.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
    # the following test uses a different set of coefficients compared to the others above
    cmin_G1 = Ocean_Night.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)

    # 11-8.6um BT Difference Test
    cmin_G2 = Ocean_Day.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
    cmin_G2 = Ocean_Night.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
    cmin_G2 = Polar_Ocean_Day.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
    cmin_G2 = Polar_Ocean_Night.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
    cmin_G1 = Ocean_Day.single_threshold_test('11-12BT_diff',
                                              viirs_data.M15.values-viirs_data.M16.values,
                                              cmin_G1)

    Ocean_Night = CloudTests(viirs_data, 'Ocean_Night', thresholds)

    cmin_G1 = Ocean_Night.single_threshold_test('11um_test', viirs_data.M15.values, cmin_G1)
    cmin_G1 = Ocean_Night.single_threshold_test('11-12um_diff',
                                                viirs_data.M15.values-viirs_data.M16.values,
                                                cmin_G1)

    c = np.ones((9, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1]))
    cmin1, c[0, :, :], bit1 = odt.simple_test(viirs_data.M15.values, thresholds['Daytime_Ocean']['bt11'], cmin1)
    cmin1, c[1, :, :], bit2 = odt.sst_test(viirs_data.M15.values, viirs_data.M16.values,
                                           viirs_data.sensor_zenith.values, viirs_data.geos_sfct.values,
                                           thresholds['Daytime_Ocean']['sst'], cmin1)
    cmin2, c[2, :, :], bit3 = odt.simple_test(viirs_data.M14.values-viirs_data.M15.values,
                                              thresholds['Daytime_Ocean']['diff_11_86um'], cmin2)
    cmin2, c[3, :, :], bit4 = odt.test_11_12_diff(viirs_data, thresholds['Daytime_Ocean']['diff11_12um'], cmin2)
    cmin2, c[4, :, :] = odt.test_11_4_diff(viirs_data.M15.values, viirs_data.M12.values,
                                           thresholds['Daytime_Ocean']['test11_4lo'], scene_flags, cmin2)
    cmin3, c[5, :, :] = odt.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
                                          thresholds['Sun_Glint'], viirs_data, cmin3)
    cmin3, c[6, :, :] = odt.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values,
                                               thresholds, scene_flags, cmin3)
    cmin3, c[7, :, :] = odt.nir_refl_test(viirs_data.M10.values, thresholds['Daytime_Ocean'],
                                          thresholds['Sun_Glint'], viirs_data, cmin3)
#
#    confidence[0, :, :] = tests.test_11um(viirs_data.M15.values, thresholds['Daytime_Ocean'])
#    confidence[1, :, :] = tests.test_11_4diff(viirs_data.M15.values, viirs_data.M13.values,
#                                              thresholds['Daytime_Ocean'], viirs_data,
#                                              thresholds['Sun_Glint']['bounds'][3])
#
#    confidence[2, :, :] = tests.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
#                                              thresholds['Sun_Glint'], viirs_data)
#
#    # Note that here I'm using M05/M07 but the corresponding hi-res channels are I1/I2
#    # IMPORTANT: conf_test_dble() needs to be verified. I don't think it's working as intended at the moment
#    confidence[3, :, :] = tests.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values,
#                                                   thresholds['Daytime_Ocean'], thresholds['Sun_Glint'])
#
#    # This test needs to be verified, for the granule I'm running everything is zero
#    confidence[4, :, :] = tests.test_11um_var(viirs_data.M15.values, thresholds['Nighttime_Ocean'],
#                                              thresholds['Daytime_Ocean_Spatial_Variability'])

    total_bit = bit1 + bit2 + bit4
    temp_confidence = cmin1 * cmin2 * cmin3 * cmin4
    confidence = cmin1 * cmin2 * cmin3 * cmin4
    # idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['ice'] == 0) & (scene_flags['uniform'] == 1) &
    #                  (confidence <= 0.99) & (confidence >= 0.05))
    # confidence[idx] = restoral.spatial(viirs_data, thresholds['Sun_Glint'], scene_flags, confidence)[idx]
    idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
                     (scene_flags['uniform'] == 1) & (confidence <= 0.95))
    confidence[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, temp_confidence)[idx]

    temp = np.zeros((viirs_data.M01.shape[0], viirs_data.M01.shape[1]))
    temp[idx] = 1
    c[8, :, :] = temp

    np.savez('test_confidence', confidence=confidence, conf_test=c,
             lat=viirs_data.latitude.values, lon=viirs_data.longitude.values)

    return confidence

    rad1 = [[255, 260, 265, 248, 223],
            [278, 285, 270, 268, 256],
            [275, 273, 266, 254, 259]]
    rad2 = [[270, 273, 271, 268, 265],
            [277, 286, 275, 277, 269],
            [280, 281, 272, 270, 267]]
    thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'

    with open(thresh_file) as f:
        text = f.read()

    rad1 = np.array(rad1)
    rad2 = np.array(rad2)

    confidence = np.zeros((2, rad1.shape[0], rad1.shape[1]))

    thresholds = yml.safe_load(text)

    confidence[0, :, :] = tests.test_11um(rad1, thresholds['Daytime_Ocean'])
    confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds['Daytime_Ocean'])

    print(f'Confidence[0,:,:]: \n {confidence[0, :, :]}')
    print(f'Confidence[1,:,:]: \n {confidence[1, :, :]}')

    return confidence


if __name__ == "__main__":