# from netCDF4 import Dataset

import xarray as xr
import numpy as np

import ancillary_data as anc

_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


def read_ancillary_data(viirs_data):

    # Ancillary files temporarily defined here. Eventually we will find a better way to pass these
    start_time = '2014-08-01 15:48:00.000'
    anc_dir = '/ships19/hercules/pveglio/neige_data/snpp_test_input/ancillary/2014_08_01_213'
    sst_file = 'oisst.20140730'
    ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.209.hdf'
    geos_cnst = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
    geos_lnd = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20140801_1530.V01.nc4'
    geos_ocn = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20140801_1530.V01.nc4'
    geos1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1500.V01.nc4'
    geos2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1800.V01.nc4'

    sst = np.empty((3232*3200, ), dtype=np.float32)
    ndvi = np.empty((3232*3200, ), dtype=np.float32)
    eco = np.empty((3232*3200, ), dtype=np.ubyte)

    geos_data = {'tpw': np.empty((3232*3200, ), dtype=np.float32),
                 'snowfr': np.empty((3232*3200, ), dtype=np.float32),
                 'icefr': np.empty((3232*3200, ), dtype=np.float32),
                 'ocnfr': np.empty((3232*3200, ), dtype=np.float32),
                 'landicefr': np.empty((3232*3200, ), dtype=np.float32),
                 'sfct': np.empty((3232*3200, ), dtype=np.float32),
                 }

    sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )),
                                  viirs_data.longitude.values.reshape((3232*3200, )),
                                  anc_dir, sst_file, sst)

    ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((3232*3200, )),
                                      viirs_data.longitude.values.reshape((3232*3200, )),
                                      anc_dir, ndvi_file, ndvi)

    eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((3232*3200, )),
                               viirs_data.longitude.values.reshape((3232*3200, )),
                               anc_dir, eco)

    geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )),
                                viirs_data.longitude.values.reshape((3232*3200, )),
                                start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data)

    data = viirs_data
    data['sst'] = (['number_of_lines', 'number_of_pixels'], sst)
    data['ndvi'] = (['number_of_lines', 'number_of_pixels'], ndvi)
    data['eco'] = (['number_of_lines', 'number_of_pixels'], eco)
    data['geos_tpw'] = (['number_of_lines', 'number_of_pixels'], geos_data['tpw'])
    data['geos_snowfr'] = (['number_of_lines', 'number_of_pixels'], geos_data['snowfr'])
    data['geos_icefr'] = (['number_of_lines', 'number_of_pixels'], geos_data['icefr'])
    data['geos_ocnfr'] = (['number_of_lines', 'number_of_pixels'], geos_data['ocnfr'])
    data['geos_landicefr'] = (['number_of_lines', 'number_of_pixels'], geos_data['landicefr'])
    data['geos_sfct'] = (['number_of_lines', 'number_of_pixels'], geos_data['sfct'])

    return data