diff --git a/conf.py b/conf.py
index ecda20600c1b6648892f0e0783eb4f804dc86478..bc586eff3cbafaa5b2a22b3e9765d0573cf045ad 100644
--- a/conf.py
+++ b/conf.py
@@ -265,7 +265,7 @@ def conf_test_dble(rad, coeffs):
confidence[(alpha1 - gamma1 <= 0) & ((rad > gamma1) | (rad < gamma2))] = 0
confidence[(alpha1 - gamma1 <= 0) & (rad <= alpha1) & (rad >= alpha2)] = 1
- confidence = np.clip(confidence, 0, 1)
+ # confidence = np.clip(confidence, 0, 1)
return confidence.reshape(radshape)
diff --git a/hires_module.py b/hires_module.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f2193b85e830def28ba0f32ae43481d2b3014b6
--- /dev/null
+++ b/hires_module.py
@@ -0,0 +1,171 @@
+import numpy as np
+# import xarray as xr
+import netCDF4 as nc
+
+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
+
+# #################################################################### #
+# 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]
+
+# MVCM output
+_cld_mask = f'{_datapath}/CLDMSK_L2_VIIRS_SNPP.A2022173.1454.001.2022174035130.nc'
+
+# 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,
+ cloud_mask_file=_cld_mask) -> None:
+
+ f = nc.Dataset(cloud_mask_file)
+
+ cm = f['geophysical_data/Cloud_Mask'][:]
+
+ f.close()
+
+ with open(threshold_file) as f:
+ text = f.read()
+ thresholds = yml.safe_load(text)
+
+ 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'
+ }
+
+ sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
+
+ viirs_data = rd.get_data(file_names, sunglint_angle, hires=True)
+
+ 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, 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['M10'] = viirs_data.I03
+ viirs_data['M12'] = viirs_data.I04
+ viirs_data['M15'] = viirs_data.I05
+
+ 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('M15', cmin_G1, test_name='11um_Test')
+
+ cmin_G2, bit2 = MyScene.oceanic_stratus_11_4um_test('I05-I04', cmin_G2,
+ test_name='11-4um_Oceanic_Stratus_Test')
+
+ cmin_G3, bit3 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
+ cmin_G3, bit4 = MyScene.vis_nir_ratio_test('I02-I01ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
+ cmin_G3, bit5 = 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.ones(confidence.shape)
+
+ 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_bit = bit1
+
+ total_confidence = np.sqrt(confidence * mvcm_conf.repeat(2, 0).repeat(2, 1))
+
+ date = mod02.split('.')[1]
+ time = mod02.split('.')[2]
+
+ np.savez(f'{data_path}/output/pyMVCM_{date}.{time}',
+ confidence=total_confidence, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values,
+ ocean_day_mask=viirs_data.Ocean_Day.values)
+
+ # 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_fcn()
diff --git a/main.py b/main.py
index 2e0268febd8f5863d7b6c203d39f72868252bf68..db6d6e21dfcd39cfdcd50132ac327952aadb7f07 100644
--- a/main.py
+++ b/main.py
@@ -8,7 +8,10 @@ import numpy as np
from glob import glob
import read_data as rd
+import spectral_tests as tst
from tests import CloudTests
+import scene as scn
+import restoral
# #################################################################### #
# TEST CASE
@@ -114,6 +117,62 @@ def main(*,
cmin_G4 = np.ones(viirs_data.M01.shape)
cmin_G5 = np.ones(viirs_data.M01.shape)
+ ##########################################################
+
+ """
+ for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert_Coast',
+ 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean',
+ 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
+ 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Land_Night',
+ 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']:
+ """
+ for scene_name in ['Ocean_Day']:
+ MyScene = tst.CloudTests(viirs_data, scene_name, thresholds)
+
+ cmin_G1, bit1 = MyScene.test_11um('M15', cmin_G1, test_name='11um_Test')
+ cmin_G1, bit2 = MyScene.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
+
+ cmin_G2, bit3 = MyScene.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test')
+ cmin_G2, bit4 = MyScene.test_11_12um_diff('M15-M16', cmin_G2, test_name='11-12um_Cirrus_Test')
+ cmin_G2, bit5 = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2,
+ test_name='11-4um_Oceanic_Stratus_Test')
+
+ cmin_G3, bit6 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
+ cmin_G3, bit7 = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
+ cmin_G3, bit8 = MyScene.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
+ cmin_G3, bit9 = MyScene.visible_reflectance_test('M128', cmin_G3, test_name='Visible_Reflectance_Test')
+ cmin_G3, bit10 = MyScene.gemi_test('GEMI', cmin_G3, test_name='GEMI_Test')
+
+ cmin_G4, bit11 = MyScene.test_1_38um_high_clouds('M09', cmin_G4, test_name='1.38um_High_Cloud_Test')
+
+ cmin_G5, bit12 = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5,
+ test_name='4-12um_BTD_Thin_Cirrus_Test')
+
+ cmin = cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5
+
+ if scene_name in ['Ocean_Day']:
+ # total_bit = np.full(viirs_data.M01.shape, 4)
+ total_bit = bit1 + bit2 + bit4 + bit11
+ 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) & (cmin <= 0.99) & (cmin >= 0.05))
+ # cmin[idx] = restoral.spatial(viirs_data, thresholds['Sun_Glint'], scene_flags, cmin)[idx]
+
+ idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
+ (scene_flags['uniform'] == 1) & (cmin <= 0.95))
+ cmin[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, cmin)[idx]
+
+ # MVCM = tst.ComputeTests(viirs_data, scene_name, thresholds)
+ # cmin = MVCM.clear_sky_restoral(cmin)
+
+ # return cmin
+ np.savez('test_ams', confidence=cmin, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values)
+
+ ##########################################################
+
+ """
# ------------------- #
# ### Testing Setup ###
# ------------------- #
@@ -283,6 +342,7 @@ def main(*,
cmin_total = cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5
return cmin_total
+ """
'''
Land_Day = CloudTests(viirs_data, 'Land_Day', thresholds)
@@ -465,4 +525,4 @@ def test_main():
if __name__ == "__main__":
- test_main()
+ main()
diff --git a/plot_data.py b/plot_data.py
index 80f6107cc2dd3a874b39cee63e2df70d9e3a1f33..31358cd950759b03bc3c17a17980a155f8133f39 100644
--- a/plot_data.py
+++ b/plot_data.py
@@ -1,4 +1,5 @@
import matplotlib.pyplot as plt
+import matplotlib
import cartopy.crs as ccrs
import numpy as np
import netCDF4 as nc
@@ -8,6 +9,8 @@ import os
from scene import test_scene
+matplotlib.use('Agg')
+
def main(scene_flag):
scn = test_scene()
@@ -78,6 +81,11 @@ def plot_confidence(numpy_file='test_confidence.npz'):
lon = data['lon']
# confidence = data['confidence']
confidence = data['confidence']
+ cm = np.zeros(confidence.shape)
+ cm[confidence > 0.66] = 1
+ cm[confidence > 0.95] = 2
+ cm[confidence > 0.99] = 3
+
plt.figure(figsize=[15, 10])
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([np.min(lon)-5, np.max(lon)+5,
@@ -85,26 +93,32 @@ def plot_confidence(numpy_file='test_confidence.npz'):
ax.coastlines(resolution='50m', color='black', linewidth=1)
- plt.pcolormesh(lon, lat, confidence, vmin=0, vmax=1,
+ plt.pcolormesh(lon, lat, cm, vmin=0, vmax=3,
transform=ccrs.PlateCarree())
- plt.title('Confidence')
- plt.colorbar()
+ # plt.title('Confidence')
+ # plt.colorbar()
- plt.savefig('figures/total_confidence.png')
+ plt.title('Cloud Mask')
+ cb = plt.colorbar(ticks=range(4))
+ cb.set_ticklabels(['confident cloudy', 'probably cloudy', 'probably clear', 'confident clear'])
+ plt.savefig('figures/python_cloudmask.png', dpi=300)
def plot_cloud_mask():
- f = nc.Dataset('/ships19/hercules/pveglio/mvcm_viirs_hires/CLDMSK_L2_VIIRS_SNPP.A2022173.1454.001.2022174035130.nc')
+ # f = nc.Dataset('/ships19/hercules/pveglio/mvcm_viirs_hires/CLDMSK_L2_VIIRS_SNPP.A2022173.1454.001.2022174035130.nc')
+ f = nc.Dataset('/ships19/hercules/pveglio/mvcm_viirs_hires/outputs/CLDMSK_L2_MODIS_SNPP.A2022173.1454.001_viirs_cris.nc')
lat = f['geolocation_data/latitude'][:]
lon = f['geolocation_data/longitude'][:]
# icm = f['geophysical_data/Integer_Cloud_Mask'][:]
cm = f['geophysical_data/Cloud_Mask'][:]
qa = f['geophysical_data/Quality_Assurance'][:]
+ sst = f['geolocation_data/SST'][:]
f.close()
# cm_flag = np.array((np.array(cm[0, :, :], np.byte) &
- # (pow(2, 1) + pow(2, 2)))/pow(2, 1))
- cm_flag = read_bits(qa[:, :, 2], 2) & ~read_bits(cm[2, :, :], 2)
+ # (pow(2, 1) + pow(2, 2)))/pow(2, 1))
+ # cm_flag = read_bits(qa[:, :, 2], 2) & ~read_bits(cm[2, :, :], 2)
+ cm_flag = read_bits(qa[:, :, 2], 0) & ~read_bits(cm[2, :, :], 0)
plt.figure(figsize=[15, 10])
ax = plt.axes(projection=ccrs.PlateCarree())
@@ -119,8 +133,8 @@ def plot_cloud_mask():
# cb.set_ticklabels(['confident cloudy', 'probably cloudy', 'probably clear', 'confident clear'])
cb = plt.colorbar(ticks=range(2))
cb.set_ticklabels(['No', 'Yes'])
- # plt.savefig('figures/ref_cloudmask.png')
- plt.savefig('figures/high_cloud_flag_11_12.png')
+ # plt.savefig('figures/ref_cloudmask.png', dpi=300)
+ plt.savefig('figures/1_38um_high_cloud_test.png', dpi=300)
def read_bits(byte, bits):
@@ -142,6 +156,6 @@ def read_bits(byte, bits):
if __name__ == "__main__":
# main(sys.argv[1])
- plot_tests(sys.argv[1])
- # plot_confidence(sys.argv[1])
+ # plot_tests(sys.argv[1])
+ plot_confidence(sys.argv[1])
# plot_cloud_mask()
diff --git a/plot_mask.py b/plot_mask.py
index a7a8fd67544e393701608e28d62fc24c60b54e67..ff5647678b77967bb8861a4b3a1d2712ce2cad99 100644
--- a/plot_mask.py
+++ b/plot_mask.py
@@ -256,7 +256,7 @@ def read_bits(byte, bits):
parser = argparse.ArgumentParser()
parser.add_argument('input')
-#parser.add_argument('test_name')
+parser.add_argument('test_name')
parser.add_argument('-c1', default='mediumblue')
parser.add_argument('-c2', default='skyblue')
parser.add_argument('-c3', default='yellowgreen')
@@ -287,7 +287,7 @@ if __name__ == "__main__":
# print(n)
main(args.input,
# n,
-# args.test_name,
+ args.test_name,
c1=args.c1,
c2=args.c2,
c3=args.c3,
diff --git a/read_data.py b/read_data.py
index e43faf37a18357100110ca1b97796a3e7fbf7a8b..1d89d4cc5959511890626a11a1f4fd76cb6bfa78 100644
--- a/read_data.py
+++ b/read_data.py
@@ -140,33 +140,51 @@ def read_ancillary_data(file_names: str,
def get_data(file_names: Dict[str, str],
- sunglint_angle: float) -> xr.Dataset:
+ sunglint_angle: float,
+ hires: bool = False) -> xr.Dataset:
mod02 = file_names['MOD02']
mod03 = file_names['MOD03']
+ img02 = file_names['IMG02']
+ img03 = file_names['IMG03']
- viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
- viirs_data = read_ancillary_data(file_names, viirs_data)
+ if hires is False:
+ viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
+ viirs_data = read_ancillary_data(file_names, viirs_data)
- m01 = viirs_data.M05.values
- m02 = viirs_data.M07.values
- r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01)
- r2 = m02 + m01 + 0.5
- r3 = r1 / r2
- gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01))
+ m01 = viirs_data.M05.values
+ m02 = viirs_data.M07.values
+ r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01)
+ r2 = m02 + m01 + 0.5
+ r3 = r1 / r2
+ gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01))
# Compute channel differences and ratios that are used in the tests
- viirs_data['M14-M15'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M14.values - viirs_data.M15.values)
- viirs_data['M15-M16'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M15.values - viirs_data.M16.values)
- viirs_data['M15-M13'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M15.values - viirs_data.M13.values)
- viirs_data['M13-M16'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M13.values - viirs_data.M16.values)
- viirs_data['M07-M05ratio'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M07.values / viirs_data.M05.values)
- viirs_data['GEMI'] = (('number_of_lines', 'number_of_pixels'), gemi)
+ viirs_data['M14-M15'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M14.values - viirs_data.M15.values)
+ viirs_data['M15-M16'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M15.values - viirs_data.M16.values)
+ viirs_data['M15-M12'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M15.values - viirs_data.M12.values)
+ viirs_data['M13-M16'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M13.values - viirs_data.M16.values)
+ viirs_data['M07-M05ratio'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M07.values / viirs_data.M05.values)
+ viirs_data['GEMI'] = (('number_of_lines', 'number_of_pixels'), gemi)
+
+ # temp value to force the code to work
+ viirs_data['M128'] = (('number_of_lines', 'number_of_pixels'), np.zeros(viirs_data.M01.shape))
+
+ else:
+ viirs_data = read_data('viirs', f'{img02}', f'{img03}')
+ viirs_data['M05'] = viirs_data.I01
+ viirs_data['M07'] = viirs_data.I02
+ viirs_data = read_ancillary_data(file_names, viirs_data)
+
+ viirs_data['I05-I04'] = (('number_of_lines_2', 'number_of_pixels_2'),
+ viirs_data.I05.values - viirs_data.I04.values)
+ viirs_data['I02-I01ratio'] = (('number_of_lines_2', 'number_of_pixels_2'),
+ viirs_data.I02.values / viirs_data.I01.values)
scene_flags = scn.find_scene(viirs_data, sunglint_angle)
scene = scn.scene_id(scene_flags)
diff --git a/restoral.py b/restoral.py
index 4438e8d7a60162c87e10da0921c1f2e2646de384..f01d7e88d1462f301c2c5d4834f7154978006bdb 100644
--- a/restoral.py
+++ b/restoral.py
@@ -24,7 +24,7 @@ def sunglint(viirs_data, threshold, bit, conf):
m02 = viirs_data.M07.values
irclr = np.ones(viirs_data.M02.shape)
- irclr[bit < 3] = 0
+ # irclr[bit < 1] = 0
var = spatial_var(m31, 0.4)
reg_var_mean = sliding_window_view(np.pad(m02, [1, 1], mode='constant'),
diff --git a/spectral_tests.py b/spectral_tests.py
index 2bbb412cf704479088a68dc99170a030e95bdec3..4c2962b35545a8bc70ec0c5b015c4a92a275f567 100644
--- a/spectral_tests.py
+++ b/spectral_tests.py
@@ -45,7 +45,8 @@ class CloudTests(object):
if test_name not in self.thresholds[self.scene_name]:
# print('Not running test for this scene')
# returns cmin. This could be changed into a keyworded argument for readability
- return args[-1]
+ test_bit = np.zeros(args[-1].shape)
+ return args[-1], test_bit
return func(self, *args, test_name, **kwargs)
return wrapper
@@ -71,7 +72,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, test_bit
@run_if_test_exists_for_scene
def surface_temperature_test(self,
@@ -89,9 +90,12 @@ class CloudTests(object):
test_name: str = 'SST_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band31].shape)
+ test_bit = np.zeros(self.data[band31].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
m31 = self.data[band31].values - 273.16
bt_diff = self.data[band31].values - self.data[band32].values
sst = self.data.geos_sfct.values - 273.16
@@ -99,15 +103,18 @@ class CloudTests(object):
c = threshold['coeffs']
- modsst = 273.16 + c[0] + c[1]*m31 + c[2]*bt_diff*sst + c[3]*bt_diff*((1/cosvza) - 1)
+ modsst = 273.16 + c[0] + c[1]*m31 + c[2]*bt_diff*sst + c[3]*bt_diff*((1.0/cosvza) - 1.0)
sfcdif = self.data.geos_sfct.values - modsst
+ idx = np.nonzero((sfcdif < threshold['thr'][1]) &
+ (self.data[self.scene_name] == 1))
+ test_bit[idx] = 1
print(f'Testing "{self.scene_name}"\n')
confidence[self.scene_idx] = conf.conf_test_new(sfcdif[self.scene_idx], threshold['thr'])
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def bt_diff_86_11um(self,
@@ -116,16 +123,22 @@ class CloudTests(object):
test_name: str = '8.6-11um_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ idx = np.nonzero((self.data[band].values < threshold['thr'][1]) &
+ (self.data[self.scene_name] == 1))
+ test_bit[idx] = 1
confidence[self.scene_idx] = conf.conf_test_new(rad, threshold['thr'])
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def test_11_12um_diff(self,
@@ -134,17 +147,25 @@ class CloudTests(object):
test_name: str = '11-12um_Cirrus_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
thr = preproc.thresholds_11_12um(self.data, threshold, self.scene_name, self.scene_idx)
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ idx = np.nonzero((rad <= thr[1, :]) &
+ (self.data[self.scene_name].values[self.scene_idx] == 1))
+ tmp_bit = test_bit[self.scene_idx]
+ tmp_bit[idx] = 1
+ test_bit[self.scene_idx] = tmp_bit
confidence[self.scene_idx] = conf.conf_test_new(rad, thr)
cmin = np.fmin(cmin, confidence)
-
- return cmin
+ test_bit = test_bit.reshape(self.data[band].shape)
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def oceanic_stratus_11_4um_test(self,
@@ -153,9 +174,12 @@ class CloudTests(object):
test_name: str = '11-4um_Oceanic_Stratus_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
thr = threshold['thr']
@@ -168,11 +192,18 @@ class CloudTests(object):
'Day_Snow', 'Night_Snow', 'Antarctic_Day']:
pass
else:
- confidence[self.scene_idx] = conf.conf_test_new(rad, thr)
+ scene_flag = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3])
+ idx = np.nonzero((self.data[band].values >= threshold['thr'][1]) &
+ (self.data[self.scene_name] == 1) &
+ (scene_flag['sunglint'] == 0))
+ test_bit[idx] = 1
+ idx = np.nonzero((self.data[self.scene_name] == 1) &
+ (scene_flag['sunglint'] == 0))
+ confidence[idx] = conf.conf_test_new(self.data[band].values[idx], thr)
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def nir_reflectance_test(self,
@@ -181,18 +212,26 @@ class CloudTests(object):
test_name: str = 'NIR_Reflectance_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
thr = preproc.thresholds_NIR(self.data, self.thresholds, self.scene_name,
test_name, self.scene_idx)
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ idx = np.nonzero((rad <= thr[1, :]) &
+ (self.data[self.scene_name].values[self.scene_idx] == 1))
+ tmp_bit = test_bit[self.scene_idx]
+ tmp_bit[idx] = 1
+ test_bit[self.scene_idx] = tmp_bit
confidence[self.scene_idx] = conf.conf_test_new(rad, thr)
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def vis_nir_ratio_test(self,
@@ -201,6 +240,8 @@ class CloudTests(object):
test_name: str = 'Vis/NIR_Reflectance_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
@@ -208,27 +249,34 @@ class CloudTests(object):
# confidence array and I don't want to think of a solution at the moment. I will need
# to figure out the logic to make it work once I'm at the stage where I want to optimize
# the code
+ qa_bit[self.scene_idx] = 1
rad = self.data[band].values # [self.scene_idx]
solar_zenith = self.data.solar_zenith.values # [self.scene_idx]
- sensor_zenith = self.data.sensor_zenith.values # [self.scene_idx]
- rel_azimuth = self.data.relative_azimuth.values # [self.scene_idx]
- sunglint_angle = self.data.sunglint_angle.values # [self.scene_idx]
-
- cos_refang = np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * \
- np.cos(rel_azimuth*_DTR) + np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR)
- refang = np.arccos(cos_refang) * _RTD
+ sunglint = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3])['sunglint']
- idx = np.nonzero((solar_zenith <= 85) & (refang <= sunglint_angle))
+ idx = np.nonzero((solar_zenith <= 85) & (sunglint == 1))
thr_no_sunglint = np.array([threshold['thr'][i] for i in range(8)])
thr_sunglint = np.array([self.thresholds['Sun_Glint']['snglnt'][i] for i in range(8)])
-
- confidence = conf.conf_test_dble(rad, thr_no_sunglint)
+ print(thr_no_sunglint)
+ print(thr_sunglint)
+ # tmp = self.thresholds['Sun_Glint']['snglnt']
+ # thr_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1])
+ confidence[sunglint == 0] = conf.conf_test_dble(rad[sunglint == 0], thr_no_sunglint)
confidence[idx] = conf.conf_test_dble(rad[idx], thr_sunglint)
+ idx = np.nonzero(((rad < thr_no_sunglint[1]) | (rad > thr_no_sunglint[4])) &
+ (self.data[self.scene_name].values == 1) & (sunglint == 0))
+ test_bit[tuple(idx)] = 1
+
+ idx = np.nonzero(((rad < thr_sunglint[1]) | (rad > thr_sunglint[4])) &
+ ((self.data[self.scene_name].values == 1) &
+ (solar_zenith <= 85) & (sunglint == 1)))
+ test_bit[tuple(idx)] = 1
+
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def test_16_21um_Reflectance(self,
@@ -237,18 +285,26 @@ class CloudTests(object):
test_name: str = '1.6_2.1um_NIR_Reflectance_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
thr = preproc.thresholds_NIR(self.data, self.thresholds, self.scene_name,
test_name, self.scene_idx)
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ idx = np.nonzero((rad <= thr[1, :]) &
+ (self.data[self.scene_name].values[self.scene_idx] == 1))
+ tmp_bit = test_bit[self.scene_idx]
+ tmp_bit[idx] = 1
+ test_bit[self.scene_idx] = tmp_bit
confidence[self.scene_idx] = conf.conf_test_new(rad, thr)
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def visible_reflectance_test(self,
@@ -257,6 +313,7 @@ class CloudTests(object):
test_name: str = 'Visible_Reflectance_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
@@ -267,7 +324,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, test_bit
@run_if_test_exists_for_scene
def gemi_test(self,
@@ -276,6 +333,7 @@ class CloudTests(object):
test_name: str = 'GEMI_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
@@ -286,7 +344,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, test_bit
@run_if_test_exists_for_scene
def test_1_38um_high_clouds(self,
@@ -295,22 +353,31 @@ class CloudTests(object):
test_name: str = '1.38um_High_Cloud_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
if self.scene_name in ['Ocean_Day', 'Polar_Day_Ocean']:
thr = preproc.thresholds_1_38um_test(self.data, self.thresholds, self.scene_name,
self.scene_idx)
else:
+ return cmin, test_bit
thr = threshold['thr']
print(f'Testing "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ idx = np.nonzero((rad <= thr[1, :]) &
+ (self.data[self.scene_name].values[self.scene_idx] == 1))
+ tmp_bit = test_bit[self.scene_idx]
+ tmp_bit[idx] = 1
+ test_bit[self.scene_idx] = tmp_bit
confidence[self.scene_idx] = conf.conf_test_new(rad, thr)
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def thin_cirrus_4_12um_BTD_test(self,
@@ -319,6 +386,7 @@ class CloudTests(object):
test_name: str = '4-12um_BTD_Thin_Cirrus_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
+ test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
if (threshold['perform'] is True and self.pixels_in_scene is True):
@@ -329,7 +397,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin
+ return cmin, test_bit
def single_threshold_test(self, test_name, band, cmin):
@@ -475,7 +543,7 @@ class ComputeTests(CloudTests):
# Group 3
cmin_G3 = self.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
- cmin_G3 = self.vis_nir_ratio_test('M07M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
+ cmin_G3 = self.vis_nir_ratio_test('M07-M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
cmin_G3 = self.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
cmin_G3 = self.visible_reflectance_test('M128', cmin_G3, test_name='Visible_Reflectance_Test')
cmin_G3 = self.gemi_test('GEMI', cmin_G3, test_name='GEMI_Test')
@@ -506,6 +574,7 @@ class ComputeTests(CloudTests):
(scene_flags['uniform'] == 1) & (cmin <= 0.95))
cmin[idx] = restoral.sunglint(self.data, self.thresholds['Sun_Glint'], total_bit, cmin)[idx]
+ """
idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) &
(scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) &
(cmin <= 0.95))
@@ -515,6 +584,7 @@ class ComputeTests(CloudTests):
(scene_flags['coast'] == 1) & (scene_flags['snow'] == 0) &
(scene_flags['ice'] == 0))
cmin[idx] = restoral.coast(self.data, self.thresholds, scene_flags, cmin)[idx]
+ """
return cmin
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 7c2565cddef859b779b751fcb488c3c673820bbc..5f5817e059e597aeae0a67b8f54210d3febe261d 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -148,7 +148,7 @@ Ocean_Day:
perform: True
SST_Test:
thr: [3.000, 2.500, 2.000, 1.0, 1.0]
- coeffs: [1.886, 0.938, 0.128, 1.094]
+ coeffs: [1.8860, 0.9380, 0.1280, 1.094]
perform: True
8.6-11um_Test:
thr: [-0.50, -1.00, -1.50, 1.0, 1.0]