added get_b1_threshold(). needs to be tested

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