# 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) in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi) in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint) 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)) sunglint_angle = np.arccos(cossna) * _RTD return sunglint_angle def correct_reflectances(): pass