30 #include <sys/types.h>
40 #include <gsl/gsl_math.h>
41 #include <gsl/gsl_spline.h>
42 #include <gsl/gsl_sort.h>
43 #include <gsl/gsl_statistics.h>
51 #define MBAND_NUM_DETECTORS 1
55 static int32_t prevScan = -1;
58 static int geoGeolocationGrp;
59 static int l1bScanLineGrp;
60 static int l1bObservationGrp;
61 static int lonId, latId, senzId, senaId, solzId, solaId, HAMSideId, scanQualityId;
64 static int extract_pixel_start = 0;
66 static short *tmpShort;
67 static unsigned char *tmpByte;
69 static size_t num_scans, num_pixels;
70 static size_t num_blue_bands, num_red_bands, num_SWIR_bands, tot_num_bands;
72 static int firstCall = 1;
73 static double starttime;
74 static double lastvalidtime;
75 static int lastvalidscan = 0;
76 static double time_interval;
78 static float latGeoFillValue = -999.9;
79 static float lonGeoFillValue = -999.9;
80 static short senzGeoFillValue = -32768;
81 static short senaGeoFillValue = -32768;
82 static short solzGeoFillValue = -32768;
83 static short solaGeoFillValue = -32768;
85 static float latL2FillValue = -999.0;
86 static float lonL2FillValue = -999.0;
87 static float senzL2FillValue = -32767;
88 static float senaL2FillValue = -32767;
89 static float solzL2FillValue = -32767;
90 static float solaL2FillValue = -32767;
93 static int Lt_blue_varid, Lt_red_varid, Lt_SWIR_varid;
94 static double **Lt_blue;
95 static double **Lt_red;
96 static double **Lt_SWIR;
116 int ncid_L1B, dimid,
status;
121 printf(
"Opening oci l1b file\n");
122 status = nc_open(
file->name, NC_NOWRITE, &ncid_L1B);
124 fprintf(
stderr,
"-E- %s line %d: nc_open(%s) failed.\n",
125 __FILE__, __LINE__,
file->name);
131 status = nc_inq_dimid(ncid_L1B,
"number_of_scans", &dimid);
133 fprintf(
stderr,
"-E- Error reading num_scans.\n");
136 nc_inq_dimlen(ncid_L1B, dimid, &num_scans);
138 status = nc_inq_dimid(ncid_L1B,
"blue_bands", &dimid);
140 fprintf(
stderr,
"-E- Error reading num_blue_bands.\n");
143 nc_inq_dimlen(ncid_L1B, dimid, &num_blue_bands);
145 status = nc_inq_dimid(ncid_L1B,
"red_bands", &dimid);
147 fprintf(
stderr,
"-E- Error reading num_red_bands.\n");
150 nc_inq_dimlen(ncid_L1B, dimid, &num_red_bands);
152 status = nc_inq_dimid(ncid_L1B,
"SWIR_bands", &dimid);
154 fprintf(
stderr,
"-E- Error reading num_SWIR_bands.\n");
157 nc_inq_dimlen(ncid_L1B, dimid, &num_SWIR_bands);
159 status = nc_inq_dimid(ncid_L1B,
"ccd_pixels", &dimid);
161 fprintf(
stderr,
"-E- Error reading num_pixels.\n");
164 nc_inq_dimlen(ncid_L1B, dimid, &num_pixels);
175 if ((nc_inq_grp_ncid(ncid_L1B,
"scan_line_attributes", &l1bScanLineGrp)) == NC_NOERR) {
177 fprintf(
stderr,
"-E- Error finding scan_line_attributes.\n");
182 status = nc_inq_varid(l1bScanLineGrp,
"HAM_side", &HAMSideId);
185 status = nc_inq_varid(l1bScanLineGrp,
"scan_quality_flags", &scanQualityId);
211 printf(
"OCI L1B Npix :%d Nlines:%d\n", (
int)num_pixels, nline);
215 status = nc_inq_attlen(ncid_L1B, NC_GLOBAL,
"time_coverage_start", &att_len);
219 fltime = (
char *) malloc(att_len + 1);
222 status = nc_get_att_text(ncid_L1B, NC_GLOBAL,
"time_coverage_start", fltime);
224 fltime[att_len] =
'\0';
231 status = nc_inq_attlen(ncid_L1B, NC_GLOBAL,
"time_coverage_end", &att_len);
235 fltime = (
char *) malloc(att_len + 1);
238 status = nc_get_att_text(ncid_L1B, NC_GLOBAL,
"time_coverage_end", fltime);
240 fltime[att_len] =
'\0';
247 time_interval = (
stoptime - starttime)/(num_scans-1);
249 if ((nc_inq_att(ncid_L1B, NC_GLOBAL,
"orbit_number", &vr_type, &vr_len) == NC_NOERR)) {
250 status = nc_get_att_int(ncid_L1B, NC_GLOBAL,
"orbit_number", &orbit_number);
259 if ((nc_inq_grp_ncid(ncid_L1B,
"observation_data", &l1bObservationGrp)) == NC_NOERR) {
261 status = nc_inq_varid(l1bObservationGrp,
"Lt_blue", &Lt_blue_varid);
263 status = nc_inq_varid(l1bObservationGrp,
"Lt_red", &Lt_red_varid);
265 status = nc_inq_varid(l1bObservationGrp,
"Lt_SWIR", &Lt_SWIR_varid);
268 fprintf(
stderr,
"-E- Error finding observation_data.\n");
272 file->sd_id = ncid_L1B;
273 file->nbands = tot_num_bands;
274 file->npix = num_pixels;
277 file->terrain_corrected = 1;
278 file->orbit_number = orbit_number;
282 "Fobar", (
void **) &Fobar);
285 printf(
"file->nbands = %d\n", (
int)
file->nbands);
288 status = nc_inq_grp_ncid(ncid_L1B,
"geolocation_data", &geoGeolocationGrp);
290 status = nc_inq_varid(geoGeolocationGrp,
"longitude", &lonId);
292 status = nc_inq_var_fill(geoGeolocationGrp, lonId,
NULL, &lonGeoFillValue);
294 status = nc_inq_varid(geoGeolocationGrp,
"latitude", &latId);
296 status = nc_inq_var_fill(geoGeolocationGrp, latId,
NULL, &latGeoFillValue);
298 status = nc_inq_varid(geoGeolocationGrp,
"sensor_zenith", &senzId);
300 status = nc_inq_var_fill(geoGeolocationGrp, senzId,
NULL, &senzGeoFillValue);
302 status = nc_inq_varid(geoGeolocationGrp,
"sensor_azimuth", &senaId);
304 status = nc_inq_var_fill(geoGeolocationGrp, senaId,
NULL, &senaGeoFillValue);
306 status = nc_inq_varid(geoGeolocationGrp,
"solar_zenith", &solzId);
308 status = nc_inq_var_fill(geoGeolocationGrp, solzId,
NULL, &solzGeoFillValue);
310 status = nc_inq_varid(geoGeolocationGrp,
"solar_azimuth", &solaId);
312 status = nc_inq_var_fill(geoGeolocationGrp, solaId,
NULL, &solaGeoFillValue);
339 latL2FillValue = info->fillValue;
342 lonL2FillValue = info->fillValue;
345 senaL2FillValue = info->fillValue;
348 senzL2FillValue = info->fillValue;
351 solaL2FillValue = info->fillValue;
354 solzL2FillValue = info->fillValue;
379 int16_t scan_year, scan_day;
382 size_t start[] = { 0, 0, 0 };
383 size_t count[] = { 1, 1, 1 };
388 int band_num, pixel_num;
389 int blue_band_num, red_band_num, SWIR_band_num;
390 int16_t scan_month, scan_dom;
391 double scan_delta_time_ms = 0;
394 for (ip = 0; ip < num_pixels; ip++) {
395 l1rec->pixnum[ip] = ip + extract_pixel_start;
403 tmpShort = (
short *) malloc(num_pixels *
sizeof(
short));
404 tmpByte = (
unsigned char *) malloc(num_pixels);
431 if (scan_delta_time_ms > 0) {
434 lastvalidtime = starttime + scan_delta_time_ms/1000;
436 lastvalidtime = lastvalidtime + (time_interval * (
line - lastvalidscan));
439 unix2yds(lastvalidtime, &scan_year, &scan_day, &scan_sec);
441 l1rec->scantime = lastvalidtime;
443 yd2md(scan_year, scan_day, &scan_month, &scan_dom);
450 count[2] = num_pixels;
451 for (blue_band_num=0; blue_band_num<num_blue_bands; blue_band_num++) {
452 start[0] = blue_band_num;
454 status = nc_get_vara_double(l1bObservationGrp, Lt_blue_varid,
start,
count, &Lt_blue[blue_band_num][0]);
456 fprintf(
stderr,
"-E- Error reading Lt_blue.\n");
466 count[2] = num_pixels;
467 for (red_band_num=0; red_band_num<num_red_bands; red_band_num++) {
468 start[0] = red_band_num;
470 status = nc_get_vara_double(l1bObservationGrp, Lt_red_varid,
start,
count, &Lt_red[red_band_num][0]);
472 fprintf(
stderr,
"-E- Error reading Lt_red.\n");
481 count[2] = num_pixels;
482 for (SWIR_band_num=0; SWIR_band_num<num_SWIR_bands; SWIR_band_num++) {
483 start[0] = SWIR_band_num;
485 status = nc_get_vara_double(l1bObservationGrp, Lt_SWIR_varid,
start,
count, &Lt_SWIR[SWIR_band_num][0]);
487 fprintf(
stderr,
"-E- Error reading Lt_SWIR.\n");
493 for (pixel_num = 0; pixel_num < num_pixels; pixel_num++) {
497 for (band_num = 0; band_num < tot_num_bands; band_num++) {
501 l1rec->Lt[pixel_num*tot_num_bands + band_num] = (
float)Lt_blue[blue_band_num][pixel_num];
502 blue_band_num = blue_band_num + 2;
504 if ((band_num >= 30) && (band_num < 60)) {
505 l1rec->Lt[pixel_num*tot_num_bands + band_num] = (
float)Lt_red[red_band_num][pixel_num];
506 red_band_num = red_band_num + 2;
508 if (band_num >= 60) {
509 l1rec->Lt[pixel_num*tot_num_bands + band_num] = (
float)Lt_SWIR[SWIR_band_num][pixel_num];
510 SWIR_band_num = SWIR_band_num + 2;
530 short scanQualityWarnMask = 2 | 8 | 128;
531 short scanQualityFailMask = 4;
532 static short scanQualityFlag = 0;
534 if (
scan != prevScan) {
536 status = nc_get_var1_short(l1bScanLineGrp, scanQualityId,
start, &scanQualityFlag);
539 if (scanQualityFlag & scanQualityFailMask) {
540 for (ip = 0; ip < num_pixels; ip++)
544 if (scanQualityFlag & scanQualityWarnMask) {
545 for (ip = 0; ip < num_pixels; ip++)
550 static unsigned char HAMSideVal = 0;
551 if (
scan != prevScan) {
553 status = nc_get_var1_uchar(l1bScanLineGrp, HAMSideId,
start, &HAMSideVal);
556 l1rec->mside = HAMSideVal;
562 count[1] = num_pixels;
566 for (
i = 0;
i < num_pixels;
i++)
567 if (
l1rec->lat[
i] == latGeoFillValue)
568 l1rec->lat[
i] = latL2FillValue;
572 for (
i = 0;
i < num_pixels;
i++)
573 if (
l1rec->lon[
i] == lonGeoFillValue)
574 l1rec->lon[
i] = lonL2FillValue;
576 status = nc_get_vara_short(geoGeolocationGrp, solzId,
start,
count, tmpShort);
578 for (
i = 0;
i < num_pixels;
i++)
579 if (tmpShort[
i] == solzGeoFillValue)
580 l1rec->solz[
i] = solzL2FillValue;
582 l1rec->solz[
i] = tmpShort[
i] * 0.01;
584 status = nc_get_vara_short(geoGeolocationGrp, solaId,
start,
count, tmpShort);
586 for (
i = 0;
i < num_pixels;
i++)
587 if (tmpShort[
i] == solaGeoFillValue)
588 l1rec->sola[
i] = solaL2FillValue;
590 l1rec->sola[
i] = tmpShort[
i] * 0.01;
592 status = nc_get_vara_short(geoGeolocationGrp, senzId,
start,
count, tmpShort);
594 for (
i = 0;
i < num_pixels;
i++)
595 if (tmpShort[
i] == senzGeoFillValue)
596 l1rec->senz[
i] = senzL2FillValue;
598 l1rec->senz[
i] = tmpShort[
i] * 0.01;
600 status = nc_get_vara_short(geoGeolocationGrp, senaId,
start,
count, tmpShort);
602 for (
i = 0;
i < num_pixels;
i++)
603 if (tmpShort[
i] == senaGeoFillValue)
604 l1rec->sena[
i] = senaL2FillValue;
606 l1rec->sena[
i] = tmpShort[
i] * 0.01;
638 unsigned char qualityFailMask = 1 | 2;
639 unsigned char qualityWarnMask = 4;
640 for (
i = 0;
i < num_pixels;
i++) {
641 if (tmpByte[
i] & qualityFailMask)
643 if (tmpByte[
i] & qualityWarnMask)
648 static double esdist = -999.9;
649 if (
scan != prevScan) {
653 int32_t yr = (int32_t)
year;
654 int32_t dy = (int32_t)
day;
655 int32_t
msec = (int32_t) (dsec * 1000.0);
777 printf(
"Closing oci l1b file\n");
783 if (tmpShort) free(tmpShort);
784 if (tmpByte) free(tmpByte);