NASA Ocean Color

ocssw V2022

 /*
  * get_avw.c
  *
  * This algorithm returns the weighted harmonic mean of the visible-range (400 – 700 nm) remote sensing reflectance (Rrs) wavelengths,
  * outputting the Apparent Visible Wavelength (AVW) in units of nanometers.
  *  The AVW is an optical water classification index, representing a one-dimensional geophysical metric that is inherently correlated
  *  to Rrs spectral shape (Vandermeulen et al. 2020).
  *
  *
  *Vandermeulen, R. A., Mannino, A., Craig, S.E., Werdell, P.J., 2020: “150 shades of green: Using the full spectrum of remote sensing reflectance
  *   to elucidate color shifts in the ocean,” Remote Sensing of Environment, 247, 111900, https://doi.org/10.1016/j.rse.2020.111900
  *
  *  Created on: Sep 30, 2020
  *      Author: M. Zhang
  */
  
 #include "l12_proto.h"
  
  
 float avw_cal_hypspectral(float *Rrs, float *wave, int nwave){
  
     float avw=BAD_FLT;
     int32_t ib;
  
     float dfirst=0,dlast=0;
     float *Secondderiv=(float *)malloc(nwave*sizeof(float));
  
     int32_t nwave_out=700-400+1;
     float *wave_out=(float *)malloc(nwave_out*sizeof(float));
     float *Rrs_out =(float *)malloc(nwave_out*sizeof(float));
     double nominator=0., denominator=0.;
  
     for(ib=0;ib<nwave_out;ib++)
         wave_out[ib]=400+ib;
  
     dfirst=first_deriv(wave,Rrs,0);
     dlast=first_deriv(wave,Rrs,nwave);
  
     spline(wave,Rrs,nwave,dfirst,dlast,Secondderiv);
  
     for(ib=0;ib<nwave_out;ib++)
         splint(wave,Rrs,Secondderiv,nwave,wave_out[ib],&Rrs_out[ib]);
  
     for(ib=0;ib<nwave_out;ib++){
         nominator+=Rrs_out[ib];
         denominator+=(Rrs_out[ib]/wave_out[ib]);
     }
     avw=nominator/denominator;
  
     free(wave_out);free(Rrs_out);
     free(Secondderiv);
     return avw;
  
 }
  
  
 float avw_cal_multispectral(float *Rrs, float *wave, int nwave){
  
     float avw=BAD_FLT;
     int32_t ib;
     float *avw_coef=input->avw_coef;
  
     double nominator=0., denominator=0.,temp=0.;
  
     for(ib=0;ib<nwave; ib++){
         nominator+=Rrs[ib];
         denominator+=(Rrs[ib]/wave[ib]);
     }
     avw=nominator/denominator;
  
     for(ib=0;ib<6;ib++)
         temp+=avw_coef[ib]*pow(avw,ib);
  
     avw=temp;
  
     return avw;
  
 }
  
 void get_avw(l2str *l2rec, float avw[]){
     filehandle* l1file = l2rec->l1rec->l1file;
  
     int32_t ip,ib;
     int32_t ipb;
     
     int32_t sensorID = l1file->sensorID;
     int32_t nbands = l1file->nbands;
     int32_t npix = l2rec->l1rec->npix;
     static float fsol;
     float *wave = l1file->fwave;
     float *Rrs;
     static float *Rrs_avw, *wave_avw;
     static int firstcall=1;
     int32_t nwave_avw;
     int32_t negative=0;
     static int ifhyper=0;
     static  int ib400=0,ib700=0;
  
     if(firstcall){
         firstcall=0;
         Rrs_avw =(float *)malloc(nbands*sizeof(float));
         wave_avw=(float *)malloc(nbands*sizeof(float));
         switch (sensorID){
         case HICO:
         case PRISM:
         case OCI:
         case OCIS:
             ifhyper=1;
             break;
         }
         if(ifhyper){
             ib400=windex(400, wave,nbands);
             ib700=windex(700, wave,nbands);
  
             if(ib400>0 && wave[ib400]>=400)
                 ib400--;
             if(ib700<nbands-1 && wave[ib700]<700)
                 ib700++;
         }
         fsol=l2rec->l1rec->fsol;
     }
  
     /*                                                      */
     /* Compute desired products at each pixel               */
     /*                                                      */
     for (ip = 0; ip < npix; ip++) {
         avw[ip]=BAD_FLT;
  
         if(l2rec->l1rec->mask[ip] || l2rec->l1rec->solz[ip] >= SOLZNIGHT)
             continue;
  
         ipb=ip*nbands;
         Rrs=&l2rec->Rrs[ipb];
  
         nwave_avw=0;
         negative=0;
         if(ifhyper){
             for(ib=ib400;ib<=ib700;ib++){
                 if(Rrs[ib]!=BAD_FLT){
                     Rrs_avw [nwave_avw]=Rrs[ib] * (l1file->Fonom[ib]*fsol/l2rec->l1rec->Fo[ib])/l2rec->outband_correction[ipb+ib];
                     wave_avw[nwave_avw]=wave[ib];
                     nwave_avw++;
                 }
                 else if(Rrs[ib]<0)
                     negative=1;
             }
         }
         else{
             for(ib=0;ib<nbands;ib++){
                 if(wave[ib]>=400 && wave[ib]<=700){
                     if(Rrs[ib]!=BAD_FLT){
                         Rrs_avw [nwave_avw]=Rrs[ib] * (l1file->Fonom[ib]*fsol/l2rec->l1rec->Fo[ib])/l2rec->outband_correction[ipb+ib];
                         wave_avw[nwave_avw]=wave[ib];
                         nwave_avw++;
                     }
                     else if(Rrs[ib]<0)
                         negative=1;
                 }
             }
         }
  
         if(negative)
             l2rec->l1rec->flags[ip] |=PRODWARN;
  
         if(ifhyper)
             avw[ip]=avw_cal_hypspectral(Rrs_avw,wave_avw,nwave_avw);
         else
             avw[ip]=avw_cal_multispectral(Rrs_avw,wave_avw,nwave_avw);
  
     }
 }
 void get_Rrs_brightness(l2str *l2rec, float Rrs_brightness[]){
  
     int32_t ip,ib;
     int32_t ipb;
  
     int32_t nbands = l2rec->l1rec->l1file->nbands;
     int32_t npix = l2rec->l1rec->npix;
     float *wave=l2rec->l1rec->l1file->fwave;
     float *Rrs;
     static float *Rrs_avw, *wave_avw;
     static int firstcall=1;
     int32_t nwave_avw;
     static  int ib400=0,ib700=0;
  
     if(firstcall){
         firstcall=0;
         Rrs_avw =(float *)malloc(nbands*sizeof(float));
         wave_avw=(float *)malloc(nbands*sizeof(float));
  
         ib400=windex(400, wave,nbands);
         ib700=windex(700, wave,nbands);
     }
  
     /*                                                      */
     /* Compute desired products at each pixel               */
     /*                                                      */
     for (ip = 0; ip < npix; ip++) {
         Rrs_brightness[ip]=BAD_FLT;
  
         if(l2rec->l1rec->mask[ip] || l2rec->l1rec->solz[ip] >= SOLZNIGHT)
             continue;
  
         ipb=ip*nbands;
         Rrs=&l2rec->Rrs[ipb];
         nwave_avw=0;
  
         for(ib=ib400;ib<=ib700;ib++){
                 if(Rrs[ib]!=BAD_FLT){
                     Rrs_avw [nwave_avw]=Rrs[ib];
                     wave_avw[nwave_avw]=wave[ib];
                     nwave_avw++;
                 }
                 else if (Rrs[ib]<0)
                     l2rec->l1rec->flags[ip]|=PRODWARN;
         }
  
         if(nwave_avw<2)
             continue;
  
         Rrs_brightness[ip]=0;
         for(ib=0;ib<nwave_avw-1;ib++){
             Rrs_brightness[ip]+=(Rrs_avw[ib]+Rrs_avw[ib+1])/2.0*(wave_avw[ib+1]-wave_avw[ib]);
         }
     }
 }
  
 void get_lambda_max(l2str *l2rec, float lambda_max[]){
  
     int32_t ip,ib;
     int32_t ipb;
  
     int32_t nbands = l2rec->l1rec->l1file->nbands;
     int32_t npix = l2rec->l1rec->npix;
     float *wave=l2rec->l1rec->l1file->fwave;
     float *Rrs;
     float Rrs_temp;
     static int firstcall=1;
     static int ib400, ib700;
     int negative=0;
  
     if(firstcall){
         ib400=windex(400, wave,nbands);
         ib700=windex(700, wave,nbands);
         firstcall=0;
     }
  
     /*                                                      */
     /* Compute desired products at each pixel               */
     /*                                                      */
     for (ip = 0; ip < npix; ip++) {
         negative=0;
         lambda_max[ip]=BAD_FLT;
  
         if(l2rec->l1rec->mask[ip] || l2rec->l1rec->solz[ip] >= SOLZNIGHT)
             continue;
  
         ipb=ip*nbands;
         Rrs=&l2rec->Rrs[ipb];
         Rrs_temp=Rrs[ib400];
  
         if(Rrs_temp>=0)
             lambda_max[ip]=wave[ib400];
         else
             negative=1;
  
         for(ib=ib400+1;ib<=ib700;ib++){
             if(Rrs[ib]<0)
                 negative=1;
             if(Rrs[ib]>Rrs_temp){
                 Rrs_temp=Rrs[ib];
                 lambda_max[ip]=wave[ib];
             }
         }
         if(negative)
             l2rec->l1rec->flags[ip] |=PRODWARN;
     }
 }