Skip to content
Snippets Groups Projects
restoral.py 2.13 KiB
Newer Older
import numpy as np

from numpy.lib.stride_tricks import sliding_window_view

_bad_data = -999.0


def spatial_var(rad, threshold):

    test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))
    test = np.abs(test.reshape(rad.shape[0], rad.shape[1], 9))

    var = np.ones(test.shape)
    var[test < threshold] = 0

    return var.sum(axis=2)


def sunglint(viirs_data, threshold, bit, conf):

    m09 = viirs_data.M02.values
    m20 = viirs_data.M12.values
    m31 = viirs_data.M15.values
    m02 = viirs_data.M07.values
    irclr = np.ones(viirs_data.M02.shape)

    irclr[bit < 3] = 0

    var = spatial_var(m31, 0.4)
    reg_var_mean = sliding_window_view(np.pad(m02, [1, 1], mode='constant'),
                                       (3, 3)).reshape(m02.shape[0], m02.shape[1], 9).mean(axis=2)
    reg_std = sliding_window_view(np.pad(m02, [1, 1], mode='constant'),
                                  (3, 3)).reshape(m02.shape[0], m02.shape[1], 9).std(axis=2)
    d37_11 = m20 - m31

    idx = np.nonzero((var == 0) & (m09 != _bad_data) & (m20 != _bad_data) & (m31 != _bad_data) & (irclr == 1) &
                     (m02 != _bad_data) & (m09 < threshold['sngm09']) & (d37_11 >= threshold['sg_tbdfl']))
    conf[idx] = 0.67

    idx = np.nonzero((var == 0) & (m09 != _bad_data) & (m20 != _bad_data) & (m31 != _bad_data) & (irclr == 1) &
                     (m02 != _bad_data) & (m09 < threshold['sngm09']) & (d37_11 >= threshold['sg_tbdfl']) &
                     (reg_var_mean != _bad_data) & (reg_var_mean*reg_std < threshold['sngm02vm']))
    conf[idx] = 0.96

    return conf


def spatial(viirs_data, threshold, scene, conf):

    m02 = viirs_data.M02.values
    m31 = viirs_data.M15.values

    var_m31 = spatial_var(m31, 0.40)
    var_m02 = spatial_var(m02, 0.0020)

    idx = np.nonzero((conf > 0.95) & (scene['day'] == 1) & (var_m31 == 0) & (var_m02 == 0))
    conf[idx] = 1

    idx = np.nonzero((conf > 0.66) & (conf <= 0.95) & (scene['day'] == 1) & (var_m31 == 0) & (var_m02 == 0))
    conf[idx] = 0.96

    idx = np.nonzero((conf <= 0.66) & (scene['day'] == 1) & (var_m31 == 0))
    conf[idx] = 0.67

    return conf