Skip to content
Snippets Groups Projects
hires_module.py 8.74 KiB
Newer Older
Paolo Veglio's avatar
Paolo Veglio committed
import numpy as np
import xarray as xr
import netCDF4 as nc
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 = ''

# MVCM output
_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))
    # 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['M07'] = viirs_data.I02
    viirs_data['M10'] = viirs_data.I03
    viirs_data['M12'] = viirs_data.I04
    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))
    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')

    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')
    # cmin_G2, bit_x = MyScene.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test')

    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')
    # 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)
    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),
                    )
                )
    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])