import ruamel_yaml as yml import numpy as np # import xarray as xr import read_data as rd import tests import ancillary_data as anc def main(): 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' anc_dir = '/ships19/hercules/pveglio/neige_data/snpp_test_input/ancillary/2014_08_01_213' sst_file = 'oisst.20140730' ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.209.hdf' geos_cnst = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4' geos_lnd = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20140801_1530.V01.nc4' geos_ocn = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20140801_1530.V01.nc4' geos1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1500.V01.nc4' geos2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1800.V01.nc4' start_time = '2014-08-01 15:48:00.000' viirs_data = rd.read_data('viirs', f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}') sst = np.empty((3232*3200, ), dtype=np.float32) ndvi = np.empty((3232*3200, ), dtype=np.float32) eco = np.empty((3232*3200, ), dtype=np.ubyte) geos_data = {'tpw': np.empty((3232*3200, ), dtype=np.float32), 'snowfr': np.empty((3232*3200, ), dtype=np.float32), 'icefr': np.empty((3232*3200, ), dtype=np.float32), 'ocnfr': np.empty((3232*3200, ), dtype=np.float32), 'landicefr': np.empty((3232*3200, ), dtype=np.float32), 'sfct': np.empty((3232*3200, ), dtype=np.float32), } tpw = np.ones((3232*3200, ), dtype=np.float32) snowfr = np.empty((3232*3200, ), dtype=np.float32) icefr = np.empty((3232*3200, ), dtype=np.float32) ocnfr = np.empty((3232*3200, ), dtype=np.float32) landicefr = np.empty((3232*3200, ), dtype=np.float32) sfct = np.empty((3232*3200, ), dtype=np.float32) sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )), viirs_data.longitude.values.reshape((3232*3200, )), anc_dir, sst_file, sst) ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((3232*3200, )), viirs_data.longitude.values.reshape((3232*3200, )), anc_dir, ndvi_file, ndvi) eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((3232*3200, )), viirs_data.longitude.values.reshape((3232*3200, )), anc_dir, eco) tpw = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )), viirs_data.longitude.values.reshape((3232*3200, )), start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst, tpw, snowfr, icefr, ocnfr, landicefr, sfct) print(tpw[:10]) #print(geos_data['tpw'][:10]) #print(geos_data['snowfr'][:10]) #print(geos_data['icefr'][:10]) #print(geos_data['ocnfr'][:10]) #print(geos_data['landicefr'][:10]) #print(geos_data['sfct'][:10]) with open(thresh_file) as f: text = f.read() thresholds = yml.safe_load(text) confidence = np.zeros((5, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1])) confidence[0, :, :] = tests.test_11um(viirs_data.M15.values, thresholds['Daytime_Ocean']) confidence[1, :, :] = tests.test_11_4diff(viirs_data.M15.values, viirs_data.M13.values, thresholds['Daytime_Ocean'], viirs_data, thresholds['Sun_Glint']['bounds'][3]) confidence[2, :, :] = tests.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'], thresholds['Sun_Glint'], viirs_data) # Note that here I'm using M05/M07 but the corresponding hi-res channelsa re I1/I2 # IMPORTANT: conf_test_dble() needs to be verified. I don't think it's working as intended at the moment confidence[3, :, :] = tests.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values, thresholds['Daytime_Ocean'], thresholds['Sun_Glint']) # This test needs to be verified, for the granule I'm running everything is zero confidence[4, :, :] = tests.test_11um_var(viirs_data.M15.values, thresholds['Nighttime_Ocean'], thresholds['Daytime_Ocean_Spatial_Variability']) return confidence def test_main(): rad1 = [[255, 260, 265, 248, 223], [278, 285, 270, 268, 256], [275, 273, 266, 254, 259]] rad2 = [[270, 273, 271, 268, 265], [277, 286, 275, 277, 269], [280, 281, 272, 270, 267]] thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml' with open(thresh_file) as f: text = f.read() rad1 = np.array(rad1) rad2 = np.array(rad2) confidence = np.zeros((2, rad1.shape[0], rad1.shape[1])) thresholds = yml.safe_load(text) confidence[0, :, :] = tests.test_11um(rad1, thresholds['Daytime_Ocean']) confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds['Daytime_Ocean']) print(f'Confidence[0,:,:]: \n {confidence[0, :, :]}') print(f'Confidence[1,:,:]: \n {confidence[1, :, :]}') return confidence if __name__ == "__main__": test_main()