get_dataday.cpp

Go to the documentation of this file.

/*
 * Given a L2 file, determine the dataday(s) it contains
 */
 
#include <cstdio>
#include <cstdlib>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <iostream>
#include <netcdf>
#include <stdlib.h>
#include "timeutils.h"
#include "version.h"
#include "sensorInfo.h"
#include "sensorDefs.h"
 
using namespace std;
using namespace netCDF;
using namespace netCDF::exceptions;
namespace bg = boost::geometry;
 
typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree>> Point_t;
typedef bg::model::linestring<Point_t> Linestring_t;
typedef bg::model::polygon<Point_t> Polygon_t;
typedef bg::model::multi_point<Point_t> MultiPoint_t;
 
enum class DL {
    NOT_CROSSED, 
    CROSSED, 
    CROSSED_NORTH_POLE, 
    CROSSED_SOUTH_POLE, 
    CROSSED_IRREGULAR};
 
Polygon_t gRing_to_gPolygon(float* gRingLats, float* gRingLons, size_t length){
    Polygon_t gPolygon;
 
    for (size_t i = 0; i < length; i++){
        double lat = boost::lexical_cast<double>(gRingLats[i]);
        double lon = boost::lexical_cast<double>(gRingLons[i]); 
        gPolygon.outer().push_back(Point_t{lon,lat});
    }
    gPolygon.outer().push_back(Point_t{boost::lexical_cast<double>(gRingLons[0]),boost::lexical_cast<double>(gRingLats[0])});
 
    return gPolygon;
}
 
void printUsage (int32_t exitStatus) {
    std::string softwareVersion;
    softwareVersion += std::to_string(VERSION_MAJOR);
    softwareVersion += ".";
    softwareVersion += std::to_string(VERSION_MINOR);
    softwareVersion += ".";
    softwareVersion += std::to_string(VERSION_PATCH);
    softwareVersion += "-";
    softwareVersion += GITSHA;
 
    cout << "get_dataday " << softwareVersion << endl;
    cout << "\nUsage: get_dataday ifile" << endl;
    cout << "where:\n\tifile is an L2 file generated by l2gen\n" << endl;
    cout << "Options:\n\t-h, --help: get this clever little usage statement" << endl;
    exit(exitStatus);
}
 
void get_datadays(
        time_t starttime, /* (in)  swath start time
                     (seconds since 1-Jan-1970 00:00:00 GMT) */
        float equatorialCrossingTime, /* (in)  sensor's nominal local equator
                              crossing time (expressed in hours) */
        DL dateLineCrossed, /* (in)  indicates whether swath crosses dateline */
        float west, /* (in)  westernmost longitude of swath */
        float east, /* (in)  easternmost longitude of swath */
        int32_t *dataday0, /* (out) dataday of swath (days since 1-Jan-1970) */
        int32_t *dataday1 /* (out) 2nd dataday for datline-spanning swaths */
        ) {
    time_t referenceTime;
    int referenceDay;
    float referenceHour;
 
    referenceTime = (time_t) (starttime + (12 - (double) equatorialCrossingTime)*3600);
    referenceDay = referenceTime / 86400;
    referenceHour = (referenceTime % 86400) / 3600.0;
 
    if (dateLineCrossed == DL::NOT_CROSSED) {
        *dataday1 = *dataday0 = referenceDay;
    } else if (referenceHour < 6) {
        *dataday0 = referenceDay - 1;
        *dataday1 = referenceDay;
    } else if (referenceHour > 18) {
        *dataday0 = referenceDay;
        *dataday1 = referenceDay + 1;
    } else if (dateLineCrossed == DL::CROSSED_NORTH_POLE || dateLineCrossed == DL::CROSSED_SOUTH_POLE) {
        *dataday1 = *dataday0 = referenceDay;
    } else {
        float westOfDateline = 180 - west; /* number of degrees west of dateline */
        float eastOfDateline = east + 180; /* number of degrees east of dateline */
 
        if (westOfDateline > eastOfDateline) {
            *dataday0 = referenceDay - 1;
            *dataday1 = referenceDay;
        } else {
            *dataday0 = referenceDay;
            *dataday1 = referenceDay + 1;
        }
    }
}
 
void getEquatorCrossingTime(int32_t sensorID, bool dayNight, time_t starttime,
        float *equatorialCrossingTime, int32_t *plusDay) {
 
    switch (sensorID) {
        case HICO:
            *equatorialCrossingTime = 12.0; // Made this up - RJH
            break;
        case MODISA:
        case VIIRSN:
        case VIIRSJ1:
            *equatorialCrossingTime = 13.5;
            break;
        case MODIST:
        case MERIS:
        case OCTS:
            *equatorialCrossingTime = 10.5;
            break;
        case OLCIS3A:
        case OLCIS3B:
            *equatorialCrossingTime = 10.0;
            break;
        case CZCS:
            *equatorialCrossingTime = 12.0;
            break;
        case SEAWIFS:
            *equatorialCrossingTime = 12.0;
            int16_t year;
            int16_t day;
            double secs;
            unix2yds(starttime, &year, &day, &secs);
            if (year > 2002) {
                // The constant 10957 makes d the number of days since 1 Jan. 2000
                int32_t d = starttime / 86400.0 - 10957.0;
                /*
                 * On 10 July 2010 (d=3843) OrbView-2/SeaWiFS was nearing the
                 * end of its orbit-raising maneuvers.  Before these maneuvers
                 * the node-crossing time was drifting further into the afternoon
                 * (first equation below).  After the orbit raising, the node-crossing
                 * time was drifting back towards noon (second equation).
 
                 * Correction equations provided by Fred Patt.
                 */
                double deg;
                if (d < 3843) {
 
                    deg = 7.7517951e-10 * d * d * d
                            - 2.1692192e-06 * d * d
                            + 0.0070669241 * d
                            - 4.1300585;
                } else {
                    /*
                     * This equation may need to be replaced with a polynomial
                     * once we get more orbit data under our belts. (16-Jul-2010 N.Kuring)
                     * ...note...the above comment was written before the demise of SeaWiFS
                     * on 11-Dec-2010 ...no further modifications necessary...
                     */
                    deg = -0.024285181 * d + 128.86093;
                }
 
                //The above polynomials yield degrees; convert to hours.
                float hours = (float) (deg / 15.);
 
                *equatorialCrossingTime = 12.0 + hours;
            }
 
            break;
        case HAWKEYE:
            *equatorialCrossingTime = 12.0;
            break;
        default:
            cerr << "-W- Unknown equator crossing time for sensorID=" << sensorID << endl;
            cerr << "-W- Assuming local noon crossing" << endl;
            *equatorialCrossingTime = 12.0;
            break;
 
    }
    // Shift the equatorial crossing time by 12 hours if doing night time binning
    // This is done so that if crossing the dateline the reference hour is being used
    // to move into a dataday in the same direction in function get_datadays
    if (dayNight) {
        *equatorialCrossingTime = *equatorialCrossingTime - 12;
        if (*equatorialCrossingTime < 0) {
            *equatorialCrossingTime = *equatorialCrossingTime + 24;
            *plusDay = 1; // if we do this, we're effectively going back
            // a day, and we so we need to add a day to the
            // output of get_datadays...
        }
    }
}
 
std::string dataday_to_isodate(int32_t dataday) {
    double datatime = dataday * 86400.0;
    std::string isoDate(unix2isodate(datatime,'G'));
    return isoDate;
}
 
int main(int argc, char** argv) {
    //double const earth_radius = 6371.229;
 
    // Define the dateline, north and south poles
    Linestring_t dateline;
    bg::append(dateline, Point_t(180.0, 89.9999));
    bg::append(dateline, Point_t(180.0, -89.9999));
    Point_t northPole = {0, 90.0};
    Point_t southPole = {0, -90.0};
 
    if (argc < 2) {
        printUsage(0);
    }
 
    string inputFile = argv[1];
    if (inputFile == "-h" || inputFile == "--help" ){
        printUsage(0);
    }
    // Open the file for read access
    
    NcFile* nc_input;
    try {
        nc_input = new NcFile(inputFile, NcFile::read );
    }
    catch( NcException& e) {
        e.what();
        cerr << "\nFailure opening  input file: " + inputFile + "\n" << endl;
        printUsage(NC2_ERR);
    }   
    
    multimap <string, NcGroupAtt> attributes;
    // Get global attributes
    char instrument[255];
    char platform[255];
    char day_night_flag[6];
    float maxEast;
    float maxWest;
    char sceneStartTimeISO[30];
 
    attributes = nc_input->getAtts();
    (attributes.find("geospatial_lon_max")->second).getValues(&maxEast);
    (attributes.find("geospatial_lon_min")->second).getValues(&maxWest);
    (attributes.find("time_coverage_start")->second).getValues(sceneStartTimeISO);
    (attributes.find("instrument")->second).getValues(instrument);
    (attributes.find("platform")->second).getValues(platform);
    (attributes.find("day_night_flag")->second).getValues(day_night_flag);
 
    int32_t sensorID = instrumentPlatform2SensorId(instrument, platform);
 
    double sceneStartTime = isodate2unix(sceneStartTimeISO);
   // Get navigation_data group attributes
    NcGroup navGroup = nc_input->getGroup("navigation_data");
 
    float* gRingLats;
    float* gRingLons;
    attributes = navGroup.getAtts();
    size_t length = (attributes.find("gringpointlatitude")->second).getAttLength();
 
    gRingLats = (float*) calloc(length,sizeof(float));
    gRingLons = (float*) calloc(length,sizeof(float));
    
    (attributes.find("gringpointlatitude")->second).getValues(gRingLats);
    (attributes.find("gringpointlongitude")->second).getValues(gRingLons);
 
    DL dateLineCrossed = DL::NOT_CROSSED;
 
    Polygon_t gPolygon = gRing_to_gPolygon(gRingLats,gRingLons,length);
    free(gRingLats);
    free(gRingLons);
    nc_input->close();
    std::vector<Point_t> result;
    bg::intersection(dateline,gPolygon,result);
    if (result.size()) {
        if (bg::within(northPole, gPolygon)) {
            dateLineCrossed = DL::CROSSED_NORTH_POLE;
            cout << "# Scene crossed North pole" << endl;
        } else if (bg::within(southPole, gPolygon)) {
            cout << "# Scene crossed South pole" << endl;
            dateLineCrossed = DL::CROSSED_SOUTH_POLE;
        } else {
            cout << "# Scene crossed dateline" << endl;
            dateLineCrossed = DL::CROSSED;
        }
    }
 
    int32_t dataday0, dataday1;
    int32_t plusDay = 0;
    float equatorialCrossingTime;
    bool nightScene = false;
    if (std::string(day_night_flag) == "Night"){
        nightScene = true;
    }
    getEquatorCrossingTime(sensorID, nightScene, sceneStartTime, &equatorialCrossingTime, &plusDay);
 
    get_datadays(sceneStartTime, equatorialCrossingTime, dateLineCrossed,
            maxWest, maxEast, &dataday0, &dataday1);
 
    std::string startDataDay = dataday_to_isodate(dataday0 + plusDay);
    std::string stopDataDay = dataday_to_isodate(dataday1 + plusDay);
    // Print out day/night flag
    cout << "Day_Night_Flag=" << day_night_flag << endl;
    // Print out the derived dataday(s)
    if (dataday0 == dataday1) {
        cout << "DataDay0=" << startDataDay.substr(0, 10) << endl;
 
    } else {
        cout << "DataDay0=" << startDataDay.substr(0, 10) << "\nDataDay1=" << stopDataDay.substr(0, 10) << endl;
    }
 
    return 0;
}