preprocess_thresholds.py

import numpy as np
import xarray as xr

import ancillary_data as anc

_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'])

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

    return xr.Dataset(data_vars=thr_dict)


def preproc(data, thresholds):
    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)
    hicut = midpt - thr_xr.adj
    # this below is for the method 2 of computing hicut
    # hicut = midpt - (thr_xr.adj * midpt)

    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


# 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