Skip to content
Snippets Groups Projects
preprocess_thresholds.py 8.29 KiB
import numpy as np
import xarray as xr

import ancillary_data as anc
import utils

_dtr = np.pi/180


# this case is written for the 11-12um Cirrus Test for scenes that follow pattern 1 (see note below)
def prepare_thresholds(data, thresholds):

    coeff_values = np.empty((data.M01.shape[0], data.M01.shape[1], 2))
    coeff_values[:, :, 0] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][0])
    coeff_values[:, :, 1] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][1])
    cmult_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['cmult'])
    adj_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['adj'])
    if 'bt1' in list(thresholds['11-12um_Cirrus_Test']):
        bt1 = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['bt1'])
    else:
        bt1 = np.full(data.M01.shape, -999)
    if 'lat' in list(thresholds['11-12um_Cirrus_Test']):
        lat = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['lat'])
    else:
        lat = np.full(data.M01.shape, -999)
    thr_dict = {'coeffs': (['number_of_lines', 'number_of_pixels', 'z'], coeff_values),
                'cmult': (['number_of_lines', 'number_of_pixels'], cmult_values),
                'adj': (['number_of_lines', 'number_of_pixels'], adj_values),
                'bt1': (['number_of_lines', 'number_of_pixels'], bt1),
                'lat': (['number_of_lines', 'number_of_pixels'], lat),
                }

    return xr.Dataset(data_vars=thr_dict)
# 3. Others
# - Land_Night
# - Polar_Night_Land
# - Polar_Night_Snow
# - Day_Snow
# - Night_Snow


def preproc(data, thresholds, scene):
    cosvza = np.cos(data.sensor_zenith * _dtr)
    schi = (1/cosvza).where(cosvza > 0, 99.0)

    schi = schi.values.reshape(np.prod(schi.shape))
    m15 = data.M15.values.reshape(np.prod(data.M15.shape))
    thr = anc.py_cithr(1, schi, m15)
    thr = thr.reshape(data.M15.shape)
    schi = schi.reshape(data.M15.shape)

    # thr_xr = xr.Dataset(np.full(data.sensor_zenith.shape, thresholds['coeffs']),
    #                     dims=('number_of_lines', 'number_of_pixels'))
    thr_xr = prepare_thresholds(data, thresholds)

    midpt = thr_xr.coeffs[:, :, 0].where((thr < 0.1) | (np.abs(schi-99) < 0.0001), thr)
    locut = midpt + (thr_xr.cmult * midpt)
    if scene in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
                 'Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
        hicut = midpt - thr_xr.adj
    elif scene in ['Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
                   'Polar_Day_Desert_Coast', 'Polar_Day_Snow']:
        hicut = midpt - (thr_xr.adj * midpt)
    elif scene in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']:
        _coeffs = {'Land_Night': 0.3, 'Polar_Night_Land': 0.3, 'Polar_Night_Snow': 0.3,
                   'Day_Snow': 0.0, 'Night_Snow': 0.3}
        midpt = midpt - (_coeffs[scene] * locut)
        if scene in ['Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
            hicut = (midpt - (0.2 * locut)).where(data.M15 < thr_xr.bt1, midpt - 1.25)
        elif scene in ['Land_Night']:
            hicut = -0.1 - np.power(90.0 - np.abs(data.latitude)/60, 4) * 1.15
            hicut = hicut.where((data.M15 < thr_xr.bt1) & (data.latitude > thr_xr.lat), 1.25)
        elif scene in ['Day_Snow']:
            hicut = locut - (thr_xr.cmult * locut)
        else:
            print('Scene not recognized\n')
    else:
        print('Scene not recognized\n')

    thr_out = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))),
                           dims=('number_of_lines', 'number_of_pixels', 'z'))
    return thr_out
    # return locut, hicut, midpt


def preproc_sst(data, thresholds):
    m31c = data.M15 - 273.16
    m32c = data.M16 - 273.16
    m31c_m32c = m31c - m32c
    sstc = data.geos_sfct - 273.16
    cosvza = np.cos(data.sensor_zenith*_dtr)

    a = thresholds['coeffs']

    modsst = 273.16 + a[0] + a[1]*m31c + a[2]*m31c_m32c*sstc + a[3]*m31c_m32c*((1/cosvza) - 1)
    sfcdif = data.geos_sfct - modsst

    return sfcdif


def preproc_nir(data, thresholds, scene):

    sza = data.solar_zenith
    band_n = 2
    # NOTE: the visud condition in the C code is equivalent to having sza <= 85
    # For the time being the visud filtering is not implemented

    a = np.array(thresholds[scene]['NIR_Reflectance_Test']['coeffs'])
    vzcpow = thresholds['VZA_correction']['vzcpow'][0]
    refang = data.sunglint_angle.values.reshape(np.prod(data.sunglint_angle.shape))
    sunglint_thresholds = thresholds['Sun_Glint']
    sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(refang.shape)
    nir_thresh = thresholds[scene]['NIR_Reflectance_Test']

    hicut0 = a[0] + a[1]*sza + a[2]*np.power(sza, 2) + a[3]*np.power(sza, 3)
    hicut0 = (hicut0 * 0.01) + nir_thresh['adj']
    hicut0 = (hicut0 + nir_thresh['bias']).values.reshape(refang.shape)
    midpt0 = hicut0 + (nir_thresh['midpt_coeff'] * nir_thresh['bias'])
    locut0 = midpt0 + (nir_thresh['locut_coeff'] * nir_thresh['bias'])
    thr = np.array([locut0, midpt0, hicut0, nir_thresh['thr'][3]*np.ones(refang.shape)])

    cosvza = np.cos(data.sensor_zenith*_dtr).values.reshape(refang.shape)
    corr_thr = np.zeros((4, refang.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]

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

    corr_thr = np.transpose(corr_thr.reshape((4, sza.shape[0], sza.shape[1])), (1, 2, 0))

    return corr_thr


def preproc_surf_temp(data, thresholds):
    thr_sfc1 = thresholds['Surface_Temperature_Test_1']
    thr_sfc2 = thresholds['Surface_Temperature_Test_2']
    thr_df1 = thresholds['Surface_Temperature_Test_df1']
    thr_df2 = thresholds['Surface_Temperature_Test_df2']
    max_vza = 70.13  # This values is set based on sensor. Check mask_processing_constants.h for MODIS value

    rs = np.prod(data.M15.shape)
    df1 = (data.M15 - data.M16).values.reshape(rs)
    df2 = (data.M15 - data.M13).values.reshape(rs)
    desert_flag = data.Desert.values.reshape(rs)
    thresh = np.ones((rs, )) * thr_sfc1

    idx = np.where((df1 >= thr_df1[0]) | ((df1 < thr_df1[0]) & ((df2 <= thr_df2[0]) | (df2 >= thr_df2[1]))))
    thresh[idx] = thr_sfc2
    idx = np.where(desert_flag == 1)
    thresh[idx] == thr_sfc1

    midpt = thresh
    idx = np.where(df1 >= thr_df1[1])
    midpt[idx] = thresh[idx] + 2.0*df1[idx]

    corr = np.power(data.sensor_zenith.values/max_vza, 4) * 3.0
    midpt = midpt.reshape(corr.shape) + corr
    locut = midpt + 2.0
    hicut = midpt - 2.0

    thr_out = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))),
                           dims=('number_of_lines', 'number_of_pixels', 'z'))

    return thr_out


def get_b1_thresholds():
    # fill_ndvi[0] is for fill_ndvi1
    # fill_ndvi[1] is for fill_ndvi2
    pass


# NOTE: 11-12um Cirrus Test
# hicut is computed in different ways depending on the scene
# 1. midpt - adj
# - Land_Day
# - Land_Day_Coast
# - Land_Day_Desert
# - Land_Day_Desert_Coast
# - Ocean_Day
# - Ocean_Night
# - Polar_Day_Ocean
# - Polar_Night_Ocean
#
# 2. midpt - (btd_thr * adj)
# - Polar_Day_Land
# - Polar_Day_Coast
# - Polar_Day_Desert
# - Polar_Day_Desert_Coast
# - Polar_Day_Snow
#
# 3. Others
# - Land_Night
# - Polar_Night_Land
# - Polar_Night_Snow
# - Day_Snow
# - Night_Snow

# NOTE: 1.38um High Cloud Test
# thresholds are not always computed the same way. In group 1 there's no preprocessing required,
# in group 2 some calcuations are needed
# 1.
# - 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
# - Polar_Day_Snow
# - Day_Snow
#
# 2.
# - Ocean_Day
# - Polar_Ocean_Day