From f78bdf68558d4e32eb6f4c9791008d652832f2f9 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Mon, 7 Nov 2022 21:36:11 +0000
Subject: [PATCH] more tests implemented

---
 main.py                          |  7 +++---
 preprocess_thresholds.py         | 41 ++++++++++++++++++++++++++++++++
 tests.py                         | 12 +++++++++-
 thresholds.mvcm.snpp.v0.0.1.yaml | 20 +++++++++++-----
 4 files changed, 70 insertions(+), 10 deletions(-)

diff --git a/main.py b/main.py
index 1133a78..3c55283 100644
--- a/main.py
+++ b/main.py
@@ -99,7 +99,7 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
                'CO2 High Clouds Test': False,
                'Water Vapor High Clouds Test': False,
                'Surface Temperature Test': False,
-               'SST Test': True,
+               'SST Test': False,
                '8.6-11um BT Difference Test': False,
                '11-12um BTD Thin Cirrus Test': False,
                '11-4um BT Difference Test': False,
@@ -107,7 +107,7 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
                'Water Vapor Cloud Test': False,
                '11um BT Variability Test': False,
                '11-4um BTD Oceanic Stratus': False,
-               'NIR Reflectance Test': False,
+               'NIR Reflectance Test': True,
                'Vis/NIR Ratio Test': False,
                '1.6um or 2.1um NIR Reflectance Test': False,
                'Visible Reflectance Test': False,
@@ -213,7 +213,8 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
 
     if perform['NIR Reflectance Test'] is True:
         for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
-            pass
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G3 = SceneType.single_thredhold_test('NIR_Reflectance_Test', 'M07', cmin_G3)
 
     if perform['Vis/NIR Ratio Test'] is True:
         for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
index 9f7bc6a..6c95d1d 100644
--- a/preprocess_thresholds.py
+++ b/preprocess_thresholds.py
@@ -2,6 +2,7 @@ import numpy as np
 import xarray as xr
 
 import ancillary_data as anc
+import utils
 
 _dtr = np.pi/180
 
@@ -64,6 +65,46 @@ def preproc_sst(data, thresholds):
     return sfcdif
 
 
+def preproc_nir(data, thresholds, scene):
+
+    sza = data.solar_zenith
+    band_n = 2
+    # NOTE: the visud condition in the C code is equivalent to having sza <= 85
+    # For the time being the visud filtering is not implemented
+
+    a = np.array(thresholds[scene]['NIR_Reflectance_Test']['coeffs'])
+    vzcpow = thresholds['VZA_correction']['vzcpow'][0]
+    refang = data.sunglint_angle.values.reshape(np.prod(data.sunglint_angle.shape))
+    sunglint_thresholds = thresholds['Sun_Glint']
+    sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(refang.shape)
+    nir_thresh = thresholds[scene]['NIR_Reflectance_Test']
+
+    hicut0 = a[0] + a[1]*sza + a[2]*np.power(sza, 2) + a[3]*np.power(sza, 3)
+    hicut0 = (hicut0 * 0.01) + nir_thresh['adj']
+    hicut0 = (hicut0 + nir_thresh['bias']).values.reshape(refang.shape)
+    midpt0 = hicut0 + (nir_thresh['midpt_coeff'] * nir_thresh['bias'])
+    locut0 = midpt0 + (nir_thresh['locut_coeff'] * nir_thresh['bias'])
+    thr = np.array([locut0, midpt0, hicut0, nir_thresh['thr'][3]*np.ones(refang.shape)])
+
+    cosvza = np.cos(data.sensor_zenith*_dtr).values.reshape(refang.shape)
+    corr_thr = np.zeros((4, refang.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:
+            dosgref = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
+            corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * \
+                (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
+            corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
+
+    corr_thr = np.transpose(corr_thr.reshape((4, sza.shape[0], sza.shape[1])), (1, 2, 0))
+
+    return corr_thr
+
+
 # NOTE: 11-12um Cirrus Test
 # hicut is computed in different ways depending on the scene
 # 1. midpt - adj
diff --git a/tests.py b/tests.py
index d2c1f4e..0508626 100644
--- a/tests.py
+++ b/tests.py
@@ -254,9 +254,14 @@ class CloudTests:
         if band == 'bad_data':
             return cmin
         print(f'Running test "{test_name}" for "{self.scene_name}"')
+
         # preproc_thresholds()
-        thr = np.array(self.thresholds[self.scene_name][test_name]['thr'])
+        if 'thr' in self.thresholds[self.scene_name][test_name]:
+            thr = np.array(self.thresholds[self.scene_name][test_name]['thr'])
+        else:
+            thr = np.array(self.thresholds[self.scene_name][test_name])
         thr_xr = xr.Dataset()
+
         if test_name == '11-12um_Cirrus_Test':
             thr_xr['threshold'] = pt.preproc(self.data, self.thresholds[self.scene_name])
             thr = np.ones((5,))  # This is only temporary to force the logic of the code
@@ -265,6 +270,9 @@ class CloudTests:
         elif test_name == 'SST_Test':
             thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'),
                                    np.ones((self.data[band].shape[0], self.data[band].shape[1], 5))*thr)
+        elif test_name == 'NIR_Reflectance_Test':
+            corr_thr = pt.preproc_nir(self.data, self.thresholds, self.scene_name)
+            thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'), corr_thr)
         else:
             thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'),
                                    np.ones((self.data[band].shape[0], self.data[band].shape[1], 5))*thr)
@@ -274,6 +282,8 @@ class CloudTests:
             data['sfcdif'] = (('number_of_lines', 'number_of_pixels'),
                               pt.preproc_sst(data, self.thresholds[self.scene_name][test_name]).values)
             band = 'sfcdif'
+        if test_name == 'NIR_Reflectance_Test':
+            pass
 
         if thr[4] == 1:
             print('test running...')
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 9fd834a..ee93cca 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -112,7 +112,13 @@ Ocean_Day:
   do_b7mid       : 0.0070
   do_b7lo        : 0.0100
   do_b7pfm       : 1.0
-  ref2           : [0.062, 0.043, 0.029, 1.0, 1.0]
+  NIR_Reflectance_Test:
+    thr: [0.062, 0.043, 0.029, 1.0, 1.0]
+    coeffs: [1.7291, 0.0715, -0.0026, 0.000025889]
+    adj: 0.0050
+    bias: 0.96
+    midpt_coeff: 0.0200
+    locut_coeff: 0.0100
   #  vnir_ratio_hi  : [1.001, 1.154, 1.205, 1.0, 1.0]
   #  vnir_ratio_lo  : [0.934, 0.889, 0.837, 1.0]
   vis_nir_ratio : [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0] # This replace dovrathi and dovratlo
@@ -125,15 +131,10 @@ Ocean_Day:
   SST_Test:
     thr: [3.000, 2.500, 2.000, 1.0, 1.0]
     coeffs: [1.886, 0.938, 0.128, 1.094]
-  b2coeffs       : [1.7291, 0.0715, -0.0026, 0.000025889]
-  b2adj          : 0.0050
-  b2mid          : 0.0200
-  b2lo           : 0.0100
   do_b26coeffs   : [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204]
   do_b26adj      : 0.0010
   do_b26_szafac  : 4.0
   do_b26pfm      : 1.0
-  b2bias_adj     : 0.96
   do_b6bias_adj  : 0.98
   do_b7bias_adj  : 0.97
 
@@ -295,6 +296,13 @@ Polar_Ocean_Day:
     coeffs: [3.0, 1.0]
     cmult: 0.3
     adj: 1.25
+  NIR_Reflectance_Test:
+    thr: [0.062, 0.043, 0.029, 1.0, 1.0]
+    coeffs: [1.7291, 0.0715, -0.0026, 0.000025889]
+    adj: 0.0050
+    bias: 0.96
+    midpt_coeff: 0.0200
+    locut_coeff: 0.0100
   pdo11_4lo      : [-11.0, -9.0,  -7.0,  1.0, 1.0]
   11um_Test: [267.0, 270.0, 273.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-- 
GitLab