From f9e22fc39e2233eddb1b4c1ef32ab6c55a69ef43 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Mon, 7 Jul 2025 19:45:39 +0000
Subject: [PATCH] reduced array size of thresholds by dropping power and binary
flag
---
mvcm/ancillary.pyx | 4 +-
mvcm/conf.py | 26 ++++----
mvcm/preprocess_thresholds.py | 62 +++++++++---------
mvcm/spectral_tests.py | 70 +++++++++++++--------
mvcm/utils.py | 8 +--
thresholds/thresholds.mvcm.snpp.v0.0.1.yaml | 31 ++++++++-
6 files changed, 122 insertions(+), 79 deletions(-)
diff --git a/mvcm/ancillary.pyx b/mvcm/ancillary.pyx
index 0062038..075a861 100644
--- a/mvcm/ancillary.pyx
+++ b/mvcm/ancillary.pyx
@@ -139,7 +139,7 @@ def py_snow_mask(cnp.ndarray[float, ndim=1] latitude, cnp.ndarray[unsigned char,
def py_conf_test(cnp.ndarray[float, ndim=1] rad,
cnp.ndarray[float, ndim=1] locut,
cnp.ndarray[float, ndim=1] hicut,
- cnp.ndarray[float, ndim=1] power,
+ power,
cnp.ndarray[float, ndim=1] midpt
):
@@ -148,7 +148,7 @@ def py_conf_test(cnp.ndarray[float, ndim=1] rad,
# cdef float[::1] c
for i in range(len(rad)):
- conf[i] = conf_test(rad[i], locut[i], hicut[i], power[i], midpt[i])
+ conf[i] = conf_test(rad[i], locut[i], hicut[i], power, midpt[i])
return conf
diff --git a/mvcm/conf.py b/mvcm/conf.py
index 0de0624..858a179 100644
--- a/mvcm/conf.py
+++ b/mvcm/conf.py
@@ -3,11 +3,13 @@
import logging
try:
- from memory_profiler import profile # real decorator
-except ImportError: # no-op fallback
- def profile(func): # type: ignore[misc]
+ from memory_profiler import profile # real decorator
+except ImportError: # no-op fallback
+
+ def profile(func): # type: ignore[misc]
return func
+
import numpy as np
import xarray as xr
@@ -120,7 +122,7 @@ def conf_test_xr(rad: xr.DataArray, thr: np.ndarray) -> xr.DataArray:
# @profile
-def conf_test_new(rad: np.ndarray, thr: np.ndarray) -> np.ndarray:
+def conf_test_new(rad: np.ndarray, thr: np.ndarray, power: float) -> np.ndarray:
"""Compute confidence based on thresholds."""
"""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).
@@ -147,11 +149,10 @@ def conf_test_new(rad: np.ndarray, thr: np.ndarray) -> np.ndarray:
if thr.ndim == 1:
thr = np.full((rad.shape[0], 4), thr[:4]).T
- coeff = np.power(2, (thr[3, :] - 1))
+ coeff = np.power(2, (power - 1))
locut = thr[0, :]
beta = thr[1, :]
hicut = thr[2, :]
- power = thr[3, :]
confidence = np.zeros(rad.shape)
alpha, gamma = np.empty(rad.shape), np.empty(rad.shape)
@@ -168,17 +169,17 @@ def conf_test_new(rad: np.ndarray, thr: np.ndarray) -> np.ndarray:
range_ = 2.0 * (beta - alpha)
s1 = (rad - alpha) / range_
idx = np.nonzero((rad <= beta) & (flipped == 0))
- confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+ confidence[idx] = coeff * np.power(s1[idx], power)
idx = np.nonzero((rad <= beta) & (flipped == 1))
- confidence[idx] = 1.0 - (coeff[idx] * np.power(s1[idx], power[idx]))
+ confidence[idx] = 1.0 - (coeff * np.power(s1[idx], power))
# Rad between beta and gamma
range_ = 2.0 * (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]))
+ confidence[idx] = 1.0 - (coeff * np.power(s1[idx], power))
idx = np.nonzero((rad > beta) & (flipped == 1))
- confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+ confidence[idx] = coeff * np.power(s1[idx], power)
# Rad outside alpha-gamma interval
confidence[(rad > gamma) & (flipped is False)] = 1
@@ -385,7 +386,7 @@ class Confidence:
if single_thresholds is not None:
if single_thresholds.ndim == 1:
- thresholds = np.full((self.rad.shape[0], 4), single_thresholds).T
+ thresholds = np.full((self.rad.shape[0], 3), single_thresholds).T
else:
thresholds = single_thresholds
self.alpha = np.where(
@@ -395,7 +396,8 @@ class Confidence:
self.gamma = np.where(
thresholds[2, :] > thresholds[0, :], thresholds[2, :], thresholds[0, :]
)
- self.power = thresholds[3, :]
+ # self.power = thresholds[3, :]
+ self.power = single_thresholds[3]
self.flip = np.where(thresholds[2, :] > thresholds[0, :], 0, 1)
self.single_conf, self.double_conf = True, False
if double_thresholds is not None:
diff --git a/mvcm/preprocess_thresholds.py b/mvcm/preprocess_thresholds.py
index 3c1289e..10fd380 100644
--- a/mvcm/preprocess_thresholds.py
+++ b/mvcm/preprocess_thresholds.py
@@ -3,11 +3,13 @@
import logging # noqa
try:
- from memory_profiler import profile # real decorator
-except ImportError: # no-op fallback
- def profile(func): # type: ignore[misc]
+ from memory_profiler import profile # real decorator
+except ImportError: # no-op fallback
+
+ def profile(func): # type: ignore[misc]
return func
+
import ancillary_data as anc
import numpy as np
import xarray as xr
@@ -169,8 +171,6 @@ def xr_thresholds_11_12um(data: xr.Dataset, thresholds: dict, scene: str) -> np.
thr_out[..., 0] = locut.reshape(shape)
thr_out[..., 1] = midpt.reshape(shape)
thr_out[..., 2] = hicut.reshape(shape)
- thr_out[..., 3] = 1.0
- thr_out[..., 4] = 1.0
return thr_out
@@ -265,9 +265,7 @@ def thresholds_11_12um(
thr_out = np.empty((locut.shape[0], 5))
thr_out[:, 0] = locut
thr_out[:, 1] = midpt
- thr_out[:, 2] = hicut
- thr_out[:, 3] = 1.0
- thr_out[:, 4] = 1.0
+ thr_out[:, 2]
if thr_out.ndim == 1:
thr_out = np.reshape(thr_out, (1, thr_out.shape[0]))
@@ -312,22 +310,22 @@ def thresholds_NIR(data, thresholds, scene, test_name, scene_idx):
hicut = hicut * nir_thresh["bias"]
midpt = hicut + (nir_thresh["midpt_coeff"] * nir_thresh["bias"])
locut = midpt + (nir_thresh["locut_coeff"] * nir_thresh["bias"])
- thr = np.array([locut, midpt, hicut, nir_thresh["thr"][3] * np.ones(refang.shape)])
+ thr = np.array([locut, midpt, hicut])
cosvza = np.cos(data.sensor_zenith.values[scene_idx] * _DTR).ravel()
- corr_thr = np.zeros((4, refang.shape[0]))
+ corr_thr = np.zeros((3, refang.shape[0]))
corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (
1.0 / np.power(cosvza[sunglint_flag == 0], vzcpow)
)
- corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
- for flag in range(1, 4):
+ # corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
+ for flag in range(1, 3):
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] * (
+ corr_thr[:, sunglint_flag == flag] = dosgref[:, sunglint_flag == flag] * (
1.0 / np.power(cosvza[sunglint_flag == flag], vzcpow)
)
- corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
+ # corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
return corr_thr
@@ -384,8 +382,6 @@ def thresholds_surface_temperature(data, thresholds, scene_idx):
thr_out[:, 0] = locut
thr_out[:, 1] = midpt
thr_out[:, 2] = hicut
- thr_out[:, 3] = 1.0
- thr_out[:, 4] = 1.0
if thr_out.ndim == 1:
thr_out = np.reshape(thr_out, (1, thr_out.shape[0]))
@@ -567,11 +563,11 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
locut = thresholds["thr"][0] * np.ones((len(scene_idx[0]),))
midpt = thresholds["thr"][1] * np.ones((len(scene_idx[0]),))
hicut = thresholds["thr"][2] * np.ones((len(scene_idx[0]),))
- power = thresholds["thr"][3] * np.ones((len(scene_idx[0]),))
+ # power = thresholds["thr"][3] * np.ones((len(scene_idx[0]),))
# out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut,
# np.ones(data.ndvi.shape), power)),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
- out_thr = np.dstack((locut, midpt, hicut, np.ones(locut.shape), power))
+ out_thr = np.dstack((locut, midpt, hicut))
return np.squeeze(out_thr.T)
@@ -579,7 +575,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
bt_bounds = thresholds["bt11_bounds"]
locut, midpt = np.empty(rad.shape), np.empty(rad.shape)
- hicut, power = np.empty(rad.shape), np.empty(rad.shape)
+ hicut = np.empty(rad.shape)
lo, hi = np.empty(rad.shape), np.empty(rad.shape)
lo_thr, hi_thr = np.empty(rad.shape), np.empty(rad.shape)
conf_range = np.empty(rad.shape)
@@ -588,13 +584,13 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
locut[idx] = thresholds["low"][0]
midpt[idx] = thresholds["low"][1]
hicut[idx] = thresholds["low"][2]
- power[idx] = thresholds["low"][3]
+ # power[idx] = thresholds["low"][3]
idx = np.nonzero(rad > bt_bounds[3])
locut[idx] = thresholds["high"][0]
midpt[idx] = thresholds["high"][1]
hicut[idx] = thresholds["high"][2]
- power[idx] = thresholds["high"][3]
+ # power[idx] = thresholds["high"][3]
# # # # #
idx = np.nonzero(
@@ -604,7 +600,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
hi[idx] = thresholds["bt11_bounds"][3]
lo_thr[idx] = thresholds["mid1"][0]
hi_thr[idx] = thresholds["mid1"][1]
- power[idx] = thresholds["mid1"][3]
+ # power[idx] = thresholds["mid1"][3]
conf_range[idx] = thresholds["mid1"][2]
idx = np.nonzero((rad >= bt_bounds[0]) & (rad < bt_bounds[1]))
@@ -612,7 +608,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
hi[idx] = thresholds["bt11_bounds"][1]
lo_thr[idx] = thresholds["mid1"][0]
hi_thr[idx] = thresholds["mid1"][1]
- power[idx] = thresholds["mid1"][3]
+ # power[idx] = thresholds["mid1"][3]
conf_range[idx] = thresholds["mid1"][2]
idx = np.nonzero((rad >= bt_bounds[1]) & (rad < bt_bounds[2]))
@@ -620,7 +616,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
hi[idx] = thresholds["bt11_bounds"][2]
lo_thr[idx] = thresholds["mid2"][0]
hi_thr[idx] = thresholds["mid2"][1]
- power[idx] = thresholds["mid2"][3]
+ # power[idx] = thresholds["mid2"][3]
conf_range[idx] = thresholds["mid2"][2]
idx = np.nonzero((rad >= bt_bounds[2]) & (rad < bt_bounds[3]))
@@ -628,7 +624,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
hi[idx] = thresholds["bt11_bounds"][3]
lo_thr[idx] = thresholds["mid3"][0]
hi_thr[idx] = thresholds["mid3"][1]
- power[idx] = thresholds["mid3"][3]
+ # power[idx] = thresholds["mid3"][3]
conf_range[idx] = thresholds["mid3"][2]
idx = np.nonzero(
@@ -648,7 +644,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
# out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut,
# np.ones(data.ndvi.shape), power)),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
- out_thr = np.squeeze(np.dstack((locut, midpt, hicut, np.ones(locut.shape), power)))
+ out_thr = np.squeeze(np.dstack((locut, midpt, hicut)))
if out_thr.ndim == 1:
out_thr = np.reshape(out_thr, (1, out_thr.shape[0]))
@@ -802,7 +798,7 @@ def vis_refl_thresholds(data, thresholds, scene, scene_idx):
# dims=('number_of_lines', 'number_of_pixels', 'z'))
# out_rad = xr.DataArray(data=m128.reshape(data.M01.shape),
# dims=('number_of_lines', 'number_of_pixels'))
- out_thr = np.squeeze(np.dstack((b1_locut, b1_midpt, b1_hicut, b1_power)))
+ out_thr = np.squeeze(np.dstack((b1_locut, b1_midpt, b1_hicut)))
out_rad = m128
@@ -819,14 +815,14 @@ def vis_refl_thresholds(data, thresholds, scene, scene_idx):
def gemi_thresholds(data, thresholds, scene_name, scene_idx):
"""Compute thresholds for GEMI test."""
ndvi = data.ndvi.values[scene_idx].ravel()
- gemi_thr = np.ones((ndvi.shape[0], 5))
+ gemi_thr = np.ones((ndvi.shape[0], 3))
idx = np.nonzero(ndvi < 0.1)
- gemi_thr[idx, :3] = thresholds["gemi0"][:3]
+ gemi_thr[idx, :] = thresholds["gemi0"][:3]
idx = np.nonzero((ndvi >= 0.1) & (ndvi < 0.2))
- gemi_thr[idx, :3] = thresholds["gemi1"][:3]
+ gemi_thr[idx, :] = thresholds["gemi1"][:3]
idx = np.nonzero((ndvi >= 0.2) & (ndvi < 0.3))
- gemi_thr[idx, :3] = thresholds["gemi2"][:3]
+ gemi_thr[idx, :] = thresholds["gemi2"][:3]
# thr_out = xr.DataArray(data=np.dstack((gemi_thr[:, :, 0],
# gemi_thr[:, :, 1], gemi_thr[:, :, 2],
@@ -903,9 +899,7 @@ def thresholds_1_38um_test(data, thresholds, scene_name, scene_idx):
# np.ones(data.ndvi.shape),
# np.ones(data.ndvi.shape))),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
- out_thr = np.squeeze(
- np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)))
- )
+ out_thr = np.squeeze(np.dstack((locut, midpt, hicut)))
if out_thr.ndim == 1:
out_thr = np.reshape(out_thr, (1, out_thr.shape[0]))
diff --git a/mvcm/spectral_tests.py b/mvcm/spectral_tests.py
index 946d8b9..c77b696 100644
--- a/mvcm/spectral_tests.py
+++ b/mvcm/spectral_tests.py
@@ -6,11 +6,13 @@ import logging
import os
try:
- from memory_profiler import profile # real decorator
-except ImportError: # no-op fallback
- def profile(func): # type: ignore[misc]
+ from memory_profiler import profile # real decorator
+except ImportError: # no-op fallback
+
+ def profile(func): # type: ignore[misc]
return func
+
import ancillary_data as anc
import numpy as np
import psutil
@@ -237,12 +239,15 @@ class CloudTests:
locut = np.full(rad.shape, threshold["thr"][0], dtype=np.float32)
midpt = np.full(rad.shape, threshold["thr"][1], dtype=np.float32)
hicut = np.full(rad.shape, threshold["thr"][2], dtype=np.float32)
- power = np.full(rad.shape, threshold["thr"][3], dtype=np.float32)
+ power = threshold["thr"][3]
kwargs["confidence"][self.scene_idx] = anc.py_conf_test(
f_rad, locut, hicut, power, midpt
)
else:
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, threshold["thr"])
+ power = threshold["thr"][3]
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(
+ rad, threshold["thr"], power
+ )
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -282,7 +287,10 @@ class CloudTests:
bits["test"][scene_idx] = self.compute_test_bits(
surface_diff_temp, sfc_thresholds[1, :], "<", scene_idx=scene_idx
)
- kwargs["confidence"][scene_idx] = conf.conf_test_new(surface_diff_temp, sfc_thresholds)
+ power = 1.0
+ kwargs["confidence"][scene_idx] = conf.conf_test_new(
+ surface_diff_temp, sfc_thresholds, power
+ )
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -303,6 +311,7 @@ class CloudTests:
bt_diff = m31 - m32
sst = self.data.geos_sst.values - 273.16
cosvza = np.cos(self.data.sensor_zenith.values * _DTR)
+ power = threshold["thr"][3]
c = threshold["coeffs"]
@@ -320,7 +329,7 @@ class CloudTests:
sfcdif[self.scene_idx], threshold["thr"][1], "<"
)
kwargs["confidence"][self.scene_idx] = conf.conf_test_new(
- sfcdif[self.scene_idx], threshold["thr"]
+ sfcdif[self.scene_idx], threshold["thr"], power
)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -342,8 +351,9 @@ class CloudTests:
# rad = self.data[band].values[self.scene_idx]
rad = self.data.M14.values[self.scene_idx] - self.data.M15.values[self.scene_idx]
+ power = threshold["thr"][3]
bits["test"][self.scene_idx] = self.compute_test_bits(rad, threshold["thr"][1], "<")
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, threshold["thr"])
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, threshold["thr"], power)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -389,7 +399,7 @@ class CloudTests:
locut = np.array(thr[0, :], dtype=np.float32)
midpt = np.array(thr[1, :], dtype=np.float32)
hicut = np.array(thr[2, :], dtype=np.float32)
- power = np.array(thr[3, :], dtype=np.float32)
+ power = threshold["power"]
logger.info("Step 4")
logger.info(f"spectral_tests: Memory usage #7: {proc.memory_info().rss / 1e6} MB")
# temp_confidence = anc.py_conf_test( f_rad, locut, hicut, power, midpt)
@@ -714,7 +724,7 @@ class CloudTests:
thr = preproc.polar_night_thresholds(
self.data, threshold, self.scene_name, test_name, self.scene_idx
)
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad11_4, thr)
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad11_4, thr, 1.0)
if self.scene_name in ["Polar_Night_Snow"]:
# rad11_12 = self.data[band1].values
@@ -731,9 +741,9 @@ class CloudTests:
self.data, threshold["bottom"], self.scene_name, test_name, idx_top
)
- kwargs["confidence"][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top)
+ kwargs["confidence"][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top, 1.0)
kwargs["confidence"][idx_bottom] = conf.conf_test_new(
- rad11_4[idx_bottom], thr_bottom
+ rad11_4[idx_bottom], thr_bottom, 1.0
)
if self.scene_name in ["Night_Snow"]:
@@ -749,9 +759,9 @@ class CloudTests:
self.data, threshold["bottom"], self.scene_name, test_name, idx_top
)
- kwargs["confidence"][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top)
+ kwargs["confidence"][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top, 1.0)
kwargs["confidence"][idx_bottom] = conf.conf_test_new(
- rad11_4[idx_bottom], thr_bottom
+ rad11_4[idx_bottom], thr_bottom, 1.0
)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -794,8 +804,11 @@ class CloudTests:
) # subtracting 1 because I don't want to count the center pixel
thr = np.array(threshold["thr"])
+ power = threshold["thr"][3]
bits["test"][self.scene_idx] = self.compute_test_bits(var[self.scene_idx], thr[1], ">")
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(var[self.scene_idx], thr)
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(
+ var[self.scene_idx], thr, power
+ )
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -844,6 +857,7 @@ class CloudTests:
)
else:
thr = threshold["thr"]
+ power = threshold["thr"][3]
logger.info(f'Running test 11-4um_Oceanic_Stratus_Test on "{self.scene_name}"\n')
# rad = self.data[band].values
@@ -865,7 +879,7 @@ class CloudTests:
idx = tuple(
np.nonzero((self.data[self.scene_name] == 1) & (sunglint_flag.values == 0))
)
- kwargs["confidence"][idx] = conf.conf_test_new(rad.values[idx], thr)
+ kwargs["confidence"][idx] = conf.conf_test_new(rad.values[idx], thr, power)
elif self.scene_name in ["Land_Day", "Land_Day_Coast"]:
idx = tuple(
@@ -873,7 +887,7 @@ class CloudTests:
)
bits["test"][idx] = 1
idx = tuple(np.nonzero(self.data[self.scene_name] == 1))
- kwargs["confidence"][idx] = conf.conf_test_new(rad.values[idx], thr)
+ kwargs["confidence"][idx] = conf.conf_test_new(rad.values[idx], thr, power)
elif self.scene_name in [
"Land_Day_Desert",
@@ -904,12 +918,11 @@ class CloudTests:
locut = np.array(thr[0, :], dtype=np.float32)
midpt = np.array(thr[1, :], dtype=np.float32)
hicut = np.array(thr[2, :], dtype=np.float32)
- power = np.array(thr[3, :], dtype=np.float32)
bits["test"][idx] = 1
kwargs["confidence"][idx] = anc.py_conf_test(f_rad, locut, hicut, power, midpt)
else:
kwargs["confidence"][self.scene_idx] = conf.conf_test_new(
- rad.values[self.scene_idx], thr
+ rad.values[self.scene_idx], thr, power
)
# these scenes have not been implemented yet... and it might not be
@@ -939,12 +952,13 @@ class CloudTests:
"NIR_Reflectance_Test",
self.scene_idx,
)
+ power = threshold["thr"][3]
logger.info(f'Running test NIR_Reflectance_Test on "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
if np.ndim(thr) == 1:
thr = thr.reshape(thr.shape[0], 1)
bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=")
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr)
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr, power)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -1028,8 +1042,9 @@ class CloudTests:
)
logger.info(f'Running test 1.6_2.1um_NIR_Reflectance_Test on "{self.scene_name}"\n')
rad = self.data[band].values[self.scene_idx]
+ power = threshold["thr"][3]
bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=")
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr)
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr, power)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -1063,7 +1078,7 @@ class CloudTests:
locut = np.array(thr[0, :], dtype=np.float32)
midpt = np.array(thr[1, :], dtype=np.float32)
hicut = np.array(thr[2, :], dtype=np.float32)
- power = np.array(thr[3, :], dtype=np.float32)
+ power = threshold["thr"][3]
kwargs["confidence"][self.scene_idx] = anc.py_conf_test(
f_rad, locut, hicut, power, midpt
)
@@ -1087,6 +1102,7 @@ class CloudTests:
bits["qa"][scene_idx] = 1
thr = preproc.gemi_thresholds(self.data, threshold, self.scene_name, scene_idx)
rad = self.data[band].values[scene_idx]
+ power = threshold["gemi0"][3]
logger.info(f'Running test GEMI_Test on "{self.scene_name}"\n')
if np.ndim(thr) == 1:
@@ -1094,7 +1110,7 @@ class CloudTests:
bits["test"][scene_idx] = self.compute_test_bits(
rad, thr[1, :], ">", scene_idx=scene_idx
)
- kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr)
+ kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr, power)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -1130,8 +1146,9 @@ class CloudTests:
)
)
bits["qa"][scene_idx] = 1
- thr = np.full((len(scene_idx[0]), 4), threshold["thr"][:4]).T
+ thr = np.full((len(scene_idx[0]), 3), threshold["thr"][:3]).T
+ power = threshold["power"]
logger.info(f'Running test 1.38um_High_Cloud_Test on "{self.scene_name}"\n')
rad = self.data[band].values[scene_idx]
if np.ndim(thr) == 1:
@@ -1139,7 +1156,7 @@ class CloudTests:
bits["test"][scene_idx] = self.compute_test_bits(
rad, thr[1, :], "<=", scene_idx=scene_idx
)
- kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr)
+ kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr, power)
cmin = np.fmin(cmin, kwargs["confidence"])
@@ -1164,9 +1181,10 @@ class CloudTests:
# rad = self.data[band].values[self.scene_idx]
thr = threshold["thr"]
rad = self.data.M13.values[self.scene_idx] - self.data.M16.values[self.scene_idx]
+ power = threshold["thr"][3]
logger.info(f'Running test 4-12um_BTD_Thin_Cirrus_Test on "{self.scene_name}"\n')
bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1], "<=")
- kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr)
+ kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr, power)
cmin = np.fmin(cmin, kwargs["confidence"])
diff --git a/mvcm/utils.py b/mvcm/utils.py
index 6c8551c..e9593a2 100644
--- a/mvcm/utils.py
+++ b/mvcm/utils.py
@@ -22,7 +22,7 @@ def sunglint_scene(refang, sunglint_thr):
def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
band = f"band{band_n}"
- sunglint_thr = np.zeros((4, len(refang)))
+ sunglint_thr = np.zeros((3, len(refang)))
# if refang > thresholds['bounds'][3]:
# sunglint = sunglint
# # dosgref[2] = hicnf
@@ -42,7 +42,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
hi_ang = thresholds["bounds"][2]
lo_ang_val = thresholds[f"{band}_10deg"][0]
hi_ang_val = thresholds[f"{band}_10deg"][1]
- power = thresholds[f"{band}_10deg"][3]
+ # power = thresholds[f"{band}_10deg"][3]
conf_range = thresholds[f"{band}_10deg"][2]
# elif (refang > thresholds['bounds'][2] and refang <= thresholds['bounds'][3]):
@@ -51,7 +51,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
hi_ang = thresholds["bounds"][3]
lo_ang_val = thresholds[f"{band}_20deg"][0]
hi_ang_val = thr[1]
- power = thresholds[f"{band}_20deg"][3]
+ # power = thresholds[f"{band}_20deg"][3]
conf_range = thresholds[f"{band}_20deg"][2]
else:
raise ValueError("Wrong sunglint flag")
@@ -61,7 +61,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
sunglint_thr[1, :] = midpt
sunglint_thr[2, :] = midpt - conf_range
sunglint_thr[0, :] = midpt + conf_range
- sunglint_thr[3, :] = power
+ # sunglint_thr[3, :] = power
return sunglint_thr
diff --git a/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
index 848bfc4..ecf5651 100644
--- a/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -8,6 +8,7 @@ Land_Day:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
+ power: 1.0
perform: True
11-4um_Oceanic_Stratus_Test:
thr: [-16.0, -14.0, -12.0, 1.0, 1.0]
@@ -21,6 +22,7 @@ Land_Day:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
dlvrat: [1.80, 1.85, 1.90, 1.0, 1.0]
@@ -32,6 +34,7 @@ Land_Day:
Land_Night:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.0
adj: 0
bt1: 270.0 # WILL FIND A BETTER NAME AT SOME POINT, MAYBE
@@ -81,6 +84,7 @@ Land_Night:
Land_Day_Coast:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.25
perform: True
@@ -98,6 +102,7 @@ Land_Day_Coast:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
# dl11_12lcmult_t2 : 0.3
@@ -107,6 +112,7 @@ Land_Day_Desert:
# lds11_12hi : [3.5, 1.0]
11-12um_Cirrus_Test:
coeffs: [3.5, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.25
perform: True
@@ -145,6 +151,7 @@ Land_Day_Desert:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
GEMI_Test:
@@ -160,6 +167,7 @@ Land_Day_Desert_Coast:
# lds11_12hi_c : [3.5, 1.0]
11-12um_Cirrus_Test:
coeffs: [3.5, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.25
perform: True
@@ -184,6 +192,7 @@ Land_Day_Desert_Coast:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
ldsbt1_c: 320.0
@@ -203,6 +212,7 @@ Ocean_Day:
perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.25
perform: True
@@ -253,6 +263,7 @@ Ocean_Day:
1.38um_High_Cloud_Test:
coeffs:
[-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204]
+ power: 1.0
adj: 0.0010
szafac: 4.0
perform: True
@@ -272,6 +283,7 @@ Ocean_Night:
perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.0
perform: True
@@ -307,6 +319,7 @@ Ocean_Night:
Polar_Day_Land:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 0.3
perform: False
@@ -320,13 +333,15 @@ Polar_Day_Land:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
Polar_Night_Land:
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
- 11-12um_Cirrus_Test:
+ 1single_1-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
@@ -392,12 +407,14 @@ Polar_Day_Coast:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
Polar_Day_Desert:
11-12um_Cirrus_Test:
coeffs: [3.5, 1.0]
+ power: 1.0
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
perform: False
@@ -419,6 +436,7 @@ Polar_Day_Desert:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
GEMI_Test:
@@ -430,6 +448,7 @@ Polar_Day_Desert:
Polar_Day_Desert_Coast:
11-12um_Cirrus_Test:
coeffs: [3.5, 1.0]
+ power: 1.0
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
perform: False
@@ -451,18 +470,21 @@ Polar_Day_Desert_Coast:
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
+ power: 1.0
tpw: -10.00
perform: True
Polar_Day_Snow:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
perform: False
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.38um_High_Cloud_Test:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
+ power: 1.0
tpw: 0.75
perform: True
# 11-4um_BT_Difference_Test_Land:
@@ -489,6 +511,7 @@ Polar_Day_Snow:
Polar_Night_Snow:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
@@ -546,6 +569,7 @@ Polar_Day_Ocean:
perform: False
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.25
perform: False
@@ -589,6 +613,7 @@ Polar_Day_Ocean:
1.38um_High_Cloud_Test:
coeffs:
[-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204]
+ power: 1.0
adj: 0.0010
szafac: 4.0
perform: True
@@ -606,6 +631,7 @@ Polar_Night_Ocean:
perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 1.0
perform: True
@@ -636,6 +662,7 @@ Polar_Night_Ocean:
Day_Snow:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0.3
adj: 0.8
perform: True
@@ -643,6 +670,7 @@ Day_Snow:
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.38um_High_Cloud_Test:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
+ power: 1.0
tpw: 0.75
perform: True
# 11-4um_BT_Difference_Test_Land:
@@ -678,6 +706,7 @@ Day_Snow:
Night_Snow:
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
+ power: 1.0
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
--
GitLab