import ruamel_yaml as yml import numpy as np # import xarray as xr from glob import glob import read_data as rd import scene as scn # import tests import ocean_day_tests as odt import restoral # #################################################################### # # 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] # thresholds: _threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.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' mod_bands = ['M01', 'M02', 'M03', 'M04', 'M06', 'M08', 'M09', 'M11', 'M13', 'M14', 'M15', 'M16'] # #################################################################### # def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03, img02=_fname_img02, img03=_fname_img03, threshold_file=_threshold_file, geos_atm_1=_geos_atm_1, geos_atm_2=_geos_atm_2, geos_land=_geos_land, geos_ocean=_geos_ocean, geos_constants=_geos_constants, ndvi_file=_ndvi_file, sst_file=_sst_file): # datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input' # fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc' # fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc' # thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml' ancillary_file_names = {'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' } viirs_data = rd.read_data('viirs', f'{mod02}', f'{mod03}') viirs_hires = rd.read_data('viirs', f'{img02}', f'{img03}') for bnd in mod_bands: viirs_hires[bnd] = (['number_of_lines', 'number_of_pixels'], viirs_data[bnd].values.repeat(2, 0).repeat(2, 1)) viirs_hires = viirs_hires.rename({'I01': 'M05', 'I02': 'M07', 'I03': 'M10', 'I04': 'M12'}) print('VIIRS data read') viirs_hires = rd.read_ancillary_data(ancillary_file_names, viirs_hires, resolution=2) print('ancillary data read') viirs_data = viirs_hires with open(threshold_file) as f: text = f.read() thresholds = yml.safe_load(text) sunglint_angle = thresholds['Sun_Glint']['bounds'][3] scene_flags = scn.find_scene(viirs_data, sunglint_angle) cmin1 = np.ones(viirs_data.M01.shape) cmin2 = np.ones(viirs_data.M01.shape) cmin3 = np.ones(viirs_data.M01.shape) cmin4 = np.ones(viirs_data.M01.shape) c = np.ones((9, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1])) cmin1, c[0, :, :], bit1 = odt.simple_test(viirs_data.M15.values, thresholds['Daytime_Ocean']['bt11'], cmin1) # For this test I changed the sst to [7.000, 6.500, 6.000, 1.0, 1.0] cmin1, c[1, :, :], bit2 = odt.sst_test(viirs_data.M15.values, viirs_data.M16.values, viirs_data.sensor_zenith.values, viirs_data.geos_sfct.values, thresholds['Daytime_Ocean']['sst'], cmin1) cmin2, c[2, :, :], bit3 = odt.simple_test(viirs_data.M14.values-viirs_data.M15.values, thresholds['Daytime_Ocean']['diff_11_86um'], cmin2) cmin2, c[3, :, :], bit4 = odt.test_11_12_diff(viirs_data, thresholds['Daytime_Ocean']['diff11_12um'], cmin2) cmin2, c[4, :, :] = odt.test_11_4_diff(viirs_data.I05.values, viirs_data.M12.values, thresholds['Daytime_Ocean']['test11_4lo'], scene_flags, cmin2) cmin3, c[5, :, :] = odt.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'], thresholds['Sun_Glint'], viirs_data, cmin3) cmin3, c[6, :, :] = odt.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values, thresholds, scene_flags, cmin3) cmin3, c[7, :, :] = odt.nir_refl_test(viirs_data.M10.values, thresholds['Daytime_Ocean'], thresholds['Sun_Glint'], viirs_data, cmin3) total_bit = bit1 + bit2 + bit4 temp_confidence = cmin1 * cmin2 * cmin3 * cmin4 confidence = cmin1 * cmin2 * cmin3 * cmin4 # 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, confidence)[idx] 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] temp = np.zeros((viirs_data.M01.shape[0], viirs_data.M01.shape[1])) temp[idx] = 1 c[8, :, :] = temp np.savez('test_confidence', confidence=confidence, conf_test=c, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values) return confidence if __name__ == "__main__": main()