add s1 indices

src/cross-level/cite-cl.c

@@ -400,6 +400,12 @@ cite_t _cite_me_[_CITE_LENGTH_] = {
     "reflectance at any given time. Remote Sensing of Environment, 162, 67-83."
     "http://dx.doi.org/10.1016/j.rse.2015.02.009",
     false }
+  { "RVI",
+    "Agapiou, A. (2020). Estimating Proportion of Vegetation Cover at the Vicinity "
+    "of Archaeological Sites Using Sentinel-1 and -2 Data, Supplemented by Crowdsourced "
+    "OpenStreetMap Geodata. Applied Sciences, 10, 4764."
+    "https://doi.org/10.3390/app10144764",
+    false }
 };
 
 

src/cross-level/cite-cl.h

@@ -60,7 +60,7 @@ enum { _CITE_FORCE_,       _CITE_L2PS_,       _CITE_ATMVAL_,
        _CITE_NDVIre2_,     _CITE_NDVIre3_,    _CITE_NDVIre1n_,
        _CITE_NDVIre2n_,    _CITE_NDVIre3n_,   _CITE_MSRre_,
        _CITE_MSRren_,      _CITE_CCI_,        _CITE_EV2_,
-       _CITE_HARMONIC_,    _CITE_LENGTH_ };
+       _CITE_HARMONIC_,    _CITE_RVI_,        _CITE_LENGTH_ };
 
 typedef struct {
   char description[NPOW_10];

src/cross-level/enum-cl.c

@@ -79,7 +79,8 @@ const tagged_enum_t _TAGGED_ENUM_IDX_[_IDX_LENGTH_] = {
   { _IDX_CRE_, "CIre"   }, { _IDX_NR1_, "NDVIre1"   }, { _IDX_NR2_, "NDVIre2"      },
   { _IDX_NR3_, "NDVIre3"}, { _IDX_N1n_, "NDVIre1n"  }, { _IDX_N2n_, "NDVIre2n"     },
   { _IDX_N3n_, "NDVIre3n"},{ _IDX_Mre_, "MSRre"     }, { _IDX_Mrn_, "MSRren"       },
-  { _IDX_CCI_, "CCI"    }, { _IDX_EV2_, "EVI2"      }, { _IDX_CSW_, "CRemoveSWIR1" }};
+  { _IDX_CCI_, "CCI"    }, { _IDX_EV2_, "EVI2"      }, { _IDX_BCR_, "BCR"          },
+  { _IDX_RVI_, "RVI"    }, { _IDX_CSW_, "CRemoveSWIR1" }};
 
 const tagged_enum_t _TAGGED_ENUM_INT_[_INT_LENGTH_] = {
   { _INT_NONE_,     "NONE"    }, { _INT_LINEAR_, "LINEAR" },

src/cross-level/enum-cl.h

@@ -157,7 +157,7 @@ enum { _IDX_BLU_, _IDX_GRN_, _IDX_RED_, _IDX_NIR_, _IDX_SW1_, _IDX_SW2_,
        _IDX_SMA_, _IDX_BVV_, _IDX_BVH_, _IDX_NDT_, _IDX_NDM_, _IDX_SW0_,
        _IDX_KNV_, _IDX_ND1_, _IDX_ND2_, _IDX_CRE_, _IDX_NR1_, _IDX_NR2_,
        _IDX_NR3_, _IDX_N1n_, _IDX_N2n_, _IDX_N3n_, _IDX_Mre_, _IDX_Mrn_,
-       _IDX_CCI_, _IDX_EV2_, _IDX_CSW_, _IDX_LENGTH_};
+       _IDX_CCI_, _IDX_EV2_, _IDX_CSW_, _IDX_BCR_, _IDX_RVI_, _IDX_LENGTH_};
 
 // standardization
 enum { _STD_NONE_, _STD_NORMAL_, _STD_CENTER_, _STD_LENGTH_ };

src/higher-level/index-hl.c

@@ -31,6 +31,8 @@ This file contains functions for computing spectral indices
 #include "gsl/gsl_blas.h"
 #include "gsl/gsl_linalg.h"
 
+#include <math.h>
+
 
 enum { TCB, TCG, TCW, TCD};
 
@@ -242,6 +244,114 @@ float ind, scale = 1000.0;
 }
 
 
+/** This function computes a spectral index time series, with a simple difference,
++++ e.g. BCR: VH[dB]-VV[dB]
+--- ard:    ARD
+--- mask_:  mask image
+--- ts:     pointer to instantly useable TSA image arrays
+--- b1:     band 1
+--- b2:     band 2
+--- nc:     number of cells
+--- nt:     number of ARD products over time
+--- nodata: nodata value
++++ Return: void
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
+void index_simple_difference(ard_t *ard, small *mask_, tsa_t *ts, int b1, int b2, int nc, int nt, short nodata){
+int p, t;
+float ind;
+
+
+  #pragma omp parallel private(t,ind) shared(ard,mask_,ts,b1,b2,nc,nt,nodata) default(none)
+  {
+
+    #pragma omp for
+    for (p=0; p<nc; p++){
+
+      if (mask_ != NULL && !mask_[p]){
+        for (t=0; t<nt; t++) ts->tss_[t][p] = nodata;
+        continue;
+      }
+
+      for (t=0; t<nt; t++){
+
+        if (!ard[t].msk[p]){
+          ts->tss_[t][p] = nodata;
+        } else {
+          ind = (ard[t].dat[b1][p] - ard[t].dat[b2][p]);
+          if (ind > SHRT_MAX || ind < SHRT_MIN){
+            ts->tss_[t][p] = nodata;
+          } else {
+            ts->tss_[t][p] = (short)(ind);
+          }
+        }
+
+      }
+
+    }
+  }
+
+  return;
+}
+
+
+/** This function computes an RVI-like time series,
++++ RVI: sqrt(1-(VV/(VV+VH)))*(4*VH/(VH+VV))
+--- ard:    ARD
+--- mask_:  mask image
+--- ts:     pointer to instantly useable TSA image arrays
+--- b1:     band 1
+--- b2:     band 2
+--- nc:     number of cells
+--- nt:     number of ARD products over time
+--- nodata: nodata value
++++ Return: void
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
+void index_rvi(ard_t *ard, small *mask_, tsa_t *ts, int b1, int b2, int nc, int nt, short nodata){
+int p, t;
+float vv_db, vh_db, vv_lin, vh_lin, dop, ind, scale = 10000.0;
+
+
+  #pragma omp parallel private(t,vv_db,vh_db,vv_lin,vh_lin,dop,ind) shared(ard,mask_,ts,b1,b2,nc,nt,nodata,scale) default(none)
+  {
+
+    #pragma omp for
+    for (p=0; p<nc; p++){
+
+      if (mask_ != NULL && !mask_[p]){
+        for (t=0; t<nt; t++) ts->tss_[t][p] = nodata;
+        continue;
+      }
+
+      for (t=0; t<nt; t++){
+
+        if (!ard[t].msk[p]){
+          ts->tss_[t][p] = nodata;
+        } else {
+          if (ard[t].dat[b2][p] != 0){
+            vv_db = (float)ard[t].dat[b1][p] / 100;
+            vh_db = (float)ard[t].dat[b2][p] / 100;
+            vv_lin = powf(10, vv_db / 10);
+            vh_lin = powf(10, vh_db / 10);
+            dop = vv_lin/(vv_lin+vh_lin);
+            ind = sqrt(1-dop)*(4*vh_lin/(vh_lin+vv_lin));
+            if (ind*scale > SHRT_MAX || ind*scale < SHRT_MIN){
+              ts->tss_[t][p] = nodata;
+            } else {
+              ts->tss_[t][p] = (short)(ind*scale);
+            }
+          } else {
+            ts->tss_[t][p] = nodata;
+          }
+        }
+
+      }
+
+    }
+  }
+
+  return;
+}
+
 /** This function computes an MSRre-like time series,
 +++ MSRre: ((b1/b2)-1)/sqrt((b1/b2)+1)
 --- ard:    ARD
@@ -1045,6 +1155,13 @@ int tsa_spectral_index(ard_t *ard, tsa_t *ts, small *mask_, int nc, int nt, int
       index_cont_remove(ard, mask_, ts, sen->swir1, sen->nir, sen->swir2, 
                         sen->w_swir1, sen->w_nir, sen->w_swir2, nc, nt, nodata);
       break;
+    case _IDX_BCR_:
+      index_simple_difference(ard, mask_, ts, sen->vh, sen->vv, nc, nt, nodata);
+      break;
+    case _IDX_RVI_:
+      cite_me(_CITE_RVI_);
+      index_rvi(ard, mask_, ts, sen->vv, sen->vh, nc, nt, nodata);
+      break;
     default:
       printf("unknown INDEX\n");
       break;      

src/higher-level/param-hl.c

@@ -690,6 +690,14 @@ int *band_ptr[_WVL_LENGTH_] = {
         v[_WVL_NIR_] = v[_WVL_SWIR1_] = v[_WVL_SWIR2_] = true;
         copy_string(tsa->index_name[idx], NPOW_02, "CSW");
         break;
+      case _IDX_BCR_:
+        v[_WVL_VV_] = v[_WVL_VH_] = true;
+        copy_string(tsa->index_name[idx], NPOW_02, "BCR");
+        break;
+      case _IDX_RVI_:
+        v[_WVL_VV_] = v[_WVL_VH_] = true;
+        copy_string(tsa->index_name[idx], NPOW_02, "RVI");
+        break;
       default:
         printf("unknown INDEX\n");
         break;