ocean_day_tests.py

import numpy as np

from numpy.lib.stride_tricks import sliding_window_view

import utils
import conf

import ancillary_data as c_tools


# ############## GROUP 1 TESTS ############## #

# 11 micron brightness temperature threshold test
def simple_test(rad, threshold, cmin):

    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))

    thr = np.array(threshold)
    confidence = np.ones(rad.shape)
    bit = np.zeros(rad.shape)

    if thr[4] == 1:
        print("simple test running")
        # the C code has the line below that I don't quite understand the purpose of.
        # It seems to be setting the bit to 0 if the BT value is greater than the midpoint
        #
        # if (m31 >= dobt11[1]) (void) set_bit(13, pxout.testbits);

        # confidence = utils.conf_test(rad, thr)
        confidence = conf.conf_test(rad, thr)
        bit[rad >= thr[1]] = 1

    return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape)


def sst_test(rad1, rad2, vza, surf_temp, threshold, cmin):

    a1 = 1.8860
    a2 = 0.9380
    a3 = 0.1280
    a4 = 1.0940

    radshape = rad1.shape
    b31 = rad1.reshape(np.prod(radshape)) - 273.16
    b32 = rad2.reshape(np.prod(radshape)) - 273.16

    thr = np.array(threshold)
    confidence = np.ones(b31.shape)
    bit = np.zeros(b31.shape)

    rad_diff = b31 - b32
    sstc = surf_temp.reshape(np.prod(radshape)) - 273.16
    mu = np.cos(vza.reshape(np.prod(radshape)) * np.pi/180.0)

    modsst = 273.16 + a1 + a2*b31 + a3*rad_diff*sstc + a4*rad_diff*((1/mu)-1)
    sfcdif = surf_temp.reshape(np.prod(radshape)) - modsst

    if thr[4] == 1:
        print('SST test running')
        confidence = conf.conf_test(sfcdif, thr)
        bit[sfcdif < thr[1]] = 1

    return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape)


def test_11_12_diff(data, threshold, cmin):
    radshape = data.M15.shape
    b31 = data.M15.values.reshape(np.prod(radshape))
    b32 = data.M15.values.reshape(np.prod(radshape))
    vza = data.sensor_zenith.values.reshape(np.prod(radshape))
    thr = np.array(threshold)

    confidence = np.ones(b31.shape)
    bit = np.zeros(b31.shape)
    rad_diff = b31 - b32

    # Get secant of viewing zenith angle
    dtr = np.pi/180
    cosvza = np.cos(vza * dtr)
    schi = np.full(cosvza.shape, 99.0)
    schi[cosvza > 0.0] = 1.0/cosvza[cosvza > 0.0]

    # Need to define this in cython
    btd_thr = c_tools.py_cithr(1, np.array(schi, dtype=np.float32), np.array(b31, dtype=np.float32))
    idx = np.nonzero((btd_thr < 0.1) | (np.abs(schi-99.0) < 0.0001))
    btd_thr[idx] = thr[0]
    locut = btd_thr + 0.3*btd_thr
    hicut = btd_thr - 1.25
    corr_thr = np.array([locut, btd_thr, hicut, np.ones(locut.shape)], dtype=np.float)

    if thr[1] == 1:
        print('11-12um diff test running')
        bit[rad_diff < thr[1]] = 1
        confidence = conf.conf_test(rad_diff, corr_thr)

    return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape)


def test_11_4_diff(rad1, rad2, threshold, scene_flags, cmin):
    radshape = rad1.shape
    b31 = rad1.reshape(np.prod(radshape))
    b20 = rad2.reshape(np.prod(radshape))
    thr = np.array(threshold)
    sunglint = scene_flags['sunglint'].reshape(np.prod(radshape))

    confidence = np.ones(b31.shape)

    if thr[4] == 1:
        print('11-4um diff test running')
        confidence[sunglint == 0] = conf.conf_test((b31-b20)[sunglint == 0], thr)

    return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape)


def vis_nir_ratio_test(rad1, rad2, threshold, scene, cmin):
    if threshold['Daytime_Ocean']['vis_nir_ratio'][6] == 1:
        print("NIR-Visible ratio test running")

        radshape = rad1.shape
        rad1 = rad1.reshape(np.prod(radshape))
        rad2 = rad2.reshape(np.prod(radshape))
        sunglint = scene['sunglint'].reshape(np.prod(radshape))
        vrat = rad2/rad1

        confidence = np.ones(rad1.shape)
        tmp = threshold['Daytime_Ocean']['vis_nir_ratio']
        thr_no_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1, 1])
        tmp = threshold['Sun_Glint']['snglnt']
        thr_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1])
        # thr_no_sunglint = np.array(threshold['Daytime_Ocean']['vis_nir_ratio'])
        # thr_sunglint = np.array(threshold['Sun_Glint']['snglnt'])
        # thr_sunglint = np.append(thr_sunglint, 1)
        # temp value to avoid linter bitching at me
        # eventually we would have the test run in two blocks as:
        # confidence[sunglint == 1] = conf.conf_test_dble(vrat[sunglint == 1], sg_threshold['snglnt'])
        # confidence[sunglint == 0] = conf.conf_test_dble(vrat[sunglint == 0], threshold['vis_nir_ratio'])
        # sunglint needs to be defined somewhere
        # thr = np.full((rad.shape[0], 4), thr[:4]).T

        # thresh = np.full((rad1.shape[0], thr_no_sunglint.shape[0]), thr_no_sunglint)
        # thresh[sunglint == 1, :6] = thr_sunglint

        confidence[sunglint == 0] = conf.conf_test_dble(vrat, thr_no_sunglint)[sunglint == 0]
        confidence[sunglint == 1] = conf.conf_test_dble(vrat, thr_sunglint)[sunglint == 1]
        # confidence = conf.conf_test_dble(vrat, thresh.T)

        return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape)


def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data, cmin):

    print("NIR reflectance test running")
    sza = viirs_data.solar_zenith.values
    refang = viirs_data.sunglint_angle.values
    vza = viirs_data.sensor_zenith.values
    dtr = np.pi/180
    # Keep in mind that band_n uses MODIS band numbers (i.e. 2=0.86um and 7=2.1um)
    # For VIIRS would be 2=M07 (0.865um) and 7=M11 (2.25um)
    band_n = 2
    vzcpow = 0.75  # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE

    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))
    confidence = np.ones(rad.shape)
    sza = sza.reshape(rad.shape)
    vza = vza.reshape(rad.shape)
    refang = refang.reshape(rad.shape)
    sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(rad.shape)

    # ref2 [5]
    # b2coeffs [4]
    # b2mid [1]
    # b2bias_adj [1]
    # b2lo [1]
    # vzcpow [3] (in different place)

    cosvza = np.cos(vza*dtr)
    coeffs = threshold['b2coeffs']
    hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3))
    hicut0 = (hicut0 * 0.01) + threshold['b2adj']
    hicut0 = hicut0 * threshold['b2bias_adj']
    midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj'])
    locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj'])
    thr = np.array([locut0, midpt0, hicut0, threshold['ref2'][3]*np.ones(rad.shape)])
    corr_thr = np.zeros((4, rad.shape[0]))

    corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0], vzcpow))
    corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
    # corr_thr[:3, :] = thr[:3, :] * (1./np.power(cosvza[:], vzcpow))
    # corr_thr[3, :] = thr[3, :]

    for flag in range(1, 4):
        if len(refang[sunglint_flag == flag]) > 0:
            sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
            corr_thr[:3, sunglint_flag == flag] = sunglint_thr[:3, sunglint_flag == flag] * (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
            corr_thr[3, sunglint_flag == flag] = sunglint_thr[3, sunglint_flag == flag]

    confidence = conf.conf_test(rad, corr_thr)

    return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape)


def nir_high_cloud_test():
    pass











def test_11um_var(rad, threshold, var_threshold):

    print("11um variability test running")
    thr = np.array(threshold['11um_var'])

    radshape = rad.shape
    var = np.zeros((radshape[0], radshape[1], 9))

    # chk_spatial2() need to figure out what this is
    # np = rg_var.num_small_diffs * 1.0
    test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))

    var[np.abs(test).reshape(radshape[0], radshape[1], 9) < var_threshold['dovar11']] = 1
    var = var.sum(axis=2).reshape(np.prod(radshape))

    rad = rad.reshape(np.prod(radshape))
    confidence = np.zeros(rad.shape)

    confidence[var == 9] = conf.conf_test(rad[var == 9], thr)

    return confidence.reshape(radshape)


def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh):

    print("11um - 4um difference test running")
    radshape = rad1.shape
    raddiff = (rad1 - rad2).reshape(np.prod(radshape))

    day = np.zeros(radshape)
    day[viirs_data.solar_zenith <= 85] = 1
    day = day.reshape(raddiff.shape)
    sunglint = np.zeros(rad1.shape)
    sunglint[viirs_data.sunglint_angle <= sg_thresh] = 1
    sunglint = sunglint.reshape(raddiff.shape)
    thr = np.array(threshold['test11_4lo'])
    confidence = np.zeros(raddiff.shape)

    # confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
    confidence[(day == 1) & (sunglint == 0)] = conf.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)

    return confidence.reshape(radshape)


def vir_refl_test(rad, threshold, viirs_data):

    print('Visible reflectance test running')

    thr = threshold['vis_refl_test']

    radshape = rad.shape()
    rad = rad.reshape(np.prod(radshape))
    confidence = np.zeros(radshape)
    vzcpow = 0.75  # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE

    vza = viirs_data.sensor_zenith.values
    dtr = np.pi/180
    cosvza = np.cos(vza*dtr)

    coeffs = utils.get_b1_thresholds()
    coeffs[:, :3] = coeffs[:, :3] * threshold['b1_bias_adj']

    # this quantity is the return of get_b1_thresholds() in the C code
    # it's defined here to keep a consistent logic with the original source, for now
    irtn = 0

    if irtn != 0:
        coeffs = thr

    coeffs[:, :3] = coeffs[:, :3] * 1/np.power(cosvza, vzcpow)

    confidence = conf.conf_test(rad, coeffs)

    return confidence.reshape(radshape)



class CloudMaskTests(object):

    def __init__(self, scene, radiance, coefficients):
        self.scene = scene
        self.coefficients = coefficients

    def select_coefficients(self):
        pass

    def test_G1(self):
        pass

    def test_G2(self):
        pass

    def test_G3(self):
        pass

    def test_G4(self):
        pass

    def overall_confidence(self):
        pass


def test():
    rad = np.random.randint(50, size=[4, 8])
    # coeffs = [5, 42, 20, 28, 15, 35, 1]
    # coeffs = [20, 28, 5, 42, 15, 35, 1]
    coeffs = [35, 15, 20, 1, 1]
    # confidence = conf_test_dble(rad, coeffs)
    confidence = test_11um(rad, coeffs)
    print(rad)
    print('\n')
    print(confidence)


if __name__ == "__main__":
    test()