Skip to content
Snippets Groups Projects
restoral.py 3.91 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
    idx = np.nonzero((m02 > 1.3) & (m01 <= 2.0))
    m02[idx] = 1.3
    idx = np.nonzero((m02 < -99) | (m02 > 1.3))
    m02[idx] = _bad_data

    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):

    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


def land(viirs_data, threshold, scene, conf):
    m04 = viirs_data.M04.values.ravel()
    m05 = viirs_data.M08.values.ravel()
    m20 = viirs_data.M12.values.ravel()
    m22 = viirs_data.M13.values.ravel()
    m31 = viirs_data.M15.values.ravel()
    eco = viirs_data.eco.values.ravel()
    desert = scene['desert'].ravel()
    conf = conf.ravel()
    tbadj = 0
    ldsbt11bd = np.array(threshold['Land_Restoral']['ldsbt11bd'])
    ldsbt11 = np.array(threshold['Land_Restoral']['ldsbt11bd'])

    irclr = 1
    hds11 = np.ones((eco.shape[0], 3)) * (ldsbt11 - tbadj)
    hds11[eco == 8, :] = ldsbt11bd - tbadj

    if irclr == 1:
        conf[m31 > hds11[:, 2]] = 1
        conf[(m31 > hds11[:, 1]) & (m31 <= hds11[:, 2])] = 0.96
        conf[m31 <= hds11[:, 1]] = 0.5

    m5_4_thr = np.full(eco.shape, threshold['Land_Restoral']['ldr5_4_thr'])
    m5_4_thr[desert == 1] = threshold['Land_Restoral']['ldsr5_4_thr']

    m5_4 = m05/m04
    md1 = m20 - m22
    md2 = m22 - m31

    idx = np.nonzero((md1 < threshold['Land_Restoral']['ld20m22']) &
                     (md2 < threshold['Land_Restoral']['ld22m31']) &
                     (m5_4 > m5_4_thr) &
                     (conf <= 0.95))
    conf[idx] = 0.96

    conf = conf.reshape(viirs_data.M01.shape)
    return conf


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

    m01 = viirs_data.M05.values
    m02 = viirs_data.M07.values
    coast_ndvi = threshold['Coastal_NDVI_Thresholds']['coast_ndvi']

    irclr = 1

    if irclr == 1:
        ndvi = (m02 - m01)/(m01 + m02)

        idx = np.nonzero((ndvi <= coast_ndvi[0]) | (ndvi >= coast_ndvi[1]))
        conf[idx] = 1

    return conf