import numpy as np
# import xarray as xr
import netCDF4 as nc
from glob import glob
try:
from ruamel import yaml as yml
except ImportError:
import ruamel_yaml as yml
import spectral_tests as tst
import restoral
import read_data as rd
import scene as scn
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0]
_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1454.001.*.uwssec.nc')[0]
# MVCM output
_cld_mask = f'{_datapath}/CLDMSK_L2_VIIRS_SNPP.A2022173.1454.001.2022174035130.nc'
# thresholds:
_threshold_file = '/home/pveglio/mvcm/thresholds.mvcm.snpp.v0.0.1.yaml'
# ancillary files:
_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4'
_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4'
_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1430.V01.nc4'
_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.V01.nc4'
_geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
_ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
_sst_file = 'oisst.20220622'
_eco_file = 'goge1_2_img.v1'
# #################################################################### #
def main(*,
data_path: str = _datapath,
mod02: str = _fname_mod02,
mod03: str = _fname_mod03,
img02: str = _fname_img02,
img03: str = _fname_img03,
threshold_file: str = _threshold_file,
geos_atm_1: str = _geos_atm_1,
geos_atm_2: str = _geos_atm_2,
geos_land: str = _geos_land,
geos_ocean: str = _geos_ocean,
geos_constants: str = _geos_constants,
ndvi_file: str = _ndvi_file,
sst_file: str = _sst_file,
cloud_mask_file=_cld_mask) -> None:
f = nc.Dataset(cloud_mask_file)
cm = f['geophysical_data/Cloud_Mask'][:]
f.close()
with open(threshold_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
file_names = {'MOD02': f'{mod02}',
'MOD03': f'{mod03}',
'IMG02': f'{img02}',
'IMG03': f'{img03}',
'GEOS_atm_1': f'{geos_atm_1}',
'GEOS_atm_2': f'{geos_atm_2}',
'GEOS_land': f'{geos_land}',
'GEOS_ocean': f'{geos_ocean}',
'GEOS_constants': f'{geos_constants}',
'NDVI': f'{ndvi_file}',
'SST': f'{sst_file}',
'ANC_DIR': f'{data_path}/ancillary'
}
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
viirs_data = rd.get_data(file_names, sunglint_angle, hires=True)
cm_flag = np.array((np.array(cm[0, :, :], np.byte) &
(pow(2, 1) + pow(2, 2)))/pow(2, 1))
mvcm_conf = np.ones(cm_flag.shape)
# conf[cm_flag == 3] = 1
mvcm_conf[cm_flag == 2] = 0.95
mvcm_conf[cm_flag == 1] = 0.66
mvcm_conf[cm_flag == 0] = 0
viirs_lores = rd.read_data('viirs', mod02, mod03)
viirs_data['M02'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M02.values.repeat(2, 0).repeat(2, 1))
viirs_data['M07'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M07.values.repeat(2, 0).repeat(2, 1))
viirs_data['M01'] = (('number_of_lines_2', 'number_of_pixels_2'), np.zeros(viirs_data.I01.shape))
viirs_data['M10'] = viirs_data.I03
viirs_data['M12'] = viirs_data.I04
viirs_data['M15'] = viirs_data.I05
cmin_G1 = np.ones(viirs_data.M15.shape)
cmin_G2 = np.ones(viirs_data.M15.shape)
cmin_G3 = np.ones(viirs_data.M15.shape)
MyScene = tst.CloudTests(viirs_data, 'Ocean_Day', thresholds)
cmin_G1, bit1 = MyScene.test_11um('M15', cmin_G1, test_name='11um_Test')
cmin_G2, bit2 = MyScene.oceanic_stratus_11_4um_test('I05-I04', cmin_G2,
test_name='11-4um_Oceanic_Stratus_Test')
cmin_G3, bit3 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
cmin_G3, bit4 = MyScene.vis_nir_ratio_test('I02-I01ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
cmin_G3, bit5 = MyScene.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
confidence = np.power(cmin_G1 * cmin_G2 * cmin_G3, 1/3)
temp_confidence = np.power(cmin_G1 * cmin_G2 * cmin_G3, 1/3)
total_bit = np.ones(confidence.shape)
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
scene_flags = scn.find_scene(viirs_data, sunglint_angle)
idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['ice'] == 0) &
(scene_flags['uniform'] == 1) & (confidence <= 0.99) & (confidence >= 0.05))
confidence[idx] = restoral.spatial(viirs_data, thresholds['Sun_Glint'], scene_flags, temp_confidence)[idx]
temp_confidence = confidence + 0
idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
(scene_flags['uniform'] == 1) & (confidence <= 0.95))
confidence[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, temp_confidence)[idx]
# total_bit = bit1
total_confidence = np.sqrt(confidence * mvcm_conf.repeat(2, 0).repeat(2, 1))
date = mod02.split('.')[1]
time = mod02.split('.')[2]
np.savez(f'{data_path}/output/pyMVCM_{date}.{time}',
confidence=total_confidence, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values,
ocean_day_mask=viirs_data.Ocean_Day.values)
# return confidence # , conf, total_bit
def read_bits(byte, bits):
# orig_shape = byte.shape
# if byte.ndim > 1:
# byte = byte.reshape((byte.shape[0]*byte.shape[1],))
if type(bits) is int:
flag = np.array((np.array(byte, np.byte) & pow(2, bits))/pow(2, bits),
dtype='bool')
else:
flag = (np.array(byte, np.byte) & pow(2, bits[0]) +
pow(2, bits[1]))/pow(2, bits[0])
flag = np.array(np.trunc(1.0-flag/2.0), dtype='bool')
# flag = flag.reshape((orig_shape))
return flag
if __name__ == "__main__":
main_fcn()