diff --git a/mvcm/__init__.py b/mvcm/__init__.py
index 3a4bfdbbff296a6f6bca5031a18989fb4c1d2301..6725e40baebbe4cf7366b09429cafc5dc37c29b7 100644
--- a/mvcm/__init__.py
+++ b/mvcm/__init__.py
@@ -1,2 +1,2 @@
"""Set up version."""
-__version__ = "0.2.1"
+__version__ = "0.2.2"
diff --git a/mvcm/c_tools/get_b1_thresholds.c b/mvcm/c_tools/get_b1_thresholds.c
index 694345a9eac3bfed60a6c182bd8d0e4dd2f37c60..6255d44738c4715c5015c9167610501932b97362 100644
--- a/mvcm/c_tools/get_b1_thresholds.c
+++ b/mvcm/c_tools/get_b1_thresholds.c
@@ -32,13 +32,13 @@ calls:
#include <stdio.h>
#include <math.h>
-#include "thresholds.h"
-#include "pixel.h"
+// #include "thresholds.h"
+// #include "pixel.h"
/**********************************************************************/
-int get_b1_thresholds(float *locut, float *hicut, float *midpt,
+int get_b1_thresholds(float ndvibk, float sctang, float *locut, float *midpt, float *hicut,
float *power)
@@ -53,62 +53,118 @@ int get_b1_thresholds(float *locut, float *hicut, float *midpt,
int interp=0;
int a, i, j, k, indvi, icnf, irtn;
+ float des_ndvi = 0.3;
static float delta_ndvi_bin = 0.1;
float x, x1, x2, y1, y2;
float coeff[10][3][4], thr[3], thr_adj[3];
+ float thr_adj_fac_land = 0.0;
+ float thr_adj_fac_desert = 0.1;
+
+ float coeffs_land[10][3][4] = {
+ {
+ {32.00000000, 0.00000000, 0.00000000, 0.00000000},
+ {42.00000000, 0.00000000, 0.00000000, 0.00000000},
+ {52.00000000, 0.00000000, 0.00000000, 0.00000000}
+ },
+ {
+ {24.00000000, 0.00000000, 0.00000000, 0.00000000},
+ {28.00000000, 0.00000000, 0.00000000, 0.00000000},
+ {32.00000000, 0.00000000, 0.00000000, 0.00000000}
+ },
+ {
+ {99.13076923, -2.00907925, 0.01492075, -0.00003531},
+ {122.19090909, -2.32652292, 0.01659848, -0.00003681},
+ {142.66293706, -2.57860528, 0.01773252, -0.00003685}
+ },
+ {
+ {85.07902098, -1.59413364, 0.01123310, -0.00002556},
+ {144.56573427, -2.81054779, 0.01967366, -0.00004324},
+ {204.35454545, -4.03411810, 0.02816667, -0.00006103}
+ },
+ {
+ {85.03846154, -1.50831391, 0.01006760, -0.00002199},
+ {165.15314685, -3.24716783, 0.02255594, -0.00004965},
+ {242.06363636, -4.90912587, 0.03445455, -0.00007587}
+ },
+ {
+ {81.00979021, -1.37731935, 0.00881294, -0.00001859},
+ {220.36783217, -4.44111888, 0.03087762, -0.00006888},
+ {359.72587413, -7.50491841, 0.05294231, -0.00011917}
+ },
+ {
+ {76.94055944, -1.35441725, 0.00896096, -0.00001952},
+ {172.36783217, -3.33144911, 0.02242308, -0.00004810},
+ {267.90909091, -5.31620047, 0.03597727, -0.00007698}
+ },
+ {
+ {85.83006993, -1.55480575, 0.01025932, -0.00002216},
+ {160.73706294, -3.07291375, 0.02041900, -0.00004330},
+ {237.37622378, -4.63444833, 0.03091317, -0.00006525}
+ },
+ {
+ {105.02447552, -1.98017094, 0.01319522, -0.00002877},
+ {135.50699301, -2.59097902, 0.01749301, -0.00003811},
+ {165.33006993, -3.18872183, 0.02171387, -0.00004734}
+ },
+ {
+ {105.02447552, -1.98017094, 0.01319522, -0.00002877},
+ {135.50699301, -2.59097902, 0.01749301, -0.00003811},
+ {165.33006993, -3.18872183, 0.02171387, -0.00004734}
+ }
+ };
+
+ float coeffs_desert[3][3][4] = {
+ {
+ {282.74916084, -5.42869658, 0.03660781, -0.00008092},
+ {344.87048951, -6.80660839, 0.04770163, -0.00010991},
+ {407.83734266, -8.20392385, 0.05894114, -0.00013926},
+ },
+ {
+ {229.50727273, -4.31606061, 0.02868182, -0.00006212},
+ {316.38517483, -6.41910256, 0.04603089, -0.00010787},
+ {403.78188811, -8.52743395, 0.06335781, -0.00015340},
+ },
+ {
+ {239.32391608, -4.49928127, 0.02977214, -0.00006436},
+ {270.31741259, -5.33625680, 0.03757168, -0.00008609},
+ {301.76013986, -6.18557498, 0.04546562, -0.00010806},
+ }
+ };
-// Transfer all coefficients to one array.
- i = 0;
+// Transfer all coefficients to one array.
+ for (i=0; i<10; i++) {
for (j=0; j<3; j++) {
- for( k=0; k<4; k++) {
-
- if(pxin.ndvibk < des_ndvi[0]) {
-
- coeff[0][j][k] = b8ndvi0[i];
- coeff[1][j][k] = b8ndvi1[i];
- coeff[2][j][k] = b8ndvi2[i];
-// printf("coeffs8: %d %d %d %f %f\n",i, j, k, coeff[2][j][k], b8ndvi2[i]);
-
- }
- else {
-
- coeff[0][j][k] = b1ndvi0[i];
- coeff[1][j][k] = b1ndvi1[i];
- coeff[2][j][k] = b1ndvi2[i];
- coeff[3][j][k] = b1ndvi3[i];
- coeff[4][j][k] = b1ndvi4[i];
- coeff[5][j][k] = b1ndvi5[i];
- coeff[6][j][k] = b1ndvi6[i];
- coeff[7][j][k] = b1ndvi7[i];
- coeff[8][j][k] = b1ndvi8[i];
- coeff[9][j][k] = b1ndvi9[i];
-// printf("coeffs1: %d %d %d %f %f\n",i, j, k, coeff[2][j][k], b1ndvi2[i]);
-
- }
-
- i++;
+ for (k=0; k<4; k++) {
+ if(ndvibk < des_ndvi) {
+ coeff[i][j][k] = coeffs_desert[i][j][k];
+ if (i > 3) {
+ break;
+ }
+ }
+ else {
+ coeff[i][j][k] = coeffs_land[i][j][k];
}
}
-
-// printf("\nband 1, 8 coefficients initialized \n");
+ }
+ }
/* Compute thresholds for current pixel as functions of background
NDVI and scattering angle.
*/
// Check for special cases and fill values ('ndvibk'=32.767).
- if(pxin.ndvibk < fill_ndvi1[0] && pxin.ndvibk > fill_ndvi2[0]) {
+ if(ndvibk < 32.767 && ndvibk > 0) {
- if(pxin.ndvibk < ndvi_bnd1[0]) {
+ if(ndvibk <0.05) {
indvi = 0;
x = 0.0;
y2 = 0.0;
interp = 0;
}
- else if(pxin.ndvibk >= ndvi_bnd2[0]) {
+ else if(ndvibk >= 0.95) {
indvi = 9;
x = 0.0;
y2 = 0.0;
@@ -119,18 +175,17 @@ int get_b1_thresholds(float *locut, float *hicut, float *midpt,
}
if(interp) {
- a = ((pxin.ndvibk / delta_ndvi_bin) - 0.5);
+ a = ((ndvibk / delta_ndvi_bin) - 0.5);
if(a < 0) {
indvi = 0;
}
else {
indvi = a;
}
-// indvi = min(nbin_ndvi, max(0, int((pxin.ndvibk / delta_ndvi_bin) - 0.5)) );
x1 = delta_ndvi_bin*indvi + delta_ndvi_bin/2.0;
x2 = x1 + delta_ndvi_bin;
- x = (pxin.ndvibk - x1) / ( x2 - x1 );
+ x = (ndvibk - x1) / ( x2 - x1 );
if(x < 0.0) x = 0.0;
if(x > 1.0) x = 1.0;
@@ -139,32 +194,24 @@ int get_b1_thresholds(float *locut, float *hicut, float *midpt,
for( icnf=0; icnf<3; icnf++) {
y1 = coeff[indvi][icnf][0] +
- (coeff[indvi][icnf][1] * pxin.sctang) +
- (coeff[indvi][icnf][2] * powf(pxin.sctang,2)) +
- (coeff[indvi][icnf][3] * powf(pxin.sctang,3));
+ (coeff[indvi][icnf][1] * sctang) +
+ (coeff[indvi][icnf][2] * powf(sctang,2)) +
+ (coeff[indvi][icnf][3] * powf(sctang,3));
if(interp) {
y2 = coeff[indvi+1][icnf][0] +
- (coeff[indvi+1][icnf][1] * pxin.sctang) +
- (coeff[indvi+1][icnf][2] * powf(pxin.sctang,2)) +
- (coeff[indvi+1][icnf][3] * powf(pxin.sctang,3));
+ (coeff[indvi+1][icnf][1] * sctang) +
+ (coeff[indvi+1][icnf][2] * powf(sctang,2)) +
+ (coeff[indvi+1][icnf][3] * powf(sctang,3));
}
thr[icnf] = ((1.0 - x) * y1) + (x * y2);
- if(pxin.ndvibk < des_ndvi[0]) {
- thr_adj[icnf] = thr[icnf] * thr_adj_fac_desert[0];
+ if(ndvibk < des_ndvi) {
+ thr_adj[icnf] = thr[icnf] * thr_adj_fac_desert;
}
else {
- thr_adj[icnf] = thr[icnf] * thr_adj_fac_land[0];
+ thr_adj[icnf] = thr[icnf] * thr_adj_fac_land;
}
-
-// printf("ndvi thr: %f %f %f %f %f %f %f\n", des_ndvi[0], fill_ndvi1[0], fill_ndvi2[0],
-// ndvi_bnd1[0], ndvi_bnd2[0], thr_adj_fac_desert[0], thr_adj_fac_land[0]);
-// printf("coeffs: %f %f %f %f %f %f %f %f\n", coeff[indvi][icnf][0], coeff[indvi][icnf][1],
-// coeff[indvi][icnf][2], coeff[indvi][icnf][3], coeff[indvi+1][icnf][0], coeff[indvi+1][icnf][1],
-// coeff[indvi+1][icnf][2], coeff[indvi+1][icnf][3]);
-// printf("interpolation: %d %d %f %f %f %f %f\n", indvi, icnf, pxin.sctang, y1, y2, thr[icnf], thr_adj[icnf]);
-
}
*(hicut) = (thr[0] + thr_adj[0]) / 100.0;
@@ -182,10 +229,6 @@ int get_b1_thresholds(float *locut, float *hicut, float *midpt,
}
-// printf("\nSubroutine get_b1_thresholds\n");
-// printf("%d %d %d %f %f %f %f %f %f %f %f \n", irtn, interp, indvi, pxin.ndvibk, pxin.sctang,
-// x1, x2, x, *(locut), *(midpt), *(hicut));
-
return irtn;
}