A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
get_dataday.cpp
Go to the documentation of this file.
1 /*
2  * Given a L2 file, determine the dataday(s) it contains
3  */
4 
5 #include "get_dataday.h"
6 #include "get_dataday.hpp"
7 #include <boost/geometry.hpp>
8 #include <boost/geometry/geometries/point_xy.hpp>
9 #include <boost/geometry/geometries/polygon.hpp>
10 #include <boost/optional.hpp>
11 #include <get_geospatial.hpp>
12 #include <cstdlib>
13 #include <iostream>
14 #include <netcdf>
15 
16 #include "sensorDefs.h"
17 #include "sensorInfo.h"
18 #include "timeutils.h"
19 #include "version.h"
20 #define _PRINT_(...) if(ddvals::verbosity) printf(__VA_ARGS__)
21 
22 namespace ddvals {
23 boost::optional<float> equatorialCrossingTime;
24 int32_t plus_day{0};
25 float deltaeqcross{0.0};
26 bool verbosity = false;
27 } // namespace ddvals
28 using namespace std;
29 using namespace netCDF;
30 using namespace netCDF::exceptions;
31 namespace bg = boost::geometry;
32 
33 typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree>> Point_t;
34 typedef bg::model::linestring<Point_t> Linestring_t;
35 typedef bg::model::polygon<Point_t> Polygon_t;
36 typedef bg::model::multi_point<Point_t> MultiPoint_t;
37 
38 //void print(std::ostream& stream, const char* format) {
39 // if (!ddvals::verbosity)
40 // return;
41 // stream << format;
42 //}
43 // h
49 Polygon_t parse_wkr_polygon(const std::string& wkt_string) {
50  Polygon_t poly;
51  boost::geometry::read_wkt(wkt_string, poly);
52  return poly;
53 }
54 
55 //template <typename T, typename... Targs>
56 //void print(std::ostream& stream, const char* format, T&& value, Targs&&... fargs) {
57 // if (!ddvals::verbosity)
58 // return;
59 // for (; *format != '\0'; format++) {
60 // if (*format == '%') {
61 // stream << value;
62 // print(stream, format + 1, fargs...);
63 // return;
64 // }
65 // stream << *format;
66 // }
67 //}
68 //
69 //void print(const char* format) {
70 // if (!ddvals::verbosity)
71 // return;
72 // cout << format;
73 //}
74 //
75 //template <typename T, typename... Targs>
76 //void print(const char* format, T&& value, Targs&&... fargs) {
77 // print(std::cout, format, value, fargs...);
78 //}
79 
80 enum class DL { NOT_CROSSED, CROSSED, CROSSED_NORTH_POLE, CROSSED_SOUTH_POLE, CROSSED_IRREGULAR };
81 
82 Polygon_t gRing_to_gPolygon(float* gRingLats, float* gRingLons, size_t length) {
83  Polygon_t gPolygon;
84 
85  for (size_t i = 0; i < length; i++) {
86  double lat = boost::lexical_cast<double>(gRingLats[i]);
87  double lon = boost::lexical_cast<double>(gRingLons[i]);
88  gPolygon.outer().push_back(Point_t{lon, lat});
89  }
90  gPolygon.outer().push_back(
91  Point_t{boost::lexical_cast<double>(gRingLons[0]), boost::lexical_cast<double>(gRingLats[0])});
92 
93  return gPolygon;
94 }
95 
96 void printUsage(int32_t exitStatus) {
97  std::string softwareVersion;
98  softwareVersion += std::to_string(VERSION_MAJOR);
99  softwareVersion += ".";
100  softwareVersion += std::to_string(VERSION_MINOR);
101  softwareVersion += ".";
102  softwareVersion += std::to_string(VERSION_PATCH);
103  softwareVersion += "-";
104  softwareVersion += GITSHA;
105 
106  cout << "get_dataday " << softwareVersion << endl;
107  cout << "\nUsage: get_dataday ifile" << endl;
108  cout << "where:\n\tifile is an L2 file generated by l2gen\n" << endl;
109  cout << "Options:\n\t-h, --help: get this clever little usage statement" << endl;
110  exit(exitStatus);
111 }
112 
113 void get_datadays(time_t starttime, /* (in) swath start time
114  (seconds since 1-Jan-1970 00:00:00 GMT) */
115  float equatorialCrossingTime, /* (in) sensor's nominal local equator
116  crossing time (expressed in hours) */
117  DL dateLineCrossed, /* (in) indicates whether swath crosses dateline */
118  float west, /* (in) westernmost longitude of swath */
119  float east, /* (in) easternmost longitude of swath */
120  int32_t* dataday0, /* (out) dataday of swath (days since 1-Jan-1970) */
121  int32_t* dataday1 /* (out) 2nd dataday for datline-spanning swaths */
122 ) {
123  time_t referenceTime;
124  int referenceDay;
125  float referenceHour;
126 
127  referenceTime = (time_t)(starttime + (12 - (double)equatorialCrossingTime) * 3600);
128  referenceDay = referenceTime / 86400;
129  referenceHour = (referenceTime % 86400) / 3600.0;
130 
131  if (dateLineCrossed == DL::NOT_CROSSED) {
132  *dataday1 = *dataday0 = referenceDay;
133  } else if (referenceHour < 6) {
134  *dataday0 = referenceDay - 1;
135  *dataday1 = referenceDay;
136  } else if (referenceHour > 18) {
137  *dataday0 = referenceDay;
138  *dataday1 = referenceDay + 1;
139  } else if (dateLineCrossed == DL::CROSSED_NORTH_POLE || dateLineCrossed == DL::CROSSED_SOUTH_POLE) {
140  *dataday1 = *dataday0 = referenceDay;
141  } else {
142  float westOfDateline = 180 - west; /* number of degrees west of dateline */
143  float eastOfDateline = east + 180; /* number of degrees east of dateline */
144 
145  if (westOfDateline > eastOfDateline) {
146  *dataday0 = referenceDay - 1;
147  *dataday1 = referenceDay;
148  } else {
149  *dataday0 = referenceDay;
150  *dataday1 = referenceDay + 1;
151  }
152  }
153 }
154 
155 void getEquatorCrossingTime(int32_t sensorID, bool dayNight, time_t starttime, float* equatorialCrossingTime,
156  int32_t* plusDay) {
157  if(std::abs(ddvals::deltaeqcross) > 1e-5)
158  *equatorialCrossingTime = 12.0 + ddvals::deltaeqcross/ 60.0;
159  else
160  {
161  switch (sensorID) {
162  case MERIS:
163  case OLCIS3A:
164  case OLCIS3B:
165  *equatorialCrossingTime = 10.0;
166  break;
167  case MODIST:
168  case OCTS:
169  *equatorialCrossingTime = 10.5;
170  break;
171  case HAWKEYE:
172  case CZCS:
173  case HICO:
174  *equatorialCrossingTime = 12.0;
175  break;
176  case OCI:
177  case OCIS:
178  case SPEXONE:
179  case HARP2:
180  *equatorialCrossingTime = 13.0;
181  break;
182  case MODISA:
183  case VIIRSN:
184  case VIIRSJ1:
185  case VIIRSJ2:
186  *equatorialCrossingTime = 13.5;
187  break;
188  case SEAWIFS:
189  *equatorialCrossingTime = 12.0;
190  int16_t year;
191  int16_t day;
192  double secs;
193  unix2yds(starttime, &year, &day, &secs);
194  if (year > 2002) {
195  // The constant 10957 makes d the number of days since 1 Jan.
196  // 2000
197  int32_t d = starttime / 86400.0 - 10957.0;
198  /*
199  * On 10 July 2010 (d=3843) OrbView-2/SeaWiFS was nearing the
200  * end of its orbit-raising maneuvers. Before these maneuvers
201  * the node-crossing time was drifting further into the
202  afternoon
203  * (first equation below). After the orbit raising, the
204  node-crossing
205  * time was drifting back towards noon (second equation).
206 
207  * Correction equations provided by Fred Patt.
208  */
209  double deg;
210  if (d < 3843) {
211  deg = 7.7517951e-10 * d * d * d - 2.1692192e-06 * d * d + 0.0070669241 * d - 4.1300585;
212  } else {
213  /*
214  * This equation may need to be replaced with a polynomial
215  * once we get more orbit data under our belts. (16-Jul-2010
216  * N.Kuring)
217  * ...note...the above comment was written before the demise
218  * of SeaWiFS on 11-Dec-2010 ...no further modifications
219  * necessary...
220  */
221  deg = -0.024285181 * d + 128.86093;
222  }
223 
224  // The above polynomials yield degrees; convert to hours.
225  float hours = (float)(deg / 15.);
226 
227  *equatorialCrossingTime = 12.0 + hours;
228  }
229 
230  break;
231  default:
232  cerr << "-W- Unknown equator crossing time for sensorID=" << sensorID << endl;
233  cerr << "-W- Assuming local noon crossing" << endl;
234  *equatorialCrossingTime = 12.0;
235  break;
236  }
237  }
238  // Shift the equatorial crossing time by 12 hours if doing night time
239  // binning This is done so that if crossing the dateline the reference hour
240  // is being used to move into a dataday in the same direction in function
241  // get_datadays
242  if (dayNight) {
243  *equatorialCrossingTime = *equatorialCrossingTime - 12;
244  if (*equatorialCrossingTime < 0) {
245  *equatorialCrossingTime = *equatorialCrossingTime + 24;
246  *plusDay = 1; // if we do this, we're effectively going back
247  // a day, and we so we need to add a day to the
248  // output of get_datadays...
249  }
250  }
251  // setting equatorialCrossingTime
252  ddvals::equatorialCrossingTime = *equatorialCrossingTime;
253  ddvals::plus_day = *plusDay;
254 }
255 
256 std::string dataday_to_isodate(int32_t dataday) {
257  double datatime = dataday * 86400.0;
258  std::string isoDate(unix2isodate(datatime, 'G'));
259  return isoDate;
260 }
261 
262 int32_t get_plusday()
263 {
264  return ddvals::plus_day;
265 }
266 
267 float get_equatorial_crossingTime() {
268  if (ddvals::equatorialCrossingTime == boost::none) {
269  std::cerr << "Equatorial Crossing Time is not set. Exiting ...";
270  exit(EXIT_FAILURE);
271  } else {
272  return ddvals::equatorialCrossingTime.get();
273  }
274 }
275 
276 std::pair<int32_t,int32_t> get_datadays(const std::string &filepath)
277 {
278  std::pair<int32_t,int32_t> res;
279  try {
280  NcFile nc_input(filepath, NcFile::read);
281  res = get_datadays(nc_input);
282  nc_input.close();
283 
284  } catch (NcException& e) {
285  e.what();
286  cerr << "\nFailure opening input file: " + filepath + "\n" << endl;
287  printUsage(NC2_ERR);
288  }
289  return res;
290 }
291 
292 std::pair<int32_t, int32_t> get_datadays(const NcFile& nc_input, float deltaeqcross, int night_flag) {
293  // double const earth_radius = 6371.229;
294 
295  // Define the dateline, north and south poles
296  Linestring_t dateline;
297  bg::append(dateline, Point_t(180.0, 89.9999));
298  bg::append(dateline, Point_t(180.0, -89.9999));
299  Point_t northPole = {0, 90.0};
300  Point_t southPole = {0, -90.0};
301  std::string instrument;
302  std::string platform;
303  std::string day_night_flag;
304  float maxEast;
305  float maxWest;
306  std::string sceneStartTimeISO;
307  Polygon_t gPolygon;
308  int32_t sensorID;
309  double sceneStartTime;
310  // starting reading the file
311 
312 
313  Geospatialbounds geo_bounds(nc_input);
314  maxEast = geo_bounds.get_max_lon();
315  maxWest = geo_bounds.get_min_lon();
316  sceneStartTimeISO = geo_bounds.get_time_coverage_start();
317  platform = geo_bounds.get_platform();
318  instrument = geo_bounds.get_instrument();
319  if (platform.empty())
320  std::cerr << "Warning: platform is not found" << std::endl;
321  if (instrument.empty())
322  std::cerr << "Warning: instrument is not found" << std::endl;
323  if (sceneStartTimeISO.empty()) {
324  std::cerr << "Error: start time is not found. Exiting" << std::endl;
325  exit(EXIT_FAILURE);
326  }
327  sensorID = instrumentPlatform2SensorId(instrument.c_str(), platform.c_str());
328  sceneStartTime = isodate2unix(sceneStartTimeISO.c_str());
329  std::string geospatial_bounds = geo_bounds.get_bounds_wkt();
330  gPolygon = parse_wkr_polygon(geospatial_bounds);
331  day_night_flag = geo_bounds.get_day_night_flag();
332  ddvals::deltaeqcross = deltaeqcross;
333 
335  DL dateLineCrossed = DL::NOT_CROSSED;
336  std::vector<Point_t> result;
337  bg::intersection(dateline, gPolygon, result);
338  if (result.size()) {
339  if (bg::within(northPole, gPolygon)) {
340  dateLineCrossed = DL::CROSSED_NORTH_POLE;
341  _PRINT_("# Scene crossed North pole\n");
342  } else if (bg::within(southPole, gPolygon)) {
343  _PRINT_("# Scene crossed South pole\n");
344  dateLineCrossed = DL::CROSSED_SOUTH_POLE;
345  } else {
346  _PRINT_("# Scene crossed dateline\n");
347  dateLineCrossed = DL::CROSSED;
348  }
349  }
350 
351  int32_t dataday0, dataday1;
352  int32_t plusDay = 0;
353  float equatorialCrossingTime;
354  bool nightScene = night_flag;
355  if (day_night_flag.c_str() == std::string("Night")) {
356  nightScene = true;
357  }
358  getEquatorCrossingTime(sensorID, nightScene, sceneStartTime, &equatorialCrossingTime, &plusDay);
359 
360  get_datadays(sceneStartTime, equatorialCrossingTime, dateLineCrossed, maxWest, maxEast, &dataday0,
361  &dataday1);
362 
363  std::string startDataDay = dataday_to_isodate(dataday0 + plusDay);
364  std::string stopDataDay = dataday_to_isodate(dataday1 + plusDay);
365  // Print out day/night flag
366  _PRINT_("Day_Night_Flag=%s\n", day_night_flag.c_str());
367  // Print out the derived dataday(s)
368  if (dataday0 == dataday1) {
369  _PRINT_("DataDay0=%s\n", startDataDay.substr(0, 10).c_str());
370 
371  } else {
372  _PRINT_("DataDay0=%s\nDataDay1=%s\n", startDataDay.substr(0, 10).c_str(), stopDataDay.substr(0, 10).c_str());
373  }
374  return {dataday0 + plusDay, dataday1 + plusDay};
375 }
376 
377 int get_datadays(const char* path, int32_t* day0, int32_t* day1) {
378  const auto data = get_datadays(path);
379  *day0 = data.first;
380  *day1 = data.second;
381  return EXIT_SUCCESS;
382 }
383 
384 void set_verbosity(int val) {
385  if (val != 0) {
386  ddvals::verbosity = true;
387  }
388 }
get_equatorial_crossingTime
float get_equatorial_crossingTime()
Get the equatorial crossing time.
Definition: get_dataday.cpp:267
OLCIS3A
#define OLCIS3A
Definition: sensorDefs.h:32
ddvals::plus_day
int32_t plus_day
Definition: get_dataday.cpp:24
_PRINT_
#define _PRINT_(...)
Definition: get_dataday.cpp:20
get_dataday.hpp
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: GEO_basic.h:72
day
int32_t day
Definition: atrem_corl1v2.h:308
DL::NOT_CROSSED
@ NOT_CROSSED
SPEXONE
#define SPEXONE
Definition: sensorDefs.h:46
instrumentPlatform2SensorId
int instrumentPlatform2SensorId(const char *instrument, const char *platform)
Definition: sensorInfo.c:405
Geospatialbounds::get_platform
const std::string & get_platform() const
Definition: get_geospatial.hpp:54
get_plusday
int32_t get_plusday()
Get the plus day.
Definition: get_dataday.cpp:262
dataday_to_isodate
std::string dataday_to_isodate(int32_t dataday)
Definition: get_dataday.cpp:198
OCI
#define OCI
Definition: sensorDefs.h:42
Linestring_t
bg::model::linestring< Point_t > Linestring_t
Definition: get_dataday.cpp:24
HARP2
#define HARP2
Definition: sensorDefs.h:47
VERSION_MINOR
#define VERSION_MINOR
Definition: version.h:2
hico.value_locate.res
def res
Definition: value_locate.py:66
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
Geospatialbounds::get_bounds_wkt
std::string get_bounds_wkt()
Definition: get_geospatial.cpp:384
DL::CROSSED_NORTH_POLE
@ CROSSED_NORTH_POLE
Geospatialbounds::get_time_coverage_start
const std::string & get_time_coverage_start() const
Definition: get_geospatial.hpp:66
DL::CROSSED_SOUTH_POLE
@ CROSSED_SOUTH_POLE
ddvals::deltaeqcross
float deltaeqcross
Definition: get_dataday.cpp:25
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
run_local.lon
lon
Definition: run_local.py:69
color_dtdb.path
string path
Definition: color_dtdb.py:221
set_verbosity
void set_verbosity(int val)
Set the verbosity.
Definition: get_dataday.cpp:384
DL
DL
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
Polygon_t
bg::model::polygon< Point_t > Polygon_t
Definition: get_dataday.cpp:25
HAWKEYE
#define HAWKEYE
Definition: sensorDefs.h:39
OCIS
#define OCIS
Definition: sensorDefs.h:43
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)
Geospatialbounds::get_instrument
const std::string & get_instrument() const
Definition: get_geospatial.hpp:58
parse_wkr_polygon
Polygon_t parse_wkr_polygon(const std::string &wkt_string)
ttps://www.vertica.com/docs/11.0.x/HTML/Content/Authoring/AnalyzingData/Geospatial/Spatial_Definition...
Definition: get_dataday.cpp:49
get_geospatial.hpp
color_dtdb.result
result
Definition: color_dtdb.py:197
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
color_dtdb.filepath
string filepath
Definition: color_dtdb.py:207
get_dataday.h
DL::CROSSED_IRREGULAR
@ CROSSED_IRREGULAR
color_dtdb.bg
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
VIIRSJ2
#define VIIRSJ2
Definition: sensorDefs.h:44
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
ddvals::verbosity
bool verbosity
Definition: get_dataday.cpp:26
SEAWIFS
#define SEAWIFS
Definition: sensorDefs.h:12
abs
#define abs(a)
Definition: misc.h:90
i
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
isodate2unix
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster