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 <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
 #include <gsl/gsl_sort_double.h>
 #include <timeutils.h>
  
 #define FILL -32767
 #define BIP 0
 #define BIL 1
 #define SENSOR "aviris"
 #define CDLFILE "AVIRIS_Level-1B_Data_Structure.cdl"
 #define NBANDS 200 // number of output AVIRIS bands
  
  
 using namespace std;
  
  
 #define VERSION "1.00"
  
 size_t interpLt(size_t startJ, size_t npixels, size_t k, size_t nbMax,
         float ncdfWavelength, 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 <= ncdfWavelength && maxWave >= ncdfWavelength) break;
         if (minWave > ncdfWavelength || n < 3)
             i0--;
         if (maxWave < ncdfWavelength || 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 (ncdfWavelength > 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, ncdfWavelength, acc);
             Lt[iwave * npixels + k] = yi;
         } else {
             Lt[iwave * npixels + k] = -9999;
         }
         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;
     int interleave = BIL;
     double scantime_first, scantime_last;
     static int first_good_scantime = 0;
     char isodatetime_first[30], isodatetime_last[30];
  
     cout << "avirisbil2ncdf " << VERSION << " ("
             << __DATE__ << " " << __TIME__ << ")" << endl;
  
     if (argc < 3) {
         cout << endl <<
                 "avirisbil2ncdf input_AVIRIS_file output_AVIRIS_file <sline (default=0)> <eline (default=last row,0=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);
  
     aviris4ocia_t **data, *temp;
     float fwhm[NBANDS];
  
     if (!checkAvProcessFile(argv[1], hdrfile, imgfile, navfile, gainfile, FILENAME_MAX))
         strncpy(hdrfile, argv[1], FILENAME_MAX);
  
     data = (aviris4ocia_t **) malloc(sizeof (aviris4ocia_t *));
     temp = open_aviris(hdrfile, imgfile, navfile, gainfile, data);
     size_t ip, nscans, npixels = temp->npix;
  
     if (eline <= 0) eline = temp->nscans;
  
     nscans = eline - sline;
  
     int ncid, grpID;
     //char** dim_names, uint32_t* dim_size, size_t n_dims
     const char *dimNames[] = {"number_of_lines", "pixels_per_line", "number_of_bands"};
     size_t dimSize[] = {nscans, npixels, NBANDS};
     size_t n_dims = 3;
     // Create ORCA output file
     static ncdfFile ncdffile;
     ncdffile.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 = ncdffile.getNcid();
  
     size_t nbands_av = temp->numBands;
     size_t nbands_aviris = NBANDS;
  
     size_t n = 0, iwave = 0;
     istringstream istr;
     string str;
     static int *ibndx, *ipbndx;
     int ipb_av, ib, ibl1, ipb;
  
  
     // Allocate space for AVIRIS Lt (short) values
     gsl_interp_accel *acc = gsl_interp_accel_alloc();
  
     float *Lt = (float *) calloc(npixels*nbands_aviris, sizeof (float));
     float *ncdfWavelength = (float *) calloc(nbands_aviris, sizeof (float));
     double* utc = (double*) (calloc(nscans, sizeof (double)));
     size_t j;
  
     ipbndx = (int*) malloc(npixels * nbands_aviris * sizeof (int));
     if (!ipbndx) {
         fprintf(stderr, "-E- %s line %d: Out of Memory allocating buffer array in readl1_aviris.\n",
                 __FILE__, __LINE__);
         exit(1);
  
     }
     ibndx = (int*) malloc(nbands_aviris * sizeof (int));
     if (!ibndx) {
         fprintf(stderr, "-E- %s line %d: Out of Memory allocating buffer array in readl1_aviris.\n",
                 __FILE__, __LINE__);
         exit(1);
  
     }
  
     //use only the selected wavelengths
     ibl1 = 0;
     for (ib = 0; ib < 31; ib++) {
         ibndx[ibl1] = ib;
         ibl1++;
     }
     for (ib = 33; ib < 95; ib++) {
         ibndx[ibl1] = ib;
         ibl1++;
     }
     for (ib = 97; ib < 159; ib++) {
         ibndx[ibl1] = ib;
         ibl1++;
     }
     for (ib = 161; ib < 206; ib++) {
         ibndx[ibl1] = ib;
         ibl1++;
     }
  
     // create an index array to map the Lt's to the wavelength indexes actually used
     for (ip = 0; ip < npixels; ip++) {
         ibl1 = 0;
         for (ib = 0; ib < 31; ib++) {
             ipb = ip * nbands_aviris + ibl1;
             switch (interleave) {
             case BIP:
                 ipb_av = ip * nbands_av + ib;
                 break;
             case BIL:
                 ipb_av = ib * npixels + ip;
                 break;
             }
             ipbndx[ipb] = ipb_av;
             ibl1++;
         }
         for (ib = 33; ib < 95; ib++) {
             ipb = ip * nbands_aviris + ibl1;
             switch (interleave) {
             case BIP:
                 ipb_av = ip * nbands_av + ib;
                 break;
             case BIL:
                 ipb_av = ib * npixels + ip;
                 break;
             }
             ipbndx[ipb] = ipb_av;
             ibl1++;
         }
         for (ib = 97; ib < 159; ib++) {
             ipb = ip * nbands_aviris + ibl1;
             switch (interleave) {
             case BIP:
                 ipb_av = ip * nbands_av + ib;
                 break;
             case BIL:
                 ipb_av = ib * npixels + ip;
                 break;
             }
             ipbndx[ipb] = ipb_av;
             ibl1++;
         }
         for (ib = 161; ib < 206; ib++) {
             ipb = ip * nbands_aviris + ibl1;
             switch (interleave) {
             case BIP:
                 ipb_av = ip * nbands_av + ib;
                 break;
             case BIL:
                 ipb_av = ib * npixels + ip;
                 break;
             }
             ipbndx[ipb] = ipb_av;
             ibl1++;
         }
     }
  
     grpID = ncdffile.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_aviris; 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 < nbands_aviris; j++) {
                 ipb = k * nbands_aviris + j;
                 if (temp->Lt[ipbndx[ipb]] > 0) {
                     Lt[iwave * npixels + k] = temp->Lt[ipbndx[ipb]];
                 } else {
                     Lt[iwave * npixels + k] = FILL;
                 }
                 ncdfWavelength[iwave] = temp->wave[ibndx[j]];
                 fwhm[iwave] = temp->fwhm[ibndx[j]];
  
                 iwave++;
             }
  
  
         } // 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_aviris; 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_aviris, 1, npixels};
         status = nc_put_vara_float(grpID, varID, start, count,
                 &Lt[0]);
         check_err(status, __LINE__, __FILE__);
         //    }
  
         utc[iscan] = temp->scantime;
         //if (utc[iscan]>0) printf("utc=%f\n",utc[iscan]);
  
         if (temp->scantime > 0) {
             scantime_last = temp->scantime;
             if (!first_good_scantime) {
                 first_good_scantime = 1;
                 scantime_first = temp->scantime;
             }
         }
  
     } // scan loop
  
     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 *navfieldsNCDF[] = {"sena", "senz", "sola", "solz", "lat", "lon"};
  
     grpID = ncdffile.getGid("navigation_data");
     int varID;
  
     size_t startNCDF[2] = {0, 0};
     size_t countNCDF[2] = {nscans, npixels};
     size_t i;
  
     for (i = 0; i < 6; i++) {
         cout << "Copying " << navfieldsAVIRIS[i] << " to " <<
                 navfieldsNCDF[i] << endl;
  
         status = nc_inq_varid(grpID, navfieldsNCDF[i], &varID);
         check_err(status, __LINE__, __FILE__);
  
         switch (i) {
         case 0:
             status = nc_put_vara_float(grpID, varID, startNCDF, countNCDF, &temp->sena[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 1:
             status = nc_put_vara_float(grpID, varID, startNCDF, countNCDF, &temp->senz[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 2:
             status = nc_put_vara_float(grpID, varID, startNCDF, countNCDF, &temp->sola[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 3:
             status = nc_put_vara_float(grpID, varID, startNCDF, countNCDF, &temp->solz[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 4:
             status = nc_put_vara_double(grpID, varID, startNCDF, countNCDF, &temp->lat[sline * npixels]);
             check_err(status, __LINE__, __FILE__);
             break;
         case 5:
             status = nc_put_vara_double(grpID, varID, startNCDF, countNCDF, &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_aviris, 0};
  
     grpID = ncdffile.getGid("sensor_band_parameters");
  
     status = nc_inq_varid(grpID, "wavelength", &varID);
  
     status = nc_put_vara_float(grpID, varID, startw, countw, ncdfWavelength);
     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", "alt"};
     const char *nscanfieldsNCDF[] = {"scan_start_time", "altitude", "range"};
     short* alt = (short*) (calloc(nscans, sizeof (short)));
  
     grpID = ncdffile.getGid("scan_line_attributes");
     countNCDF[1] = 0;
     for (i = 0; i < 3; i++) {
         cout << "Copying " << nscanfieldsAVIRIS[i] << " to " <<
                 nscanfieldsNCDF[i] << endl;
  
         status = nc_inq_varid(grpID, nscanfieldsNCDF[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:
             for (i = 0; i < nscans; i++) alt[i] = temp->alt[sline + i];
             status = nc_put_vara_short(grpID, varID, starts, counts, alt);
             check_err(status, __LINE__, __FILE__);
             break;
         case 2:
             for (i = 0; i < nscans; i++) alt[i] = temp->alt[sline + i]*1000;
             status = nc_put_vara_short(grpID, varID, starts, counts, alt);
             check_err(status, __LINE__, __FILE__);
             break;
         default:
             break;
         }
     }
  
     //  outfile.close();
  
     status = close_aviris(temp);
  
     return 0;
 }
  

l1mapgen.stderr

stderr

Definition: l1mapgen.py:204

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

CDLFILE

#define CDLFILE

Definition: avirisbil2ncdf.cpp:27

NULL

#define NULL

Definition: decode_rs.h:63

main

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

Definition: avirisbil2ncdf.cpp:99

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

FILL

#define FILL

Definition: avirisbil2ncdf.cpp:23

cdl_utils.h

BIL

#define BIL

Definition: avirisbil2ncdf.cpp:25

hico.value_locate.x

list x

Definition: value_locate.py:64

color_dtdb.y

Definition: color_dtdb.py:430

BIP

#define BIP

Definition: avirisbil2ncdf.cpp:24

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

VERSION

#define VERSION

Definition: avirisbil2ncdf.cpp:34

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

interpLt

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

Definition: avirisbil2ncdf.cpp:36

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)

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

int i

Definition: decode_rs.h:71

strcpy

How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")

ncdfFile

Definition: cdl_utils.h:9

int k

Definition: decode_rs.h:73