Newer
Older
import ruamel_yaml as yml
import numpy as np
import xarray as xr
from glob import glob
import scene as scn
from tests import CloudTests
# import tests
import ocean_day_tests as odt
Paolo Veglio
committed
import restoral
Paolo Veglio
committed
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
Paolo Veglio
committed
_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]
Paolo Veglio
committed
_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]
Paolo Veglio
committed
# thresholds:
_threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
Paolo Veglio
committed
# 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'
Paolo Veglio
committed
_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'
Paolo Veglio
committed
# #################################################################### #
Paolo Veglio
committed
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'
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'
}
Paolo Veglio
committed
with open(threshold_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
viirs_data = rd.get_data(file_names, sunglint_angle)
# scene_xr = xr.Dataset()
# for s in scn._scene_list:
# scene_xr[s] = (('number_of_lines', 'number_of_pixels'), scn.scene_id[s])
# scene_xr['latitude'] = viirs_xr.latitude
# scene_xr['longitude'] = viirs_xr.longitude
#
# viirs_data = xr.Dataset(viirs_xr, coords=scene_xr)
# viirs_data.drop_vars(['latitude', 'longitude'])
cmin_G1 = np.ones(viirs_data.M01.shape)
cmin_G2 = np.ones(viirs_data.M01.shape)
# cmin_test = {'Ocean_Day': np.ones(viirs_data.M01.shape),
# 'Polar_Ocean_Day': np.ones(viirs_data.M01.shape),
# 'Polar_Ocean_Night': 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)
Land_Day = CloudTests(viirs_data, 'Land_Day', thresholds)
Land_Night = CloudTests(viirs_data, 'Land_Night', thresholds)
Land_Day_Coast = CloudTests(viirs_data, 'Land_Day_Coast', thresholds)
Land_Day_Desert = CloudTests(viirs_data, 'Land_Day_Desert', thresholds)
Land_Day_Desert_Coast = CloudTests(viirs_data, 'Land_Day_Desert_Coast', thresholds)
Ocean_Day = CloudTests(viirs_data, 'Ocean_Day', thresholds)
Ocean_Night = CloudTests(viirs_data, 'Ocean_Night', thresholds)
Polar_Day_Land = CloudTests(viirs_data, 'Polar_Day_Land', thresholds)
Polar_Night_Land = CloudTests(viirs_data, 'Polar_Night_Land', thresholds)
Polar_Day_Coast = CloudTests(viirs_data, 'Polar_Day_Coast', thresholds)
Polar_Day_Desert = CloudTests(viirs_data, 'Polar_Day_Desert', thresholds)
Polar_Day_Desert_Coast = CloudTests(viirs_data, 'Polar_Day_Desert_Coast', thresholds)
Polar_Day_Snow = CloudTests(viirs_data, 'Polar_Day_Snow', thresholds)
Polar_Night_Snow = CloudTests(viirs_data, 'Polar_Night_Snow', thresholds)
Polar_Ocean_Day = CloudTests(viirs_data, 'Polar_Ocean_Day', thresholds)
Polar_Ocean_Night = CloudTests(viirs_data, 'Polar_Ocean_Night', thresholds)
Day_Snow = CloudTests(viirs_data, 'Day_Snow', thresholds)
Night_Snow = CloudTests(viirs_data, 'Night_Snow', thresholds)
Antarctic_Day = CloudTests(viirs_data, 'Antarctic_Day', thresholds)
# 11um BT Test
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
cmin_G1 = Ocean_Day.single_threshold_test('11um_Test', 'M15', cmin_G1)
cmin_G1 = Polar_Ocean_Day.single_threshold_test('11um_Test', 'M15', cmin_G1)
cmin_G1 = Polar_Ocean_Night.single_threshold_test('11um_Test', 'M15', cmin_G1)
# CO2 High Clouds Test
# NOTE: VIIRS doesn't have the MODIS equivalent of B35 so this test is not performed
cmin_G1 = Land_Day.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Night.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Coast.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Desert.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Desert_Coast.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Ocean_Day.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Ocean_Night.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Day_Snow.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Night_Snow.single_threshold_test('CO2_High_Clouds_Test', 'bad_data', cmin_G1)
# Water Vapor High Clouds Test
# NOTE: VIIRS doesn't have the MODIS equivalent of B27 so this test is not performed
cmin_G1 = Land_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Desert.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Land_Day_Desert_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Ocean_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Ocean_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Day_Land.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Night_Land.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Day_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Day_Desert.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Day_Desert_Coast.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Day_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Night_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Ocean_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Polar_Ocean_Night.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Day_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Night_Snow.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
cmin_G1 = Antarctic_Day.single_threshold_test('Water_Vapor_High_Clouds_Test', 'bad_data', cmin_G1)
# Surface Temperature Test
# ## NOTE: This requires some calculations for the thresholds.
# Moreover this test is using the 11um - 12um difference instead of a single channel
# Also, look at the test carefully for these two cases. Polar_Night_Land uses some hardcoded coeffs
# (i.e. *_df1 and *_df2) that might be worth moving in the thresholds file for consistency
cmin_G1 = Land_Night.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)
cmin_G1 = Polar_Night_Land.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)
# SST Test
# ## NOTE: This requires some calculations for the thresholds.
# Moreover this test is using the 11um - 12um difference instead of a single channel
cmin_G1 = Ocean_Day.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
cmin_G1 = Polar_Ocean_Day.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
cmin_G1 = Polar_Ocean_Night.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
# the following test uses a different set of coefficients compared to the others above
cmin_G1 = Ocean_Night.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
# 11-8.6um BT Difference Test
cmin_G2 = Ocean_Day.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
cmin_G2 = Ocean_Night.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
cmin_G2 = Polar_Ocean_Day.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
cmin_G2 = Polar_Ocean_Night.single_threshold_test('11-8.6um_Test', 'M15-M14', cmin_G2)
# 11-12um BT Difference
cmin_G1 = Ocean_Day.single_threshold_test('11-12BT_diff',
viirs_data.M15.values-viirs_data.M16.values,
cmin_G1)
Ocean_Night = CloudTests(viirs_data, 'Ocean_Night', thresholds)
cmin_G1 = Ocean_Night.single_threshold_test('11um_test', viirs_data.M15.values, cmin_G1)
cmin_G1 = Ocean_Night.single_threshold_test('11-12um_diff',
viirs_data.M15.values-viirs_data.M16.values,
cmin_G1)
Paolo Veglio
committed
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)
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.M15.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)
#
# 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 channels are 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'])
Paolo Veglio
committed
total_bit = bit1 + bit2 + bit4
temp_confidence = cmin1 * cmin2 * cmin3 * cmin4
confidence = cmin1 * cmin2 * cmin3 * cmin4
Paolo Veglio
committed
# 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]
Paolo Veglio
committed
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
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__":