Newer
Older
import sys
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
import importlib
importlib.reload(tst)
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
_fname_mod02 = ''
_fname_mod03 = ''
_fname_img02 = ''
_fname_img03 = ''
_cld_mask = ''
# 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,
out_fname: str = 'temp.nc') -> None:
check = nc.Dataset(mod02)
if 'M01' not in list(check['observation_data'].variables):
check.close()
return
if check.isopen():
check.close()
print(threshold_file)
with open(threshold_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
file_names = {'MOD02': f'{data_path}/{mod02}',
'MOD03': f'{data_path}/{mod03}',
'IMG02': f'{data_path}/{img02}',
'IMG03': f'{data_path}/{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)
if 'I01' not in viirs_data:
return
# 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, data_path + '/' + mod03)
viirs_data['M02'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M02.values.repeat(2, 0).repeat(2, 1))
Paolo Veglio
committed
# 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))
Paolo Veglio
committed
viirs_data['M07'] = viirs_data.I02
viirs_data['M10'] = viirs_data.I03
viirs_data['M12'] = viirs_data.I04
Paolo Veglio
committed
viirs_data['M14'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M14.values.repeat(2, 0).repeat(2, 1))
viirs_data['M15'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M15.values.repeat(2, 0).repeat(2, 1))
viirs_data['M16'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M16.values.repeat(2, 0).repeat(2, 1))
Paolo Veglio
committed
viirs_data['M15-M16'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_data.M15.values - viirs_data.M16.values)
viirs_data['M14-M15'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_data.M14.values - viirs_data.M15.values)
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('I05', cmin_G1, test_name='11um_Test')
cmin_G1, bit2 = MyScene.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
Paolo Veglio
committed
cmin_G2, bit4 = MyScene.test_11_12um_diff('M15-M16', cmin_G2, test_name='11-12um_Cirrus_Test')
cmin_G2, bit3 = MyScene.oceanic_stratus_11_4um_test('I05-I04', cmin_G2,
test_name='11-4um_Oceanic_Stratus_Test')
Paolo Veglio
committed
# cmin_G2, bit_x = MyScene.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test')
Paolo Veglio
committed
cmin_G3, bit4_1 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
cmin_G3, bit5 = MyScene.vis_nir_ratio_test('I02-I01ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
Paolo Veglio
committed
# visnir, bit5 = MyScene.vis_nir_ratio_test('I02-I01ratio', np.ones(viirs_data.M15.shape), test_name='Vis/NIR_Ratio_Test')
cmin_G3, bit6 = 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.zeros(confidence.shape)
Paolo Veglio
committed
total_bit = bit1 + bit2 + bit4
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_confidence = np.sqrt(confidence * mvcm_conf.repeat(2, 0).repeat(2, 1))
date = mod02.split('.')[1]
time = mod02.split('.')[2]
# mvcm_confidence = mvcm_conf.repeat(2, 0).repeat(2, 1)
out_xr = xr.Dataset(
data_vars=dict(
hires_confidence=(['number_of_lines', 'number_of_pixels'], confidence),
ocean_day_mask=(['number_of_lines', 'number_of_pixels'], viirs_data.Ocean_Day.values),
),
coords=dict(
latitude=(['number_of_lines', 'number_of_pixels'], viirs_data.latitude.values),
longitude=(['number_of_lines', 'number_of_pixels'], viirs_data.longitude.values),
)
)
Paolo Veglio
committed
out_path = '/ships19/hercules/pveglio/mvcm_viirs_hires/outputs'
out_xr.to_netcdf(f'{out_path}/hires_MVCM.{date}.{time}.nc')
# out_xr.to_netcdf(out_fname)
# # np.savez(f'{data_path}/pyMVCM_{date}.{time}',
# np.savez(f'{data_path}/outputs/pyMVCM_{date}.{time}',
# confidence=total_confidence, py_conf=confidence, mvcm_conf=mvcm_conf,
# lat=viirs_data.latitude.values, lon=viirs_data.longitude.values,
# ocean_day_mask=viirs_data.Ocean_Day.values)
# # c1=cmin_G1, c2=cmin_G2, c3=cmin_G3)
# 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(data_path=sys.argv[1],
mod02=sys.argv[2],
mod03=sys.argv[3],
img02=sys.argv[4],
img03=sys.argv[5],
threshold_file=sys.argv[6],
geos_atm_1=sys.argv[7],
geos_atm_2=sys.argv[8],
geos_land=sys.argv[9],
geos_ocean=sys.argv[10],
geos_constants=sys.argv[11],
ndvi_file=sys.argv[12],
sst_file=sys.argv[13],
out_fname=sys.argv[14])