read_data.py

# from netCDF4 import Dataset

import xarray as xr
import numpy as np

_DTR = np.pi/180.
_RTD = 180./np.pi


def read_data(sensor: str, l1b_filename: str, geo_filename: str):

    data = xr.open_dataset(l1b_filename, group='observation_data', decode_cf=False)
    in_data = xr.Dataset()
    if sensor.lower() == 'viirs':
        for band in list(data.variables):
            if 'reflectance' in data[band].long_name:
                if hasattr(data[band], 'VCST_scale_factor'):
                    scale_factor = data[band].VCST_scale_factor * data[band].bias_correction
                else:
                    scale_factor = data[band].radiance_scale_factor
                in_data[band] = (('number_of_lines', 'number_of_pixels'),
                                 data[band].values * scale_factor)

            elif 'radiance' in data[band].long_name:
                in_data[band] = (('number_of_lines', 'number_of_pixels'),
                                 data[f'{band}_brightness_temperature_lut'].values[data[band].values])
            else:
                pass

    data = xr.open_dataset(geo_filename, group='geolocation_data')

    in_data = in_data.merge(data)

    relazi = relative_azimuth_angle(data.sensor_azimuth.values, data.solar_azimuth.values)
    sunglint = sun_glint_angle(data.sensor_zenith.values, data.solar_zenith.values, relazi)
    scatt_angle = scattering_angle(data.solar_zenith.values, data.sensor_zenith.values, relazi)

    in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi)
    in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint)
    in_data['scattering_angle'] = (('number_of_lines', 'number_of_pixels'), scatt_angle)

    return in_data


def relative_azimuth_angle(sensor_azimuth: float, solar_azimuth: float) -> float:
    rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
    return rel_azimuth


def sun_glint_angle(sensor_zenith: float, solar_zenith: float, rel_azimuth: float) -> float:
    cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
              np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
    cossna[cossna > 1] = 1
    sunglint_angle = np.arccos(cossna) * _RTD
    return sunglint_angle


def scattering_angle(solar_zenith, sensor_zenith, relative_azimuth):
    cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) -
                             np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) *
                             np.cos(relative_azimuth*_DTR))
    scatt_angle = np.arccos(cos_scatt_angle) * _RTD
    return scatt_angle


def correct_reflectances():
    pass