15 #include <boost/algorithm/string.hpp>
16 #include <boost/geometry.hpp>
17 #include <boost/geometry/geometries/point_xy.hpp>
18 #include <boost/geometry/geometries/polygon.hpp>
19 #include <boost/geometry/geometries/box.hpp>
20 #include <boost/foreach.hpp>
42 using namespace netCDF;
43 using namespace netCDF::exceptions;
52 bool enableDebugPrint=
true;
56 #define PI 3.141592653589793
57 #define MTILT_DIMS_2 20
58 #define LTILT_DIMS_2 2
59 #define MAXALLOCPERBIN 20
60 #define EARTH_RADIUS 6378.14
62 static const int max_l3b_products = 400;
69 static int32_t nrows = -1;
71 static int64_t *basebin;
78 static float32 **data_values;
79 static double **data_areas;
80 static int16_t **file_index;
81 static uint8_t **data_quality;
82 static float64 **time_value;
83 static double time_tai93;
108 typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>>
Point_t;
110 typedef bg::model::box<Point_t>
Box_t;
112 #define VERSION "7.0.3"
113 #define PROGRAM "l2bin"
348 if (strcmp(prodname, l2.prodname[
index]) == 0)
350 if (
index == l2.nprod)
357 int32_t ibin = bin - basebin[
krow];
364 (int64_t) bin < basebin[
krow])
373 printf(
"DC: %10ld %12d %8.2f %8.2f\n",
393 if (file_index[ibin] ==
NULL) {
394 file_index[ibin] = (int16_t *)
397 data_values[ibin] = (float32 *)
399 if (data_values[ibin] == 0x0) {
401 printf(
"Allocation failed for data_values[ibin]: %d %s\n",
406 data_areas[ibin] = (
double *) calloc(
n_allocperbin,
sizeof (
double));
407 if (data_areas[ibin] == 0x0) {
409 printf(
"Allocation failed for data_areas[ibin]: %d %s\n",
414 data_quality[ibin] = (uint8_t *)
416 if (data_quality[ibin] == 0x0) {
418 printf(
"Allocation failed for data_quality[ibin]: %d %s\n",
423 time_value[ibin] = (float64 *)
425 if (time_value[ibin] == 0x0) {
427 printf(
"Allocation failed for time_value[ibin]: %d %s\n",
441 if (input.qual_prod[0] != 0) {
445 if (strncmp(input.suite,
"SST", 3) == 0
447 data_quality[ibin][
nobs[ibin]] = 3;
453 time_value[ibin][
nobs[ibin]] = time_tai93;
457 if (input.composite_prod[0] != 0) {
459 f32 = l2_str[
ifile].l2_data[idx][ipixl];
463 if (
nobs[ibin] != 0) {
464 if (strcmp(input.composite_scheme,
"max") == 0) {
477 data_areas[ibin][
nobs[ibin]] = areaFrac;
481 for (int32_t l3_iprod = 0; l3_iprod <
l3b_nprod; l3_iprod++) {
488 if (input.qual_prod[0] != 0)
489 data_quality[ibin][
nobs[ibin]] = 4;
491 if (input.composite_prod[0] != 0) {
502 (l2_str[
ifile].l2_flags[ipixl] >>
526 if (input.composite_prod[0] != 0)
nobs[ibin] = 1;
532 file_index[ibin] = (int16_t *)
533 realloc(file_index[ibin],
537 data_values[ibin] = (float32 *)
538 realloc(data_values[ibin],
541 if (data_values[ibin] == 0x0) {
543 printf(
"Reallocation failed for data_values[ibin]: %d %s\n",
548 data_areas[ibin] = (
double *) realloc(data_areas[ibin],
550 if (data_areas[ibin] == 0x0) {
552 printf(
"Reallocation failed for data_areas[ibin]: %d %s\n",
557 data_quality[ibin] = (uint8_t *)
558 realloc(data_quality[ibin],
561 if (data_quality[ibin] == 0x0) {
563 printf(
"Reallocation failed for data_quality[ibin]: %d %s\n",
568 time_value[ibin] = (float64 *)
569 realloc(time_value[ibin],
572 if (time_value[ibin] == 0x0) {
574 printf(
"Reallocation failed for time_value[ibin]: %d %s\n",
592 if(enableDebugPrint) {
593 printf(
"lat = [%f, %f, %f, %f, %f]\n",
s, n, n,
s,
s);
594 printf(
"lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
595 printf(
"plt.plot(lon, lat)\n\n");
603 if(!bg::disjoint(box, pixelBox)) {
605 if(bg::intersection(box, pixelBox,
output)) {
606 double intersectArea = bg::area(
output);
607 if(intersectArea > 0) {
609 double binArea = (n-
s) * (e-w);
610 areaFrac = intersectArea / binArea;
625 if(enableDebugPrint) {
626 printf(
"lat = [%f, %f, %f, %f, %f]\n",
s, n, n,
s,
s);
627 printf(
"lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
628 printf(
"plt.plot(lon, lat)\n\n");
636 if(!bg::disjoint(box, pixelPoly)) {
637 std::deque<Polygon_t>
output;
638 if(bg::intersection(box, pixelPoly,
output)) {
639 double binArea = (n-
s) * (e-w);
642 double intersectArea = bg::area(
p);
643 if(intersectArea > 0.0) {
645 areaFrac += intersectArea / binArea;
668 areas.emplace(bin, areaFrac);
680 areas.emplace(bin, areaFrac);
690 void getBins(int32_t
ifile, int32_t ipixl, std::map<uint64_t, double> &areas) {
694 double lat0 = l2_str[
ifile].latitude[ipixl];
695 double lon0 = l2_str[
ifile].longitude[ipixl];
703 if(enableDebugPrint) {
704 printf(
"lat=%f\n", lat0);
705 printf(
"lon=%f\n", lon0);
706 printf(
"plt.plot(lon, lat, \"ro\")\n\n");
715 double deltaLat = 180.0 / shape->
getNumRows();
718 if(input.area_weighting == 1) {
723 pixelBox.min_corner().set<0>(l2_str[
ifile].longitude[ipixl] - l2_str[
ifile].lon1[ipixl]);
724 pixelBox.min_corner().set<1>(l2_str[
ifile].latitude[ipixl] - l2_str[
ifile].lat1[ipixl]);
725 pixelBox.max_corner().set<0>(l2_str[
ifile].longitude[ipixl] + l2_str[
ifile].lon1[ipixl]);
726 pixelBox.max_corner().set<1>(l2_str[
ifile].latitude[ipixl] + l2_str[
ifile].lat1[ipixl]);
731 if(enableDebugPrint) {
732 printf(
"lat = [%f, %f, %f, %f, %f]\n",
733 pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
734 pixelBox.min_corner().y(), pixelBox.min_corner().y());
735 printf(
"lon = [%f, %f, %f, %f, %f]\n",
736 pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
737 pixelBox.max_corner().x(), pixelBox.min_corner().x());
738 printf(
"plt.plot(lon, lat)\n\n");
751 lat = lat0 - deltaLat;
759 }
else if(input.area_weighting == 2) {
779 if((lonMax - lonMin) > 180) {
782 lonMax = tmpF + 360.0;
785 pixelBox.min_corner().set<0>(lonMin);
786 pixelBox.min_corner().set<1>(latMin);
787 pixelBox.max_corner().set<0>(lonMax);
788 pixelBox.max_corner().set<1>(latMax);
794 if(enableDebugPrint) {
795 printf(
"lat = [%f, %f, %f, %f, %f]\n",
796 pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
797 pixelBox.min_corner().y(), pixelBox.min_corner().y());
798 printf(
"lon = [%f, %f, %f, %f, %f]\n",
799 pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
800 pixelBox.max_corner().x(), pixelBox.min_corner().x());
801 printf(
"plt.plot(lon, lat)\n\n");
814 lat = lat0 - deltaLat;
829 if(enableDebugPrint) {
830 printf(
"lat = [%f, %f, %f, %f, %f]\n",
831 l2_str[
ifile].lat1[ipixl],
832 l2_str[
ifile].lat2[ipixl],
833 l2_str[
ifile].lat2[ipixl+1],
834 l2_str[
ifile].lat1[ipixl+1],
835 l2_str[
ifile].lat1[ipixl]);
836 printf(
"lon = [%f, %f, %f, %f, %f]\n",
837 l2_str[
ifile].lon1[ipixl],
838 l2_str[
ifile].lon2[ipixl],
839 l2_str[
ifile].lon2[ipixl+1],
840 l2_str[
ifile].lon1[ipixl+1],
841 l2_str[
ifile].lon1[ipixl]);
842 printf(
"plt.plot(lon, lat)\n\n");
847 bg::append(pixelPoly.outer(),
Point_t(l2_str[
ifile].lon1[ipixl], l2_str[
ifile].lat1[ipixl]));
848 bg::append(pixelPoly.outer(),
Point_t(l2_str[
ifile].lon2[ipixl], l2_str[
ifile].lat2[ipixl]));
849 bg::append(pixelPoly.outer(),
Point_t(l2_str[
ifile].lon2[ipixl+1], l2_str[
ifile].lat2[ipixl+1]));
850 bg::append(pixelPoly.outer(),
Point_t(l2_str[
ifile].lon1[ipixl+1], l2_str[
ifile].lat1[ipixl+1]));
851 bg::append(pixelPoly.outer(),
Point_t(l2_str[
ifile].lon1[ipixl], l2_str[
ifile].lat1[ipixl]));
854 bg::correct(pixelPoly);
865 lat = lat0 - deltaLat;
876 int main(
int argc,
char **argv) {
883 int32_t
ifile, jsrow, ipixl;
888 int32_t n_active_files;
892 int32_t n_filled_bins;
893 int64_t total_filled_bins = 0;
898 float32 *slat =
NULL;
899 float32 *elat =
NULL;
900 float32 *clat =
NULL;
901 float32 *slon =
NULL;
902 float32 *elon =
NULL;
903 float32 *clon =
NULL;
909 int32_t igroup, ngroup;
913 static unsigned char *scan_in_rowgroup[
MAXNFILES];
918 static ptrdiff_t stride[] = { 20, 5, 1 };
921 int64_t *binnum_data;
925 time_t diffday_beg, diffday_end;
932 uint32_t flagusemask;
936 int32_t proc_day_beg, proc_day_end;
938 uint8_t *best_qual, qual_max_allowed;
943 float32 northmost = -90.0, southmost = 90.0, eastmost = -180.0, westmost = 180.0;
952 float onorth, osouth, owest, oeast, eqcross;
954 int32_t dataday0, dataday1, startdate, enddate;
955 double ddstime, ddetime, dbldate;
957 float **dolat, **dolon, *lonArray, *latArray;
958 float *dnlat[2], *dnlon[2];
960 int32_t *years, *days, *msecs;
963 vector<NcDim> dims_nc;
965 int32_t small_dims[3];
967 static meta_l2Type meta_l2;
970 static float lonRanges[4];
977 vector<string> l2_prodname;
978 vector<string> l3_prodname;
994 struct mallinfo minfo;
1014 switch(input.area_weighting) {
1027 nfiles = input.files.size();
1028 printf(
"%d input files\n", nfiles);
1037 if ((fileused = (
int *) calloc(nfiles,
sizeof(
int))) ==
NULL) {
1038 printf(
"-E- Problem allocating memory for fileused element of metadata structure\n");
1042 proc_day_beg = input.sday;
1043 proc_day_end = input.eday;
1044 syear = (int32_t) input.sday / 1000.;
1047 eyear = (int32_t) input.eday / 1000.;
1051 if (strcmp(input.resolve,
"1D") == 0) {
1054 }
else if (strcmp(input.resolve,
"HD") == 0) {
1057 }
else if (strcmp(input.resolve,
"QD") == 0) {
1060 }
else if (strcmp(input.resolve,
"36") == 0) {
1063 }
else if (strcmp(input.resolve,
"18") == 0) {
1066 }
else if (strcmp(input.resolve,
"9") == 0) {
1069 }
else if (strcmp(input.resolve,
"4") == 0) {
1072 }
else if (strcmp(input.resolve,
"2") == 0) {
1075 }
else if (strcmp(input.resolve,
"1") == 0) {
1078 }
else if (strcmp(input.resolve,
"H") == 0) {
1081 }
else if (strcmp(input.resolve,
"Q") == 0) {
1084 }
else if (strcmp(input.resolve,
"HQ") == 0) {
1087 }
else if (strcmp(input.resolve,
"HH") == 0) {
1092 printf(
"Grid resolution not defined.\n");
1095 dlat = 180. / nrows;
1096 bool is64bit = (nrows > 2160 * 16);
1098 qual_max_allowed = input.qual_max;
1100 printf(
"Resolution: %s\n", input.resolve);
1101 printf(
"Max Qual Allowed: %d\n", input.qual_max);
1105 strcpy(buf, l2_str[0].flagnames);
1107 printf(
"flagusemask: %d\n", flagusemask);
1108 printf(
"required: %d\n",
required);
1111 if (boost::iequals(input.l3bprod,
"ALL")) {
1112 vector<string> l2prods(l2_str[0].prodname, l2_str[0].prodname+l2_str[0].nprod);
1113 string joined = boost::algorithm::join(l2prods,
":");
1114 strcpy(input.l3bprod, joined.c_str());
1117 string parms = input.parms;
1118 boost::ireplace_all(parms,
"l3bprod = ALL",
"l3bprod = " + joined);
1119 strcpy(input.parms, parms.c_str());
1123 string delim1 =
":, ";
1124 string l3bprod = input.l3bprod;
1125 boost::trim_if(l3bprod, boost::is_any_of(delim1));
1126 vector<string> prodparam;
1127 boost::algorithm::split(prodparam, l3bprod,
1128 boost::is_any_of(delim1));
1131 string delim2 =
";=";
1134 for (
size_t iprod = 0; iprod < prodparam.size(); iprod++) {
1137 vector<string> thisprod;
1138 boost::algorithm::split(thisprod, prodparam[iprod],
1139 boost::is_any_of(delim2));
1142 if(l2_iprod == -1) {
1143 printf(
"-E- Product: %s was not found in the L2 file: %s\n", thisprod[0].c_str(), l2_str[0].
filename);
1146 int32_t tmpThirdDim = l2_str[0].thirdDim[l2_iprod];
1147 if (tmpThirdDim == 1) {
1148 l2_prodname.push_back(thisprod[0]);
1149 l3_prodname.push_back(thisprod[0]);
1151 if (thisprod.size() > 1) {
1152 min_value.push_back((float32) strtod(thisprod[1].c_str(),
NULL));
1157 for (int32_t ibpp=0; ibpp < tmpThirdDim; ibpp++) {
1159 string subprodname = thisprod[0];
1160 string tag = to_string(l2_str[0].
wavelength[ibpp]);
1161 subprodname.append(
"_");
1162 subprodname.append(tag);
1164 if (wave_string !=
",ALL,") {
1165 size_t found = wave_string.find(
","+tag+
",");
1166 if (found == string::npos) {
1171 l2_prodname.push_back(thisprod[0]);
1172 l3_prodname.push_back(subprodname);
1174 if (thisprod.size() > 1) {
1175 min_value.push_back((float32) strtod(thisprod[1].c_str(),
NULL));
1185 printf(
"-E- Number of output products is greater than %d\n", max_l3b_products);
1192 bscan_row[
i] =
NULL;
1193 escan_row[
i] =
NULL;
1196 scan_in_rowgroup[
i] =
NULL;
1207 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
1209 vector<string> operands;
1210 boost::algorithm::split(operands, l2_prodname[iprod],
1211 boost::is_any_of(
"/"));
1216 printf(
"L3 product: \"%s\" not found in L2 file \"%s\".\n",
1225 if (operands.size() > 1) {
1226 boost::replace_all(l3_prodname[iprod],
"/",
"_");
1229 printf(
"L3 product: \"%s\" not found in L2 file \"%s\".\n",
1239 if (input.qual_prod[0] != 0) {
1242 printf(
"Quality product: \"%s\" not found in L2 file \"%s\".\n",
1243 input.qual_prod, l2_str[
ifile].filename);
1250 if (input.composite_prod[0] != 0) {
1253 printf(
"Composite product: \"%s\" not found in L2 file \"%s\".\n",
1254 input.composite_prod, l2_str[
ifile].filename);
1265 if (input.composite_prod[0] != 0) {
1266 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
1267 if(l3_prodname[iprod] == input.composite_prod) {
1273 printf(
"Composite product: \"%s\" not found in L3 product list.\n",
1274 input.composite_prod);
1283 slat = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1284 elat = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1286 bscan_row[
ifile] = (int32_t *) calloc(l2_str[
ifile].nrec,
sizeof(int32_t));
1287 escan_row[
ifile] = (int32_t *) calloc(l2_str[
ifile].nrec,
sizeof(int32_t));
1291 NcFile nc_input(input.files[
ifile], NcFile::read);
1292 nc_input.getAtt(
"instrument").getValues(
instrument);
1293 nc_input.getAtt(
"platform").getValues(
platform);
1294 NcGroup scan_grp = nc_input.getGroup(
"scan_line_attributes");
1295 scan_grp.getVar(
"slat").getVar(slat);
1296 scan_grp.getVar(
"elat").getVar(elat);
1299 }
catch (NcException
const & e) {
1302 }
catch (exception
const & e) {
1303 cerr <<
"-X- " << __FILE__
" line " << __LINE__ <<
": "
1304 << e.what() << endl;
1313 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
1326 for (jsrow = 0; jsrow < l2_str[
ifile].nrec; jsrow++) {
1327 escan_row[
ifile][jsrow] = (int32_t) ((90 + elat[jsrow]) / dlat);
1328 bscan_row[
ifile][jsrow] = (int32_t) ((90 + slat[jsrow]) / dlat);
1330 if (escan_row[
ifile][jsrow] > bscan_row[
ifile][jsrow]) {
1331 k = escan_row[
ifile][jsrow];
1332 escan_row[
ifile][jsrow] = bscan_row[
ifile][jsrow];
1333 bscan_row[
ifile][jsrow] =
k;
1335 escan_row[
ifile][jsrow] -= 10;
1336 bscan_row[
ifile][jsrow] += 10;
1346 n_active_files = nfiles;
1348 if (input.deltaeqcross == 0.0) {
1375 int d = l2_stimes[
ifile] / 86400 - 10957;
1387 deg = 7.7517951e-10 * d * d * d
1388 - 2.1692192e-06 * d * d
1398 deg = -0.024285181 * d + 128.86093;
1402 float hours = (
float) (deg / 15.);
1404 eqcross = 12 + hours;
1412 "-Warning- Unknown equator crossing time for sensorID=%d file=%s\n ",
1414 fprintf(
stderr,
"using crossing time of %f\n",eqcross);
1419 eqcross = 12.0 + (input.deltaeqcross / 60.);
1425 eqcross = eqcross - 12;
1427 eqcross = eqcross + 24;
1434 slon = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1435 elon = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1436 clon = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1437 elat = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1438 slat = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1439 clat = (float32 *) calloc(l2_str[
ifile].nrec,
sizeof(float32));
1440 years = (int32_t *) calloc(l2_str[
ifile].nrec,
sizeof(int32_t));
1441 days = (int32_t *) calloc(l2_str[
ifile].nrec,
sizeof(int32_t));
1442 msecs = (int32_t *) calloc(l2_str[
ifile].nrec,
sizeof(int32_t));
1443 lonArray = (
float *) malloc(
sizeof(
float)*l2_str[
ifile].nrec * l2_str[
ifile].nsamp);
1444 latArray = (
float *) malloc(
sizeof(
float)*l2_str[
ifile].nrec * l2_str[
ifile].nsamp);
1447 NcFile nc_input(input.files[
ifile], NcFile::read);
1449 NcGroup scan_grp = nc_input.getGroup(
"scan_line_attributes");
1450 scan_grp.getVar(
"slon").getVar(slon);
1451 scan_grp.getVar(
"elon").getVar(elon);
1452 scan_grp.getVar(
"clon").getVar(clon);
1453 scan_grp.getVar(
"slat").getVar(slat);
1454 scan_grp.getVar(
"elat").getVar(elat);
1455 scan_grp.getVar(
"clat").getVar(clat);
1456 scan_grp.getVar(
"year").getVar(years);
1457 scan_grp.getVar(
"day").getVar(days);
1458 scan_grp.getVar(
"msec").getVar(msecs);
1460 NcGroup nav_grp = nc_input.getGroup(
"navigation_data");
1461 nav_grp.getVar(
"longitude").getVar(lonArray);
1462 nav_grp.getVar(
"latitude").getVar(latArray);
1464 dims_nc = nav_grp.getVar(
"latitude").getDims();
1465 for (
i = 0;
i < 2;
i++) {
1466 dims[
i] = (int32_t) dims_nc[
i].getSize();
1467 if (dims[
i] / stride[
i] < 5)
1469 small_dims[
i] = dims[
i] / stride[
i];
1474 }
catch (NcException
const & e) {
1477 }
catch (exception
const & e) {
1478 cerr <<
"-X- " << __FILE__
" line " << __LINE__ <<
": "
1479 << e.what() << endl;
1483 dolat = (
float **) malloc(
sizeof(
float *)*small_dims[0]);
1484 dolon = (
float **) malloc(
sizeof(
float *)*small_dims[0]);
1487 for (
i = 0;
i < small_dims[0];
i++) {
1488 dolat[
i] = (
float *) calloc(small_dims[1],
sizeof(
float));
1489 dolon[
i] = (
float *) calloc(small_dims[1],
sizeof(
float));
1494 for (
j = 0;
j < small_dims[0];
j++) {
1495 k = dims[1] * (
j * stride[0]);
1496 for (
i = 0;
i < small_dims[1];
i++) {
1497 dolat[
j][
i] = latArray[
k];
1498 dolon[
j][
i] = lonArray[
k];
1509 brk_scan[
ifile] = 0;
1513 if (strcasecmp(input.prodtype,
"regional") == 0) {
1514 printf(
"%s brk:%5d %5d %3d %6d\n",
1515 buf, brk_scan[
ifile],
1525 for (
i = 0;
i < small_dims[0];
i++) {
1542 if ((dorn = (int8_t **) malloc(small_dims[0] *
sizeof(int8_t *))) ==
NULL) {
1543 printf(
"%s -Error: Cannot allocate memory to dorn array\n",
1547 for (
i = 0;
i < small_dims[0];
i++) {
1548 if ((dorn[
i] = (int8_t *) calloc(small_dims[1],
sizeof(int8_t))) ==
NULL) {
1549 printf(
"%s -Error: Cannot allocate memory to dorn array\n",
1558 diffday_beg = date - startdate;
1559 diffday_end = date - enddate;
1568 small_dims[1], small_dims[0],
1571 n, dnlat, dnlon, dorn);
1573 if (dnStatus != 0) {
1575 "Consider QC fail for file: %s\n...look at the file "
1576 "though, as I might be lying...\n",
1578 if (ret_status == 0) {
1582 if (input.fileuse[0] != 0) {
1583 char *qcFailFile = (
char*) malloc(strlen(input.fileuse) + 10);
1584 strcpy(qcFailFile, input.fileuse);
1585 strcat(qcFailFile,
".qcfail");
1586 fp2 = fopen(qcFailFile,
"w");
1596 if (dnStatus == 0 && n[input.night] > 0) {
1598 &onorth, &osouth, &owest, &oeast, &dateline);
1602 &dataday0, &dataday1);
1603 dataday0 += plusday;
1604 dataday1 += plusday;
1605 if (dataday1 == dataday0) {
1606 if (dataday0 < startdate || dataday1 > enddate)
1607 brk_scan[
ifile] = -9999;
1609 brk_scan[
ifile] = 0;
1611 if (dataday1 < startdate)
1612 brk_scan[
ifile] = -9999;
1613 else if (dataday0 > enddate) {
1614 brk_scan[
ifile] = -9999;
1616 if (dataday1 > enddate)
1617 brk_scan[
ifile] = 1;
1619 brk_scan[
ifile] = -1;
1628 brk_scan[
ifile] = -9999;
1631 for (
i = 0;
i < small_dims[0];
i++) {
1644 if (brk_scan[
ifile] == -9999) n_active_files--;
1656 for (
i = 0;
i < small_dims[0];
i++) {
1667 if (ret_status == 120)
1674 basebin = (int64_t *) calloc(nrows + 1,
sizeof(int64_t));
1675 for (
i = 0;
i <= nrows;
i++) {
1680 printf(
"total number of bins: %ld\n", (
long int) basebin[nrows] - 1);
1685 strcpy(buf, input.ofile);
1686 status = nc_create(buf, NC_NETCDF4, &fileid_w);
1689 status = nc_def_grp(fileid_w,
"level-3_binned_data", &out_grp);
1693 ds_id.
fid = fileid_w;
1702 printf(
"-E- %s:%d Could not define binList variable in output file\n",
1703 __FILE__, __LINE__);
1707 char** prodnames = (
char**) malloc(l3_prodname.size() *
sizeof(
char*));
1708 for(
size_t i=0;
i<l3_prodname.size();
i++) {
1709 prodnames[
i] =
strdup(l3_prodname[
i].c_str());
1713 printf(
"-E- %s:%d Could not define binData variable in output file\n",
1714 __FILE__, __LINE__);
1717 for(
size_t i=0;
i<l3_prodname.size();
i++) {
1727 printf(
"-E- %s:%d Could not define binIndex variable in output file\n",
1728 __FILE__, __LINE__);
1732 if (input.qual_prod[0] != 0) {
1735 printf(
"-E- %s:%d Could not define quality variable in output file\n",
1736 __FILE__, __LINE__);
1743 beg = (int64_t *) calloc(nrows,
sizeof(int64_t));
1744 ext = (int32_t *) calloc(nrows,
sizeof(int32_t));
1753 for (
i = 0;
i < nrows;
i++) {
1757 if (i32 < 0 || i32 >= nrows) {
1758 printf(
"%d %d\n",
i, nrows);
1763 binIndex64nc[
i].
begin = beg[i32];
1767 binIndex32nc[
i].
begin = beg[i32];
1776 printf(
"row_group not defined, using 270.\n");
1784 for (
i = nrows;
i > 0;
i--) {
1785 if ((nrows %
i) == 0) {
1793 printf(
"Input row_group: %d Actual row_group: %d\n",
1800 for (
krow = 0, igroup = 0; igroup < ngroup; igroup++) {
1814 tmnow = localtime(&tnow);
1815 printf(
"krow:%6d out of %6d (%6.2f to %6.2f) ",
1817 ((float32) (
krow) / nrows) * 180 - 90,
1819 printf(
"%s", asctime(tmnow));
1830 scan_in_rowgroup[
ifile] = (
unsigned char *)
1831 calloc(l2_str[
ifile].nrec + 1,
sizeof (
unsigned char));
1833 for (jsrow = 0; jsrow < l2_str[
ifile].nrec; jsrow++) {
1834 scan_in_rowgroup[
ifile][jsrow] = 1;
1837 scan_in_rowgroup[
ifile][jsrow] = 255;
1844 for (jsrow = 0; jsrow < l2_str[
ifile].nrec; jsrow++) {
1845 if (scan_in_rowgroup[
ifile][jsrow] == 1) {
1856 if (within_flag == 0) {
1859 free(scan_in_rowgroup[
ifile]);
1886 data_values = (float32 **) calloc(
n_bins_in_group,
sizeof(float32 *));
1889 data_quality = (uint8_t **) calloc(
n_bins_in_group,
sizeof(uint8_t *));
1895 n_active_files * l2_str[0].nrec * l2_str[0].nsamp / 50000000;
1903 printf(
"%-20s:%8d\n\n",
"# allocated per bin",
n_allocperbin);
1919 if (brk_scan[
ifile] == -9999)
continue;
1923 for (jsrow = 0; jsrow < l2_str[
ifile].nrec; jsrow++) {
1924 if (scan_in_rowgroup[
ifile][jsrow] == 1) {
1928 if (jsrow == l2_str[
ifile].nrec) {
1936 ntilts = l2_str[
ifile].ntilts;
1937 for (
i = 0;
i < ntilts;
i++) {
1938 tilt_flags[
i] = l2_str[
ifile].tilt_flags[
i];
1939 tilt_ranges[0][
i] = l2_str[
ifile].tilt_ranges[0][
i];
1940 tilt_ranges[1][
i] = l2_str[
ifile].tilt_ranges[1][
i];
1947 diffday_beg = date - startdate;
1948 diffday_end = date - enddate;
1952 for (jsrow = 0; jsrow < l2_str[
ifile].nrec; jsrow++) {
1956 if (scan_in_rowgroup[
ifile][jsrow] != 1)
continue;
1961 scan_in_rowgroup[
ifile]);
1963 if (
status == 5)
continue;
1974 ((
double) eqcross - 12) * 3600.);
1976 ((
double) eqcross + 12) * 3600.);
1978 int scan_crosses_dateline = 0;
1979 if (brk_scan[
ifile] != 0) {
1980 int npixls = l2_str[
ifile].nsamp - 1;
1981 float slat = l2_str[
ifile].longitude[0];
1982 float elat = l2_str[
ifile].longitude[npixls];
1983 if (
abs(slat - elat) > 180)
1984 scan_crosses_dateline = 1;
1987 if (scan_crosses_dateline == 0
1988 && (strcasecmp(input.prodtype,
"regional") != 0)) {
1990 if (dbldate < ddstime)
continue;
1991 if (dbldate > ddetime)
continue;
1996 for (
i = 0;
i < ntilts;
i++) {
1997 if ((jsrow + 1) <= tilt_ranges[1][
i]) {
2003 if (l2_str[
ifile].nsamp == 0)
continue;
2005 if ((jsrow % 100) == 0 && input.verbose) {
2006 printf(
"ifile:%4d jsrow:%6d nsamp:%8d\n",
ifile, jsrow, l2_str[
ifile].nsamp);
2014 for (ipixl = 0; ipixl < l2_str[
ifile].nsamp; ipixl++) {
2017 if(ipixl>=510 && ipixl<512)
2018 enableDebugPrint =
true;
2020 enableDebugPrint =
false;
2025 flagcheck = (l2_str[
ifile].l2_flags[ipixl]);
2026 if ((flagcheck & flagusemask) != 0)
2033 if (scan_crosses_dateline == 1) {
2037 if (input.night == 1) {
2039 if ((brk_scan[
ifile] == -1) &&
2040 (diffday_beg == -1) &&
2043 if ((brk_scan[
ifile] == +1) &&
2044 (diffday_end == 0) &&
2049 if ((brk_scan[
ifile] == -1) &&
2050 (diffday_beg <= 0) &&
2053 if ((brk_scan[
ifile] == +1) &&
2054 (diffday_end >= 0) &&
2061 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
2063 for(
int i=0;
i<l2_str[
ifile].thirdDim[l2_iprod];
i++) {
2064 f32 = l2_str[
ifile].l2_data[l2_iprod][ipixl*l2_str[
ifile].thirdDim[l2_iprod]+
i];
2073 if (bad_value == 1)
continue;
2077 if (input.lonwest != 0.0 || input.loneast != 0.0) {
2082 if(input.area_weighting) {
2084 std::map<uint64_t, double> areas;
2086 for(
auto const&
val : areas) {
2103 if (input.meminfo) {
2104 int32_t total_alloc_space;
2107 total_alloc_space = 0;
2111 printf(
"Used space: %10d\n", minfo.uordblks);
2112 printf(
"Allo space: %10d\n", total_alloc_space * (2 +
l3b_nprod * 4));
2118 if (input.verbose) {
2120 tmnow = localtime(&tnow);
2121 printf(
"krow:%5d After data_value fill: %s\n",
krow, asctime(tmnow));
2127 if (
nobs[ibin] > 0 &&
nobs[ibin] < input.minobs)
2130 if (
nobs[ibin] != 0) n_filled_bins++;
2138 if (n_filled_bins > 0) {
2154 if (input.qual_prod[0] != 0 &&
nobs[ibin] > 0) {
2155 best_qual[ibin] = 255;
2156 for (
j = 0;
j <
nobs[ibin];
j++)
2157 if (data_quality[ibin][
j] < best_qual[ibin])
2158 best_qual[ibin] = data_quality[ibin][
j];
2161 for (
j = 0;
j <
nobs[ibin];
j++) {
2162 if ((data_quality[ibin][
j] <= best_qual[ibin]) &&
2163 (data_quality[ibin][
j] <= qual_max_allowed)) {
2165 data_areas[ibin][
k] = data_areas[ibin][
j];
2166 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
2176 if (
nobs[ibin] == 0) n_filled_bins--;
2179 if (
nobs[ibin] != 0) {
2180 bin = (uint64_t) ibin + basebin[
krow];
2197 for (
j = 0;
j <
nobs[ibin];
j++) {
2198 if (file_index[ibin][
j] ==
ifile)
2199 area += data_areas[ibin][
j];
2205 for (
j = 0;
j <
nobs[ibin];
j++)
2206 time_rec += time_value[ibin][
j];
2222 if (latbin > northmost) northmost = latbin;
2223 if (latbin < southmost) southmost = latbin;
2225 float minNeg = lonRanges[0];
2226 float maxNeg = lonRanges[1];
2227 float minPos = lonRanges[2];
2228 float maxPos = lonRanges[3];
2230 minNeg = fmin(minNeg, lonbin);
2231 maxNeg = fmax(maxNeg, lonbin);
2234 minPos = fmin(minPos, lonbin);
2235 maxPos = fmax(maxPos, lonbin);
2238 float max_angle = 90.0;
2239 int8_t lon000 = ((minPos - maxNeg) < max_angle);
2240 int8_t lon180 = ((maxPos - minNeg) > (360.0 - max_angle));
2241 int8_t polar = (lon000 && lon180) ||
2242 (90.0 - fmax(northmost, -1 * southmost) <= FLT_EPSILON);
2244 if (minNeg >= maxNeg) {
2247 }
else if (minPos >= maxPos) {
2257 }
else if (lon000) {
2260 }
else if (lon180) {
2264 lonRanges[0] = minNeg;
2265 lonRanges[1] = maxNeg;
2266 lonRanges[2] = minPos;
2267 lonRanges[3] = maxPos;
2273 if (n_filled_bins == 0)
goto freemem;
2278 printf(
"%-20s:%8d\n",
"# filled bins", n_filled_bins);
2297 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
2304 if (
nobs[ibin] == 0)
continue;
2316 float32 pixval = data_values[ibin][
j *
l3b_nprod + iprod];
2317 double pixarea = data_areas[ibin][
j];
2319 area_file = pixarea;
2320 sum_file = pixval * pixarea;
2321 sum2_file = pixval * pixarea * pixval * pixarea;
2322 prev_file = file_index[ibin][
j];
2325 for (
j = 1;
j <
nobs[ibin];
j++) {
2326 pixval = data_values[ibin][
j *
l3b_nprod + iprod];
2327 pixarea = data_areas[ibin][
j];
2329 if (file_index[ibin][
j] == prev_file) {
2331 area_file += pixarea;
2332 sum_file += pixval * pixarea;
2333 sum2_file += pixval * pixarea * pixval * pixarea;
2338 sum_bin[ibin] += (sum_file / sqrt(area_file));
2339 sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2343 area_file = pixarea;
2344 sum_file = pixval * pixarea;
2345 sum2_file = pixval * pixarea * pixval * pixarea;
2346 prev_file = file_index[ibin][
j];
2351 sum_bin[ibin] += (sum_file / sqrt(area_file));
2352 sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2358 float *productData = (
float *) calloc(2 * n_filled_bins,
sizeof(
float));
2364 if (
nobs[ibin] != 0) {
2365 productData[
i * 2] = sum_bin[ibin];
2366 productData[
i * 2 + 1] = sum2_bin[ibin];
2384 if (input.qual_prod[0] != 0) {
2385 uint8_t *qualityData = (uint8_t *) calloc(n_filled_bins, 1);
2389 if (
nobs[ibin] != 0) {
2390 memcpy(&qualityData[
i], &best_qual[ibin], 1);
2401 if (sum_bin !=
NULL) free(sum_bin);
2402 if (sum2_bin !=
NULL) free(sum2_bin);
2406 binnum_data = (int64_t *) calloc(n_filled_bins,
sizeof(int64_t));
2410 if (
nobs[ibin] != 0) {
2411 binnum_data[
i] = (int64_t) ibin + basebin[
krow];
2413 if (i < 0 || i >= n_filled_bins) {
2414 printf(
"Error: %d %d %d %d\n",
2421 get_beg_ext(n_filled_bins, binnum_data, basebin, nrows, beg,
ext);
2425 total_filled_bins += n_filled_bins;
2431 if (input.verbose) {
2433 tmnow = localtime(&tnow);
2434 printf(
"krow:%5d After bin processing: %s",
krow, asctime(tmnow));
2442 if (i32 < 0 || i32 >= nrows) {
2443 printf(
"Error: %d %d\n",
i,
krow);
2473 if (scan_in_rowgroup[
ifile] !=
NULL) free(scan_in_rowgroup[
ifile]);
2476 if (input.verbose) {
2478 tmnow = localtime(&tnow);
2479 printf(
"krow:%5d Befre free dynic mem: %s",
krow, asctime(tmnow));
2483 if (file_index[
i] !=
NULL) free(file_index[
i]);
2484 if (data_values[
i] !=
NULL) free(data_values[
i]);
2485 if (data_areas[
i] !=
NULL) free(data_areas[
i]);
2486 if (data_quality[
i] !=
NULL) free(data_quality[
i]);
2487 if (time_value[
i] !=
NULL) free(time_value[
i]);
2490 if (input.verbose) {
2492 tmnow = localtime(&tnow);
2493 printf(
"krow:%5d After free dynic mem: %s",
krow, asctime(tmnow));
2508 if (input.verbose) {
2510 tmnow = localtime(&tnow);
2511 printf(
"krow:%5d %s",
krow, asctime(tmnow));
2513 printf(
"total_filled_bins: %ld\n", (
long int) total_filled_bins);
2515 if (total_filled_bins == 0) {
2517 strcat(buf, input.ofile);
2519 printf(
"%s\n", buf);
2535 if (input.fileuse[0] != 0) {
2536 fp2 = fopen(input.fileuse,
"w");
2538 if (brk_scan[
ifile] != -9999) {
2539 fileused[
ifile] = 1;
2546 if (brk_scan[
ifile] != -9999) {
2547 fileused[
ifile] = 1;
2555 printf(
"Writing Global Attributes\n");
2561 if (meta_l2.sensor_name !=
NULL)
2566 strcat(meta_l3b.
title,
" Level-3 Binned Data");
2567 if (meta_l2.mission !=
NULL)
2571 strcat(meta_l3b.
pversion, input.pversion);
2585 double startTime =
yds2unix(2100, 1, 0);
2586 double endTime =
yds2unix(1900, 1, 0);
2589 for (
i = 0;
i < nfiles;
i++) {
2590 if (fileused[
i] == 1) {
2593 meta_l3b.
end_orb = l2_str[
i].orbit;
2600 (float64) (l2_str[
i].
smsec / 1000.0));
2601 if (tmpTime < startTime)
2602 startTime = tmpTime;
2605 (float64) (l2_str[
i].
emsec / 1000.0));
2606 if (tmpTime > endTime)
2611 strcat(meta_l3b.
infiles,
",");
2621 meta_l3b.
north = northmost;
2622 meta_l3b.
south = southmost;
2623 meta_l3b.
east = eastmost;
2624 meta_l3b.
west = westmost;
2627 for (
i = 1;
i < argc;
i++) {
2629 strcat(buf, argv[
i]);
2636 for (iprod = 0; iprod <
l3b_nprod; iprod++) {
2645 printf(
"-E- units metadata is too long. Remove some products from product list.\n");
2649 buf[strlen(buf) - 1] = 0;
2653 int64_t totalNumBins;
2654 totalNumBins = basebin[nrows] - 1;
2658 "L3B", input.suite, input.pversion);
2676 printf(
"Freeing Dynamic Memory\n");
2700 printf(
"Freeing L2 arrays\n");
2709 status = nc_close(fileid_w);