From 0b2781efcdcba876040d4859f78b8cbf7f3b0dd7 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Thu, 1 Dec 2022 22:31:48 +0000
Subject: [PATCH] first handful of tests validated. code cleanup still going
---
conf.py | 73 ++++++++++++++++++++++++++
preprocess_thresholds.py | 5 +-
read_data.py | 4 +-
tests.py | 88 +++++++++++++++++++++-----------
thresholds.mvcm.snpp.v0.0.1.yaml | 66 +++++++++++++++---------
5 files changed, 178 insertions(+), 58 deletions(-)
diff --git a/conf.py b/conf.py
index 5e4cf1a..161789b 100644
--- a/conf.py
+++ b/conf.py
@@ -21,6 +21,79 @@ def test_dble(flipped=False):
print(c)
+def conf_test_new(rad: np.ndarray,
+ thr: np.ndarray) -> np.ndarray:
+ '''
+ Assuming a linear function between min and max confidence level, the plot below shows
+ how the confidence (y axis) is computed as function of radiance (x axis).
+ This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
+ flipped.
+ gamma
+ c 1 ________
+ o | /
+ n | /
+ f | /
+ i | beta /
+ d 1/2 |....../
+ e | /
+ n | /
+ c | /
+ e 0________/
+ | alpha
+ --------- radiance ---------->
+ '''
+
+ radshape = rad.shape
+ rad = rad.ravel()
+ thr = np.array(thr)
+
+ if thr.ndim == 1:
+ thr = np.full((rad.shape[0], 4), thr[:4]).T
+
+ coeff = np.power(2, (thr[3] - 1))
+ hicut = thr[2, :]
+ beta = thr[1, :]
+ locut = thr[0, :]
+ power = thr[3, :]
+ confidence = np.zeros(rad.shape)
+
+ alpha, gamma = np.empty(rad.shape), np.empty(rad.shape)
+ flipped = np.zeros(rad.shape)
+
+ gamma[hicut > locut] = thr[2, hicut > locut]
+ alpha[hicut > locut] = thr[0, hicut > locut]
+ flipped[hicut > locut] = 0
+ gamma[hicut < locut] = thr[0, hicut < locut]
+ alpha[hicut < locut] = thr[2, hicut < locut]
+ flipped[hicut < locut] = 1
+
+ # Rad between alpha and beta
+ range_ = 2. * (beta - alpha)
+ s1 = (rad - alpha) / range_
+ idx = np.nonzero((rad <= beta) & (flipped == 0))
+ confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+ idx = np.nonzero((rad <= beta) & (flipped == 1))
+ confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])
+
+ # Rad between beta and gamma
+ range_ = 2. * (beta - gamma)
+ s1 = (rad - gamma) / range_
+ idx = np.nonzero((rad > beta) & (flipped == 0))
+ confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])
+ idx = np.nonzero((rad > beta) & (flipped == 1))
+ confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+
+ # Rad outside alpha-gamma interval
+ confidence[(rad > gamma) & (flipped is False)] = 1
+ confidence[(rad < alpha) & (flipped is False)] = 0
+ confidence[(rad > gamma) & (flipped is True)] = 0
+ confidence[(rad < alpha) & (flipped is True)] = 1
+
+ confidence = np.clip(confidence, 0, 1)
+
+ return confidence.reshape(radshape)
+
+
def conf_test(rad, thr):
'''
Assuming a linear function between min and max confidence level, the plot below shows
diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
index fd101d1..3f4ed92 100644
--- a/preprocess_thresholds.py
+++ b/preprocess_thresholds.py
@@ -7,6 +7,7 @@ import utils
from numpy.lib.stride_tricks import sliding_window_view
_dtr = np.pi/180
+_DTR = np.pi/180
# this case is written for the 11-12um Cirrus Test for scenes that follow pattern 1 (see note below)
@@ -42,7 +43,7 @@ def prepare_thresholds(data, thresholds):
def preproc(data, thresholds, scene):
- cosvza = np.cos(data.sensor_zenith * _dtr)
+ cosvza = np.cos(data.sensor_zenith * _DTR)
schi = (1/cosvza).where(cosvza > 0, 99.0)
schi = schi.values.reshape(np.prod(schi.shape))
@@ -90,7 +91,7 @@ def preproc_sst(data, thresholds):
m32c = data.M16 - 273.16
m31c_m32c = m31c - m32c
sstc = data.geos_sfct - 273.16
- cosvza = np.cos(data.sensor_zenith*_dtr)
+ cosvza = np.cos(data.sensor_zenith*_DTR)
a = thresholds['coeffs']
diff --git a/read_data.py b/read_data.py
index 166f1bd..3507048 100644
--- a/read_data.py
+++ b/read_data.py
@@ -149,8 +149,8 @@ def get_data(file_names: Dict[str, str],
viirs_data = read_ancillary_data(file_names, viirs_data)
# Compute channel differences and ratios that are used in the tests
- viirs_data['M15-M14'] = (('number_of_lines', 'number_of_pixels'),
- viirs_data.M15.values - viirs_data.M14.values)
+ 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'),
diff --git a/tests.py b/tests.py
index 82df1fd..2adb710 100644
--- a/tests.py
+++ b/tests.py
@@ -2,6 +2,7 @@ import numpy as np
import xarray as xr
from numpy.lib.stride_tricks import sliding_window_view
+from typing import Dict
import utils
import conf
@@ -10,6 +11,8 @@ import preprocess_thresholds as pt
import importlib
+_DTR = np.pi/180
+
# this is used for testing, eventually we want to remove it
importlib.reload(pt)
@@ -193,31 +196,6 @@ def vis_nir_ratio_test(rad1, rad2, threshold, sg_threshold):
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
-
-
# old class, doesn't use xarray much
class CloudTests_old:
@@ -247,18 +225,68 @@ class CloudTests_old:
# new class to try to use xarray more extensively
-class CloudTests:
+class CloudTests(object):
- def __init__(self, data, scene_name, thresholds):
+ def __init__(self,
+ data: xr.Dataset,
+ scene_name: str,
+ thresholds: Dict) -> None:
self.data = data
self.scene_name = scene_name
self.thresholds = thresholds
+ self.scene_idx = tuple(np.nonzero(data[scene_name] == 1))
+
+ def test_11um(self,
+ band: str,
+ cmin: np.ndarray) -> np.ndarray:
+ confidence = np.ones(self.data.M01.shape)
+ threshold = self.thresholds[self.scene_name]['11um_Test']
+
+ if threshold['perform'] is True:
+ print(f'Testing "{self.scene_name}"\n')
+ rad = self.data[band].values[self.scene_idx]
+ confidence[self.scene_idx] = conf.conf_test_new(rad, threshold['thr'])
+
+ cmin = np.fmin(cmin, confidence)
+
+ return cmin
+
+ def sst_test(self,
+ band: str,
+ cmin: np.ndarray) -> np.ndarray:
+
+ confidence = np.ones(self.data.M01.shape)
+ threshold = self.thresholds[self.scene_name]['SST_Test']
+
+ if threshold['perform'] is True:
+ m31 = self.data.M15.values - 273.16
+ bt_diff = self.data['M15-M16'].values
+ sst = self.data.geos_sfct.values - 273.16
+ cosvza = np.cos(self.data.sensor_zenith.values*_DTR)
+
+ c = threshold['coeffs']
+
+ modsst = 273.16 + c[0] + c[1]*m31 + c[2]*bt_diff*sst + c[3]*bt_diff*((1/cosvza) - 1)
+ sfcdif = self.data.geos_sfct.values - modsst
+
+ 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
+
+ def bt_diff_86_11um(self,
+ band: str,
+ cmin: np.ndarray) -> np.ndarray:
- def test_11um(self, band):
confidence = np.ones(self.data.M01.shape)
- if self.thresholds[self.scene_name]['11um_Test']['perform'] is True:
+ threshold = self.thresholds[self.scene_name]['8.6-11um_Test']
- confidence = conf.conf_test(data, band)
+ if threshold['perform'] is True:
+ print(f'Testing "{self.scene_name}"\n')
+ rad = self.data[band].values[self.scene_idx]
+ confidence[self.scene_idx] = conf.conf_test_new(rad, threshold['thr'])
cmin = np.fmin(cmin, confidence)
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 4ec8be0..0cffd7a 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -120,14 +120,21 @@ Land_Day_Desert_Coast:
# lds11_12hcadj_c : 1.25
Ocean_Day:
+ 11um_Test:
+ thr: [267.0, 270.0, 273.0, 1.0, 1.0]
+ perform: True
+ SST_Test:
+ thr: [3.000, 2.500, 2.000, 1.0, 1.0]
+ coeffs: [1.886, 0.938, 0.128, 1.094]
+ perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
- 11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.0, 1.0, 1.0]
- 11um_Test:
- thr: [267.0, 270.0, 273.0, 1.0, 1.0]
+ 8.6-11um_Test:
+ thr: [-0.50, -1.00, -1.50, 1.0, 1.0]
perform: True
+ 11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.0, 1.0, 1.0]
CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.6_2.1um_NIR_Reflectance_Test:
@@ -152,10 +159,6 @@ Ocean_Day:
# Hi-Mid-Lo--Lo-Mid-Hi
# Last two values are power and switch
# to turn the test on/off
- 8.6-11um_Test: [-0.50, -1.00, -1.50, 1.0, 1.0]
- SST_Test:
- thr: [3.000, 2.500, 2.000, 1.0, 1.0]
- coeffs: [1.886, 0.938, 0.128, 1.094]
1.38um_High_Cloud_Test:
coeffs: [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204]
adj: 0.0010
@@ -165,20 +168,25 @@ Ocean_Day:
do_b7bias_adj : 0.97
Ocean_Night:
+ 11um_Test:
+ thr: [267.0, 270.0, 273.0, 2.0, 1.0]
+ perform: True
+ SST_Test:
+ thr: [3.000, 2.500, 2.000, 1.0, 1.0]
+ coeffs: [1.8860, 0.9380, 0.1280, 1.094]
+ perform: True
+ 8.6-11um_Test:
+ thr: [-0.5, -1.0, -1.5, 1.0, 1.0]
+ perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.0
11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 1.0, 1.0]
- 11um_Test: [267.0, 270.0, 273.0, 2.0, 1.0]
CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
no86_73 : [14.0, 15.0, 16.0, 1.0, 1.0]
11um_Variability_Test: [3.0, 6.0, 7.0, 1.0, 1.0]
- 8.6-11um_Test: [-0.5, -1.0, -1.5, 1.0, 1.0]
- SST_Test:
- thr: [3.000, 2.500, 2.000, 1.0, 1.0]
- coeffs: [1.8860, 0.9380, 0.1280, 1.094]
11-4um_BT_Difference_Test:
coeffs: [-0.7093, 0.1128, -0.1567]
corr: 0.0 # no_intadj
@@ -339,7 +347,17 @@ Polar_Night_Snow:
pnbt3 : 230.0
pntpw : 0.2
-Polar_Ocean_Day:
+Polar_Day_Ocean:
+ 11um_Test:
+ thr: [267.0, 270.0, 273.0, 1.0, 1.0]
+ perform: True
+ SST_Test:
+ thr: [3.000, 2.500, 2.000, 1.0, 1.0]
+ coeffs: [1.886, 0.938, 0.128, 1.094]
+ perform: True
+ 8.6-11um_Test:
+ thr: [-0.50, -1.00, -1.50, 1.0, 1.0]
+ perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
@@ -359,13 +377,8 @@ Polar_Ocean_Day:
bias: 0.98
perform_test: True
11-4um_Oceanic_Stratus_Test: [-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]
Vis/NIR_Ratio_Test: [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0]
- 8.6-11um_Test: [-0.50, -1.00, -1.50, 1.0, 1.0]
- SST_Test:
- thr: [3.000, 2.500, 2.000, 1.0, 1.0]
- coeffs: [1.886, 0.938, 0.128, 1.094]
pdo_b26pfm : 1.0
pdo_b2bias_adj : 0.96
pdo_b7bias_adj : 0.97
@@ -375,21 +388,26 @@ Polar_Ocean_Day:
szafac: 4.0
perform: True
-Polar_Ocean_Night:
+Polar_Night_Ocean:
+ 11um_Test:
+ thr: [267.0, 270.0, 273.0, 1.0, 1.0]
+ perform: True
+ SST_Test:
+ thr: [3.000, 2.500, 2.000, 1.0, 1.0]
+ coeffs: [1.8860, 0.9380, 0.1280, 1.094]
+ perform: True
+ 8.6-11um_Test:
+ thr: [-0.5, -1.0, -1.5, 1.0, 1.0]
+ perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.0
pnobt1 : 280.0
11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 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]
pno86_73 : [14.0, 15.0, 16.0, 1.0, 1.0]
11um_Variability_Test: [3.0, 6.0, 7.0, 1.0, 1.0]
- 8.6-11um_Test: [-0.5, -1.0, -1.5, 1.0, 1.0]
- SST_Test:
- thr: [3.000, 2.500, 2.000, 1.0, 1.0]
- coeffs: [1.8860, 0.9380, 0.1280, 1.094]
11-4um_BT_Difference_Test:
coeffs: [-0.7093, 0.1128, -0.1567]
corr: 0.0 # pno_intadj
--
GitLab