l1c_str.cpp

Go to the documentation of this file.

//
//  l1c_str.cpp
//  
//
//  Created by Martin Montes on 7/17/2021
#include "l1c_str.h"  
#include <iostream>
#include <string>
#include "l1c_filehandle.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <netcdf.h>
#include <nc4utils.h>
#include "libnav.h"
#include <timeutils.h>
#include <genutils.h>
 
 
#include "allocate4d.h"
#include <allocate3d.h>
#include <allocate2d.h>
 
//static float *Fobar; // reflectance to radiance conversion factors
//static int extract_pixel_start = 0;
 
static short *tmpShort;
 
// whole file stuff
static size_t expected_num_blue_bands = 60; //120;
static size_t expected_num_red_bands = 60; //120;
static size_t expected_num_SWIR_bands = 9;
 
static size_t nviews;
//static size_t nbands;
static size_t num_scans, num_pixels;
static size_t num_blue_bands, num_red_bands, num_SWIR_bands;
//static size_t tot_num_bands = 123; //238;
static int ncid_L1B;
 
// scan line attributes
static double *scan_time; // seconds of day
static int32_t scan_time_year, scan_time_month, scan_time_day;
//static int32_t scan_time_hour,scan_time_min;
//static double scan_time_secs;
 
//static uint_8 *scan_quality;
 
// geolocation data
static int geolocationGrp;
static int lonId, latId, senzId, senaId, solzId, solaId;
static float latFillValue = -999.0;
static float lonFillValue = -999.0;
static short senzFillValue = -32768;
static float senzScale = 0.01;
static float senzOffset = 0.0;
static short senaFillValue = -32768;
static float senaScale = 0.01;
static float senaOffset = 0.0;
static short solzFillValue = -32768;
static float solzScale = 0.01;
static float solzOffset = 0.0;
static short solaFillValue = -32768;
static float solaScale = 0.01;
static float solaOffset = 0.0;
 
// Observation data
static int observationGrp;
static int Lt_blueId, Lt_redId, Lt_SWIRId;
//static float **Lt_blue;                //[num_blue_bands][num_pixels], This scan
//static float **Lt_red;                 //[num_red_bands][num_pixels], This scan
//static float **Lt_SWIR;                //[num_SWIR_bands][num_pixels], This scan
 
//static float ****Lt_all;
 
static float Lt_blueFillValue = -999.0;
static float Lt_redFillValue = -999.0;
static float Lt_SWIRFillValue = -999.0;
 
 
//Navigation data
//static int navGrp;
//static int ovId,opId,attID;
//static float *orbv;//for each scan we have xpixels dim
//static float *orbp;
//static float *attang;
//static double *timeArr;
 
static int status,dimid;
 
using namespace std;
 
 
 
 
 
 
 
 
 
 
namespace l1c {
 
     l1c_str::l1c_str() :  att_ang{-999.0,-999.0,-999.0}, orb_pos{-999.0,-999.0,-999.0}, orb_vel{-999.0,-999.0,-999.0} {
  
    //init constructor
       //global attributes
     npix=-1;
     iscan=0;
     nscan=-1;
     //scan line attributes--
     ev_mid_time=nullptr;
     scan_quality_flag=nullptr;
     spix=-1;
     epix=-1;
     dpix=-1;
     
     // structure--pointers to data arrays
     Lt=nullptr; //dim depends between sensors--OCI is bands x pixels --
     Lt_blue=nullptr;                //[num_views][num_pol][num_blue_bands][num_pixels]
     Lt_red=nullptr;
     Lt_SWIR=nullptr;
     
     //sensor/sun geometry
     senz=nullptr;
     sena=nullptr;
     solz=nullptr;
     sola=nullptr;
     delphi=nullptr;
     scattang=nullptr;     
 
     //geolocation--
     senazpix=nullptr;
     latpix=nullptr;
     lonpix=nullptr;
     latpix2=nullptr;
     lonpix2=nullptr;
     latnad=nullptr;
     lonnad=nullptr;
     lonershift=nullptr;
     terr_height=nullptr;
     cloud_height=nullptr;
     
     
     l1cfile=nullptr;
   
    //ancill info?       
}
 
 
l1c_str::~l1c_str() {
}
 
int32_t l1c_str::openl1b_ocis_l1c(l1c_str *l1cstr,l1c_filehandle *l1cfile,int16_t *file_id){
 
    std::string str;
    const char *ptstr;
 
 // Open the netcdf4 input file    
    str=l1cfile->l1b_name; 
    ptstr=str.c_str();
    printf("Opening OCIS L1B file\n");
    cout<<str<<endl;
 
    status = nc_open(ptstr, NC_NOWRITE, &ncid_L1B);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- %s line %d: nc_open(%s) failed.\n",
                __FILE__, __LINE__, ptstr);
        return (1);
    }
 
    // num_scans
    status = nc_inq_dimid(ncid_L1B, "number_of_scans", &dimid);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- Error reading number_of_scans.\n");
        exit(EXIT_FAILURE);
    }
    nc_inq_dimlen(ncid_L1B, dimid, &num_scans);
    l1cstr->nscan=num_scans;
 
    // num_pixels
    status = nc_inq_dimid(ncid_L1B, "ccd_pixels", &dimid);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- Error reading num_pixels.\n");
        exit(EXIT_FAILURE);
    }
    nc_inq_dimlen(ncid_L1B, dimid, &num_pixels);
    l1cstr->npix=num_pixels;
 
    // num_blue_bands
    status = nc_inq_dimid(ncid_L1B, "blue_bands", &dimid);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- Error reading num_blue_bands.\n");
        exit(EXIT_FAILURE);
    }
    nc_inq_dimlen(ncid_L1B, dimid, &num_blue_bands);
    if(num_blue_bands < expected_num_blue_bands) {
        fprintf(stderr, "-E- Not enough blue bands, expecting %d, found %d.\n",
                (int)expected_num_blue_bands, (int)num_blue_bands);
        exit(EXIT_FAILURE);
    }
 
    // num_red_bands
    status = nc_inq_dimid(ncid_L1B, "red_bands", &dimid);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- Error reading num_red_bands.\n");
        exit(EXIT_FAILURE);
    };
    nc_inq_dimlen(ncid_L1B, dimid, &num_red_bands);
    if(num_red_bands < expected_num_red_bands) {
        fprintf(stderr, "-E- Not enough red bands, expecting %d, found %d.\n",
                (int)expected_num_red_bands, (int)num_red_bands);
        exit(EXIT_FAILURE);
    }
 
    // num_SWIR_bands
    status = nc_inq_dimid(ncid_L1B, "SWIR_bands", &dimid);
    if (status != NC_NOERR) {
        fprintf(stderr, "-E- Error reading num_SWIR_bands.\n");
        exit(EXIT_FAILURE);
    };
    nc_inq_dimlen(ncid_L1B, dimid, &num_SWIR_bands);
    if(num_SWIR_bands < expected_num_SWIR_bands) {
        fprintf(stderr, "-E- Not enough SWIR bands, expecting %d, found %d.\n",
                (int)expected_num_SWIR_bands, (int)num_SWIR_bands);
        exit(EXIT_FAILURE);
    }
 
//    if (want_verbose) {
        printf("OCI L1B Npix  :%d Nlines:%d\n", (int)num_pixels, (int)num_scans);
//    } // want_verbose
 
 
// allocate all of the data
    tmpShort = (short*) calloc(num_pixels, sizeof(short));
    scan_time = (double*) calloc(num_scans, sizeof(double));
 
    l1cstr->senazpix=(float*)calloc(num_pixels,sizeof(float));
    l1cstr->latpix=(float*)calloc(num_pixels,sizeof(float));
    l1cstr->lonpix=(float*)calloc(num_pixels,sizeof(float));
    l1cstr->latpix2=(float*)calloc(num_pixels,sizeof(float));
    l1cstr->lonpix2=(float*)calloc(num_pixels,sizeof(float));
 
 
    l1cstr->Lt_blue = allocate2d_float(num_blue_bands, num_pixels);
    l1cstr->Lt_red = allocate2d_float(num_red_bands, num_pixels);
    l1cstr->Lt_SWIR = allocate2d_float(num_SWIR_bands, num_pixels);
 
    // Get group id from L1B file for GROUP scan_line_attributes.
    int groupid;
    if ((nc_inq_grp_ncid(ncid_L1B, "scan_line_attributes", &groupid)) == NC_NOERR) {
    } else {
        fprintf(stderr, "-E- Error finding scan_line_attributes.\n");
        exit(EXIT_FAILURE);
    }
    int varid;
    double scan_timeFillValue = -999.9;
    status = nc_inq_varid(groupid, "time", &varid);
    if(status == NC_NOERR) {
        status = nc_inq_var_fill(groupid, varid, NULL, &scan_timeFillValue);
        check_err(status, __LINE__, __FILE__);
        status = nc_get_var_double(groupid, varid, scan_time);
        check_err(status, __LINE__, __FILE__);
        status = nc_get_att_int(groupid, varid, "year", &scan_time_year);
        check_err(status, __LINE__, __FILE__);
        status = nc_get_att_int(groupid, varid, "month", &scan_time_month);
        check_err(status, __LINE__, __FILE__);
        status = nc_get_att_int(groupid, varid, "day", &scan_time_day);
        check_err(status, __LINE__, __FILE__);
    } else {
        status = nc_inq_varid(groupid, "ev_mid_time", &varid);
        check_err(status, __LINE__, __FILE__);
        status = nc_inq_var_fill(groupid, varid, NULL, &scan_timeFillValue);
        check_err(status, __LINE__, __FILE__);
        status = nc_get_var_double(groupid, varid, scan_time);
        check_err(status, __LINE__, __FILE__);
 
    // get start time
        size_t att_len;
    status = nc_inq_attlen(ncid_L1B, NC_GLOBAL, "time_coverage_start", &att_len);
    check_err(status, __LINE__, __FILE__);
 
        // allocate required space before retrieving values 
        char* time_str = (char *) malloc(att_len + 1); // + 1 for trailing null 
 
        // get attribute values 
        status = nc_get_att_text(ncid_L1B, NC_GLOBAL, "time_coverage_start", time_str);
    check_err(status, __LINE__, __FILE__);
        time_str[att_len] = '\0';
 
        double start_time = isodate2unix(time_str);
        int16_t syear, smon, sday;
        double secs;
        unix2ymds(start_time, &syear, &smon, &sday, &secs);
        scan_time_year = syear;
        scan_time_month = smon;
        scan_time_day = sday;
 
    }
 
    for(size_t i=0; i<num_scans; i++) {
        if(scan_time[i] == scan_timeFillValue)
            scan_time[i] = BAD_FLT;
    }
 
    // read the orbit#
    int orbit_number;
    status = nc_get_att_int(ncid_L1B, NC_GLOBAL, "orbit_number", &orbit_number);
    check_err(status, __LINE__, __FILE__);
 
    // Setup geofile pointers
    status = nc_inq_grp_ncid(ncid_L1B, "geolocation_data", &geolocationGrp);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_varid(geolocationGrp, "longitude", &lonId);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_var_fill(geolocationGrp, lonId, NULL, &lonFillValue);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_varid(geolocationGrp, "latitude", &latId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(geolocationGrp, latId, NULL, &latFillValue);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(geolocationGrp, "sensor_zenith", &senzId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(geolocationGrp, senzId, NULL, &senzFillValue);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, senzId, "scale_factor", &senzScale);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, senzId, "add_offset", &senzOffset);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(geolocationGrp, "sensor_azimuth", &senaId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(geolocationGrp, senaId, NULL, &senaFillValue);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, senaId, "scale_factor", &senaScale);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, senaId, "add_offset", &senaOffset);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(geolocationGrp, "solar_zenith", &solzId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(geolocationGrp, solzId, NULL, &solzFillValue);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, solzId, "scale_factor", &solzScale);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, solzId, "add_offset", &solzOffset);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(geolocationGrp, "solar_azimuth", &solaId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(geolocationGrp, solaId, NULL, &solaFillValue);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, solaId, "scale_factor", &solaScale);
    check_err(status, __LINE__, __FILE__);
    status = nc_get_att_float(geolocationGrp, solaId, "add_offset", &solaOffset);
    check_err(status, __LINE__, __FILE__);
 
 
    // get IDs for the observations
    status = nc_inq_grp_ncid(ncid_L1B, "observation_data", &observationGrp);
    check_err(status, __LINE__, __FILE__);
 
    // Get varids for each of the Lt_*
    status = nc_inq_varid(observationGrp, "Lt_blue", &Lt_blueId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(observationGrp, Lt_blueId, NULL, &Lt_blueFillValue);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(observationGrp, "Lt_red", &Lt_redId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(observationGrp, Lt_redId, NULL, &Lt_redFillValue);
    check_err(status, __LINE__, __FILE__);
 
    status = nc_inq_varid(observationGrp, "Lt_SWIR", &Lt_SWIRId);
    check_err(status, __LINE__, __FILE__);
    status = nc_inq_var_fill(observationGrp, Lt_SWIRId, NULL, &Lt_SWIRFillValue);
    check_err(status, __LINE__, __FILE__);
 
 
    l1cfile->sd_id = ncid_L1B;
    l1cfile->nband_blue=num_blue_bands;
    l1cfile->nband_red=num_red_bands;
    l1cfile->nband_swir=num_SWIR_bands;
    l1cfile->nbands=num_blue_bands+num_red_bands+num_SWIR_bands;
    nviews=2;
    l1cfile->n_views=nviews;
 
    cout<<"file->nbands_b ="<<l1cfile->nband_blue<<endl;
    cout<<"file->nbands_r = "<<l1cfile->nband_red<<endl;
    cout<<"file->nbands_swir = "<<l1cfile->nband_swir<<endl;
 
    l1cfile->npix = num_pixels;
    l1cfile->nadpix=(num_pixels-1)/2;//nadir pixel index
    l1cfile->nscan = num_scans;
    l1cfile->ndets = 1;
    l1cfile->terrain_corrected = 1; // presumed.
    l1cfile->orbit_number = orbit_number;
 
//    rdsensorinfo(l1cfile->sensorID, input->evalmask,
//            "Fobar", (void **) &Fobar);
 
//    if (want_verbose)
//        printf("file->nbands = %d\n", (int) file->nbands);
 
return 0;
}
 
 
int32_t l1c_str::readl1b_ocis_l1c(l1c_str *l1cstr,l1c_filehandle *l1cfile, int16_t *file_id,int32_t sline){
  size_t oneline;
  size_t start[] = {0,0},start2[]={0,0},start3[] = {0,0,0};
  size_t count[] = {1,1},count2[] = {1,1,1};
  string str;
 
 
  str=l1cfile->l1b_name;
 
  printf("Reading OCIS L1B file\n");
  cout<<str<<endl;
 
   //assigning mem for Lt of different bands **********************
 
//GEOLOCATION
  l1cstr->iscan=sline;
 
  start[0] = sline;
  start[1] = 0;
  count[0] = 1;
  count[1] = num_pixels;               // 1 line at a time
 
  status = nc_get_vara_float(geolocationGrp, senaId, start, count, l1cstr->senazpix);//sensor azimuth -1800 to 1800 degrees, so must be divdied by 100
  status = nc_get_vara_float(geolocationGrp, latId, start, count, l1cstr->latpix);
  check_err(status, __LINE__, __FILE__);
  status = nc_get_vara_float(geolocationGrp, lonId, start, count, l1cstr->lonpix);
  check_err(status, __LINE__, __FILE__);
 
 
  oneline=sline;
 
  if(oneline<num_scans-1)
    start2[0] = sline+1;
  else
    start2[0] = sline;
 
  start2[1] = 0;
 
 
  status = nc_get_vara_float(geolocationGrp, latId, start2, count, l1cstr->latpix2);
  check_err(status, __LINE__, __FILE__);
  status = nc_get_vara_float(geolocationGrp, lonId, start2, count, l1cstr->lonpix2);
  check_err(status, __LINE__, __FILE__);
 
 
//RADIANCES
  start3[0] = 0;
  start3[1] = sline;
  start3[2] = 0;
  count2[0] = num_blue_bands;
  count2[1] = 1;               // 1 line at a time
  count2[2] = num_pixels;
 
  status = nc_get_vara_float(observationGrp, Lt_blueId, start3, count2, l1cstr->Lt_blue[0]);//these are 2-D arrays------------ #bands x #pixels
  check_err(status, __LINE__, __FILE__);
 
  count2[0] = num_red_bands;
  status = nc_get_vara_float(observationGrp, Lt_redId, start3, count2, l1cstr->Lt_red[0]);
  check_err(status, __LINE__, __FILE__);
 
  count2[0] = num_SWIR_bands;
  status = nc_get_vara_float(observationGrp, Lt_SWIRId, start3, count2, l1cstr->Lt_SWIR[0]);
  check_err(status, __LINE__, __FILE__);
 
return 0;
}
 
 
int32_t l1c_str::closel1b_ocis_l1c(l1c_str *l1cstr,l1c_filehandle *l1cfile){
     string str;
 
     str=l1cfile->l1b_name;
     printf("Closing ocis l1b file\n");
     cout<<str<<endl;
        status = nc_close(l1cfile->sd_id);
          check_err(status, __LINE__, __FILE__);
 
 
    // Free memory
    // From openl1b_ocis
     if (l1cstr->senazpix != nullptr)
       delete [] (l1cstr->senazpix);
     if (l1cstr->latpix != nullptr)
       delete [] (l1cstr->latpix);
     if (l1cstr->lonpix != nullptr) 
        delete [] (l1cstr->lonpix);  
     if (l1cstr->latpix2 != nullptr)
       delete [] (l1cstr->latpix2);
     if (l1cstr->lonpix2 != nullptr)
        delete [] (l1cstr->lonpix2);
 
    if (tmpShort) free(tmpShort);
    if (scan_time) free(scan_time);
    if (l1cstr->Lt_blue) free2d_float(l1cstr->Lt_blue);
    if (l1cstr->Lt_red) free2d_float(l1cstr->Lt_red);
    if (l1cstr->Lt_SWIR) free2d_float(l1cstr->Lt_SWIR);
 
    l1cstr->latpix=nullptr;
    l1cstr->latpix2=nullptr;
    l1cstr->lonpix=nullptr;
    l1cstr->lonpix2=nullptr;
 
    l1cstr->Lt_blue=nullptr;
    l1cstr->Lt_red=nullptr;
    l1cstr->Lt_SWIR=nullptr;
 
return 0;
}
 
 
} // namespace l1c