NASA Ocean Color

ocssw V2022

 /*---------------------------------------------------------------------*/
 /* get_ndvi.c - vegetation index classification for MSl12.             */
 /*                                                                     */
 /* Inputs:                                                             */
 /*     l2rec - level-2 structure containing one complete scan after    */
 /*             atmospheric correction.                                 */
 /* Outputs:                                                            */
 /*     ndvi  - vegetation index for land, 1 value per pixel.           */
 /*                                                                     */
 /* Written by: Bryan Franz, SAIC-GSC, February 2000                    */
 /*                                                                     */
 /*---------------------------------------------------------------------*/
  
 #include <stdlib.h>
 #include <math.h>
 #include "l12_proto.h"
  
 static float undef = -2.0;
 static int ibblue = -1;
 static int ibred = -1;
 static int ibnir = -1;
  
 //static float minval = -2.0;
 //static float maxval =  2.0;
  
 // as per C. Tucker (11/2014), no cut-off at -2
 static float minval = -1000.0;
 static float maxval = 1000.0;
 static int32_t mask = LAND;
  
 void get_ndvi(l1str *l1rec, float ndvi[]) {
     int32_t ip, ipb;
     float red;
     float nir;
  
     if (ibred < 0 || ibnir < 0) {
         printf("NDVI requires bands near 670 and 865nm\n");
         exit(1);
     }
  
     for (ip = 0; ip < l1rec->npix; ip++) {
         ndvi[ip] = undef;
  
         ipb = l1rec->l1file->nbands*ip;
         red = l1rec->rhos[ipb + ibred];
         nir = l1rec->rhos[ipb + ibnir];
  
         if (l1rec->dem[ip] < -500 ||
                 (l1rec->flags[ip] & mask) == 0 ||
                 red <= 0.0 || nir <= 0.0 ||
                 red > 1.0 || nir > 1.0) {
             l1rec->flags[ip] |= PRODFAIL;
             continue;
         } else {
             ndvi[ip] = (nir - red)
                     / (nir + red);
  
             ndvi[ip] = MIN(MAX(ndvi[ip], minval), maxval);
         }
     }
 }
  
 void get_evi(l1str *l1rec, float evi[]) {
     static float L = 1.0, c1 = 6.0, c2 = 7.5;
  
     int32_t ip, ipb;
     float blu;
     float red;
     float nir;
     double val;
  
     if (ibblue < 0 || ibred < 0 || ibnir < 0) {
         printf("EVI requires bands near 412, 670, and 865nm\n");
         exit(1);
     }
  
     for (ip = 0; ip < l1rec->npix; ip++) {
         evi[ip] = undef;
  
         ipb = l1rec->l1file->nbands*ip;
         blu = l1rec->rhos[ipb + ibblue];
         red = l1rec->rhos[ipb + ibred];
         nir = l1rec->rhos[ipb + ibnir];
  
         if (l1rec->dem[ip] < -500 ||
                 (l1rec->flags[ip] & mask) == 0 ||
                 red <= 0.0 || nir <= 0.0 ||
                 red > 1.0 || nir > 1.0) {
             l1rec->flags[ip] |= PRODFAIL;
             continue;
  
         } else {
  
             if (blu > 0.0 && (blu <= red || red <= nir)) {
  
                 /* Most cases - EVI formula */
  
                 if ((val = L + nir + c1 * red - c2 * blu) == 0)
                     continue;
                 else
                     evi[ip] = 2.5 * (nir - red) / val;
  
             } else {
  
                 /* Backup - SAVI formula */
  
                 if ((val = 0.5 + nir + red) == 0)
                     continue;
                 else
                     evi[ip] = 1.5 * (nir - red) / val;
             }
             evi[ip] = MIN(MAX(evi[ip], minval), maxval);
  
         }
     }
 }
  
 // From Compton Tucker 07/30/2016
 //EVI3=2.5*(nir-red)/(nir+6*red-7.5*blue+1)                 
 //EVI2=2.5*(NIR-Red)/(NIR+2.4*Red+1)
  
 void get_evi2(l1str *l1rec, float evi2[]) {
     int32_t ip, ipb;
     float red;
     float nir;
  
     if (ibred < 0 || ibnir < 0) {
         printf("EVI2 requires bands near 670, and 865nm\n");
         exit(1);
     }
  
     for (ip = 0; ip < l1rec->npix; ip++) {
         evi2[ip] = undef;
  
         ipb = l1rec->l1file->nbands*ip;
         red = l1rec->rhos[ipb + ibred];
         nir = l1rec->rhos[ipb + ibnir];
  
         if (l1rec->dem[ip] < -500 ||
                 (l1rec->flags[ip] & mask) == 0 ||
                 red <= 0.0 || nir <= 0.0 ||
                 red > 1.0 || nir > 1.0) {
             l1rec->flags[ip] |= PRODFAIL;
             continue;
  
         } else {
  
             if (red <= nir) {
                 evi2[ip] = 2.5 * (nir - red) / (nir + 2.4 * red + 1);
             }
             evi2[ip] = MIN(MAX(evi2[ip], minval), maxval);
         }
     }
 }
  
 void get_evi3(l1str *l1rec, float evi3[]) {
     static float L = 1.0, c1 = 6.0, c2 = 7.5;
  
     int32_t ip, ipb;
     float blu;
     float red;
     float nir;
     double val;
  
     if (ibblue < 0 || ibred < 0 || ibnir < 0) {
         printf("EVI requires bands near 412, 670, and 865nm\n");
         exit(1);
     }
  
     for (ip = 0; ip < l1rec->npix; ip++) {
         evi3[ip] = undef;
  
         ipb = l1rec->l1file->nbands*ip;
         blu = l1rec->rhos[ipb + ibblue];
         red = l1rec->rhos[ipb + ibred];
         nir = l1rec->rhos[ipb + ibnir];
  
         if (l1rec->dem[ip] < -500 ||
                 (l1rec->flags[ip] & mask) == 0 ||
                 red <= 0.0 || nir <= 0.0 ||
                 red > 1.0 || nir > 1.0) {
             l1rec->flags[ip] |= PRODFAIL;
             continue;
  
         } else {
  
  
             if (blu > 0.0 && (blu <= red || red <= nir)) {
  
                 /* Most cases - EVI formula */
  
                 if ((val = L + nir + c1 * red - c2 * blu) == 0)
                     continue;
                 else
                     evi3[ip] = 2.5 * (nir - red) / val;
  
             }
             evi3[ip] = MIN(MAX(evi3[ip], minval), maxval);
         }
     }
 }
  
 void get_ndvi_evi(l1str *l1rec, int prodnum, float prod[]) {
  
     static int firstCall = 1;
     if (firstCall) {
         ibblue = bindex_get(412);
         ibred = bindex_get(670);
         ibnir = bindex_get(865);
  
         if (l1rec->l1file->sensorID == MODISA ||
                 l1rec->l1file->sensorID == MODIST) {
             ibred = bindex_get(645);
             ibnir = bindex_get(859);
         }
     }
  
     switch (prodnum) {
     case CAT_ndvi:
         get_ndvi(l1rec, prod);
         break;
     case CAT_evi:
         get_evi(l1rec, prod);
         break;
     case CAT_evi2:
         get_evi2(l1rec, prod);
         break;
     case CAT_evi3:
         get_evi3(l1rec, prod);
         break;
     default:
         printf("Error: %s : Unknown product specifier: %d\n", __FILE__, prodnum);
         exit(FATAL_ERROR);
         break;
     }
  
 }

MAX

#define MAX(A, B)

Definition: swl0_utils.h:26

MIN

#define MIN(x, y)

Definition: rice.h:169

red

void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf, int ic1, int jc1, int jcf, int kc, float ***c, float **s)

#define L(lambda, T)

Definition: PreprocessP.h:185

CAT_evi2

#define CAT_evi2

Definition: l2prod.h:73

l1rec

read l1rec

Definition: HOWTO_Add_a_sensor.txt:72

get_ndvi

void get_ndvi(l1str *l1rec, float ndvi[])

Definition: get_ndvi.c:31

MODIST

#define MODIST

Definition: sensorDefs.h:18

PRODFAIL

#define PRODFAIL

Definition: l2_flags.h:41

get_evi

void get_evi(l1str *l1rec, float evi[])

Definition: get_ndvi.c:63

l12_proto.h

bindex_get

int bindex_get(int32_t wave)

Definition: windex.c:45

get_ndvi_evi

void get_ndvi_evi(l1str *l1rec, int prodnum, float prod[])

Definition: get_ndvi.c:203

CAT_evi

#define CAT_evi

Definition: l2prod.h:45

CAT_evi3

#define CAT_evi3

Definition: l2prod.h:74

FATAL_ERROR

#define FATAL_ERROR

Definition: swl0_parms.h:5

mask

a context in which it is NOT documented to do so subscript which cannot be easily calculated when extracting TONS attitude data from the Terra L0 files Corrected several defects in extraction of entrained ephemeris and and as HDF file for both the L1A and Geolocation enabling retrieval of South Polar DEM data Resolved Bug by changing to opent the geolocation file only after a successful read of the L1A and also by checking for fatal errors from not restoring C5 and to report how many of those high resolution values were water in the new WaterPresent SDS Added valid_range attribute to Land SeaMask Changed to bilinearly interpolate the geoid_height to remove artifacts at one degree lines Made corrections to const qualification of pointers allowed by new version of M API library Removed casts that are no longer for same not the geoid Corrected off by one error in calculation of high resolution offsets Corrected parsing of maneuver list configuration parameter Corrected to set Height SDS to fill values when geolocation when for elevation and land water mask

Definition: HISTORY.txt:114

LAND

#define LAND

Definition: l2_flags.h:12

get_evi3

void get_evi3(l1str *l1rec, float evi3[])

Definition: get_ndvi.c:156

CAT_ndvi

#define CAT_ndvi

Definition: l2prod.h:39

val

msiBandIdx val

Definition: l1c_msi.cpp:34

MODISA

#define MODISA

Definition: sensorDefs.h:19

get_evi2

void get_evi2(l1str *l1rec, float evi2[])

Definition: get_ndvi.c:122