DDProcess.cpp

Go to the documentation of this file.

/*******************************************************************************
 *
 * NAME: DDProcess.cpp
 *
 * DESCRIPTION: Process granule using selected algorithm.
 *
 *  Created on: Aug 25, 2020
 *      Author: Sam Anderson
 *
 *******************************************************************************/
 
#include <pugixml.hpp>
 
#include <DDataset.hpp>
#include <DDProcess.h>
#include <DDSensor.h>
#include <DDOptions.h>
#include <DDAncillary.h>
#include "darktarget/DtAlgorithm.h"
#include "deepblue/DbAlgorithm.h"
 
using namespace std;
using namespace pugi;
 
const vector<string> DDProcess::out_groups = { "geolocation", "ancillary",
        "observations", "geophysical_data", "statistics" };
 
const map<string,dtype> DDProcess::input_names = {
        {"cloud_mask", dtype::UBYTE},
        {"land_water", dtype::UBYTE},
        {"relative_azimuth", dtype::FLOAT},
        {"solar_zenith", dtype::FLOAT},
        {"sensor_zenith", dtype::FLOAT},
        {"scattering_angle", dtype::FLOAT},
        {"elevation", dtype::FLOAT},
        {"wind_speed", dtype::FLOAT},
        {"wind_direction", dtype::FLOAT},
        {"water_vapor", dtype::FLOAT},
        {"ozone", dtype::FLOAT},
        {"surface_pressure", dtype::FLOAT},
        {"chlorophyll", dtype::FLOAT},
};
 
const map<string,dtype> DDProcess::dtdb_names = {
        {"quality_flag", dtype::SHORT},
        {"aerosol_type", dtype::SHORT},
        {"l2_flags", dtype::UINT},
        {"fmf_550", dtype::FLOAT},
        {"angstrom", dtype::FLOAT},
        {"ndvi", dtype::FLOAT},
        {"residual_error", dtype::FLOAT},
};
 
const map<string,rhot_band> DDProcess::rhot_band_names = {
        {"rhot_410", rhot_band::W410},           // m01
        {"rhot_445", rhot_band::W445},           // m02
        {"rhot_490", rhot_band::W490},           // m03
        {"rhot_550", rhot_band::W550},           // m04
        {"rhot_670", rhot_band::W670},           // m05
        {"rhot_865", rhot_band::W865},           // m07
        {"rhot_1240", rhot_band::W1240},         // m08
        {"rhot_1380", rhot_band::W1380},         // m09
        {"rhot_1610", rhot_band::W1610},         // m10
        {"rhot_2250", rhot_band::W2250},         // m11
};
 
const map<string,aot_band> DDProcess::aot_band_names = {
        {"aot_410", aot_band::W410},           // m01
        {"aot_490", aot_band::W490},           // m03
        {"aot_550", aot_band::W550},           // m04
        {"aot_670", aot_band::W670},           // m05
        {"aot_865", aot_band::W865},           // m07
        {"aot_1240", aot_band::W1240},         // m08
        {"aot_1610", aot_band::W1610},         // m10
        {"aot_2250", aot_band::W2250},         // m11
};
 
const map<string,srf_band> DDProcess::srf_band_names = {
        {"surface_410", srf_band::W410},           // m01
        {"surface_490", srf_band::W490},           // m03
        {"surface_670", srf_band::W670},           // m05
        {"surface_2250", srf_band::W2250},         // m11
};
 
/**************************************************************************
 * NAME: DDProcess()
 *
 * DESCRIPTION: Class Constructor
 *
 *************************************************************************/
DDProcess::DDProcess()
{
    lines_ = 0.0;
    pixels_ = 0.0;
    alg_ = algorithm::NOTHING;
    instrument_ = SENSOR::NOTHING;
    psensor_ = nullptr;
    pancillary_ = nullptr;
    pl_ = nullptr;
    po_ = nullptr;
    bday_ = true;
    lprw_ = 1;
    bglintmask_=true;
    bcloudmask_=true;
    bgascorrect_=true;
}
 
/**************************************************************************
 * NAME: ~DDProcess()
 *
 * DESCRIPTION: Class Destructor
 *
 *************************************************************************/
DDProcess::~DDProcess()
{
    if (psensor_ != 0) {
        delete psensor_;
        psensor_ = 0L;
    }
    if (pancillary_ != 0) {
        delete pancillary_;
        pancillary_ = 0L;
    }
}
 
/**************************************************************************
 * NAME: initialize()
 *
 * DESCRIPTION: Initializes data and object classes for granule
 *
 *************************************************************************/
int DDProcess::initialize()
{
    int status = DTDB_SUCCESS;
 
    string salg = get_option(ALGORITHM);
    if (salg == "darktarget") {
        std::cerr << "DDProcess:: Algorithm Dark Target" << std::endl;
        alg_ = algorithm::DARKTARGET;
    } else if (salg == "deepblue") {
        std::cerr << "DDProcess:: Algorithm Deep Blue" << std::endl;
        alg_ = algorithm::DEEPBLUE;
    } else {
        std::cerr << "DDProcess:: Invalid algorithm, exiting ..." << std::endl;
        return DTDB_FAIL;
    }
    status = query_granule();
    if (status != DTDB_SUCCESS) {
        cerr << "DDProcess:: Query granule failure" << endl;
        return DTDB_FAIL;
    } else if (!bday_) {
        cerr << "DDProcess:: Night granule detected" << endl;
        return DTDB_FAIL;
    }
    switch (instrument_) {
    case SENSOR::POCI:
        psensor_ = new POCI();
        break;
    case SENSOR::VIIRS:
        psensor_ = new VIIRS();
        break;
    case SENSOR::NOTHING:
    default:
        cerr << "DDProcess:: Invalid sensor" << endl;
        return DTDB_FAIL;
        break;
    }
    pancillary_ = new DDAncillary();
    if (alg_== algorithm::DARKTARGET) {
        pl_ = new DtAlgLand();
        po_ = new DtAlgOcean();
    } else if (alg_== algorithm::DEEPBLUE) {
        pl_ = new DbAlgLand();
        po_ = new DbAlgOcean();
    }
    lprw_ = (size_t) get_option_int(INPUT_LINES_PER_RW);
    lprw_ = (lprw_<3) ? 3 : lprw_;
 
    vector<string> products = pl_->get_products();
    out_products.insert(out_products.end(), products.begin(), products.end());
    products = po_->get_products();
    out_products.insert(out_products.end(), products.begin(), products.end());
 
    return status;
}
 
/**************************************************************************
 * NAME: query_granule()
 *
 * DESCRIPTION: Read granule data essential for initialization
 *
 *************************************************************************/
int DDProcess::query_granule()
{
    int status = DTDB_SUCCESS;
    string filepath = get_option(INPUT_L1B);
 
    if (!filepath.empty()) {
        NcFile* nc_input;
        try {
            nc_input = new NcFile(filepath, NcFile::read );
        }
        catch( NcException& e) {
            e.what();
            cout << "DDProcess:: Failure opening L1B file" << endl;
            return DTDB_FAIL;
        }
        string str;
        multimap <string, NcGroupAtt> attributes = nc_input->getAtts(NcGroup::Current);
        multimap <string, NcGroupAtt>::iterator it;
        if ((it = attributes.find("instrument")) != attributes.end()) {
            (it->second).getValues(str);
        }
        if (str.find("OCI") != string::npos) {
            instrument_ =  SENSOR::POCI;
        } else if (str.find("VIIRS") != string::npos) {
            instrument_ = SENSOR::VIIRS;
        } else {
            cout << "DDProcess:: Instrument not supported" << endl;
            return DTDB_FAIL;
        }
        switch (instrument_) {
        case SENSOR::VIIRS:
            if ((it = attributes.find("DayNightFlag")) != attributes.end()) {
                (it->second).getValues(str);
            } else if((it = attributes.find("day_night_flag")) != attributes.end()) {
                    (it->second).getValues(str);
            }
            if (str != "Night") {
                bday_ = true;
            } else {
                bday_ = false;
            }
            lines_ = nc_input->getDim("number_of_lines").getSize();
            pixels_ = nc_input->getDim("number_of_pixels").getSize();
            title_ = "VIIRS Level-2 Data";
            break;
        case SENSOR::POCI:
            bday_ = true;
            lines_ = nc_input->getDim("number_of_scans").getSize();
            pixels_ = nc_input->getDim("ccd_pixels").getSize();
            title_ = "OCI Level-2 Data";
            break;
        default:
            cout << "DDProcess:: Failure reading dimensions" << endl;
            return DTDB_FAIL;
        }
        delete nc_input;
    }  else {
        cout << "Invalid L1B filename ..." << endl;
        status = DTDB_FAIL;
    }
    return status;
}
 
/**************************************************************************
 * NAME: process()
 *
 * DESCRIPTION: Processes data and produces output values for granule
 *
 *************************************************************************/
int DDProcess::process()
{
    int status = DTDB_SUCCESS;
 
    map<string,ddata*> imap = psensor_->create( {0, 0}, {lines_, pixels_} );
    if (static_cast<ddval<int>*>(imap["status"])->val != DTDB_SUCCESS) {
        cerr << "DDProcess:: Create sensor failure at line " << endl;
        return DTDB_FAIL;
    }
 
    ddstr* dstr = new ddstr(title_);
    imap.insert({"title", dstr});
    ddma<unsigned char,2>* plw = static_cast<ddma<unsigned char,2>*>(imap["land_water"]);
    string parfilepath = get_option(INPUT_PAR);
    dstr = new ddstr(parfilepath);
    imap.insert({"config_file", dstr});
    ddval<unsigned char>* bval = new ddval<unsigned char>(dtype::UBYTE, get_bool(BOOL_GAS_CORRECTION));
    imap.insert({"bgascorrect", bval});
    bgascorrect_ = bval->val;
    bval = new ddval<unsigned char>(dtype::UBYTE, get_bool(BOOL_MASK_GLINT));
    imap.insert({"bglintmask", bval});
    bglintmask_ = bval->val;
    bval = new ddval<unsigned char>(dtype::UBYTE, get_bool(BOOL_MASK_CLOUD));
    imap.insert({"bcloudmask", bval});
    bcloudmask_ = bval->val;
 
    status = pancillary_->initialize( (static_cast<ddval<int>*>(imap["start_hour"]))->val,
                                      (static_cast<ddval<int>*>(imap["start_minute"]))->val);
    if (status != DTDB_SUCCESS) {
        cerr << "DDProcess:: Ancillary initialization failure" << endl;
        return status;
    }
    cerr << "DDProcess:: Initializing Land Process" << endl;
    status = pl_->initialize(imap);
    if (status != DTDB_SUCCESS) {
        cerr << "DDProcess:: Land initialization failure" << endl;
        return status;
    }
    cerr << "DDProcess:: Initializing Ocean Process" << endl;
    status = po_->initialize(imap);
    if (status != DTDB_SUCCESS) {
        cerr << "DDProcess:: Ocean initialization failure" << endl;
        return status;
    }
 
    cerr << "DDProcess:: Processing Granule" << endl;
    status = create_nc4(imap);
    if (status != DTDB_SUCCESS) {
        cerr << "DDProcess:: Output file creation failure" << endl;
        return status;
    }
 
    for (auto &it : imap) {
        string name = (string) it.first;
        if (name != "land_water") {
            delete imap[name];
        }
    }
 
    map<string, ddma<float, 2>*> flmap;
    map<string, ddma<unsigned int, 2>*> uimap;
    map<string, ddma<short, 2>*> shmap;
    map<string, ddma<unsigned char, 2>*> ubmap;
    map<string, ddata*> smap;
    map<string, ddata*> amap;
    map<string, ddata*> pmap;
 
    for (auto &it : dtdb_names) {
        string name = (string) it.first;
        if ( it.second == dtype::FLOAT ) {
            ddma<float, 2>* ddfl =
                new ddma<float, 2>(dtype::FLOAT, {0,0}, {lprw_, pixels_});
            flmap.insert({ name, ddfl });
        } else if ( it.second == dtype::UINT) {
            ddma<unsigned int, 2>* ddin =
                new ddma<unsigned int, 2>(dtype::UINT, {0,0}, {lprw_, pixels_});
            uimap.insert({ name, ddin });
        } else if ( it.second == dtype::SHORT ) {
            ddma<short, 2>* ddsh =
                new ddma<short, 2>(dtype::SHORT, {0,0}, {lprw_, pixels_});
            shmap.insert({ name, ddsh });
        } else if ( it.second == dtype::UBYTE ) {
            ddma<unsigned char, 2>* ddub =
                new ddma<unsigned char, 2>(dtype::UBYTE, {0,0}, {lprw_, pixels_});
            ubmap.insert({ name, ddub });
        }
    }
    ddma<unsigned char, 2>* ddcm =
        new ddma<unsigned char, 2>(dtype::UBYTE, {0,0}, {lprw_, pixels_});
    ubmap.insert({ "cloud_mask", ddcm });
 
    for (auto &it : rhot_band_names) {
        string name = (string) it.first;
        ddma<float, 2>* ddfl =
                new ddma<float, 2>(dtype::FLOAT, {0,0}, {lprw_, pixels_});
        string oname = "rhot_" + name.substr(1);
        flmap.insert({ oname, ddfl });
    }
    for (auto &it : aot_band_names) {
        string name = (string) it.first;
        ddma<float, 2>* ddfl =
                new ddma<float, 2>(dtype::FLOAT, {0,0}, {lprw_, pixels_});
        flmap.insert({ name, ddfl });
    }
    for (auto &it : srf_band_names) {
        string name = (string) it.first;
        ddma<float, 2>* ddfl =
                new ddma<float, 2>(dtype::FLOAT, {0,0}, {lprw_, pixels_});
        flmap.insert({ name, ddfl });
    }
    ddma<unsigned char, 2>* ddlw =
            new ddma<unsigned char, 2>(dtype::UBYTE, {0,0}, {lprw_, pixels_});
 
// for special aot_380
    ddma<float, 2>* ddfl =
            new ddma<float, 2>(dtype::FLOAT, {0,0}, {lprw_, pixels_});
    flmap.insert({ "aot_380", ddfl });
 
    size_t cline = 0;
    while ( cline < lines_ ) {
        imap.clear();
        size_t sy = cline;
        size_t cy = lprw_;
        size_t left = lines_-cline;
        if (left < lprw_) {
            cy = left;
        }
        for (auto &it : smap) {
            delete smap[(string) it.first];
        }
        smap.clear();
        smap = psensor_->read( {sy, 0}, {cy, pixels_} );
        if (static_cast<ddval<int>*>(smap["status"])->val != DTDB_SUCCESS) {
            cerr << "DDProcess:: Read sensor failure at line " << sy << endl;
        }
        delete smap["status"];
        smap.erase("status");
        imap.insert(smap.begin(), smap.end());
        for (auto &it : amap) {
            delete amap[(string) it.first];
        }
        amap.clear();
        amap = pancillary_->read(imap);
        if (static_cast<ddval<int>*>(amap["status"])->val != DTDB_SUCCESS) {
            cerr << "DDProcess:: Read ancillary failure at line " << sy << endl;
        }
        delete amap["status"];
        amap.erase("status");
        imap.insert(amap.begin(), amap.end());
        for ( size_t iy = 0; iy < cy; iy++ ) {
            for ( size_t ix = 0; ix < pixels_; ix++ ) {
                if (ix==5 && iy+sy==1714) {
//                  bool testmode = true;
                }
                if (plw->pts[sy+iy][ix] == 0) {
                    pmap = po_->process({iy,ix}, {cy,pixels_}, imap);
                    ddlw->pts[iy][ix] = 0;
                    if (static_cast<ddval<int>*>(pmap["status"])->val != DTDB_SUCCESS) {
                        cerr << "DDProcess:: Ocean process failure at line " << iy+sy << " pixel " << ix << endl;
                    }
                } else if (plw->pts[sy+iy][ix] == 1) {
                    pmap = pl_->process({iy,ix}, {cy,pixels_}, imap);
                    ddlw->pts[iy][ix] = 1;
                    if (static_cast<ddval<int>*>(pmap["status"])->val != DTDB_SUCCESS) {
                        cerr << "DDProcess:: Land process failure at line " << iy+sy << " pixel " << ix << endl;
                    }
                }
                for (auto &it : flmap) {
                    string name = (string) it.first;
                    if (pmap.find(name) != pmap.end()) {
                        flmap[name]->pts[iy][ix] =
                            (static_cast<ddval<float>*>(pmap[name]))->val;
                    }
                }
                for (auto &it : uimap) {
                    string name = (string) it.first;
                    if (pmap.find(name) != pmap.end()) {
                        uimap[name]->pts[iy][ix] =
                            (static_cast<ddval<unsigned int>*>(pmap[name]))->val;
                    }
                }
                for (auto &it : shmap) {
                    string name = (string) it.first;
                    if (pmap.find(name) != pmap.end()) {
                        shmap[name]->pts[iy][ix] =
                            (static_cast<ddval<short>*>(pmap[name]))->val;
                    }
                }
                for (auto &it : ubmap) {
                    string name = (string) it.first;
                    if (pmap.find(name) != pmap.end()) {
                        ubmap[name]->pts[iy][ix] =
                            (static_cast<ddval<short>*>(pmap[name]))->val;
                    }
                }
                if (pmap.find("cloud_mask") != pmap.end()) {
                    ubmap["cloud_mask"]->pts[iy][ix] =
                            (static_cast<ddval<unsigned char>*>(pmap["cloud_mask"]))->val;
                }
                for (auto &it : pmap) {
                    delete pmap[(string) it.first];
                }
                pmap.clear();
            }
        }
        delete amap["cloud_mask"];
        amap.erase("cloud_mask");
        imap.erase("cloud_mask");
        for (auto &it : flmap) {
            string name = (string) it.first;
            flmap[name]->start = {cline, 0};
            flmap[name]->count = {cy,pixels_};
            imap.insert({name, flmap[name]});
        }
        for (auto &it : uimap) {
            string name = (string) it.first;
            uimap[name]->start = {cline, 0};
            uimap[name]->count = {cy,pixels_};
            imap.insert({name, uimap[name]});
        }
        for (auto &it : shmap) {
            string name = (string) it.first;
            shmap[name]->start = {cline, 0};
            shmap[name]->count = {cy,pixels_};
            imap.insert({name, shmap[name]});
        }
        for (auto &it : ubmap) {
            string name = (string) it.first;
            ubmap[name]->start = {cline, 0};
            ubmap[name]->count = {cy,pixels_};
            imap.insert({name, ubmap[name]});
        }
        ddlw->start = {cline,0};
        ddlw->count = {cy,pixels_};
        imap.insert({"land_water", ddlw});
        cerr << ".";
        status = write_nc4( imap );
        cline += lprw_;
    }
    for (auto &it : imap) {
        delete imap[(string) it.first];
    }
    imap.clear();
 
    delete plw;
    delete pl_;
    delete po_;
 
    cerr << " DONE! " << endl;
 
    if (get_bool(BOOL_STATISTICS)) {
        vector<string> dnames = {"fmf_550","angstrom","aot_410","aot_490",\
                "aot_550","aot_670","aot_865","aot_1240","aot_1610","aot_2250"};
        size_t window = 8;
 
        cerr << "DDProcess:: Generating statistics on decimated samples" << endl;
        status = write_decimated(dnames, window);
        if (status != DTDB_SUCCESS) {
            cerr << "DDProcess:: Write decimated statistics failure" << endl;
        }
        cerr << "ALL DONE! " << endl;
    }
 
    return status;
}
 
/**************************************************************************
 * NAME: create_nc4()
 *
 * DESCRIPTION: Create NetCDF4 product file
 *
 *************************************************************************/
int DDProcess::create_nc4( map<string, ddata*> imap )
{
    int status = DTDB_SUCCESS;
    map<string, ddata*> omap;
 
    bool bshort = (bool) get_bool(BOOL_SHORT_FORMAT);
 
    vector<size_t> sdims, cdims;
    sdims = {0,0};
    cdims = {lines_, pixels_};
 
    for (auto &it : DDProcess::input_names) {
        string name = (string) it.first;
        ddata* odat = new ddata(it.second, sdims, cdims, nullptr);
        omap.insert({name, odat});
    }
    for (auto &it : DDProcess::dtdb_names) {
        string name = (string) it.first;
        auto rslt = find(out_products.begin(), out_products.end(), name);
        if( !get_bool(name) || (rslt == out_products.end())) {
            continue;
        }
        ddata* odat = new ddata(it.second, sdims, cdims, nullptr);
        omap.insert({name, odat});
    }
    for (auto &it : DDProcess::rhot_band_names) {
        string name = (string) it.first;
        ddata* odat = new ddata(dtype::FLOAT, sdims, cdims, nullptr);
        omap.insert({name, odat});
    }
    for (auto &it : DDProcess::aot_band_names) {
        string name = (string) it.first;
        auto rslt = find(out_products.begin(), out_products.end(), name);
        if( !get_bool(name) || (rslt == out_products.end())) {
            continue;
        }
        ddata* odat = new ddata(dtype::FLOAT, sdims, cdims, nullptr);
        omap.insert({name, odat});
    }
    for (auto &it : DDProcess::srf_band_names) {
        string name = (string) it.first;
        auto rslt = find(out_products.begin(), out_products.end(), name);
        if( !get_bool(name) || (rslt == out_products.end())) {
            continue;
        }
        ddata* odat = new ddata(dtype::FLOAT, sdims, cdims, nullptr);
        omap.insert({name, odat});
    }
 
// for special aot_380
    auto rslt = find(out_products.begin(), out_products.end(), "aot_380");
    if( get_bool("aot_380") && (rslt != out_products.end())) {
        ddata* odat = new ddata(dtype::FLOAT, sdims, cdims, nullptr);
        omap.insert({"aot_380", odat});
    }
 
// Create output file
    string output_filepath = get_option(OUTPUT_NC4);
    if (output_filepath.empty()) {
        cerr << "\nDDProcess:: Failure locating output file path.\n" << endl;
        return DTDB_FAIL;
    }
    NcFile* nc_output;
    try {
        nc_output = new NcFile( output_filepath, NcFile::replace );
    }
    catch( NcException& e) {
        e.what();
        cerr << "\nDDProcess:: Failure creating product file: " + output_filepath + "\n" << endl;
        return DTDB_FAIL;
    }
    try {
        string xfilepath = get_option(PRODUCT_XML);
        xml_document doc;
        xml_parse_result result = doc.load_file(xfilepath.c_str());
        if ( result.status !=  status_ok ) {
            cerr << "DDProcess:: Failure opening XML product configuration file " +
                    xfilepath << endl;
            return DTDB_FAIL;
        }
        xml_node products = doc.child("products");
 
        for (xml_node gatt: products.children()) {
            string sname = string(gatt.name());
            if (sname == "group" || sname == "metadata") {
                continue;
            }
            nc_output->putAtt(sname, gatt.child_value());
        }
 
        for (xml_node meta: products.children("metadata"))
        {
            string mname = meta.attribute("name").value();
            if (imap.find(mname) != imap.end()) {
                ddata* dsp = imap[mname];
                dsp->putAtt(nc_output, mname);
            }
        }
 
        NcDim lines_dim = nc_output->addDim( "number_of_lines", lines_ );
        NcDim pixels_dim = nc_output->addDim( "pixels_per_line", pixels_ );
        NcDim all_bands_dim = nc_output->addDim( "number_of_bands", NTWL );
        dmap_.insert({lines_dim.getSize(), lines_dim});
        dmap_.insert({pixels_dim.getSize(), pixels_dim});
        dmap_.insert({all_bands_dim.getSize(), all_bands_dim});
        vector<NcDim> pdims;
 
        for (xml_node group: products.children("group"))
        {
            string gname = group.attribute("name").value();
            if (!get_bool(gname)) {
                continue;
            }
            bool boutput = false;
            for(string str : out_groups) {
                if (str == gname) boutput = true;
            }
            NcGroup grp = nc_output->addGroup(gname);
            for (xml_node product: group.children("product")) {
                string pname = product.attribute("name").value();
                auto rslt = find(out_products.begin(), out_products.end(), pname);
                if(gname==BOOL_DATA && (!get_bool(pname) || (rslt == out_products.end()))) {
                    continue;
                }
                ddata* dsp;
                if (omap.find(pname) != omap.end()) {
                    dsp = omap[pname];
                } else if (imap.find(pname) != imap.end()) {
                    dsp = imap[pname];
                } else {
                    cout << "  --" << pname << " not found" << endl;
                    continue;
                }
                bool bnotshort = (gname == "sensor_band_parameters" ||
                                  gname == "scan_line_attributes"   ||
                                  gname == "navigation_data");
                if (!bnotshort) {
                    dsp->setFormat(bshort);
                }
                pdims.clear();
                for (size_t i=0; i<dsp->count.size(); i++) {
                    NcDim tdim = dmap_[dsp->count[i]];
                    pdims.push_back(tdim);
                }
                NcVar var;
                if(bshort && dsp->type==dtype::FLOAT &&
                        dsp->count.size()==2 && !bnotshort) {
                    var = grp.addVar(pname, NcType(NC_SHORT), pdims);
                } else {
                    var = grp.addVar(pname, NcType((nc_type)dsp->type), pdims);
                }
                vector<size_t> chunksizes = {lprw_, pixels_};
                var.setCompression(bShuffleFilter, bDeflateFilter, deflateLevel);
                var.setChunking(var.nc_CHUNKED, chunksizes);
 
                float vmax,vmin;
                for (xml_node patt: product.children()) {
                    string attr_name = string(patt.name());
                    if (attr_name == "type") {
                        continue;
                    }
                    if (attr_name == "_FillValue") {
                        string pval_str = string(patt.child_value());
                        dsp->setFill(var, pval_str);
                    } else if (attr_name == "valid_min") {
                        vmin = boost::lexical_cast<float> (patt.child_value());
                        dsp->putAtt(var, "valid_min", patt.child_value());
                    } else if (attr_name == "valid_max") {
                        vmax = boost::lexical_cast<float> (patt.child_value());
                        dsp->putAtt(var, "valid_max", patt.child_value());
                    } else if (attr_name == "valid_logmin" ) {
                        if (!bcloudmask_) {
                            vmin = boost::lexical_cast<float> (patt.child_value());
                            dsp->putAtt(var, "valid_min", patt.child_value());
                        }
                    } else if (attr_name == "valid_logmax" ) {
                        if (!bcloudmask_) {
                            vmax = boost::lexical_cast<float> (patt.child_value());
                            dsp->putAtt(var, "valid_max", patt.child_value());
                        }
                    } else {
                        var.putAtt(attr_name, patt.child_value());
                    }
                }
                if(bshort && dsp->type==dtype::FLOAT &&
                        dsp->count.size()==2 && !bnotshort) {
                    float scale = (vmax-vmin)/2.0/MAX_SHORT;
                    float offset = (vmax+vmin)/2.0;
                    var.putAtt("scale_factor", ncFloat, scale);
                    var.putAtt("add_offset", ncFloat, offset);
                    var.putAtt("valid_min", ncShort, -MAX_SHORT);
                    var.putAtt("valid_max", ncShort, MAX_SHORT);
                }
                if (!boutput) {
                    dsp->putVar(var);
                }
            }
        }
    }
    catch( NcException& e) {
        e.what();
        cerr << "\nDDProcess:: Failure initializing product file: " +
                output_filepath + "\n" << endl;
        status = DTDB_FAIL;
    }
 
    for (auto &it : omap) {
        delete omap[(string) it.first];
    }
    omap.clear();
    dmap_.clear();
 
    delete nc_output;
 
    return status;
}
 
/**************************************************************************
 * NAME: write_nc4()
 *
 * DESCRIPTION: Write a single line to output file.
 *
 *************************************************************************/
int  DDProcess::write_nc4( map<string, ddata*> omap )
{
    int status = DTDB_SUCCESS;
 
    bool bshort = (bool) get_bool(BOOL_SHORT_FORMAT);
 
    string output_filepath = get_option(OUTPUT_NC4);
    if (output_filepath.empty()) {
        cerr << "\nDDProcess:: Failure locating output file path.\n" << endl;
        return DTDB_FAIL;
    }
    NcFile* nc_output;
    try {
        nc_output = new NcFile( output_filepath, NcFile::write );
    }
    catch( NcException& e) {
        e.what();
        cerr << "\nDDProcess:: Failure writing to product file: " +
                output_filepath + "\n" << endl;
        return DTDB_FAIL;
    }
 
    string xfilepath = get_option(PRODUCT_XML);
    xml_document doc;
    xml_parse_result result = doc.load_file(xfilepath.c_str());
    if ( result.status !=  status_ok ) {
        cerr << "DDProcess:: Failure opening XML product configuration file " +
                xfilepath << endl;
        return DTDB_FAIL;
    }
    xml_node products = doc.child("products");
 
    for (string gname : out_groups) {
        try {
            if (!get_bool(gname)) {
                continue;
            }
            xml_node group =
                 products.find_child_by_attribute("group", "name", gname.c_str());
            string xname = group.attribute("name").value();
            NcGroup grp = nc_output->getGroup(gname);
            for (xml_node product: group.children("product")) {
                string pname = product.attribute("name").value();
                auto rslt = find(out_products.begin(), out_products.end(), pname);
                if(gname==BOOL_DATA && (!get_bool(pname) || (rslt == out_products.end()))) {
                    continue;
                }
                if (omap.find(pname) != omap.end()) {
                    ddata* dsp = omap[pname];
                    if (!(pname=="latitude" || pname=="longitude")) {
                        dsp->setFormat(bshort);
                    }
                    NcVar var = grp.getVar(pname);
                    if (!var.isNull()) {
                        dsp->putVar(var);
                    }
                }
            }
        }
        catch( NcException& e) {
            e.what();
            cerr << "\nDDProcess:: Failure writing line " + output_filepath + "\n" << endl;
            status = DTDB_FAIL;
        }
    }
    delete nc_output;
 
    return status;
}
 
/**************************************************************************
 * NAME: write_decimated()
 *
 * DESCRIPTION: Compute and write decimated data sets.
 *
 *************************************************************************/
int DDProcess::write_decimated(vector<string> names, size_t& nwin)
{
    int status = DTDB_SUCCESS;
 
    string filepath = get_option(OUTPUT_NC4);
 
    if (filepath.empty()) {
        cout << "DDAlgorithm:: Missing L2 file name for decimation" << endl;
        status = DTDB_FAIL;
        return status;
    }
    NcFile* nc_output;
    try {
        nc_output = new NcFile(filepath, NcFile::write);
    }
    catch( NcException& e) {
        e.what();
        cout << "DDAlgorithm:: Failure opening L2 file for decimation" << endl;
        return DTDB_FAIL;
    }
    size_t lines = nc_output->getDim("number_of_lines").getSize();
    size_t pixels = nc_output->getDim("pixels_per_line").getSize();
    vector<size_t> startp = {0,0};
    vector<size_t> countp = {lines, pixels};
    size_t dlines = lines/nwin+1;
    size_t dpixels = pixels/nwin+1;
    vector<size_t> startd = {0,0,0};
    vector<size_t> countd = {3, dlines, dpixels};
 
    NcGroup dgrp = nc_output->addGroup("statistics");
    NcDim sets_dim = nc_output->addDim( "number_of_views", countd[0] );
    NcDim lines_dim = nc_output->addDim( "number_of_decimated_lines", countd[1] );
    NcDim samples_dim = nc_output->addDim( "Decimated_samples_per_line", countd[2] );
    vector<NcDim> dimsd = {sets_dim, lines_dim, samples_dim};
    ddma<float,3>* ddec = new ddma<float,3>(dtype::FLOAT, startd, countd);
 
    for (auto &it : names) {
        string name = (string) it;
        NcGroup pgrp = nc_output->getGroup("geophysical_data");
        ddma<float,2>* ddat =
                new ddma<float,2>(pgrp, name, dtype::FLOAT, startp, countp);
        if (ddat->ptr == nullptr) {
            delete ddat;
            continue;
        }
 
        fill(ddec->pts.origin(), ddec->pts.origin()+ddec->pts.num_elements(),0);
        string strdec = name +"_"+ to_string(nwin) +"x"+ to_string(nwin);
        NcVar vard = dgrp.addVar(strdec, NcType(NC_FLOAT), dimsd);
        if (vard.isNull()) {
            cout << "DDProcess::write_decimated() Failure to create netcdf vars" << endl;
            return DTDB_FAIL;
        }
        string lname = strdec + " Mean, Standard deviation, and sample counts on grid decimated by " + to_string(nwin);
        vard.putAtt("long_name", lname);
        ddec->setFill(vard, to_string(DFILL_FLOAT));
        vector<size_t> chunksizes = {countd[0], countd[1], countd[2]};
        vard.setCompression(bShuffleFilter, bDeflateFilter, deflateLevel);
        vard.setChunking(vard.nc_CHUNKED, chunksizes);
 
        for (size_t j=0; j<lines; j++) {
            for (size_t i=0; i<pixels; i++) {
                if (ddat->pts[j][i] < DFILL_TEST) continue;
                size_t iidx = i/nwin;
                size_t jidx = j/nwin;
                ddec->pts[0][jidx][iidx] += ddat->pts[j][i];
                ddec->pts[2][jidx][iidx]++;
            }
        }
        for (size_t j=0; j<dlines; j++) {
            for (size_t i=0; i<dpixels; i++) {
                if (ddec->pts[2][j][i] >9)  {
                    ddec->pts[0][j][i] /= ddec->pts[2][j][i];
                } else {
                    ddec->pts[0][j][i] = DFILL_FLOAT;
                }
            }
        }
        for (size_t j=0; j<lines; j++) {
            for (size_t i=0; i<pixels; i++) {
                if (ddat->pts[j][i] < DFILL_TEST) continue;
                size_t iidx = i/nwin;
                size_t jidx = j/nwin;
                ddec->pts[1][jidx][iidx] +=
                        pow(ddat->pts[j][i]-ddec->pts[0][jidx][iidx],2.0);
            }
        }
        for (size_t j=0; j<dlines; j++) {
            for (size_t i=0; i<dpixels; i++) {
                if (ddec->pts[2][j][i] >9)  {
                    ddec->pts[1][j][i] /= (ddec->pts[2][j][i]-1);
                    ddec->pts[1][j][i] = sqrt(ddec->pts[1][j][i]);
                } else {
                    ddec->pts[1][j][i] = DFILL_FLOAT;
                }
            }
        }
        ddec->putVar(vard);
 
        delete ddat;
    }
    delete ddec;
    delete nc_output;
 
    return status;
};