NASA Ocean Color

ocssw V2022

 /*
  * 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;
 }
  

convert2ocrvc.instrument

instrument

Definition: convert2ocrvc.py:29

OLCIS3A

#define OLCIS3A

Definition: sensorDefs.h:32

MDN.parameters.float

float

Definition: parameters.py:23

day

int32_t day

Definition: atrem_corl1v2.h:308

DL::NOT_CROSSED

@ NOT_CROSSED

instrumentPlatform2SensorId

int instrumentPlatform2SensorId(const char *instrument, const char *platform)

Definition: sensorInfo.c:405

dataday_to_isodate

std::string dataday_to_isodate(int32_t dataday)

Definition: get_dataday.cpp:198

Linestring_t

bg::model::linestring< Point_t > Linestring_t

Definition: get_dataday.cpp:24

VERSION_MINOR

#define VERSION_MINOR

Definition: version.h:2

get_datadays

void get_datadays(time_t starttime, float equatorialCrossingTime, DL dateLineCrossed, float west, float east, int32_t *dataday0, int32_t *dataday1)

Definition: get_dataday.cpp:65

VIIRSN

#define VIIRSN

Definition: sensorDefs.h:23

MERIS

#define MERIS

Definition: sensorDefs.h:22

GITSHA

#define GITSHA

Definition: version.h:4

MODIST

#define MODIST

Definition: sensorDefs.h:18

VERSION_PATCH

#define VERSION_PATCH

Definition: version.h:3

DL::CROSSED_NORTH_POLE

@ CROSSED_NORTH_POLE

DL::CROSSED_SOUTH_POLE

@ CROSSED_SOUTH_POLE

string

@ string

Definition: GEO_read_param_file.c:64

VERSION_MAJOR

#define VERSION_MAJOR

Definition: version.h:1

printUsage

void printUsage(int32_t exitStatus)

prints error code

Definition: get_dataday.cpp:48

DL::CROSSED

@ CROSSED

run_local.lat

lat

Definition: run_local.py:69

main

int main(int argc, char **argv)

Definition: get_dataday.cpp:204

run_local.lon

lon

Definition: run_local.py:69

Definition: get_dataday.cpp:28

gRing_to_gPolygon

Polygon_t gRing_to_gPolygon(float *gRingLats, float *gRingLons, size_t length)

Definition: get_dataday.cpp:35

oci_hdr_strip.argc

argc

Definition: oci_hdr_strip.py:7

Polygon_t

bg::model::polygon< Point_t > Polygon_t

Definition: get_dataday.cpp:25

HAWKEYE

#define HAWKEYE

Definition: sensorDefs.h:39

Point_t

bg::model::point< double, 2, bg::cs::spherical_equatorial< bg::degree > > Point_t

Definition: get_dataday.cpp:23

unix2yds

void unix2yds(double usec, short *year, short *day, double *secs)

color_dtdb.result

result

Definition: color_dtdb.py:197

convert2ocrvc.platform

platform

Definition: convert2ocrvc.py:30

DL::CROSSED_IRREGULAR

@ CROSSED_IRREGULAR

color_dtdb.bg

Definition: color_dtdb.py:451

OCTS

#define OCTS

Definition: sensorDefs.h:14

geometry

subroutine geometry

Definition: atrem_app_refl_f90_cubeio.f:1329

getEquatorCrossingTime

void getEquatorCrossingTime(int32_t sensorID, bool dayNight, time_t starttime, float *equatorialCrossingTime, int32_t *plusDay)

Definition: get_dataday.cpp:108

unix2isodate

char * unix2isodate(double dtime, char zone)

Definition: unix2isodate.c:10

MultiPoint_t

bg::model::multi_point< Point_t > MultiPoint_t

Definition: get_dataday.cpp:26

CZCS

#define CZCS

Definition: sensorDefs.h:17

version.h

timeutils.h

seadasutils.ancDB.starttime

starttime

Definition: ancDB.py:314

OLCIS3B

#define OLCIS3B

Definition: sensorDefs.h:41

ddvals::equatorialCrossingTime

boost::optional< float > equatorialCrossingTime

Definition: get_dataday.cpp:23

HICO

#define HICO

Definition: sensorDefs.h:25

SEAWIFS

#define SEAWIFS

Definition: sensorDefs.h:12

int i

Definition: decode_rs.h:71

sensorID

int32_t sensorID[MAXNFILES]

Definition: l2bin.cpp:91

MODISA

#define MODISA

Definition: sensorDefs.h:19

VIIRSJ1

#define VIIRSJ1

Definition: sensorDefs.h:37

sensorInfo.h

sensorDefs.h

attributes

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 attributes

Definition: HISTORY.txt:65

isodate2unix

double isodate2unix(const char *isodate)

Definition: unix2isodate.c:61