NASA Ocean Color

ocssw V2022

 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 #include <unistd.h>
  
 #include <sstream>
 #include <iomanip> 
 #include <stdint.h>
 #include <iostream>
 #include <fstream>
  
 #include "cdl_utils.h"
 #include "netcdf.h"
 #include "nc4utils.h"
 #include "nc_gridutils.h"
 #include "aviris.h"
 #include <timeutils.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
 #include <gsl/gsl_sort_double.h>
  
 #define SENSOR "ocia"
 #define CDLFILE "OCI_Level-1B_Data_Structure.cdl"
 #define NBANDS 61 // number of OCI bands
 #define FILL -32767
  
 extern "C" void handle_error(int status);
  
 using namespace std;
  
  
 #define VERSION "1.00"
  
 size_t interpLt(size_t startJ, size_t npixels, size_t k, size_t nbMax,
         float ociaWavelength, size_t iwave, aviris4ocia_t* data,
         gsl_interp_accel* acc, float* Lt) {
     size_t j, n;
     size_t i0 = startJ, i1 = startJ + 5;
     float minWave, maxWave;
  
     double* x = (double*) (calloc(nbMax, sizeof (double)));
     double* y = (double*) (calloc(nbMax, sizeof (double)));
     size_t iter = 0;
     while (i0 >= 0 && i1 < nbMax) {
         n = 0;
         minWave = 99999;
         maxWave = -1;
         for (j = i0; j < i1; j++) {
             if (data->Lt[j * npixels + k] > 0) {
                 x[n] = data->wave[j];
                 y[n] = data->Lt[j * npixels + k];
                 n++;
                 if (minWave > data->wave[j]) minWave = data->wave[j];
                 if (maxWave < data->wave[j]) maxWave = data->wave[j];
             }
         }
  
         if (n > 2 && minWave <= ociaWavelength && maxWave >= ociaWavelength) break;
         if (minWave > ociaWavelength || n < 3)
             i0--;
         if (maxWave < ociaWavelength || n < 3)
             i1++;
         iter++;
     }
     if (n > 2) {
         // Allocate space for spline
         gsl_spline* spline = gsl_spline_alloc(gsl_interp_cspline, n);
         gsl_sort2(x, 1, y, 1, n);
         if (ociaWavelength > x[0]) {
             // Sort wavelengths and corresponding Lt values
             gsl_sort2(x, 1, y, 1, n);
             // Initiate spline
             gsl_spline_init(spline, x, y, n);
             // Generate interpolated values
             double yi = gsl_spline_eval(spline, ociaWavelength, acc);
             Lt[iwave * npixels + k] = yi;
         } else {
             Lt[iwave * npixels + k] = FILL;
         }
         iwave++;
         gsl_spline_free(spline);
     }
     free(x);
     free(y);
     return iwave;
 }
  
 //    Modification history:
 //  Programmer     Organization   Date     Ver   Description of change
 //  ----------     ------------   ----     ---   ---------------------
 //  Joel Gales     FutureTech     12/18/14 0.01  Original development
 //  Joel Gales     FutureTech     01/12/15 0.10  Add support for writing
 //                                               interpolated ORCA Lt values
 //                                               and navigation
 //  Rick Healy     SAIC           06/06/16 1.00  Integrated L2gen L1 reader for Aviris
  
 int main(int argc, char* argv[]) {
     int status;
     size_t sline = 0, eline = 0;
     char filename[FILENAME_MAX], hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX],
             navfile[FILENAME_MAX], gainfile[FILENAME_MAX];
     ;
     char *filedir;
     double scantime_first, scantime_last;
     static int first_good_scantime = 0;
     cout << "avirisbil2oci " << VERSION << " ("
             << __DATE__ << " " << __TIME__ << ")" << endl;
  
     if (argc < 3) {
         cout << endl <<
                 "avirisbil2oci input_AVIRIS_file output_OCIA_file <sline (default=0)> <eline (default=last row) <img file(optional)> <nav file(optional)> <gain file(optional)>>"
                 << endl;
         return 0;
     }
  
     if (argc >= 4) {
         sscanf(argv[3], "%ld", &sline);
         if (sline < 0) sline = 0;
     }
  
     if (argc >= 5) {
         sscanf(argv[4], "%ld", &eline);
         if (eline < 0) eline = 0;
     }
     if (argc >= 6) {
         sscanf(argv[5], "%s", imgfile);
         if (eline < 0) eline = 0;
     }
     if (argc >= 7) {
         sscanf(argv[6], "%s", navfile);
         if (eline < 0) eline = 0;
     }
     if (argc >= 8) {
         sscanf(argv[7], "%s", gainfile);
         if (eline < 0) eline = 0;
     }
     //  aviris4orca_t *&data;
  
     if ((filedir = getenv("OCDATAROOT")) == NULL) {
         printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
         return (-1);
     }
     strcpy(filename, filedir);
  
     strcat(filename, "/");
     strcat(filename, SENSOR);
     strcat(filename, "/");
  
     strcat(filename, CDLFILE);
  
     if (!checkAvProcessFile(argv[1], hdrfile, imgfile, navfile, gainfile, FILENAME_MAX))
         strncpy(hdrfile, argv[1], FILENAME_MAX);
  
     aviris4ocia_t **data, *temp;
     float fwhm[NBANDS];
     data = (aviris4ocia_t **) malloc(sizeof (aviris4ocia_t *));
     temp = open_aviris(hdrfile, imgfile, navfile, gainfile, data);
     size_t nscans, ip, npixels = temp->npix;
  
     if (eline <= 0) eline = temp->nscans;
  
     if (eline == sline && sline > 0) sline--;
     nscans = eline - sline;
  
     int ncid, grpID;
     //char** dim_names, uint32_t* dim_size, size_t n_dims
     const char *dimNames[] = {"number_of_scans", "number_of_pixels", "number_of_bands"};
     size_t dimSize[] = {nscans, npixels, NBANDS};
     size_t n_dims = 3;
  
     // Create ORCA output file
     static ncdfFile ociafile;
     ociafile.cdlCreateDim(argv[2], filename, dimNames, dimSize, n_dims, nscans);
  
     //  nc_def_dim(orcafile.ncid, "number_of_lines", nscans, &dimid[0]);
     //  nc_def_dim(orcafile.ncid, "pixels_per_line", nscans, &dimid[1]);
     //
     // Copy AVIRIS nav fields to ORCA
     //
     ncid = ociafile.getNcid();
  
     size_t nbands_av = temp->numBands;
     size_t nbands_ocia = NBANDS;
  
     size_t n = 0, iwave = 0;
     istringstream istr;
     string str;
     char isodatetime_first[30], isodatetime_last[30];
  
     // Allocate space for AVIRIS Lt (short) values
     //  size_t start[3] = {0,0,0};
     //  size_t count[3] = {1,npixels,0};
  
     gsl_interp_accel *acc = gsl_interp_accel_alloc();
  
     float *Lt = (float *) calloc(npixels*nbands_ocia, sizeof (float));
     float *ociaWavelength = (float *) calloc(nbands_ocia, sizeof (float));
     double* utc = (double*) (calloc(nscans, sizeof (double)));
     float* alt = (float*) (calloc(nscans, sizeof (float)));
     size_t j;
  
     grpID = ociafile.getGid("observation_data");
  
     for (size_t iscan = sline; iscan < eline; iscan++) {
  
  
         // Read AVIRIS Lt's for this scan
         if (read_aviris(temp, iscan)) {
  
             cout << "Can't read file " << argv[1] << endl;
             exit(1);
  
         }
         for (ip = 0; ip < npixels; ip++) for (iwave = 0; iwave < nbands_ocia; iwave++) Lt[iwave * npixels + ip] = FILL;
  
         if ((iscan % 100) == 0) {
             cout << "Processing scan: " << iscan << " out of " << nscans << endl;
             //printf("year=%d doy=%d msec=%d month=%d day=%d hour=%d min=%d sec=%f\n ",temp->year, temp->doy, temp->msec,temp->month,temp->day,temp->hour,temp->min,temp->sec);
         }
         // For each pixel in scan ...
         for (size_t k = 0; k < npixels; k++) {
  
             // Get number of good AVIRIS wavelengths for each pixel
             n = 0;
             j = 0;
             for (j = 0; j < nbands_av; j++) {
                 if (temp->Lt[j * npixels + k] > 0) n++;
             }
  
             // Bail if no good AVIRIS Lt values
             if (n == 0) continue;
  
  
             // Generate x & y arrays for spline
             // x - AVIRIS wavelengths
             // y - AVIRIS Lt values
             iwave = 0;
             for (j = 0; j < 31; j++) {
                 if (temp->Lt[j * npixels + k] > 0) {
                     ociaWavelength[iwave] = temp->wave[j];
                     Lt[iwave * npixels + k] = temp->Lt[j * npixels + k];
                     fwhm[iwave] = temp->fwhm[j];
                 }
                 iwave++;
             }
             for (j = 33; j < 57; j++) {
                 if (temp->Lt[j * npixels + k] > 0) {
                     ociaWavelength[iwave] = temp->wave[j];
                     Lt[iwave * npixels + k] = temp->Lt[j * npixels + k];
                     fwhm[iwave] = temp->fwhm[j];
                 }
                 iwave++;
             }
  
             ociaWavelength[iwave] = 936.646;
             fwhm[iwave] = 32.4702967;
             Lt[iwave * npixels + k] = (11.9269624 * temp->Lt[60 * npixels + k] +
                     10.78153 * temp->Lt[61 * npixels + k] +
                     9.76180425 * temp->Lt[62 * npixels + k]) / fwhm[iwave];
             iwave++;
  
             ociaWavelength[iwave] = 1238.750;
             fwhm[iwave] = 20.3179033;
             Lt[iwave * npixels + k] = (10.1448143 * temp->Lt[92 * npixels + k] +
                     10.173089 * temp->Lt[93 * npixels + k]) / fwhm[iwave];
             iwave++;
  
             Lt[iwave * npixels + k] = temp->Lt[109 * npixels + k];
             fwhm[iwave] = temp->fwhm[109];
             ociaWavelength[iwave] = temp->wave[109];
             iwave++;
  
             ociaWavelength[iwave] = 1641.530;
             fwhm[iwave] = 53.6097915;
             Lt[iwave * npixels + k] = (10.7227049 * temp->Lt[133 * npixels + k] +
                     10.7221941 * temp->Lt[134 * npixels + k] +
                     10.7218208 * temp->Lt[135 * npixels + k] +
                     10.721585 * temp->Lt[136 * npixels + k] +
                     10.7214867 * temp->Lt[137 * npixels + k]) / fwhm[iwave];
             iwave++;
             //iwave = interpLt(133, npixels, k, nbands_av, orcaWavelength[iwave], iwave, temp, acc,Lt);
             ociaWavelength[iwave] = 2126.795;
             fwhm[iwave] = 54.1850639;
             Lt[iwave * npixels + k] = (10.8702638 * temp->Lt[184 * npixels + k] +
                     10.8539551 * temp->Lt[185 * npixels + k] +
                     10.8373296 * temp->Lt[186 * npixels + k] +
                     10.8203873 * temp->Lt[187 * npixels + k] +
                     10.8031281 * temp->Lt[188 * npixels + k]) / fwhm[iwave];
             iwave++;
             //iwave = interpLt(184, npixels, k, nbands_av, orcaWavelength[iwave], iwave, temp, acc,Lt);
             ociaWavelength[iwave] = 2246.658;
             fwhm[iwave] = 53.0639591;
             Lt[iwave * npixels + k] = (10.6536474 * temp->Lt[196 * npixels + k] +
                     10.6335365 * temp->Lt[197 * npixels + k] +
                     10.6131087 * temp->Lt[198 * npixels + k] +
                     10.592364 * temp->Lt[199 * npixels + k] +
                     10.5713025 * temp->Lt[200 * npixels + k]) / fwhm[iwave];
             iwave++;
             //iwave = interpLt(196, npixels, k, nbands_av, orcaWavelength[iwave], iwave, temp, acc,Lt);
  
  
         } // pixel loop
  
         if (NBANDS != iwave) {
             printf("Oops, something's wrong iwave (%d) != nbands (%d) \n", (int) iwave, NBANDS);
             exit(-1);
         }
         // Write Lt's
         //    for (size_t iwave=0; iwave<nbands_ocia; iwave++ ) {
         int varID;
  
         //      printf("Lambda(%d) = %d\n",iwave+1,(int)(orcaWavelength[iwave]+.5));
         status = nc_inq_varid(grpID, "Lt", &varID);
         check_err(status, __LINE__, __FILE__);
  
         size_t start[3] = {0, iscan, 0};
         size_t count[3] = {nbands_ocia, 1, npixels};
         status = nc_put_vara_float(grpID, varID, start, count,
                 &Lt[0]);
         //      &Lt[iwave*npixels]);
         check_err(status, __LINE__, __FILE__);
         //    }
         alt[iscan] = temp->alt[iscan];
         utc[iscan] = temp->scantime;
         if (temp->scantime > 0) {
             scantime_last = temp->scantime;
             if (!first_good_scantime) {
                 first_good_scantime = 1;
                 scantime_first = temp->scantime;
                 //          printf("utc=%f\n",utc[iscan]);
                 //          strncpy(isodatetime_first,unix2isodate(scantime_first, 'G'),30);
                 //          printf("%s\n",isodatetime_first);
             }
             //      strncpy(isodatetime_last,unix2isodate(scantime_last, 'G'),30);
             //      printf("Last=%s\n",isodatetime_last);
         }
  
  
     } // scan loop
     int varID;
  
     if (scantime_first < scantime_last) {
         strncpy(isodatetime_first, unix2isodate(scantime_first, 'G'), 30);
         strncpy(isodatetime_last, unix2isodate(scantime_last, 'G'), 30);
     } else {
         strncpy(isodatetime_last, unix2isodate(scantime_first, 'G'), 30);
         strncpy(isodatetime_first, unix2isodate(scantime_last, 'G'), 30);
     }
     status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_start", strlen(isodatetime_first) + 1, isodatetime_first);
     if (status != NC_NOERR) {
         printf("-E- %s %d: %s for %s\n",
                 __FILE__, __LINE__, nc_strerror(status), "time_coverage_start");
         exit(1);
     }
     status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_end", strlen(isodatetime_last) + 1, isodatetime_last);
     if (status != NC_NOERR) {
         printf("-E- %s %d: %s for %s\n",
                 __FILE__, __LINE__, nc_strerror(status), "time_coverage_end");
         exit(1);
     }
  
     const char *navfieldsAVIRIS[] = {"to-sensor_azimuth", "to-sensor_zenith", "to-sun_azimuth", "to-sun_zenith", "lat", "lon"};
     const char *navfieldsOCIA[] = {"sensor_azimuth", "sensor_zenith", "solar_azimuth", "solar_zenith", "latitude", "longitude"};
  
     grpID = ociafile.getGid("geolocation_data");
  
     size_t startOCIA[2] = {0, 0};
     size_t countOCIA[2] = {nscans, npixels};
     size_t i;
  
     for (i = 0; i < 6; i++) {
         cout << "Copying " << navfieldsAVIRIS[i] << " to " <<
                 navfieldsOCIA[i] << endl;
  
         status = nc_inq_varid(grpID, navfieldsOCIA[i], &varID);
         check_err(status, __LINE__, __FILE__);
  
         switch (i) {
         case 0:
             status = nc_put_vara_float(grpID, varID, startOCIA, countOCIA, &temp->sena[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 1:
             status = nc_put_vara_float(grpID, varID, startOCIA, countOCIA, &temp->senz[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 2:
             status = nc_put_vara_float(grpID, varID, startOCIA, countOCIA, &temp->sola[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 3:
             status = nc_put_vara_float(grpID, varID, startOCIA, countOCIA, &temp->solz[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 4:
             status = nc_put_vara_double(grpID, varID, startOCIA, countOCIA, &temp->lat[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 5:
             status = nc_put_vara_double(grpID, varID, startOCIA, countOCIA, &temp->lon[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         default:
             break;
         }
     }
  
     gsl_interp_accel_free(acc);
  
     size_t startw[2] = {0, 0};
     size_t countw[2] = {nbands_ocia, 0};
  
     grpID = ociafile.getGid("sensor_band_parameters");
  
     status = nc_inq_varid(grpID, "wavelength", &varID);
  
     status = nc_put_vara_float(grpID, varID, startw, countw, ociaWavelength);
     check_err(status, __LINE__, __FILE__);
  
     status = nc_inq_varid(grpID, "fwhm", &varID);
  
     status = nc_put_vara_float(grpID, varID, startw, countw, fwhm);
     check_err(status, __LINE__, __FILE__);
  
     size_t starts[2] = {0, 0};
     size_t counts[2] = {nscans, 0};
  
     const char *nscanfieldsAVIRIS[] = {"utc", "alt"};
     const char *nscanfieldsOCIA[] = {"scan_start_time","altitude"};
  
     for (i = 0; i < nscans; i++) {
         //utc[i] = temp->scantime;
         //alt[i] = temp->alt[sline + i]*1000;
     }
     grpID = ociafile.getGid("scan_line_attributes");
     countOCIA[1] = 0;
     for (i = 0; i < 2; i++) {
         cout << "Copying " << nscanfieldsAVIRIS[i] << " to " <<
                 nscanfieldsOCIA[i] << endl;
  
         status = nc_inq_varid(grpID, nscanfieldsOCIA[i], &varID);
         check_err(status, __LINE__, __FILE__);
  
         switch (i) {
         case 0:
             status = nc_put_vara_double(grpID, varID, starts, counts, utc);
             check_err(status, __LINE__, __FILE__);
             break;
         case 1:
             status = nc_put_vara_float(grpID, varID, starts, counts, alt);
             check_err(status, __LINE__, __FILE__);
             break;
         default:
             break;
         }
     }
  
     //  outfile.close();
  
     status = close_aviris(temp);
  
     return 0;
 }

VERSION

#define VERSION

Definition: avirisbil2oci.cpp:33

int j

Definition: decode_rs.h:73

status

int status

Definition: l1_czcs_hdf.c:32

ncdfFile::getNcid

int getNcid()

Definition: cdl_utils.cpp:746

check_err

void check_err(const int stat, const int line, const char *file)

Definition: nc4utils.c:35

close_aviris

int close_aviris(aviris4ocia_t *data)

Definition: read_aviris.c:822

main

int main(int argc, char *argv[])

Definition: avirisbil2oci.cpp:98

NULL

#define NULL

Definition: decode_rs.h:63

interpLt

size_t interpLt(size_t startJ, size_t npixels, size_t k, size_t nbMax, float ociaWavelength, size_t iwave, aviris4ocia_t *data, gsl_interp_accel *acc, float *Lt)

Definition: avirisbil2oci.cpp:35

spline

subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)

Definition: phs.f:1348

nc4utils.h

counts

uint8 * counts

Definition: l1_czcs_hdf.c:30

open_aviris

aviris4ocia_t * open_aviris(char *filename, char *imgfile, char *navfile, char *gainfile, aviris4ocia_t **data)

Definition: read_aviris.c:175

checkAvProcessFile

int checkAvProcessFile(char *filename, char *hdrfile, char *imgfile, char *navfile, char *gainfile, int itemsize)

Definition: read_aviris.c:831

fwhm

float fwhm[NBANDS]

Definition: atrem_corl1.h:144

cdl_utils.h

hico.value_locate.x

list x

Definition: value_locate.py:64

color_dtdb.y

Definition: color_dtdb.py:430

oci_hdr_strip.argc

argc

Definition: oci_hdr_strip.py:7

ncdfFile::getGid

int getGid(const char *grpName)

Definition: cdl_utils.cpp:736

read_aviris

int read_aviris(aviris4ocia_t *data, int32_t recnum)

Definition: read_aviris.c:710

filename

char filename[FILENAME_MAX]

Definition: atrem_corl1.h:122

nc_gridutils.h

seadasutils.aquarius_utils.start

start

Definition: aquarius_utils.py:27

data

no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data

Definition: HISTORY.txt:356

SENSOR

Definition: DDProcess.h:81

CDLFILE

#define CDLFILE

Definition: avirisbil2oci.cpp:24

hico.interpolateFields2hicotimes.hdrfile

string hdrfile

Definition: interpolateFields2hicotimes.py:806

iscan

int32_t iscan

Definition: l1_hmodis_hdf.c:745

unix2isodate

char * unix2isodate(double dtime, char zone)

Definition: unix2isodate.c:10

l1cextract.eline

eline

Definition: l1cextract.py:392

aviris.h

strncpy

PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)

NBANDS

#define NBANDS

Definition: avirisbil2oci.cpp:25

ncdfFile::cdlCreateDim

int cdlCreateDim(char *l1_filename, char *cdl_filename, const char **dim_names, size_t *dim_size, size_t n_dims, size_t numScans)

Definition: cdl_utils.cpp:294

timeutils.h

FILL

#define FILL

Definition: avirisbil2oci.cpp:26

int i

Definition: decode_rs.h:71