From c2c16d330aa6e1229cc9d55453b89feee8c5af34 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Mon, 18 Jul 2022 21:51:21 +0000
Subject: [PATCH] changes to main code to simulate the group testing found in
MVCM
---
main.py | 106 ++++++++++------
ocean_day_tests.py | 305 +++++++++++++++++++++++++++++++++++++++++++++
read_data.py | 2 +-
tests.py | 2 +
utils.py | 2 +-
5 files changed, 378 insertions(+), 39 deletions(-)
create mode 100644 ocean_day_tests.py
diff --git a/main.py b/main.py
index 6c7c0ed..0541ac5 100644
--- a/main.py
+++ b/main.py
@@ -2,29 +2,34 @@ import ruamel_yaml as yml
import numpy as np
# import xarray as xr
+from glob import glob
+
import read_data as rd
-import tests
+import scene as scn
+# import tests
+import ocean_day_tests as odt
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
-_fname_mod02 = 'VNP02MOD.A2022173.0336.001.2022173154147.uwssec.nc'
-_fname_mod03 = 'VNP03MOD.A2022173.0336.001.2022173155356.uwssec.nc'
-_fname_img02 = 'VNP02IMG.A2022173.0336.001.2022173154147.uwssec.nc'
-_fname_img03 = 'VNP03IMG.A2022173.0336.001.2022173155356.uwssec.nc'
+_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec.nc')[0]
+_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0]
+_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec.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_0300.V01.nc4'
-_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_0600.V01.nc4'
-_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_0330.V01.nc4'
-_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_0330.V01.nc4'
+_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'
# #################################################################### #
@@ -39,17 +44,17 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
# 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'{data_path}/ancillary/{geos_atm_1}',
- 'GEOS_atm_2': f'{data_path}/ancillary/{geos_atm_2}',
- 'GEOS_land': f'{data_path}/ancillary/{geos_land}',
- 'GEOS_ocean': f'{data_path}/ancillary/{geos_ocean}',
- 'GEOS_constants': f'{data_path}/ancillary/{geos_constants}',
- 'NDVI': f'{data_path}/ancillary/{ndvi_file}',
- 'SST': f'{data_path}/ancillary/{sst_file}',
+ 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'{data_path}/{mod02}', f'{data_path}/{mod03}')
+ viirs_data = rd.read_data('viirs', f'{mod02}', f'{mod03}')
viirs_data = rd.read_ancillary_data(ancillary_file_names, viirs_data)
@@ -57,26 +62,53 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
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 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'])
-
- return confidence
+ 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)
+
+ # confidence = np.zeros((8, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1]))
+ cmin1 = odt.simple_test(viirs_data.M15.values, thresholds['Daytime_Ocean']['bt11'], cmin1)
+ cmin1 = 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 = odt.simple_test(viirs_data.M14.values-viirs_data.M15.values,
+ thresholds['Daytime_Ocean']['diff_11_86um'], cmin2)
+ # PLACEHOLDER FOR 11-12 MICRON BTD TEST UNTIL I INCLUDE THE CITHR() FUNCTION FROM C IN CYTHON
+ cmin2 = odt.test_11_12_diff(viirs_data.M15.values, viirs_data.M16.values,
+ thresholds['Daytime_Ocean']['diff11_12um'], cmin2)
+ cmin2 = odt.test_11_4_diff(viirs_data.M15.values, viirs_data.M12.values,
+ thresholds['Daytime_Ocean']['test11_4lo'], cmin2)
+ cmin3 = odt.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
+ thresholds['Sun_Glint'], viirs_data, cmin3)
+ cmin3 = odt.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values,
+ thresholds['Daytime_Ocean'], thresholds['Sun_Glint'], cmin3)
+ cmin3 = 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'])
+
+ confidence = cmin1 * cmin2 * cmin3 * cmin4
+
+ return cmin1, cmin2, cmin3
def test_main():
diff --git a/ocean_day_tests.py b/ocean_day_tests.py
new file mode 100644
index 0000000..e0f1db5
--- /dev/null
+++ b/ocean_day_tests.py
@@ -0,0 +1,305 @@
+import numpy as np
+
+from numpy.lib.stride_tricks import sliding_window_view
+
+import utils
+import conf
+
+
+# ############## GROUP 1 TESTS ############## #
+
+# 11 micron brightness temperature threshold test
+def simple_test(rad, threshold, cmin):
+
+ radshape = rad.shape
+ rad = rad.reshape(np.prod(radshape))
+
+ thr = np.array(threshold)
+ confidence = np.ones(rad.shape)
+
+ if thr[4] == 1:
+ print("simple test running")
+ # the C code has the line below that I don't quite understand the purpose of.
+ # It seems to be setting the bit to 0 if the BT value is greater than the midpoint
+ #
+ # if (m31 >= dobt11[1]) (void) set_bit(13, pxout.testbits);
+
+ # confidence = utils.conf_test(rad, thr)
+ confidence = conf.conf_test(rad, thr)
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def sst_test(rad1, rad2, vza, surf_temp, threshold, cmin):
+
+ a1 = 1.8860
+ a2 = 0.9380
+ a3 = 0.1280
+ a4 = 1.0940
+
+ radshape = rad1.shape
+ b31 = rad1.reshape(np.prod(radshape)) - 273.16
+ b32 = rad2.reshape(np.prod(radshape)) - 273.16
+
+ thr = np.array(threshold)
+ confidence = np.ones(b31.shape)
+
+ rad_diff = b31 - b32
+ sstc = surf_temp.reshape(np.prod(radshape)) - 273.16
+ mu = np.cos(vza.reshape(np.prod(radshape)) * np.pi/180.0)
+
+ modsst = 273.16 + a1 + a2*b31 + a3*rad_diff*sstc + a4*rad_diff*((1/mu)-1)
+ sfcdif = surf_temp.reshape(np.prod(radshape)) - modsst
+
+ if thr[4] == 1:
+ print('SST test running')
+ confidence = conf.conf_test(sfcdif, thr)
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def test_11_12_diff(rad1, rad2, threshold, cmin):
+ radshape = rad1.shape
+ b31 = rad1.reshape(np.prod(radshape))
+ b32 = rad2.reshape(np.prod(radshape))
+ thr = np.array(threshold)
+
+ confidence = np.ones(b31.shape)
+ rad_diff = b31 - b32
+
+ if thr[1] == 1:
+ print('11-12um diff test running')
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def test_11_4_diff(rad1, rad2, threshold, cmin):
+ radshape = rad1.shape
+ b31 = rad1.reshape(np.prod(radshape))
+ b20 = rad2.reshape(np.prod(radshape))
+ thr = np.array(threshold)
+
+ confidence = np.ones(b31.shape)
+
+ if thr[4] == 1:
+ print('11-4um diff test running')
+ confidence = conf.conf_test(b31-b20, thr)
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def vis_nir_ratio_test(rad1, rad2, threshold, sg_threshold, cmin):
+ print("NIR-Visible ratio test running")
+ if threshold['vis_nir_ratio'][6] == 1:
+
+ radshape = rad1.shape
+ rad1 = rad1.reshape(np.prod(radshape))
+ rad2 = rad2.reshape(np.prod(radshape))
+
+ vrat = rad2/rad1
+
+ thresh = np.zeros((7,))
+
+ # temp value to avoid linter bitching at me
+ # eventually we would have the test run in two blocks as:
+ # confidence[sunglint == 1] = conf.conf_test_dble(vrat[sunglint == 1], sg_threshold['snglnt'])
+ # confidence[sunglint == 0] = conf.conf_test_dble(vrat[sunglint == 0], threshold['vis_nir_ratio'])
+ # sunglint needs to be defined somewhere
+ sunglint = 0
+ if sunglint:
+ thresh = threshold['snglnt']
+ else:
+ thresh = threshold['vis_nir_ratio']
+
+ confidence = conf.conf_test_dble(vrat, thresh)
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data, cmin):
+
+ print("NIR reflectance test running")
+ sza = viirs_data.solar_zenith.values
+ refang = viirs_data.sunglint_angle.values
+ vza = viirs_data.sensor_zenith.values
+ dtr = np.pi/180
+ # Keep in mind that band_n uses MODIS band numbers (i.e. 2=0.86um and 7=2.1um)
+ # For VIIRS would be 2=M07 (0.865um) and 7=M11 (2.25um)
+ band_n = 2
+ vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE
+
+ radshape = rad.shape
+ rad = rad.reshape(np.prod(radshape))
+ confidence = np.ones(rad.shape)
+ sza = sza.reshape(rad.shape)
+ vza = vza.reshape(rad.shape)
+ refang = refang.reshape(rad.shape)
+ sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(rad.shape)
+
+ # ref2 [5]
+ # b2coeffs [4]
+ # b2mid [1]
+ # b2bias_adj [1]
+ # b2lo [1]
+ # vzcpow [3] (in different place)
+
+ cosvza = np.cos(vza*dtr)
+ coeffs = threshold['b2coeffs']
+ hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3))
+ hicut0 = (hicut0 * 0.01) + threshold['b2adj']
+ hicut0 = hicut0 * threshold['b2bias_adj']
+ midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj'])
+ locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj'])
+ thr = np.array([locut0, midpt0, hicut0, threshold['ref2'][3]*np.ones(rad.shape)])
+ # corr_thr = np.zeros((4, 4))
+ corr_thr = np.zeros((4, rad.shape[0]))
+
+ corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0], vzcpow))
+ corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
+
+ for flag in range(1, 4):
+ if len(refang[sunglint_flag == flag]) > 0:
+ sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
+ corr_thr[:3, sunglint_flag == flag] = sunglint_thr[:3, sunglint_flag == flag] * (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
+ corr_thr[3, sunglint_flag == flag] = sunglint_thr[3, sunglint_flag == flag]
+
+ confidence = conf.conf_test(rad, corr_thr)
+
+ return np.minimum(cmin, confidence.reshape(radshape))
+
+
+def nir_high_cloud_test():
+ pass
+
+
+
+
+
+
+
+
+
+
+
+def test_11um_var(rad, threshold, var_threshold):
+
+ print("11um variability test running")
+ thr = np.array(threshold['11um_var'])
+
+ radshape = rad.shape
+ var = np.zeros((radshape[0], radshape[1], 9))
+
+ # chk_spatial2() need to figure out what this is
+ # np = rg_var.num_small_diffs * 1.0
+ test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))
+
+ var[np.abs(test).reshape(radshape[0], radshape[1], 9) < var_threshold['dovar11']] = 1
+ var = var.sum(axis=2).reshape(np.prod(radshape))
+
+ rad = rad.reshape(np.prod(radshape))
+ confidence = np.zeros(rad.shape)
+
+ confidence[var == 9] = conf.conf_test(rad[var == 9], thr)
+
+ return confidence.reshape(radshape)
+
+
+def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh):
+
+ print("11um - 4um difference test running")
+ radshape = rad1.shape
+ raddiff = (rad1 - rad2).reshape(np.prod(radshape))
+
+ day = np.zeros(radshape)
+ day[viirs_data.solar_zenith <= 85] = 1
+ day = day.reshape(raddiff.shape)
+ sunglint = np.zeros(rad1.shape)
+ sunglint[viirs_data.sunglint_angle <= sg_thresh] = 1
+ sunglint = sunglint.reshape(raddiff.shape)
+ thr = np.array(threshold['test11_4lo'])
+ confidence = np.zeros(raddiff.shape)
+
+ # confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
+ confidence[(day == 1) & (sunglint == 0)] = conf.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
+
+ return confidence.reshape(radshape)
+
+
+def vir_refl_test(rad, threshold, viirs_data):
+
+ print('Visible reflectance test running')
+
+ thr = threshold['vis_refl_test']
+
+ radshape = rad.shape()
+ rad = rad.reshape(np.prod(radshape))
+ confidence = np.zeros(radshape)
+ vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE
+
+ vza = viirs_data.sensor_zenith.values
+ dtr = np.pi/180
+ cosvza = np.cos(vza*dtr)
+
+ coeffs = utils.get_b1_thresholds()
+ coeffs[:, :3] = coeffs[:, :3] * threshold['b1_bias_adj']
+
+ # this quantity is the return of get_b1_thresholds() in the C code
+ # it's defined here to keep a consistent logic with the original source, for now
+ irtn = 0
+
+ if irtn != 0:
+ coeffs = thr
+
+ coeffs[:, :3] = coeffs[:, :3] * 1/np.power(cosvza, vzcpow)
+
+ confidence = conf.conf_test(rad, coeffs)
+
+ return confidence.reshape(radshape)
+
+
+
+class CloudMaskTests(object):
+
+ def __init__(self, scene, radiance, coefficients):
+ self.scene = scene
+ self.coefficients = coefficients
+
+ def select_coefficients(self):
+ pass
+
+ def test_G1(self):
+ pass
+
+ def test_G2(self):
+ pass
+
+ def test_G3(self):
+ pass
+
+ def test_G4(self):
+ pass
+
+ def overall_confidence(self):
+ pass
+
+
+def test():
+ rad = np.random.randint(50, size=[4, 8])
+ # coeffs = [5, 42, 20, 28, 15, 35, 1]
+ # coeffs = [20, 28, 5, 42, 15, 35, 1]
+ coeffs = [35, 15, 20, 1, 1]
+ # confidence = conf_test_dble(rad, coeffs)
+ confidence = test_11um(rad, coeffs)
+ print(rad)
+ print('\n')
+ print(confidence)
+
+
+if __name__ == "__main__":
+ test()
+
+
+
+
+
+
diff --git a/read_data.py b/read_data.py
index e0b7088..864f3d1 100644
--- a/read_data.py
+++ b/read_data.py
@@ -72,7 +72,7 @@ def correct_reflectances():
def read_ancillary_data(file_names, viirs_data):
# Ancillary files temporarily defined here. Eventually we will find a better way to pass these
- start_time = '2022-06-22 03:36:00.000'
+ start_time = '2022-06-22 14:54:00.000'
anc_dir = file_names['ANC_DIR']
sst_file = file_names['SST']
ndvi_file = file_names['NDVI']
diff --git a/tests.py b/tests.py
index 044b206..7a5edd7 100644
--- a/tests.py
+++ b/tests.py
@@ -6,6 +6,8 @@ import utils
import conf
+# ############## GROUP 1 TESTS ############## #
+
def test_11um(rad, threshold):
radshape = rad.shape
diff --git a/utils.py b/utils.py
index 7f90e49..6b45056 100644
--- a/utils.py
+++ b/utils.py
@@ -32,7 +32,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
# if refang <= thresholds['bounds'][1]:
if sunglint_flag == 3:
- sunglint_thr = np.full((4, len(refang)), thresholds[f'{band}_0deg'])
+ sunglint_thr = np.full((len(refang), 4), thresholds[f'{band}_0deg']).T
# sunglint_thr[:, :] = thresholds[f'{band}_0deg']
else:
--
GitLab