DtAlgorithm.h

Go to the documentation of this file.

/*
 * DtAlgorithm.h
 *
 *  Created on: Nov 3, 2016
 *      Author: ssander2
 */
 
#ifndef INCLUDE_DTALGORITHM_H_
#define INCLUDE_DTALGORITHM_H_
 
#include <DDAlgorithm.h>
#include "DtLutNetcdf.h"
 
using namespace std;
 
enum DTBANDS_ENUM {
    D488,           // m03
    D550,           // m04
    D670,           // m05
    D865,           // m07
    D1240,          // m08
    D1610,          // m10
    D2250,          // m11
    D412,           // m01
    D1380,          // m09
    D11000,         // m15
    MAXDBAND
};
 
enum DTLBANDS_ENUM {
    DL488,           // m03
    DL550,           // m04
    DL670,           // m05
    DL2250,          // m11
    MAXDLBAND
};
 
enum AEROSOL_TYPE {
    CONTN=1,
    BACKG,
    SMOKE,
    URBAN,
    COARSE,
    MAX_AERO
};
 
 
class DtAlgorithm : public DDAlgorithm
{
    friend class DtCloudMaskLand;
    friend class DtCloudMaskOcean;
    friend class DtSedimentMask;
public:
 
// static definitions
    static constexpr int            LINE = 1;
    static constexpr int            GRIDX = LINE;
    static constexpr int            GRIDY = LINE;
    static constexpr int            NUM_CASES_SMALL = 4;
    static constexpr int            NUM_CASES_BIG = 5;
    static constexpr int            NUM_LUTS = 4;
    static constexpr int            NUM_LATS = 180;
    static constexpr int            NUM_LONS = 360;
 
// angle limits
    static constexpr float          MINMTHET0=0.0;
    static constexpr float          MAXMTHET0=80.0;
    static constexpr float          MINMTHET=0.0;
    static constexpr float          MAXMTHET=72.0;
    static constexpr float          MINMPHI=0.0;
    static constexpr float          MAXMPHI=180.0;
 
    static constexpr float          MAXTAU=5.0;
    static constexpr float          MINTAU=-0.01;
 
    static constexpr int            NUM_RFL_BANDS = 10;
    static constexpr int            DT_OWAVES = 7;
    static constexpr int            DT_LWAVES = 3;
    static constexpr int            MODIS_WAVES = 7;
    static constexpr int            VIIRS_WAVES = 3;
    static const float              wind_[WIND_LUT_ENTRIES];
    static const float              pressure_[P_LEVELS];
 
    enum SIZES {
        ISMALL,
        IBIG
    };
 
    enum STATS {
        IBEST,
        IAVG
    };
 
    dtGasCorrectionLUT gc_lut_;
    unsigned char cmask_;
 
    float rfld_[NUM_RFL_BANDS];
 
    DtAlgorithm();
 
    virtual ~DtAlgorithm ();
 
    virtual int initialize( map<string, ddata*> imap );
 
    virtual vector<string>  get_products() { return {};};
 
    virtual map<string, ddata*> process(vector<size_t> start, vector<size_t> count,
            map<string, ddata*> imap);
 
    int compute_gas_correction();
 
    int fit_line( float x[], float y[], float sig[], int ndata,
                  float& A, float& B );
    float   cloud_fraction_;
    float   scatter_angle_;
    float   glint_angle_;
    float   glint_refl_;
    float   ndvi_;
 
protected:
 
    // LUT angle-based indices
    int   SZA_0_;
    int   SZA_1_;
    int   THE_0_;
    int   THE_1_;
    int   PHI_0_;
    int   PHI_1_;
 
    int interp_extrap( int num, float xin,
            float x[], float y[], float& yout );
    int mean_std( int n, float* data, float& mean, float& sdev );
 
    int mean_std_weighted( int n, float* data, float& mean,
            float& sdev, float* weight );
 
    int sort_index( int numPts, float array[], int index[] );
 
    int sort_inplace( int numPts, float array1[], float array2[] );
 
    int set_byte( short val, short bitPos, short& target);
 
    int compute_glint_angle();
 
    int compute_scatter_angle(float& scat_angle);
 
};
 
class DtCloudMaskOcean;
 
class DtAlgOcean : public DtAlgorithm
{
public:
 
// algorithm thresholds
    static constexpr float THRHLD550 = 0.0025;
    static constexpr float THRHLD470 = 0.4;
    static constexpr float THRHLD138_1 = 0.03;
    static constexpr float THRHLD138_2 = 0.005;
    static constexpr float THRHLD_DUST = 0.75;
    static constexpr float THRHLD_CIRRUS = 0.03;
    static constexpr float Threshold_LSQ_Error = 3.7;
 
    static constexpr float  GLINT_ANGLE_THRESHOLD = 40.0;
//    static constexpr float  GLINT_ANGLE_THRESHOLD = 30.0;
    static constexpr float  LOW_THRESHOLD = -0.05;
    static constexpr float  HIGH_THRESHOLD = 5.0;
    static constexpr float  GLINT_REFL_THRESHOLD = 0.01570796326; //0.005*PI;
 
    static constexpr int NUMCASES = 4;
    static constexpr int NUMCASEB = 5;
    static constexpr int NWAV = 7;
    static constexpr int NAOT = 6;
    static constexpr int NTH0 = 11;
    static constexpr int NTHET = 16;
    static constexpr int NPHI = 16;
    static constexpr int NUM_ARRAY_ELEMENTS = 77;
    static constexpr int NUM_SIZES = 2;
    static constexpr int NUM_STATS = 2;
 
// data from LUTs
    dtOceanAerosolLUT lut_;
    DtCloudMaskOcean* cm_;
 
// global working data
//  float scatter_angle_[NUMCASES*NUMCASEB];
 
// interpolate_rayleigh values based on reported angles
    float refl_rayl_[NWAV];
 
// interpolate_angle
    float refl_big_[NWAV][NUM_CASES_BIG][NAOT];
    float refl_small_[NWAV][NUM_CASES_SMALL][NAOT];
 
// interpolate_sza_albedo, interpolate_tau_albedo
    float albedo_R_small_tau_[NWAV][NUM_CASES_SMALL][NAOT];
    float albedo_R_big_tau_[NWAV][NUM_CASES_BIG][NAOT];
    float albedo_T_small_tau_[NWAV][NUM_CASES_SMALL][NAOT];
    float albedo_T_big_tau_[NWAV][NUM_CASES_BIG][NAOT];
 
// compute_average
    float refl_[NWAV];
    float sdev_[NWAV];
    int   numData_[NWAV];
    int   good_pixels_[NWAV];
 
// results by case
    float tau_[NWAV][NUMCASES*NUMCASEB];
    float tau_small_[NWAV][NUMCASES*NUMCASEB];
    float tau_big_[NWAV][NUMCASES*NUMCASEB];
    float backscatter_[NWAV][NUMCASES*NUMCASEB];
    float assym_[NWAV][NUMCASES*NUMCASEB];
    float refl_flux_[NWAV][NUMCASES*NUMCASEB];
    float trans_flux_[NWAV][NUMCASES*NUMCASEB];
    float angstrom_exp_[NUM_STATS][NUMCASES*NUMCASEB];
    float solution_index_[NUM_SIZES][NUMCASES*NUMCASEB];
    float ccn_[NUMCASES*NUMCASEB];
    float eff_radius_[NUMCASES*NUMCASEB];
    float eff_variance_[NUMCASES*NUMCASEB];
    float mass_con_ocean_[NUMCASES*NUMCASEB];
    float xmin_[2][NUMCASES*NUMCASEB];
    float funmin_[2][NUMCASES*NUMCASEB];
    float tau_X55_[2][NUMCASES*NUMCASEB];
    float error_[2][NWAV][NUMCASES*NUMCASEB];
 
// results best and averaged
    float refl_avg_[NWAV];
    float sdev_avg_[NWAV];
    float tau_avg_[NUM_STATS][NWAV];
    float tau_small_avg_[NUM_STATS][NWAV];
    float tau_big_avg_[NUM_STATS][NWAV];
    float backscatter_avg_[NUM_STATS][NWAV];
    float assym_avg_[NUM_STATS][NWAV];
    float refl_flux_avg_[NUM_STATS][NWAV];
    float trans_flux_avg_[NUM_STATS][NWAV];
    float angstrom_exp_avg_[NUM_STATS][NUM_SIZES];
    float solution_index_avg_[NUM_STATS][NUM_SIZES];
    float ccn_avg_[NUM_STATS];
    float eff_radius_avg_[NUM_STATS];
    float eff_variance_avg_[NUM_STATS];
    float mass_con_ocean_avg_[NUM_STATS];
    float xmin_avg_[NUM_STATS];
    float funmin_avg_[NUM_STATS];
    float tau_X55_avg_[NUM_STATS];
 
// quality control
    short qcontrol_;
    short qcontrol_exclude_;
    short qcontrol_special_;
    short qcontrol_cirrus_;
    short quality_dust_flag_glint_;
    short quality_dust_flag_off_glint_;
    short quality_to_pass_[2];
    short quality_flag_[12];
 
// Intermediate output data per cell
    float sds_refl_[NWAV];
    float sds_refl_sdev_[NWAV];
    short sds_numPixels_[NWAV];
    float sds_tau_best_[NWAV];
    float sds_tau_avg_[NWAV];
    float sds_tau_big_best_[NWAV];
    float sds_tau_big_avg_[NWAV];
    float sds_tau_small_best_[NWAV];
    float sds_tau_small_avg_[NWAV];
    float sds_assy_best_[NWAV];
    float sds_assy_avg_[NWAV];
    float sds_back_best_[NWAV];
    float sds_back_avg_[NWAV];
    float sds_reff_best_[NWAV];
    float sds_reff_avg_[NWAV];
    float sds_tranf_best_[NWAV];
    float sds_tranf_avg_[NWAV];
    float sds_Tau_Land_Ocean_;
    float sds_Tau_Land_Ocean_img_;
    float sds_Small_Weighting_[NUM_SIZES];
    float sds_Least_Error_[NUM_SIZES];
    float sds_tau_X55_[NUM_SIZES];
    float sds_EffRad_[NUM_SIZES];
    float sds_EffVar_[NUM_SIZES];
    short sds_Sol_Index_Small_[NUM_SIZES];
    short sds_Sol_Index_Large_[NUM_SIZES];
    float sds_Angs_Coeff1_[NUM_SIZES];
    float sds_Angs_Coeff2_[NUM_SIZES];
    float sds_Mass_Conc_[NUM_SIZES];
    float sds_CCN_[NUM_SIZES];
    float sds_AOT_model_[NUMCASES+NUMCASEB];
    short sds_land_ocean_quality_;
 
    DtAlgOcean();
 
    ~DtAlgOcean ();
 
    int initialize( map<string, ddata*> imap );
 
    vector<string>  get_products() { return {"cloud_mask", "quality_flag", "aerosol_type",
        "l2_flags", "scattering_angle", "fmf_550", "angstrom", "aot_380",
        "aot_490", "aot_550", "aot_670","aot_865", "aot_1240", "aot_1610", "aot_2250"};};
 
    map<string, ddata*> process(vector<size_t> start, vector<size_t> count,
            map<string, ddata*> imap);
 
protected:
 
    int index_geometry( float sza, float azim, float phi);
 
    int interpolate_rayleigh();
 
    int interpolate_angles();
 
    int compute_tau_flux( int iBig, int iSmall, int iSol );
 
    int compute_avg_refl();
 
    int compute_average_to_500m();
 
    int store_reflectance();
 
    int store_output();
 
    int assign_quality();
 
    int set_fill_out();
 
    int average_output();
 
    int compute_minimum(int iBig, int iSmall, int iSol);
 
    int compute_minimum_baseline(int iBig, int iSmall, int iSol);
 
    float fun_tau( float xmin, int iBig, int iSmall, int iSol );
};
 
class DtCloudMaskLand;
 
class DtAlgLand : public DtAlgorithm
{
public:
 
// algorithm constants
    static constexpr int NLTHET0 = 11;
    static constexpr int NLTHE = 15;
    static constexpr int NLPHI = 16;
    static constexpr int NLWAV = 4;
    static constexpr int NLTAU = 7;
    static constexpr int NLTABLE = 5;
    static constexpr int NLETA = 13;
    static constexpr int NLSIZE = 2;
    static constexpr int DTABLE = 5;
 
    static constexpr float THR213MIN_1 = 0.01;
    static constexpr float THR213MAX_1 = 0.25;
    static constexpr float THR213MIN_2 = 0.25;
    static constexpr float THR213MAX_2 = 0.25;
    static constexpr float PRESSURE_P0 = 1013.0;
    static constexpr float DLAT = 0.5;
    static constexpr float DLON = 0.5;
 
    static constexpr int QA_LAND = 5;
 
// data from LUTs
    dtLandAerosolLUT lut_;
    DtCloudMaskLand* cm_;
 
    int   season_;
    float scatter_angle_;
    int   mtable_[NLSIZE];
 
// interpolate_
    float refl_ray_nl_[NLWAV];
    float opth_nl_[NLSIZE][NLWAV][NLTAU];
    float int_nl_[NLSIZE][NLWAV][NLTAU];
    float fd_nl_[NLSIZE][NLWAV][NLTAU];
    float t_nl_[NLSIZE][NLWAV][NLTAU];
    float fdt_nl_[NLSIZE][NLWAV][NLTAU];
    float sbar_nl_[NLSIZE][NLWAV][NLTAU];
// average
    float refl_inter_[NLWAV][GRIDX*GRIDY];
    float refl_[NLWAV]; // mean reflectance
    float sdev_[NLWAV]; // reflectance standard deviation
    short error_;
// retrieve_first, retrieve_second
    float rho_star_[NLSIZE][NLWAV][NLTAU];
    float rho_star_tot_[NLWAV][NLTAU];
    float rho_S212_[NLSIZE][NLTAU];
    float rho_S212_tot_[NLTAU];
    float errwave_[NLWAV];
    float rho_S466_;
    float yint_466_;
    float slope_466_;
    float rho_S644_;
    float yint_644_;
    float slope_644_;
    short eta_flag_;
// results
    float aot_d_[NLWAV];  // tau corrected
    float aot_f_[NLWAV];  // tau fine (small)
    float aot_c_[NLWAV];  // tau coarse (big)
    float rho_sfc_[NLWAV];  // surface reflectance
    float eta_;  // dust weighting
    float err644_;  // fitting error
    float masscon_;  // mass concentration
    float angstrom_;  // angs coefficient
    float ndvi_;  // ndvi
    short iaer_;  // aerosol type
    int   good_pixels_[NLWAV]; // good pixels
    int   num_pixels_used_;  // = good_pixels_[W659]
 
// program and quality control
    short iproc_;
    short ifinish_;
    float thresh_min_;
    float thresh_max_;
    short return_quality_cirrus_;
    short success_ret_;
    short fail_ret_;
    short quality_flag_for_joint_[2];
    short quality_flag_for_retr_;
    short qcontrol_special_;
 
// Intermediate output data per cell
    short sds_qcontrol_[QA_LAND];
    short sds_aerosol_type_;
    float sds_scat_angle_;
    float sds_fitting_error_;
    float sds_mass_conc_;
    float sds_cloud_fraction_;
    float sds_angs_coeff_;
    float sds_dust_weighting_;
    float sds_ndvi_;
    float sds_numpixels_[NLWAV];
    float sds_tau_corrected_[NLWAV];
    float sds_refl_[NLWAV];
    float sds_refl_std_[NLWAV];
    float sds_tau_small_[NLWAV];
    float sds_tau_big_[NLWAV];
    float sds_surface_reflectance_[NLWAV];
 
    DtAlgLand();
 
    ~DtAlgLand ();
 
    int initialize( map<string, ddata*> imap );
 
    vector<string>  get_products() { return {"cloud_mask", "quality_flag", "aerosol_type",
        "l2_flags", "scattering_angle", "fmf_550", "angstrom", "aot_380",
        "aot_490", "aot_550", "aot_670", "aot_2250"};};
 
    map<string, ddata*> process(vector<size_t> start, vector<size_t> count,
            map<string, ddata*> imap);
 
protected:
 
    int index_geometry( float sza, float azim, float phi);
 
//  int select_lut( int iSize, int iTable );
 
    int interpolate_rayleigh();
 
    int interpolate_angle();
 
    int interpolate_elevation();
 
    int compute_average(size_t iy, size_t ix);
 
    int compute_cloudmask_ndvi(size_t iy, size_t ix,
                                int& iCldRed, int& iCldBlue);
 
    int simulate_toa();
 
    int retrieve_first();
 
    int retrieve_second();
 
    int assign_quality();
 
    int store_output();
 
    int fill_values();
 
};
 
#endif