From 472b2d9c0505c7b4b36a2ead04b025f223f896e2 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Fri, 3 Jun 2022 13:25:56 +0000
Subject: [PATCH] added get_b1_threshold(). needs to be tested
---
utils.py | 75 ++++++++++++++++++++++++++++++++++++++++++++++++--------
1 file changed, 65 insertions(+), 10 deletions(-)
diff --git a/utils.py b/utils.py
index 6734164..21fb39b 100644
--- a/utils.py
+++ b/utils.py
@@ -227,16 +227,71 @@ def conf_test_dble(rad, coeffs):
return confidence
-def get_b1_threshold(thresholds):
+def get_b1_threshold(thresholds, ndvi_thresholds):
# this needs to be read from the thresholds file
des_ndvi = 0.25
pxin_ndvibk = 0 # placeholder until I figure out how to pass this value
-
- if pxin_ndvibk < des_ndvi:
- coeffs = thresholds['Coeffs_Band8_land_thresh']
- else:
- coeffs = thresholds['Coeffs_Band1_land_thresh']
-
- # to be continued ...
-
- return coeffs
+ pxin_sctang = 0 # placeholder until I figure out how to pass this value
+ delta_ndvi_bin = 0.1 # defined as is in the original C code
+ thr = np.zeros((4, len(pxin_ndvibk)))
+ thr_adj = np.zeros((4, len(pxin_ndvibk)))
+
+ coeffs = np.full((10, 3, 4, len(pxin_ndvibk)))
+
+ indvi = np.zeros((len(pxin_ndvibk), ))
+ x = np.zeros((len(pxin_ndvibk), ))
+ y1 = np.zeros((len(pxin_ndvibk), ))
+ y2 = np.zeros((len(pxin_ndvibk), ))
+ interp = np.zeros((len(pxin_ndvibk), ))
+
+ # Check for special cases and fill values ('ndvibk'=32.767)
+ # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
+ # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
+ # (pxin_ndvibk < ndvi_thresholds['ndvi_bnd1']))
+ idx = np.nonzero(pxin_ndvibk < ndvi_thresholds['ndvi_bnd1'])
+ # all gets set to zero, so no need to do anything for this idx
+
+ # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
+ # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
+ # (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2']))
+ idx = np.nonzero(pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2'])
+ indvi[idx] = 9
+
+ # else
+ # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
+ # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
+ # (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) &
+ # (pxin_ndvibk < ndvi_thresholds['ndvi_bnd2']))
+ idx = np.nonzero((pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) &
+ (pxin_ndvibk < ndvi_thresholds['ndvi_bnd2']))
+ interp[idx] = 1 # this is not really needed anymore, maybe
+ indvi = (pxin_ndvibk[idx] / delta_ndvi_bin) - 0.5
+ indvi[indvi < 0] = 0
+ x = (pxin_ndvibk - delta_ndvi_bin*indvi + delta_ndvi_bin/2) / delta_ndvi_bin
+ x[x < 0] = 0
+ x[x > 1] = 1
+
+ for i in range(3):
+ # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
+ # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]))
+
+ y1[:, i] = (coeffs[indvi, i, 0] +
+ coeffs[indvi, i, 1]*pxin_sctang +
+ coeffs[indvi, i, 2]*pxin_sctang**2 +
+ coeffs[indvi, i, 3]*pxin_sctang**3)
+
+ # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
+ # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & interp == 1)
+ y2[idx, i] = (coeffs[indvi[idx]+1, i, 0] +
+ coeffs[indvi[idx]+1, i, 1]*pxin_sctang[idx] +
+ coeffs[indvi[idx]+1, i, 2]*pxin_sctang[idx]**2 +
+ coeffs[indvi[idx]+1, i, 3]*pxin_sctang[idx]**3)
+
+ thr[:, i] = ((1.0 - x) * y1) + (x * y2)
+
+ thr_adj[pxin_ndvibk < des_ndvi, i] = thr[pxin_ndvibk < des_ndvi, 0]*ndvi_thresholds['adj_fac_desert']
+ thr_adj[pxin_ndvibk >= des_ndvi, i] = thr[pxin_ndvibk >= des_ndvi, 0]*ndvi_thresholds['adj_fac_land']
+
+ thr[:, 3] = 2.0
+
+ return thr
--
GitLab