OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
dataday.c
Go to the documentation of this file.
1 /*
2  * dataday.c
3  *
4  * Created on: Jan 22, 2015
5  * Author: rhealy
6  * Functions to calculate dataday based on equatorial crossing
7  * and swath vertices
8  */
9 #include <math.h>
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <stdint.h>
13 #include <string.h>
14 #include "dataday.h"
15 #define PI 3.1415926535897932384626433832795029L
16 
18  time_t starttime, /* (in) swath start time
19  (seconds since 1-Jan-1970 00:00:00 GMT) */
20  float eqxhour, /* (in) sensor's nominal local equator crossing time
21  (expressed in hours) */
22  int isnight, /* (in) =0 for daytime outlines, =1 for nighttime */
23  int dateline, /* (in) indicates whether swath crosses dateline */
24  float west, /* (in) westernmost longitude of swath */
25  float east, /* (in) easternmost longitude of swath */
26  int32_t *dataday0, /* (out) dataday of swath (days since 1-Jan-1970) */
27  int32_t *dataday1 /* (out) 2nd dataday for datline-spanning swaths */
28  ) {
29  time_t reftime;
30  int refday;
31  float refhour;
32 
33  reftime = (time_t) (starttime + (12 - (double) eqxhour) * 3600);
34  refday = reftime / 86400;
35  refhour = (reftime % 86400) / 3600.0;
36 
37  if (dateline == DATELINE_NOT_CROSSED) {
38  *dataday1 = *dataday0 = refday;
39  } else if (refhour < 6) {
40  *dataday0 = refday - 1;
41  *dataday1 = refday;
42  } else if (refhour > 18) {
43  *dataday0 = refday;
44  *dataday1 = refday + 1;
45  } else if (dateline == DATELINE_NORTH_POLE || dateline == DATELINE_SOUTH_POLE) {
46  *dataday1 = *dataday0 = refday;
47  } else {
48  float wod; /* number of degrees west of dateline */
49  float eod; /* number of degrees east of dateline */
50 
51  /* The 180-degree meridian is the dateline. */
52  wod = 180 - west;
53  eod = east + 180;
54  // }
55  if (wod > eod) {
56  *dataday0 = refday - 1;
57  *dataday1 = refday;
58  } else {
59  *dataday0 = refday;
60  *dataday1 = refday + 1;
61  }
62  }
63 }
64 
66  int isnight, /* (in) =0 for daytime outlines, =1 for nighttime */
67  int n, /* (in) number of elements in outline arrays */
68  float *olat, /* (in) outline latitudes in clockwise order */
69  float *olon, /* (in) outline longitudes in clockwise order */
70  float *north, /* (out) northernmost latitude covered by swath */
71  float *south, /* (out) southernmost latitude covered by swath */
72  float *west, /* (out) westernmost longitude covered by swath */
73  float *east, /* (out) easternmost longitude covered by swath */
74  int *dateline /* (out) pole and 180-deg. meridian crossing info */
75  ) {
76 
77  float maxlat, minlat, maxlon, minlon;
78  int p;
79  int w180, e180; /* westward and eastward dateline crossings */
80  int w000, e000; /* the same for the 0-degree dateline */
81  int irreg = 0; /* flag to signal irregular swath outline */
82 
83  /*
84  Traverse the swath outline looking for 180-degree meridian crossings
85  and keeping track of coordinate extrema.
86  */
87  minlat = INVALID_COORD;
88  minlon = INVALID_COORD;
89  maxlat = -INVALID_COORD;
90  maxlon = -INVALID_COORD;
91  w180 = e180 = w000 = e000 = 0;
92  for (p = 1; p <= n; p++) {
93  float deltalon;
94  int pmodn;
95 
96  pmodn = p % n;
97 
98  /*
99  Trap cases where the outline just touches the 180-degree meridian.
100  I don't want those recorded as a dateline crossing unless the outline
101  continues on the other side of the meridian.
102  */
103  if (olon[pmodn] == -180 && olon[p - 1] > 0) olon[pmodn] = 180;
104  else if (olon[pmodn] == 180 && olon[p - 1] < 0) olon[pmodn] = -180;
105 
106  deltalon = olon[pmodn] - olon[p - 1] + 360 * (e180 - w180);
107  if (deltalon < -180) {
108  /* 180-degree meridian crossed heading eastward */
109  e180++;
110  } else if (deltalon > 180) {
111  /* 180-degree meridian crossed heading westward */
112  w180++;
113  } else if (olon[pmodn] >= 0 && olon[p - 1] - 360 * (e180 - w180) < 0) {
114  /* 0-degree meridian crossed heading eastward */
115  e000++;
116  } else if (olon[pmodn] < 0 && olon[p - 1] - 360 * (e180 - w180) >= 0) {
117  /* 0-degree meridian crossed heading westward */
118  w000++;
119  }
120  if (abs(e180 - w180) > 1) {
121  irreg = 1;
122  }
123  if (abs(e000 - w000) > 1) {
124  irreg = 1;
125  }
126  olon[pmodn] += 360 * (e180 - w180);
127 
128  if (olon[pmodn] > maxlon) maxlon = olon[pmodn];
129  if (olon[pmodn] < minlon) minlon = olon[pmodn];
130  if (olat[pmodn] > maxlat) maxlat = olat[pmodn];
131  if (olat[pmodn] < minlat) minlat = olat[pmodn];
132  }
133 
134  /*
135  Set the dateline argument. For daytime outlines, the dateline
136  is the 180-degree meridian. For nighttime outlines, the dateline
137  is the prime meridian.
138  */
139  if (isnight) {
140  if ((w000 + e000) == 0) {
141  /* The swath did not cross the 0-degree meridian. */
142  *dateline = DATELINE_NOT_CROSSED;
143  } else if ((w000 + e000) & 1) {
144  /*
145  An odd number of crossings of the 0-degree meridian means that
146  one of the poles is included in the swath. I am assuming that a
147  swath never includes both poles.
148  */
149  if (90 - maxlat < minlat + 90) {
150  /* It's the north pole. */
151  *dateline = DATELINE_NORTH_POLE;
152  maxlat = 90;
153  } else {
154  /* It's the south pole. */
155  *dateline = DATELINE_SOUTH_POLE;
156  minlat = -90;
157  }
158  minlon = -180;
159  maxlon = 180;
160  } else {
161  /* The swath crosses the 0-degree meridian but doesn't include a pole. */
162  *dateline = DATELINE_CROSSED;
163  }
164  } else {
165  /* Set the dateline argument. */
166  if ((w180 + e180) == 0) {
167  /* The swath did not cross the 180-degree meridian. */
168  *dateline = DATELINE_NOT_CROSSED;
169  } else if ((w180 + e180) & 1) {
170  /*
171  An odd number of crossings of the 180-degree meridian means that
172  one of the poles is included in the swath. Since this code is
173  for sun-synchronous daytime sensors, I am assuming that a swath
174  never includes both poles.
175  */
176  if (90 - maxlat < minlat + 90) {
177  /* It's the north pole. */
178  *dateline = DATELINE_NORTH_POLE;
179  maxlat = 90;
180  } else {
181  /* It's the south pole. */
182  *dateline = DATELINE_SOUTH_POLE;
183  minlat = -90;
184  }
185  minlon = -180;
186  maxlon = 180;
187  } else {
188  /* The swath crosses the 180-deg. meridian but doesn't include a pole. */
189  *dateline = DATELINE_CROSSED;
190  }
191  }
192  if (irreg) *dateline = DATELINE_IRREGULAR;
193 
194  /* Put the longitudes in the range [-180,180]. */
195  for (p = 0; p < n; p++) {
196  while (olon[p] < -180) olon[p] += 360;
197  while (olon[p] >= 180) olon[p] -= 360;
198  }
199  while (minlon < -180) minlon += 360;
200  while (maxlon > 180) maxlon -= 360;
201 
202  /* Copy the values to the calling function. */
203  *north = maxlat;
204  *south = minlat;
205  *west = minlon;
206  *east = maxlon;
207 }
208 
209 int findGeoLimits(int32_t npix, float *lat, int32_t *spix, int32_t *cpix, int32_t *epix) {
210  *spix = -1;
211  *cpix = -1;
212  *epix = -1;
213  // find spix
214  int32_t i;
215  for(i=0; i<npix; i++) {
216  if(lat[i] >= -90.0 && lat[i] <= 90.0) {
217  *spix = i;
218  break;
219  }
220  }
221  if(i==npix) {
222  return 0;
223  }
224 
225  // find epix
226  for(i=npix-1; i>=0; i--) {
227  if(lat[i] >= -90.0 && lat[i] <= 90.0) {
228  *epix = i;
229  break;
230  }
231  }
232 
233  *cpix = (*spix + *epix) / 2;
234  if(lat[*cpix] >= -90.0 && lat[*cpix] <= 90.0) {
235  return 1;
236  }
237  *spix = -1;
238  *cpix = -1;
239  *epix = -1;
240  return 0;
241 }
242 
243 
245  int32_t *year,
246  int32_t *dayOfYear,
247  int32_t *msecondOfDay,
248  int32_t wid, /* (in) width of lat/lon arrays */
249  int32_t hgt, /* (in) height of lat/lon arrays */
250  float **lat, /* (in) pixel latitudes */
251  float **lon, /* (in) pixel longitudes */
252  int32_t n[2], /* (out) number of vertices per outline */
253  float *olat[2], /* (out) outline latitudes (olat[2][n[i]]) */
254  float *olon[2], /* (out) outline longitudes (olon[2][n[i]]) */
255  int8_t **dorn /* (out) 2D array, dorn[hgt][wid], will hold the */
256  /* day (=0) or night (=1) state of each pixel */
257  /* in the input arrays. */
258  ) {
259 
260  struct coord * list[2] = { NULL, NULL };
261  int s, fgs, lgs, night_tally = 0;
262  int8_t isnight = -1;
263  double secondOfDay;
264  int idayOfYear, iyear;
265  int status = 0;
266  int32_t spix; /* first good pix */
267  int32_t cpix; /* center good pix */
268  int32_t epix; /* lasy good pix */
269 
270 
271  // search for first good scan
272  for (s = 0; s < hgt; s++) {
273  if(findGeoLimits(wid, lat[s], &spix, &cpix, &epix)) {
274  break;
275  }
276  }
277 
278  if (s == hgt) {
279  fprintf(stderr, "-W- %s line %d: ", __FILE__, __LINE__);
280  fprintf(stderr, "No valid latitudes found\n ");
281  status = 1;
282  return status;
283  } else {
284  fgs = s;
285  }
286  // search for last good scan
287  for (s = hgt - 1; s > fgs; s--) {
288  if(findGeoLimits(wid, lat[s], &spix, &cpix, &epix)) {
289  break;
290  }
291  }
292 
293  lgs = s;
294  if (lgs == fgs) {
295  fprintf(stderr, "-W- %s line %d: ", __FILE__, __LINE__);
296  fprintf(stderr, "Insufficient valid latitudes found\n ");
297  status = 1;
298  return status;
299  }
300  /*
301  Save some time by first checking the four corners of
302  the scene to see if it is day, night, or mixed.
303  */
304  for (s = fgs; s <= lgs; s += lgs - fgs) {
305 
306  float sv[3], sundist;
307  int p;
308 
309  /* Get the solar position vector. */
310  secondOfDay = msecondOfDay[s] / 1000.;
311  idayOfYear = dayOfYear[s];
312  iyear = year[s];
313 
314  l_sun_(&iyear, &idayOfYear, &secondOfDay, sv, &sundist);
315 
316  if(!findGeoLimits(wid, lat[s], &spix, &cpix, &epix)) {
317  fprintf(stderr, "-W- %s line %d: ", __FILE__, __LINE__);
318  fprintf(stderr, "Insufficient valid latitudes found in scan %d\n", s);
319  status = 1;
320  return status;
321  }
322 
323  for (p = spix; p <= epix; p += epix-spix) {
324 
325  float pv[3];
326  double coslat;
327 
328  /* Compute the position vector of this pixel. */
329  coslat = cos(lat[s][p] * PI / 180);
330  pv[0] = coslat * cos(lon[s][p] * PI / 180);
331  pv[1] = coslat * sin(lon[s][p] * PI / 180);
332  pv[2] = sin(lat[s][p] * PI / 180);
333 
334  if (pv[0] * sv[0] + pv[1] * sv[1] + pv[2] * sv[2] <= 0) {
335  /*
336  If the dot product of the sun vector and the pixel
337  vector is less than or equal to zero (i.e. the angle between
338  the two vectors is greater than or equal to 90 degrees) then the
339  pixel is on the night side of the planet.
340  */
341  night_tally++;
342  }
343  /*
344  I don't expect any 1-pixel-wide scenes, but, just in case,
345  the following line will prevent infinite loops.
346  */
347  if (epix-spix < 1) break;
348  }
349  /*
350  The following line prevents single-scan-line segments
351  from causing an infinite loop.
352  */
353  if (lgs - fgs < 1) break;
354  }
355 
356  if (night_tally == 0 || night_tally == 4) {
357  int dn, i, j;
358 
359  dn = night_tally / 4; /* dn==0 : daytime dn==1 : nighttime */
360 
361  /* Compute the number of vertices in the scene outline. */
362  n[dn] = 2 * (wid + lgs - fgs - 1);
363 
364  /* Allocate space for that many latitudes and longitudes. */
365  MALLOC(olat[dn], float, n[dn]);
366  MALLOC(olon[dn], float, n[dn]);
367  float tlat, tlon;
368  /* ... and fill them in. */
369  j = 0;
370  // fill in first good scan
371  for (i = 0; i < wid; i++) {
372  tlat = lat[fgs][i];
373  tlon = lon[fgs][i];
374  if ((tlat >= -90 && tlat <= 90) &&
375  (tlon >= -180 && tlon <= 180)) {
376  olat[dn][j] = tlat;
377  olon[dn][j] = tlon;
378  j++;
379  }
380 
381  }
382  // fill in end of scan line
383  for (i = fgs+1; i < lgs; i++) {
384  if(findGeoLimits(wid, lat[i], &spix, &cpix, &epix)) {
385  tlat = lat[i][epix];
386  tlon = lon[i][epix];
387  if ((tlat >= -90 && tlat <= 90) &&
388  (tlon >= -180 && tlon <= 180)) {
389  olat[dn][j] = tlat;
390  olon[dn][j] = tlon;
391  j++;
392  }
393 
394  }
395  }
396  // fill in last good scan line backwards
397  for (i = wid - 1; i >= 0; i--) {
398  tlat = lat[lgs][i];
399  tlon = lon[lgs][i];
400  if ((tlat >= -90 && tlat <= 90) &&
401  (tlon >= -180 && tlon <= 180)) {
402  olat[dn][j] = tlat;
403  olon[dn][j] = tlon;
404  j++;
405  }
406 
407  }
408  // fill in beginning of scan lines backwards
409  for (i = lgs - 1; i > fgs; i--) {
410  if(findGeoLimits(wid, lat[i], &spix, &cpix, &epix)) {
411  tlat = lat[i][spix];
412  tlon = lon[i][spix];
413  if ((tlat >= -90 && tlat <= 90) &&
414  (tlon >= -180 && tlon <= 180)) {
415  olat[dn][j] = tlat;
416  olon[dn][j] = tlon;
417  j++;
418  }
419  }
420  }
421  // set the real number of valid indices
422  n[dn] = j;
423  /* The other outline is empty in a non-mixed scene. */
424  n[!dn] = 0;
425  olat[!dn] = NULL;
426  olon[!dn] = NULL;
427 
428  } else {
429  /*
430  This is a mixed scene, so I need to determine
431  the day/night status of every pixel.
432  */
433 
434  for (s = fgs; s <= lgs; s++) {
435 
436  float sv[3], sundist, prevlat, prevlon;
437  int p;
438 
439  /* Get the solar position vector. */
440  secondOfDay = msecondOfDay[s] / 1000.;
441  idayOfYear = dayOfYear[s];
442  iyear = year[s];
443  l_sun_(&iyear, &idayOfYear, &secondOfDay, sv, &sundist);
444 
445  if(!findGeoLimits(wid, lat[s], &spix, &cpix, &epix))
446  continue;
447  for (p = spix; p <= epix; p++) {
448 
449  float pv[3];
450  double coslat;
451  int daynightchange, i;
452 
453  /* Scan the pixels boustrophedonically. */
454  if (s & 1) { /* If this an odd scan line ... */
455  i = epix + spix - p; /* then scan from right to left; */
456  } else { /* otherwise, ... */
457  i = p; /* scan from left to right. */
458  }
459 
460  if (lat[s][i] < -90 ||
461  lat[s][i] > 90 ||
462  lon[s][i] < -180 ||
463  lon[s][i] > 180) continue;
464 
465  /* Compute the position vector of this pixel. */
466  coslat = cos(lat[s][i] * PI / 180);
467  pv[0] = coslat * cos(lon[s][i] * PI / 180);
468  pv[1] = coslat * sin(lon[s][i] * PI / 180);
469  pv[2] = sin(lat[s][i] * PI / 180);
470 
471  if (pv[0] * sv[0] + pv[1] * sv[1] + pv[2] * sv[2] > 0) {
472  /*
473  If the dot product of the sun vector and the pixel
474  vector is greater than zero (i.e. the angle between
475  the two vectors is less than 90 degrees) then the
476  pixel is on the day side of the planet.
477  */
478  daynightchange = isnight != 0;
479  isnight = 0;
480  } else {
481  daynightchange = isnight != 1;
482  isnight = 1;
483  }
484 
485  /*
486  If the pixel is on the edge of the scene or on the edge of a
487  day/night transition, then add its coordinates to the appropriate
488  end (head or tail) of the appropriate outline (day or night).
489  */
490  if (s == fgs || i == epix) {
491  push(&list[isnight], tail, lat[s][i], lon[s][i]);
492  } else if (i == spix) {
493  push(&list[isnight], head, lat[s][i], lon[s][i]);
494  } else if (s == lgs) {
495  if (s & 1) { /* moving towards the left */
496  push(&list[isnight], tail, lat[s][i], lon[s][i]);
497  } else { /* moving towards the right */
498  push(&list[isnight], head, lat[s][i], lon[s][i]);
499  }
500  } else if (daynightchange) {
501  if (s & 1) { /* moving towards the left */
502  push(&list[isnight], tail, lat[s][i], lon[s][i]);
503  push(&list[!isnight], head, prevlat, prevlon);
504  } else { /* moving towards the right */
505  push(&list[isnight], head, lat[s][i], lon[s][i]);
506  push(&list[!isnight], tail, prevlat, prevlon);
507  }
508  }
509  dorn[s][i] = isnight;
510  prevlat = lat[s][i];
511  prevlon = lon[s][i];
512  }
513  }
514  for (s = 0; s < 2; s++) {
515  int count = 0;
516  struct coord *p;
517  float *lats, *lons, *latp, *lonp, prevlat = -999, prevlon = -999;
518 
519  /* Count the number of vertices in this outline. */
520  for (p = list[s]; p != NULL; p = p->next) count++;
521 
522  if (count > 0) {
523  /* Allocate space for that many latitudes and longitudes. */
524  MALLOC(lats, float, count);
525  MALLOC(lons, float, count);
526 
527  /*
528  Copy the latitudes and longitudes from the linked list to
529  the arrays omitting any duplicate coordinates in the process.
530  */
531  latp = lats;
532  lonp = lons;
533  for (p = list[s]; p != NULL; p = p->next) {
534  if (p->lat != prevlat && p->lon != prevlon) {
535  prevlat = *latp++ = p->lat;
536  prevlon = *lonp++ = p->lon;
537  } else {
538  count--;
539  }
540  }
541 
542  if (count <= 0) {
543  fprintf(stderr, "-E- %s line %d: ", __FILE__, __LINE__);
544  fprintf(stderr,
545  "Program logic error. Likely some bad navigation got in - fix the code!.\n");
546  status = 2;
547  }
548 
549  /* Deallocate the memory in the linked list. */
550  if (list[s] != NULL) {
551  struct coord *pn;
552  for (p = list[s]; p != NULL; p = pn) {
553  pn = p->next;
554  free(p);
555  p = NULL;
556  }
557  }
558 
559  /* Return the output arrays to the calling program. */
560  n[s] = count;
561  olat[s] = lats;
562  olon[s] = lons;
563  } else {
564  n[s] = 0;
565  olat[s] = NULL;
566  olon[s] = NULL;
567  }
568  }
569  }
570  return status;
571 }
572 
573 void push(
574  struct coord **vrtx,
575  enum hort hort,
576  float lat,
577  float lon
578  ) {
579  struct coord *newvrtx;
580 
581  MALLOC(newvrtx, struct coord, 1);
582 
583  newvrtx->lat = lat;
584  newvrtx->lon = lon;
585 
586  if (hort == head) {
587  newvrtx->next = *vrtx;
588  *vrtx = newvrtx;
589  } else if (hort == tail) {
590  struct coord *p;
591 
592  newvrtx->next = NULL;
593 
594  if (*vrtx == NULL) {
595  *vrtx = newvrtx;
596  } else {
597  /*
598  Go to the end of the list. Note that
599  the final semicolon is intentional.
600  */
601  for (p = *vrtx; p->next != NULL; p = p->next);
602 
603  p->next = newvrtx;
604  }
605  }
606 }
607 
#define DATELINE_NORTH_POLE
Definition: dataday.h:15
float * hgt
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
list(APPEND LIBS ${PGSTK_LIBRARIES}) add_executable(atteph_info_modis atteph_info_modis.c) target_link_libraries(atteph_info_modis $
Definition: CMakeLists.txt:7
#define NULL
Definition: decode_rs.h:63
float lat
Definition: dataday.h:41
#define INVALID_COORD
Definition: dataday.h:18
@ head
Definition: dataday.h:37
void get_datadays(time_t starttime, float eqxhour, int isnight, int dateline, float west, float east, int32_t *dataday0, int32_t *dataday1)
Definition: dataday.c:17
float * lat
#define DATELINE_IRREGULAR
Definition: dataday.h:17
#define MALLOC(ptr, typ, num)
Definition: bindepths.c:45
float lon
Definition: dataday.h:42
struct coord * next
Definition: dataday.h:43
hort
Definition: dataday.h:36
@ tail
Definition: dataday.h:37
#define DATELINE_NOT_CROSSED
Definition: dataday.h:13
#define DATELINE_SOUTH_POLE
Definition: dataday.h:16
int32 spix
Definition: l1_czcs_hdf.c:21
float * lon
data_t s[NROOTS]
Definition: decode_rs.h:75
int32 epix
Definition: l1_czcs_hdf.c:23
#define PI
Definition: dataday.c:15
Definition: dataday.h:40
void get_coord_extrema(int isnight, int n, float *olat, float *olon, float *north, float *south, float *west, float *east, int *dateline)
Definition: dataday.c:65
#define abs(a)
Definition: misc.h:90
int i
Definition: decode_rs.h:71
int findGeoLimits(int32_t npix, float *lat, int32_t *spix, int32_t *cpix, int32_t *epix)
Definition: dataday.c:209
int npix
Definition: get_cmp.c:27
void l_sun_(int *iyr, int *iday, double *sec, float *sunr, float *rs)
float p[MODELMAX]
Definition: atrem_corl1.h:131
int daynight_outlines(int32_t *year, int32_t *dayOfYear, int32_t *msecondOfDay, int32_t wid, int32_t hgt, float **lat, float **lon, int32_t n[2], float *olat[2], float *olon[2], int8_t **dorn)
Definition: dataday.c:244
void push(struct coord **vrtx, enum hort hort, float lat, float lon)
Definition: dataday.c:573
#define DATELINE_CROSSED
Definition: dataday.h:14
int count
Definition: decode_rs.h:79