/******************************************************************************
 *
 * Project:  NOAA Polar Orbiter Level 1b Dataset Reader (AVHRR)
 * Purpose:  Can read NOAA-9(F)-NOAA-17(M) AVHRR datasets
 * Author:   Andrey Kiselev, dron@ak4719.spb.edu
 *
 * Some format info at: http://www.sat.dundee.ac.uk/noaa1b.html
 *
 ******************************************************************************
 * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * ----------------------------------------------------------------------------
 * Lagrange interpolation suitable for NOAA level 1B file formats.
 * Submitted by Andrew Brooks <arb@sat.dundee.ac.uk>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_string.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
#include "ogr_srs_api.h"

#include <algorithm>

typedef enum
{                     // File formats
    L1B_NONE,         // Not a L1B format
    L1B_NOAA9,        // NOAA-9/14
    L1B_NOAA15,       // NOAA-15/METOP-2
    L1B_NOAA15_NOHDR  // NOAA-15/METOP-2 without ARS header
} L1BFileFormat;

typedef enum
{            // Spacecrafts:
    TIROSN,  // TIROS-N
             // NOAA are given a letter before launch and a number after launch
    NOAA6,   // NOAA-6(A)
    NOAAB,   // NOAA-B
    NOAA7,   // NOAA-7(C)
    NOAA8,   // NOAA-8(E)
    NOAA9_UNKNOWN,  // Some NOAA-18 and NOAA-19 HRPT are recognized like that...
    NOAA9,          // NOAA-9(F)
    NOAA10,         // NOAA-10(G)
    NOAA11,         // NOAA-11(H)
    NOAA12,         // NOAA-12(D)
    NOAA13,         // NOAA-13(I)
    NOAA14,         // NOAA-14(J)
    NOAA15,         // NOAA-15(K)
    NOAA16,         // NOAA-16(L)
    NOAA17,         // NOAA-17(M)
    NOAA18,         // NOAA-18(N)
    NOAA19,         // NOAA-19(N')
    // MetOp are given a number before launch and a letter after launch
    METOP2,  // METOP-A(2)
    METOP1,  // METOP-B(1)
    METOP3,  // METOP-C(3)
} L1BSpaceCraftdID;

typedef enum
{  // Product types
    HRPT,
    LAC,
    GAC,
    FRAC
} L1BProductType;

typedef enum
{  // Data format
    PACKED10BIT,
    UNPACKED8BIT,
    UNPACKED16BIT
} L1BDataFormat;

typedef enum
{        // Receiving stations names:
    DU,  // Dundee, Scotland, UK
    GC,  // Fairbanks, Alaska, USA (formerly Gilmore Creek)
    HO,  // Honolulu, Hawaii, USA
    MO,  // Monterey, California, USA
    WE,  // Western Europe CDA, Lannion, France
    SO,  // SOCC (Satellite Operations Control Center), Suitland, Maryland, USA
    WI,  // Wallops Island, Virginia, USA
    SV,  // Svalbard, Norway
    UNKNOWN_STATION
} L1BReceivingStation;

typedef enum
{         // Data processing centers:
    CMS,  // Centre de Meteorologie Spatiale - Lannion, France
    DSS,  // Dundee Satellite Receiving Station - Dundee, Scotland, UK
    NSS,  // NOAA/NESDIS - Suitland, Maryland, USA
    UKM,  // United Kingdom Meteorological Office - Bracknell, England, UK
    UNKNOWN_CENTER
} L1BProcessingCenter;

typedef enum
{  // AVHRR Earth location indication
    ASCEND,
    DESCEND
} L1BAscendOrDescend;

/************************************************************************/
/*                      AVHRR band widths                               */
/************************************************************************/

static const char *const apszBandDesc[] = {
    // NOAA-7 -- METOP-2 channels
    "AVHRR Channel 1:  0.58  micrometers -- 0.68 micrometers",
    "AVHRR Channel 2:  0.725 micrometers -- 1.10 micrometers",
    "AVHRR Channel 3:  3.55  micrometers -- 3.93 micrometers",
    "AVHRR Channel 4:  10.3  micrometers -- 11.3 micrometers",
    "AVHRR Channel 5:  11.5  micrometers -- 12.5 micrometers",  // not in
                                                                // NOAA-6,-8,-10
                                                                // NOAA-13
    "AVHRR Channel 5:  11.4  micrometers -- 12.4 micrometers",
    // NOAA-15 -- METOP-2
    "AVHRR Channel 3A: 1.58  micrometers -- 1.64 micrometers",
    "AVHRR Channel 3B: 3.55  micrometers -- 3.93 micrometers"};

/************************************************************************/
/*      L1B file format related constants                               */
/************************************************************************/

// Length of the string containing dataset name
#define L1B_DATASET_NAME_SIZE 42
#define L1B_NOAA9_HEADER_SIZE 122   // Terabit memory (TBM) header length
#define L1B_NOAA9_HDR_NAME_OFF 30   // Dataset name offset
#define L1B_NOAA9_HDR_SRC_OFF 70    // Receiving station name offset
#define L1B_NOAA9_HDR_CHAN_OFF 97   // Selected channels map offset
#define L1B_NOAA9_HDR_CHAN_SIZE 20  // Length of selected channels map
#define L1B_NOAA9_HDR_WORD_OFF 117  // Sensor data word size offset

// Archive Retrieval System (ARS) header
#define L1B_NOAA15_HEADER_SIZE 512
#define L1B_NOAA15_HDR_CHAN_OFF 97   // Selected channels map offset
#define L1B_NOAA15_HDR_CHAN_SIZE 20  // Length of selected channels map
#define L1B_NOAA15_HDR_WORD_OFF 117  // Sensor data word size offset

// Length of header record filled with the data
#define L1B_NOAA9_HDR_REC_SIZE 146
#define L1B_NOAA9_HDR_REC_ID_OFF 0      // Spacecraft ID offset
#define L1B_NOAA9_HDR_REC_PROD_OFF 1    // Data type offset
#define L1B_NOAA9_HDR_REC_DSTAT_OFF 34  // DACS status offset

/* See
 * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83132-2.htm */
// Length of header record filled with the data
#define L1B_NOAA15_HDR_REC_SIZE 992
#define L1B_NOAA15_HDR_REC_SITE_OFF 0  // Dataset creation site ID offset
#define L1B_NOAA15_HDR_REC_FORMAT_VERSION_OFF                                  \
    4  // NOAA Level 1b Format Version Number
#define L1B_NOAA15_HDR_REC_FORMAT_VERSION_YEAR_OFF                             \
    6  // Level 1b Format Version Year (e.g., 1999)
#define L1B_NOAA15_HDR_REC_FORMAT_VERSION_DAY_OFF                              \
    8  // Level 1b Format Version Day of Year (e.g., 365)
#define L1B_NOAA15_HDR_REC_LOGICAL_REC_LENGTH_OFF                              \
    10  // Logical Record Length of source Level 1b data set prior to processing
#define L1B_NOAA15_HDR_REC_BLOCK_SIZE_OFF                                      \
    12  // Block Size of source Level 1b data set prior to processing
#define L1B_NOAA15_HDR_REC_HDR_REC_COUNT_OFF                                   \
    14  // Count of Header Records in this Data Set
#define L1B_NOAA15_HDR_REC_NAME_OFF 22   // Dataset name
#define L1B_NOAA15_HDR_REC_ID_OFF 72     // Spacecraft ID offset
#define L1B_NOAA15_HDR_REC_PROD_OFF 76   // Data type offset
#define L1B_NOAA15_HDR_REC_STAT_OFF 116  // Instrument status offset
#define L1B_NOAA15_HDR_REC_DATA_RECORD_COUNT_OFF 128
#define L1B_NOAA15_HDR_REC_CALIBRATED_SCANLINE_COUNT_OFF 130
#define L1B_NOAA15_HDR_REC_MISSING_SCANLINE_COUNT_OFF 132
#define L1B_NOAA15_HDR_REC_SRC_OFF 154  // Receiving station name offset
#define L1B_NOAA15_HDR_REC_ELLIPSOID_OFF 328

/* This only apply if L1B_HIGH_GCP_DENSITY is explicitly set to NO */
/* otherwise we will report more GCPs */
#define DESIRED_GCPS_PER_LINE 11
#define DESIRED_LINES_OF_GCPS 20

// Fixed values used to scale GCPs coordinates in AVHRR records
#define L1B_NOAA9_GCP_SCALE 128.0
#define L1B_NOAA15_GCP_SCALE 10000.0

/************************************************************************/
/* ==================================================================== */
/*                      TimeCode (helper class)                         */
/* ==================================================================== */
/************************************************************************/

#define L1B_TIMECODE_LENGTH 100

class TimeCode
{
    long lYear;
    long lDay;
    long lMillisecond;
    char szString[L1B_TIMECODE_LENGTH];

  public:
    TimeCode() : lYear(0), lDay(0), lMillisecond(0)
    {
        memset(szString, 0, sizeof(szString));
    }

    void SetYear(long year)
    {
        lYear = year;
    }

    void SetDay(long day)
    {
        lDay = day;
    }

    void SetMillisecond(long millisecond)
    {
        lMillisecond = millisecond;
    }

    long GetYear() const
    {
        return lYear;
    }

    long GetDay() const
    {
        return lDay;
    }

    long GetMillisecond() const
    {
        return lMillisecond;
    }

    char *PrintTime()
    {
        snprintf(szString, L1B_TIMECODE_LENGTH,
                 "year: %ld, day: %ld, millisecond: %ld", lYear, lDay,
                 lMillisecond);
        return szString;
    }
};

#undef L1B_TIMECODE_LENGTH

/************************************************************************/
/* ==================================================================== */
/*                              L1BDataset                              */
/* ==================================================================== */
/************************************************************************/
class L1BGeolocDataset;
class L1BGeolocRasterBand;
class L1BSolarZenithAnglesDataset;
class L1BSolarZenithAnglesRasterBand;
class L1BNOAA15AnglesDataset;
class L1BNOAA15AnglesRasterBand;
class L1BCloudsDataset;
class L1BCloudsRasterBand;

class L1BDataset final : public GDALPamDataset
{
    friend class L1BRasterBand;
    friend class L1BMaskBand;
    friend class L1BGeolocDataset;
    friend class L1BGeolocRasterBand;
    friend class L1BSolarZenithAnglesDataset;
    friend class L1BSolarZenithAnglesRasterBand;
    friend class L1BNOAA15AnglesDataset;
    friend class L1BNOAA15AnglesRasterBand;
    friend class L1BCloudsDataset;
    friend class L1BCloudsRasterBand;

    // char        pszRevolution[6]; // Five-digit number identifying spacecraft
    // revolution
    L1BReceivingStation eSource;      // Source of data (receiving station name)
    L1BProcessingCenter eProcCenter;  // Data processing center
    TimeCode sStartTime;
    TimeCode sStopTime;

    int bHighGCPDensityStrategy;
    GDAL_GCP *pasGCPList;
    int nGCPCount;
    int iGCPOffset;
    int iGCPCodeOffset;
    int iCLAVRStart;
    int nGCPsPerLine;
    int eLocationIndicator;
    int iGCPStart;
    int iGCPStep;

    L1BFileFormat eL1BFormat;
    int nBufferSize;
    L1BSpaceCraftdID eSpacecraftID;
    L1BProductType eProductType;  // LAC, GAC, HRPT, FRAC
    L1BDataFormat iDataFormat;    // 10-bit packed or 16-bit unpacked
    int nRecordDataStart;
    int nRecordDataEnd;
    int nDataStartOffset;
    int nRecordSize;
    int nRecordSizeFromHeader;
    GUInt32 iInstrumentStatus;
    GUInt32 iChannelsMask;

    OGRSpatialReference m_oGCPSRS{};

    VSILFILE *fp;

    int bGuessDataFormat;

    int bByteSwap;

    int bExposeMaskBand;
    GDALRasterBand *poMaskBand;

    void ProcessRecordHeaders();
    int FetchGCPs(GDAL_GCP *, GByte *, int) const;
    static void FetchNOAA9TimeCode(TimeCode *, const GByte *, int *);
    void FetchNOAA15TimeCode(TimeCode *, const GByte *, int *) const;
    void FetchTimeCode(TimeCode *psTime, const void *pRecordHeader,
                       int *peLocationIndicator) const;
    CPLErr ProcessDatasetHeader(const char *pszFilename);
    int ComputeFileOffsets();

    void FetchMetadata();
    void FetchMetadataNOAA15();

    vsi_l_offset GetLineOffset(int nBlockYOff) const;

    GUInt16 GetUInt16(const void *pabyData) const;
    GInt16 GetInt16(const void *pabyData) const;
    GUInt32 GetUInt32(const void *pabyData) const;
    GInt32 GetInt32(const void *pabyData) const;

    static L1BFileFormat DetectFormat(const char *pszFilename,
                                      const GByte *pabyHeader,
                                      int nHeaderBytes);

  public:
    explicit L1BDataset(L1BFileFormat);
    virtual ~L1BDataset();

    virtual int GetGCPCount() override;
    const OGRSpatialReference *GetGCPSpatialRef() const override;
    virtual const GDAL_GCP *GetGCPs() override;

    static int Identify(GDALOpenInfo *);
    static GDALDataset *Open(GDALOpenInfo *);
};

/************************************************************************/
/* ==================================================================== */
/*                            L1BRasterBand                             */
/* ==================================================================== */
/************************************************************************/

class L1BRasterBand final : public GDALPamRasterBand
{
    friend class L1BDataset;

  public:
    L1BRasterBand(L1BDataset *, int);

    //    virtual double GetNoDataValue( int *pbSuccess = NULL );
    virtual CPLErr IReadBlock(int, int, void *) override;
    virtual GDALRasterBand *GetMaskBand() override;
    virtual int GetMaskFlags() override;
};

/************************************************************************/
/* ==================================================================== */
/*                            L1BMaskBand                               */
/* ==================================================================== */
/************************************************************************/

class L1BMaskBand final : public GDALPamRasterBand
{
    friend class L1BDataset;

  public:
    explicit L1BMaskBand(L1BDataset *);

    virtual CPLErr IReadBlock(int, int, void *) override;
};

/************************************************************************/
/*                            L1BMaskBand()                             */
/************************************************************************/

L1BMaskBand::L1BMaskBand(L1BDataset *poDSIn)
{
    CPLAssert(poDSIn->eL1BFormat == L1B_NOAA15 ||
              poDSIn->eL1BFormat == L1B_NOAA15_NOHDR);

    poDS = poDSIn;
    eDataType = GDT_Byte;

    nRasterXSize = poDS->GetRasterXSize();
    nRasterYSize = poDS->GetRasterYSize();
    nBlockXSize = poDS->GetRasterXSize();
    nBlockYSize = 1;
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr L1BMaskBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
                               void *pImage)
{
    L1BDataset *poGDS = (L1BDataset *)poDS;

    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poGDS->fp, poGDS->GetLineOffset(nBlockYOff) + 24, SEEK_SET));

    GByte abyData[4];
    CPL_IGNORE_RET_VAL(VSIFReadL(abyData, 1, 4, poGDS->fp));
    GUInt32 n32 = poGDS->GetUInt32(abyData);

    if ((n32 >> 31) != 0) /* fatal flag */
        memset(pImage, 0, nBlockXSize);
    else
        memset(pImage, 255, nBlockXSize);

    return CE_None;
}

/************************************************************************/
/*                           L1BRasterBand()                            */
/************************************************************************/

L1BRasterBand::L1BRasterBand(L1BDataset *poDSIn, int nBandIn)

{
    poDS = poDSIn;
    nBand = nBandIn;
    eDataType = GDT_UInt16;

    nBlockXSize = poDS->GetRasterXSize();
    nBlockYSize = 1;
}

/************************************************************************/
/*                           GetMaskBand()                              */
/************************************************************************/

GDALRasterBand *L1BRasterBand::GetMaskBand()
{
    L1BDataset *poGDS = (L1BDataset *)poDS;
    if (poGDS->poMaskBand)
        return poGDS->poMaskBand;
    return GDALPamRasterBand::GetMaskBand();
}

/************************************************************************/
/*                           GetMaskFlags()                             */
/************************************************************************/

int L1BRasterBand::GetMaskFlags()
{
    L1BDataset *poGDS = (L1BDataset *)poDS;
    if (poGDS->poMaskBand)
        return GMF_PER_DATASET;
    return GDALPamRasterBand::GetMaskFlags();
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr L1BRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
                                 void *pImage)
{
    L1BDataset *poGDS = (L1BDataset *)poDS;

    /* -------------------------------------------------------------------- */
    /*      Seek to data.                                                   */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poGDS->fp, poGDS->GetLineOffset(nBlockYOff), SEEK_SET));

    /* -------------------------------------------------------------------- */
    /*      Read data into the buffer.                                      */
    /* -------------------------------------------------------------------- */
    GUInt16 *iScan = nullptr;  // Unpacked 16-bit scanline buffer
    int i, j;

    switch (poGDS->iDataFormat)
    {
        case PACKED10BIT:
        {
            // Read packed scanline
            GUInt32 *iRawScan = (GUInt32 *)CPLMalloc(poGDS->nRecordSize);
            CPL_IGNORE_RET_VAL(
                VSIFReadL(iRawScan, 1, poGDS->nRecordSize, poGDS->fp));

            iScan = (GUInt16 *)CPLMalloc(poGDS->nBufferSize);
            j = 0;
            for (i = poGDS->nRecordDataStart / (int)sizeof(iRawScan[0]);
                 i < poGDS->nRecordDataEnd / (int)sizeof(iRawScan[0]); i++)
            {
                GUInt32 iWord1 = poGDS->GetUInt32(&iRawScan[i]);
                GUInt32 iWord2 = iWord1 & 0x3FF00000;

                iScan[j++] = (GUInt16)(iWord2 >> 20);
                iWord2 = iWord1 & 0x000FFC00;
                iScan[j++] = (GUInt16)(iWord2 >> 10);
                iScan[j++] = (GUInt16)(iWord1 & 0x000003FF);
            }
            CPLFree(iRawScan);
        }
        break;
        case UNPACKED16BIT:
        {
            // Read unpacked scanline
            GUInt16 *iRawScan = (GUInt16 *)CPLMalloc(poGDS->nRecordSize);
            CPL_IGNORE_RET_VAL(
                VSIFReadL(iRawScan, 1, poGDS->nRecordSize, poGDS->fp));

            iScan = (GUInt16 *)CPLMalloc(
                sizeof(GUInt16) * poGDS->GetRasterXSize() * poGDS->nBands);
            for (i = 0; i < poGDS->GetRasterXSize() * poGDS->nBands; i++)
            {
                iScan[i] =
                    poGDS->GetUInt16(&iRawScan[poGDS->nRecordDataStart /
                                                   (int)sizeof(iRawScan[0]) +
                                               i]);
            }
            CPLFree(iRawScan);
        }
        break;
        case UNPACKED8BIT:
        {
            // Read 8-bit unpacked scanline
            GByte *byRawScan = (GByte *)CPLMalloc(poGDS->nRecordSize);
            CPL_IGNORE_RET_VAL(
                VSIFReadL(byRawScan, 1, poGDS->nRecordSize, poGDS->fp));

            iScan = (GUInt16 *)CPLMalloc(
                sizeof(GUInt16) * poGDS->GetRasterXSize() * poGDS->nBands);
            for (i = 0; i < poGDS->GetRasterXSize() * poGDS->nBands; i++)
                iScan[i] = byRawScan[poGDS->nRecordDataStart /
                                         (int)sizeof(byRawScan[0]) +
                                     i];
            CPLFree(byRawScan);
        }
        break;
        default:  // NOTREACHED
            break;
    }

    int nBlockSize = nBlockXSize * nBlockYSize;
    if (poGDS->eLocationIndicator == DESCEND)
    {
        for (i = 0, j = 0; i < nBlockSize; i++)
        {
            ((GUInt16 *)pImage)[i] = iScan[j + nBand - 1];
            j += poGDS->nBands;
        }
    }
    else
    {
        for (i = nBlockSize - 1, j = 0; i >= 0; i--)
        {
            ((GUInt16 *)pImage)[i] = iScan[j + nBand - 1];
            j += poGDS->nBands;
        }
    }

    CPLFree(iScan);
    return CE_None;
}

/************************************************************************/
/*                           L1BDataset()                               */
/************************************************************************/

L1BDataset::L1BDataset(L1BFileFormat eL1BFormatIn)
    : eSource(UNKNOWN_STATION), eProcCenter(UNKNOWN_CENTER),
      // sStartTime
      // sStopTime
      bHighGCPDensityStrategy(
          CPLTestBool(CPLGetConfigOption("L1B_HIGH_GCP_DENSITY", "TRUE"))),
      pasGCPList(nullptr), nGCPCount(0), iGCPOffset(0), iGCPCodeOffset(0),
      iCLAVRStart(0), nGCPsPerLine(0),
      eLocationIndicator(DESCEND),  // XXX: should be initialized
      iGCPStart(0), iGCPStep(0), eL1BFormat(eL1BFormatIn), nBufferSize(0),
      eSpacecraftID(TIROSN), eProductType(HRPT), iDataFormat(PACKED10BIT),
      nRecordDataStart(0), nRecordDataEnd(0), nDataStartOffset(0),
      nRecordSize(0), nRecordSizeFromHeader(0), iInstrumentStatus(0),
      iChannelsMask(0), fp(nullptr), bGuessDataFormat(FALSE),
      // L1B is normally big-endian ordered, so byte-swap on little-endian CPU.
      bByteSwap(CPL_IS_LSB), bExposeMaskBand(FALSE), poMaskBand(nullptr)
{
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    m_oGCPSRS.importFromWkt(
        "GEOGCS[\"WGS 72\",DATUM[\"WGS_1972\","
        "SPHEROID[\"WGS 72\",6378135,298.26,AUTHORITY[\"EPSG\",7043]],"
        "TOWGS84[0,0,4.5,0,0,0.554,0.2263],AUTHORITY[\"EPSG\",6322]],"
        "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],"
        "UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],"
        "AUTHORITY[\"EPSG\",4322]]");
}

/************************************************************************/
/*                            ~L1BDataset()                             */
/************************************************************************/

L1BDataset::~L1BDataset()

{
    FlushCache(true);

    if (nGCPCount > 0)
    {
        GDALDeinitGCPs(nGCPCount, pasGCPList);
        CPLFree(pasGCPList);
    }
    if (fp != nullptr)
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
    delete poMaskBand;
}

/************************************************************************/
/*                          GetLineOffset()                             */
/************************************************************************/

vsi_l_offset L1BDataset::GetLineOffset(int nBlockYOff) const
{
    return (eLocationIndicator == DESCEND)
               ? nDataStartOffset + (vsi_l_offset)nBlockYOff * nRecordSize
               : nDataStartOffset +
                     (vsi_l_offset)(nRasterYSize - nBlockYOff - 1) *
                         nRecordSize;
}

/************************************************************************/
/*                            GetGCPCount()                             */
/************************************************************************/

int L1BDataset::GetGCPCount()

{
    return nGCPCount;
}

/************************************************************************/
/*                          GetGCPSpatialRef()                          */
/************************************************************************/

const OGRSpatialReference *L1BDataset::GetGCPSpatialRef() const

{
    if (nGCPCount > 0 && !m_oGCPSRS.IsEmpty())
        return &m_oGCPSRS;
    else
        return nullptr;
}

/************************************************************************/
/*                               GetGCPs()                              */
/************************************************************************/

const GDAL_GCP *L1BDataset::GetGCPs()
{
    return pasGCPList;
}

/************************************************************************/
/*      Byte swapping helpers                                           */
/************************************************************************/

GUInt16 L1BDataset::GetUInt16(const void *pabyData) const
{
    GUInt16 iTemp;
    memcpy(&iTemp, pabyData, 2);
    if (bByteSwap)
        return CPL_SWAP16(iTemp);
    return iTemp;
}

GInt16 L1BDataset::GetInt16(const void *pabyData) const
{
    GInt16 iTemp;
    memcpy(&iTemp, pabyData, 2);
    if (bByteSwap)
        return CPL_SWAP16(iTemp);
    return iTemp;
}

GUInt32 L1BDataset::GetUInt32(const void *pabyData) const
{
    GUInt32 lTemp;
    memcpy(&lTemp, pabyData, 4);
    if (bByteSwap)
        return CPL_SWAP32(lTemp);
    return lTemp;
}

GInt32 L1BDataset::GetInt32(const void *pabyData) const
{
    GInt32 lTemp;
    memcpy(&lTemp, pabyData, 4);
    if (bByteSwap)
        return CPL_SWAP32(lTemp);
    return lTemp;
}

/************************************************************************/
/*      Fetch timecode from the record header (NOAA9-NOAA14 version)    */
/************************************************************************/

void L1BDataset::FetchNOAA9TimeCode(TimeCode *psTime,
                                    const GByte *piRecordHeader,
                                    int *peLocationIndicator)
{
    GUInt32 lTemp;

    lTemp = ((piRecordHeader[2] >> 1) & 0x7F);
    psTime->SetYear((lTemp > 77)
                        ? (lTemp + 1900)
                        : (lTemp + 2000));  // Avoid `Year 2000' problem
    psTime->SetDay((GUInt32)(piRecordHeader[2] & 0x01) << 8 |
                   (GUInt32)piRecordHeader[3]);
    psTime->SetMillisecond(((GUInt32)(piRecordHeader[4] & 0x07) << 24) |
                           ((GUInt32)piRecordHeader[5] << 16) |
                           ((GUInt32)piRecordHeader[6] << 8) |
                           (GUInt32)piRecordHeader[7]);
    if (peLocationIndicator)
    {
        *peLocationIndicator =
            ((piRecordHeader[8] & 0x02) == 0) ? ASCEND : DESCEND;
    }
}

/************************************************************************/
/*      Fetch timecode from the record header (NOAA15-METOP2 version)   */
/************************************************************************/

void L1BDataset::FetchNOAA15TimeCode(TimeCode *psTime,
                                     const GByte *pabyRecordHeader,
                                     int *peLocationIndicator) const
{
    psTime->SetYear(GetUInt16(pabyRecordHeader + 2));
    psTime->SetDay(GetUInt16(pabyRecordHeader + 4));
    psTime->SetMillisecond(GetUInt32(pabyRecordHeader + 8));
    if (peLocationIndicator)
    {
        // FIXME: hemisphere
        *peLocationIndicator =
            ((GetUInt16(pabyRecordHeader + 12) & 0x8000) == 0) ? ASCEND
                                                               : DESCEND;
    }
}

/************************************************************************/
/*                          FetchTimeCode()                             */
/************************************************************************/

void L1BDataset::FetchTimeCode(TimeCode *psTime, const void *pRecordHeader,
                               int *peLocationIndicator) const
{
    if (eSpacecraftID <= NOAA14)
    {
        FetchNOAA9TimeCode(psTime, (const GByte *)pRecordHeader,
                           peLocationIndicator);
    }
    else
    {
        FetchNOAA15TimeCode(psTime, (const GByte *)pRecordHeader,
                            peLocationIndicator);
    }
}

/************************************************************************/
/*      Fetch GCPs from the individual scanlines                        */
/************************************************************************/

int L1BDataset::FetchGCPs(GDAL_GCP *pasGCPListRow, GByte *pabyRecordHeader,
                          int iLine) const
{
    // LAC and HRPT GCPs are tied to the center of pixel,
    // GAC ones are slightly displaced.
    double dfDelta = (eProductType == GAC) ? 0.9 : 0.5;
    double dfPixel = (eLocationIndicator == DESCEND)
                         ? iGCPStart + dfDelta
                         : (nRasterXSize - (iGCPStart + dfDelta));

    int nGCPs;
    if (eSpacecraftID <= NOAA14)
    {
        // NOAA9-NOAA14 records have an indicator of number of working GCPs.
        // Number of good GCPs may be smaller than the total amount of points.
        nGCPs = (*(pabyRecordHeader + iGCPCodeOffset) < nGCPsPerLine)
                    ? *(pabyRecordHeader + iGCPCodeOffset)
                    : nGCPsPerLine;
#ifdef DEBUG_VERBOSE
        CPLDebug("L1B", "iGCPCodeOffset=%d, nGCPsPerLine=%d, nGoodGCPs=%d",
                 iGCPCodeOffset, nGCPsPerLine, nGCPs);
#endif
    }
    else
        nGCPs = nGCPsPerLine;

    pabyRecordHeader += iGCPOffset;

    int nGCPCountRow = 0;
    while (nGCPs--)
    {
        if (eSpacecraftID <= NOAA14)
        {
            GInt16 nRawY = GetInt16(pabyRecordHeader);
            pabyRecordHeader += sizeof(GInt16);
            GInt16 nRawX = GetInt16(pabyRecordHeader);
            pabyRecordHeader += sizeof(GInt16);

            pasGCPListRow[nGCPCountRow].dfGCPY = nRawY / L1B_NOAA9_GCP_SCALE;
            pasGCPListRow[nGCPCountRow].dfGCPX = nRawX / L1B_NOAA9_GCP_SCALE;
        }
        else
        {
            GInt32 nRawY = GetInt32(pabyRecordHeader);
            pabyRecordHeader += sizeof(GInt32);
            GInt32 nRawX = GetInt32(pabyRecordHeader);
            pabyRecordHeader += sizeof(GInt32);

            pasGCPListRow[nGCPCountRow].dfGCPY = nRawY / L1B_NOAA15_GCP_SCALE;
            pasGCPListRow[nGCPCountRow].dfGCPX = nRawX / L1B_NOAA15_GCP_SCALE;
        }

        if (pasGCPListRow[nGCPCountRow].dfGCPX < -180 ||
            pasGCPListRow[nGCPCountRow].dfGCPX > 180 ||
            pasGCPListRow[nGCPCountRow].dfGCPY < -90 ||
            pasGCPListRow[nGCPCountRow].dfGCPY > 90)
            continue;

        pasGCPListRow[nGCPCountRow].dfGCPZ = 0.0;
        pasGCPListRow[nGCPCountRow].dfGCPPixel = dfPixel;
        dfPixel += (eLocationIndicator == DESCEND) ? iGCPStep : -iGCPStep;
        pasGCPListRow[nGCPCountRow].dfGCPLine =
            (double)((eLocationIndicator == DESCEND)
                         ? iLine
                         : nRasterYSize - iLine - 1) +
            0.5;
        nGCPCountRow++;
    }
    return nGCPCountRow;
}

/************************************************************************/
/*                      ProcessRecordHeaders()                          */
/************************************************************************/

void L1BDataset::ProcessRecordHeaders()
{
    void *pRecordHeader = CPLCalloc(1, nRecordDataStart);

    CPL_IGNORE_RET_VAL(VSIFSeekL(fp, nDataStartOffset, SEEK_SET));
    CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));

    FetchTimeCode(&sStartTime, pRecordHeader, &eLocationIndicator);

    CPL_IGNORE_RET_VAL(VSIFSeekL(
        fp, nDataStartOffset + (nRasterYSize - 1) * nRecordSize, SEEK_SET));
    CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));

    FetchTimeCode(&sStopTime, pRecordHeader, nullptr);

    /* -------------------------------------------------------------------- */
    /*      Pick a skip factor so that we will get roughly 20 lines         */
    /*      worth of GCPs.  That should give respectable coverage on all    */
    /*      but the longest swaths.                                         */
    /* -------------------------------------------------------------------- */
    int nTargetLines;
    double dfLineStep = 0.0;

    if (bHighGCPDensityStrategy)
    {
        if (nRasterYSize < nGCPsPerLine)
        {
            nTargetLines = nRasterYSize;
        }
        else
        {
            int nColStep;
            nColStep = nRasterXSize / nGCPsPerLine;
            if (nRasterYSize >= nRasterXSize)
            {
                dfLineStep = nColStep;
            }
            else
            {
                dfLineStep = nRasterYSize / nGCPsPerLine;
            }
            nTargetLines = static_cast<int>(nRasterYSize / dfLineStep);
        }
    }
    else
    {
        nTargetLines = std::min(DESIRED_LINES_OF_GCPS, nRasterYSize);
    }
    if (nTargetLines > 1)
        dfLineStep = 1.0 * (nRasterYSize - 1) / (nTargetLines - 1);

    /* -------------------------------------------------------------------- */
    /*      Initialize the GCP list.                                        */
    /* -------------------------------------------------------------------- */
    const int nExpectedGCPs = nTargetLines * nGCPsPerLine;
    if (nExpectedGCPs > 0)
    {
        pasGCPList =
            (GDAL_GCP *)VSI_CALLOC_VERBOSE(nExpectedGCPs, sizeof(GDAL_GCP));
        if (pasGCPList == nullptr)
        {
            CPLFree(pRecordHeader);
            return;
        }
        GDALInitGCPs(nExpectedGCPs, pasGCPList);
    }

    /* -------------------------------------------------------------------- */
    /*      Fetch the GCPs for each selected line.  We force the last       */
    /*      line sampled to be the last line in the dataset even if that    */
    /*      leaves a bigger than expected gap.                              */
    /* -------------------------------------------------------------------- */
    int iStep;
    int iPrevLine = -1;

    for (iStep = 0; iStep < nTargetLines; iStep++)
    {
        int iLine;

        if (iStep == nTargetLines - 1)
            iLine = nRasterYSize - 1;
        else
            iLine = (int)(dfLineStep * iStep);
        if (iLine == iPrevLine)
            continue;
        iPrevLine = iLine;

        CPL_IGNORE_RET_VAL(
            VSIFSeekL(fp, nDataStartOffset + iLine * nRecordSize, SEEK_SET));
        CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));

        int nGCPsOnThisLine =
            FetchGCPs(pasGCPList + nGCPCount, (GByte *)pRecordHeader, iLine);

        if (!bHighGCPDensityStrategy)
        {
            /* --------------------------------------------------------------------
             */
            /*      We don't really want too many GCPs per line.  Downsample to
             */
            /*      11 per line. */
            /* --------------------------------------------------------------------
             */

            const int nDesiredGCPsPerLine =
                std::min(DESIRED_GCPS_PER_LINE, nGCPsOnThisLine);
            int nGCPStep =
                (nDesiredGCPsPerLine > 1)
                    ? (nGCPsOnThisLine - 1) / (nDesiredGCPsPerLine - 1)
                    : 1;
            int iSrcGCP = nGCPCount;
            int iDstGCP = nGCPCount;

            if (nGCPStep == 0)
                nGCPStep = 1;

            for (int iGCP = 0; iGCP < nDesiredGCPsPerLine; iGCP++)
            {
                if (iGCP == nDesiredGCPsPerLine - 1)
                    iSrcGCP = nGCPCount + nGCPsOnThisLine - 1;
                else
                    iSrcGCP += nGCPStep;
                iDstGCP++;

                pasGCPList[iDstGCP].dfGCPX = pasGCPList[iSrcGCP].dfGCPX;
                pasGCPList[iDstGCP].dfGCPY = pasGCPList[iSrcGCP].dfGCPY;
                pasGCPList[iDstGCP].dfGCPPixel = pasGCPList[iSrcGCP].dfGCPPixel;
                pasGCPList[iDstGCP].dfGCPLine = pasGCPList[iSrcGCP].dfGCPLine;
            }

            nGCPCount += nDesiredGCPsPerLine;
        }
        else
        {
            nGCPCount += nGCPsOnThisLine;
        }
    }

    if (nGCPCount < nExpectedGCPs)
    {
        GDALDeinitGCPs(nExpectedGCPs - nGCPCount, pasGCPList + nGCPCount);
        if (nGCPCount == 0)
        {
            CPLFree(pasGCPList);
            pasGCPList = nullptr;
        }
    }

    CPLFree(pRecordHeader);

    /* -------------------------------------------------------------------- */
    /*      Set fetched information as metadata records                     */
    /* -------------------------------------------------------------------- */
    // Time of first scanline
    SetMetadataItem("START", sStartTime.PrintTime());
    // Time of last scanline
    SetMetadataItem("STOP", sStopTime.PrintTime());
    // AVHRR Earth location indication

    switch (eLocationIndicator)
    {
        case ASCEND:
            SetMetadataItem("LOCATION", "Ascending");
            break;
        case DESCEND:
        default:
            SetMetadataItem("LOCATION", "Descending");
            break;
    }
}

/************************************************************************/
/*                           FetchMetadata()                            */
/************************************************************************/

void L1BDataset::FetchMetadata()
{
    if (eL1BFormat != L1B_NOAA9)
    {
        FetchMetadataNOAA15();
        return;
    }

    std::string osDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", "");
    if (osDir.empty())
    {
        osDir = CPLGetPathSafe(GetDescription());
        if (osDir.empty())
            osDir = ".";
    }
    const std::string osMetadataFile =
        std::string(osDir)
            .append("/")
            .append(CPLGetFilename(GetDescription()))
            .append("_metadata.csv");
    VSILFILE *fpCSV = VSIFOpenL(osMetadataFile.c_str(), "wb");
    if (fpCSV == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Cannot create metadata file : %s", osMetadataFile.c_str());
        return;
    }

    CPL_IGNORE_RET_VAL(
        VSIFPrintfL(fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "FATAL_FLAG,TIME_ERROR,DATA_GAP,DATA_JITTER,INSUFFICIENT_DATA_"
               "FOR_CAL,NO_EARTH_LOCATION,DESCEND,P_N_STATUS,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "BIT_SYNC_STATUS,SYNC_ERROR,FRAME_SYNC_ERROR,FLYWHEELING,BIT_"
               "SLIPPAGE,C3_SBBC,C4_SBBC,C5_SBBC,"));
    CPL_IGNORE_RET_VAL(
        VSIFPrintfL(fpCSV, "TIP_PARITY_FRAME_1,TIP_PARITY_FRAME_2,TIP_PARITY_"
                           "FRAME_3,TIP_PARITY_FRAME_4,TIP_PARITY_FRAME_5,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "SYNC_ERRORS,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "CAL_SLOPE_C1,CAL_INTERCEPT_C1,CAL_SLOPE_C2,CAL_INTERCEPT_C2,"
               "CAL_SLOPE_C3,CAL_INTERCEPT_C3,CAL_SLOPE_C4,CAL_INTERCEPT_C4,"
               "CAL_SLOPE_C5,CAL_INTERCEPT_C5,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "NUM_SOLZENANGLES_EARTHLOCPNTS"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);

    for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
    {
        /* --------------------------------------------------------------------
         */
        /*      Seek to data. */
        /* --------------------------------------------------------------------
         */
        CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));

        CPL_IGNORE_RET_VAL(
            VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));

        GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);

        TimeCode timeCode;
        FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);

        CPL_IGNORE_RET_VAL(
            VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,", nScanlineNumber, nBlockYOff,
                        (int)timeCode.GetYear(), (int)timeCode.GetDay(),
                        (int)timeCode.GetMillisecond()));
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[8] >> 7) & 1,
            (pabyRecordHeader[8] >> 6) & 1, (pabyRecordHeader[8] >> 5) & 1,
            (pabyRecordHeader[8] >> 4) & 1, (pabyRecordHeader[8] >> 3) & 1,
            (pabyRecordHeader[8] >> 2) & 1, (pabyRecordHeader[8] >> 1) & 1,
            (pabyRecordHeader[8] >> 0) & 1));
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[9] >> 7) & 1,
            (pabyRecordHeader[9] >> 6) & 1, (pabyRecordHeader[9] >> 5) & 1,
            (pabyRecordHeader[9] >> 4) & 1, (pabyRecordHeader[9] >> 3) & 1,
            (pabyRecordHeader[9] >> 2) & 1, (pabyRecordHeader[9] >> 1) & 1,
            (pabyRecordHeader[9] >> 0) & 1));
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,", (pabyRecordHeader[10] >> 7) & 1,
            (pabyRecordHeader[10] >> 6) & 1, (pabyRecordHeader[10] >> 5) & 1,
            (pabyRecordHeader[10] >> 4) & 1, (pabyRecordHeader[10] >> 3) & 1));
        CPL_IGNORE_RET_VAL(
            VSIFPrintfL(fpCSV, "%d,", pabyRecordHeader[11] >> 2));
        GInt32 i32;
        for (int i = 0; i < 10; i++)
        {
            i32 = GetInt32(pabyRecordHeader + 12 + 4 * i);
            /* Scales :
             * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-3.htm
             */
            if ((i % 2) == 0)
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 30.0)));
            else
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 22.0)));
        }
        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d", pabyRecordHeader[52]));
        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
    }

    CPLFree(pabyRecordHeader);
    CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
}

/************************************************************************/
/*                         FetchMetadataNOAA15()                        */
/************************************************************************/

void L1BDataset::FetchMetadataNOAA15()
{
    int i, j;
    std::string osDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", "");
    if (osDir.empty())
    {
        osDir = CPLGetPathSafe(GetDescription());
        if (osDir.empty())
            osDir = ".";
    }
    const std::string osMetadataFile =
        std::string(osDir)
            .append("/")
            .append(CPLGetFilename(GetDescription()))
            .append("_metadata.csv");
    VSILFILE *fpCSV = VSIFOpenL(osMetadataFile.c_str(), "wb");
    if (fpCSV == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Cannot create metadata file : %s", osMetadataFile.c_str());
        return;
    }

    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,SAT_CLOCK_DRIF_DELTA,"
               "SOUTHBOUND,SCANTIME_CORRECTED,C3_SELECT,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV,
        "FATAL_FLAG,TIME_ERROR,DATA_GAP,INSUFFICIENT_DATA_FOR_CAL,"
        "NO_EARTH_LOCATION,FIRST_GOOD_TIME_AFTER_CLOCK_UPDATE,"
        "INSTRUMENT_STATUS_CHANGED,SYNC_LOCK_DROPPED,"
        "FRAME_SYNC_ERROR,FRAME_SYNC_DROPPED_LOCK,FLYWHEELING,"
        "BIT_SLIPPAGE,TIP_PARITY_ERROR,REFLECTED_SUNLIGHT_C3B,"
        "REFLECTED_SUNLIGHT_C4,REFLECTED_SUNLIGHT_C5,RESYNC,P_N_STATUS,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "BAD_TIME_CAN_BE_INFERRED,BAD_TIME_CANNOT_BE_INFERRED,"
               "TIME_DISCONTINUITY,REPEAT_SCAN_TIME,"));
    CPL_IGNORE_RET_VAL(
        VSIFPrintfL(fpCSV, "UNCALIBRATED_BAD_TIME,CALIBRATED_FEWER_SCANLINES,"
                           "UNCALIBRATED_BAD_PRT,CALIBRATED_MARGINAL_PRT,"
                           "UNCALIBRATED_CHANNELS,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "NO_EARTH_LOC_BAD_TIME,EARTH_LOC_QUESTIONABLE_TIME,"
               "EARTH_LOC_QUESTIONABLE,EARTH_LOC_VERY_QUESTIONABLE,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV,
        "C3B_UNCALIBRATED,C3B_QUESTIONABLE,C3B_ALL_BLACKBODY,"
        "C3B_ALL_SPACEVIEW,C3B_MARGINAL_BLACKBODY,C3B_MARGINAL_SPACEVIEW,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV,
        "C4_UNCALIBRATED,C4_QUESTIONABLE,C4_ALL_BLACKBODY,"
        "C4_ALL_SPACEVIEW,C4_MARGINAL_BLACKBODY,C4_MARGINAL_SPACEVIEW,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV,
        "C5_UNCALIBRATED,C5_QUESTIONABLE,C5_ALL_BLACKBODY,"
        "C5_ALL_SPACEVIEW,C5_MARGINAL_BLACKBODY,C5_MARGINAL_SPACEVIEW,"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "BIT_ERRORS,"));
    for (i = 0; i < 3; i++)
    {
        const char *pszChannel = (i == 0) ? "C1" : (i == 1) ? "C2" : "C3A";
        for (j = 0; j < 3; j++)
        {
            const char *pszType = (j == 0)   ? "OP"
                                  : (j == 1) ? "TEST"
                                             : "PRELAUNCH";
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_1,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_1,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_2,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_2,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERSECTION,",
                                           pszType, pszChannel));
        }
    }
    for (i = 0; i < 3; i++)
    {
        const char *pszChannel = (i == 0) ? "C3B" : (i == 1) ? "C4" : "C5";
        for (j = 0; j < 2; j++)
        {
            const char *pszType = (j == 0) ? "OP" : "TEST";
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_1,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_2,",
                                           pszType, pszChannel));
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_3,",
                                           pszType, pszChannel));
        }
    }
    CPL_IGNORE_RET_VAL(VSIFPrintfL(
        fpCSV, "EARTH_LOC_CORR_TIP_EULER,EARTH_LOC_IND,"
               "SPACECRAFT_ATT_CTRL,ATT_SMODE,ATT_PASSIVE_WHEEL_TEST,"
               "TIME_TIP_EULER,TIP_EULER_ROLL,TIP_EULER_PITCH,TIP_EULER_YAW,"
               "SPACECRAFT_ALT"));
    CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);
    GInt16 i16;
    GUInt16 n16;
    GInt32 i32;
    GUInt32 n32;

    for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
    {
        /* --------------------------------------------------------------------
         */
        /*      Seek to data. */
        /* --------------------------------------------------------------------
         */
        CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));

        CPL_IGNORE_RET_VAL(
            VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));

        GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);

        TimeCode timeCode;
        FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);

        /* Clock drift delta */
        i16 = GetInt16(pabyRecordHeader + 6);
        /* Scanline bit field */
        n16 = GetInt16(pabyRecordHeader + 12);

        CPL_IGNORE_RET_VAL(
            VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,", nScanlineNumber,
                        nBlockYOff, (int)timeCode.GetYear(),
                        (int)timeCode.GetDay(), (int)timeCode.GetMillisecond(),
                        i16, (n16 >> 15) & 1, (n16 >> 14) & 1, (n16)&3));

        n32 = GetUInt32(pabyRecordHeader + 24);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,",
            (n32 >> 31) & 1, (n32 >> 30) & 1, (n32 >> 29) & 1, (n32 >> 28) & 1,
            (n32 >> 27) & 1, (n32 >> 26) & 1, (n32 >> 25) & 1, (n32 >> 24) & 1,
            (n32 >> 23) & 1, (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1,
            (n32 >> 8) & 1, (n32 >> 6) & 3, (n32 >> 4) & 3, (n32 >> 2) & 3,
            (n32 >> 1) & 1, (n32 >> 0) & 1));

        n32 = GetUInt32(pabyRecordHeader + 28);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,", (n32 >> 23) & 1,
            (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1, (n32 >> 15) & 1,
            (n32 >> 14) & 1, (n32 >> 13) & 1, (n32 >> 12) & 1, (n32 >> 11) & 1,
            (n32 >> 7) & 1, (n32 >> 6) & 1, (n32 >> 5) & 1, (n32 >> 4) & 1));

        for (i = 0; i < 3; i++)
        {
            n16 = GetUInt16(pabyRecordHeader + 32 + 2 * i);
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,",
                                           (n16 >> 7) & 1, (n16 >> 6) & 1,
                                           (n16 >> 5) & 1, (n16 >> 4) & 1,
                                           (n16 >> 2) & 1, (n16 >> 1) & 1));
        }

        /* Bit errors */
        n16 = GetUInt16(pabyRecordHeader + 38);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n16));

        int nOffset = 48;
        for (i = 0; i < 3; i++)
        {
            for (j = 0; j < 3; j++)
            {
                i32 = GetInt32(pabyRecordHeader + nOffset);
                nOffset += 4;
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 7.0)));
                i32 = GetInt32(pabyRecordHeader + nOffset);
                nOffset += 4;
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
                i32 = GetInt32(pabyRecordHeader + nOffset);
                nOffset += 4;
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 7.0)));
                i32 = GetInt32(pabyRecordHeader + nOffset);
                nOffset += 4;
                CPL_IGNORE_RET_VAL(
                    VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
                i32 = GetInt32(pabyRecordHeader + nOffset);
                nOffset += 4;
                CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", i32));
            }
        }
        for (i = 0; i < 18; i++)
        {
            i32 = GetInt32(pabyRecordHeader + nOffset);
            nOffset += 4;
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
        }

        n32 = GetUInt32(pabyRecordHeader + 312);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(
            fpCSV, "%d,%d,%d,%d,%d,", (n32 >> 16) & 1, (n32 >> 12) & 15,
            (n32 >> 8) & 15, (n32 >> 4) & 15, (n32 >> 0) & 15));

        n32 = GetUInt32(pabyRecordHeader + 316);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n32));

        for (i = 0; i < 3; i++)
        {
            i16 = GetUInt16(pabyRecordHeader + 320 + 2 * i);
            CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i16 / pow(10.0, 3.0)));
        }

        n16 = GetUInt16(pabyRecordHeader + 326);
        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f", n16 / pow(10.0, 1.0)));

        CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
    }

    CPLFree(pabyRecordHeader);
    CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
}

/************************************************************************/
/*                           EBCDICToASCII                              */
/************************************************************************/

constexpr GByte EBCDICToASCII[] = {
    0x00, 0x01, 0x02, 0x03, 0x9C, 0x09, 0x86, 0x7F, 0x97, 0x8D, 0x8E, 0x0B,
    0x0C, 0x0D, 0x0E, 0x0F, 0x10, 0x11, 0x12, 0x13, 0x9D, 0x85, 0x08, 0x87,
    0x18, 0x19, 0x92, 0x8F, 0x1C, 0x1D, 0x1E, 0x1F, 0x80, 0x81, 0x82, 0x83,
    0x84, 0x0A, 0x17, 0x1B, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x05, 0x06, 0x07,
    0x90, 0x91, 0x16, 0x93, 0x94, 0x95, 0x96, 0x04, 0x98, 0x99, 0x9A, 0x9B,
    0x14, 0x15, 0x9E, 0x1A, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0xA2, 0x2E, 0x3C, 0x28, 0x2B, 0x7C, 0x26, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x21, 0x24, 0x2A, 0x29, 0x3B, 0xAC,
    0x2D, 0x2F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xA6, 0x2C,
    0x25, 0x5F, 0x3E, 0x3F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x60, 0x3A, 0x23, 0x40, 0x27, 0x3D, 0x22, 0x00, 0x61, 0x62, 0x63,
    0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x6A, 0x6B, 0x6C, 0x6D, 0x6E, 0x6F, 0x70, 0x71, 0x72, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x7E, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78,
    0x79, 0x7A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x7B, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x7D, 0x4A, 0x4B, 0x4C, 0x4D, 0x4E, 0x4F, 0x50,
    0x51, 0x52, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x5C, 0x00, 0x53, 0x54,
    0x55, 0x56, 0x57, 0x58, 0x59, 0x5A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x9F,
};

/************************************************************************/
/*                      ProcessDatasetHeader()                          */
/************************************************************************/

CPLErr L1BDataset::ProcessDatasetHeader(const char *pszFilename)
{
    char szDatasetName[L1B_DATASET_NAME_SIZE + 1];

    if (eL1BFormat == L1B_NOAA9)
    {
        GByte abyTBMHeader[L1B_NOAA9_HEADER_SIZE];

        if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
            VSIFReadL(abyTBMHeader, 1, L1B_NOAA9_HEADER_SIZE, fp) <
                L1B_NOAA9_HEADER_SIZE)
        {
            CPLDebug("L1B", "Can't read NOAA-9/14 TBM header.");
            return CE_Failure;
        }

        // If dataset name in EBCDIC, decode it in ASCII
        if (*(abyTBMHeader + 8 + 25) == 'K' &&
            *(abyTBMHeader + 8 + 30) == 'K' &&
            *(abyTBMHeader + 8 + 33) == 'K' &&
            *(abyTBMHeader + 8 + 40) == 'K' &&
            *(abyTBMHeader + 8 + 46) == 'K' &&
            *(abyTBMHeader + 8 + 52) == 'K' && *(abyTBMHeader + 8 + 61) == 'K')
        {
            for (int i = 0; i < L1B_DATASET_NAME_SIZE; i++)
                abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i] =
                    EBCDICToASCII[abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i]];
        }

        // Fetch dataset name. NOAA-9/14 datasets contain the names in TBM
        // header only, so read it there.
        memcpy(szDatasetName, abyTBMHeader + L1B_NOAA9_HDR_NAME_OFF,
               L1B_DATASET_NAME_SIZE);
        szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';

        // Deal with a few NOAA <= 9 datasets with no dataset name in TBM header
        if (memcmp(szDatasetName,
                   "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0"
                   "\0\0\0\0\0\0\0\0\0\0\0\0\0",
                   L1B_DATASET_NAME_SIZE) == 0 &&
            strlen(pszFilename) == L1B_DATASET_NAME_SIZE)
        {
            strcpy(szDatasetName, pszFilename);
        }

        // Determine processing center where the dataset was created
        if (STARTS_WITH_CI(szDatasetName, "CMS"))
            eProcCenter = CMS;
        else if (STARTS_WITH_CI(szDatasetName, "DSS"))
            eProcCenter = DSS;
        else if (STARTS_WITH_CI(szDatasetName, "NSS"))
            eProcCenter = NSS;
        else if (STARTS_WITH_CI(szDatasetName, "UKM"))
            eProcCenter = UKM;
        else
            eProcCenter = UNKNOWN_CENTER;

        // Determine number of bands
        int i;
        for (i = 0; i < L1B_NOAA9_HDR_CHAN_SIZE; i++)
        {
            if (abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 1 ||
                abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 'Y')
            {
                nBands++;
                iChannelsMask |= (1 << i);
            }
        }
        if (nBands == 0 || nBands > 5)
        {
            nBands = 5;
            iChannelsMask = 0x1F;
        }

        // Determine data format (10-bit packed or 8/16-bit unpacked)
        if (STARTS_WITH_CI((const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF,
                           "10"))
            iDataFormat = PACKED10BIT;
        else if (STARTS_WITH_CI(
                     (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "16"))
            iDataFormat = UNPACKED16BIT;
        else if (STARTS_WITH_CI(
                     (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "08"))
            iDataFormat = UNPACKED8BIT;
        else if (STARTS_WITH_CI((const char *)abyTBMHeader +
                                    L1B_NOAA9_HDR_WORD_OFF,
                                "  ") ||
                 abyTBMHeader[L1B_NOAA9_HDR_WORD_OFF] == '\0')
            /* Empty string can be found in the following samples :
                http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/franh.1b (10
               bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/frang.1b (10
               bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/calfilel.1b
               (16 bit)
                http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/rapnzg.1b (16
               bit)
                ftp://ftp.sat.dundee.ac.uk/misc/testdata/noaa12/hrptnoaa1b.dat
               (10 bit)
            */
            bGuessDataFormat = TRUE;
        else
        {
#ifdef DEBUG
            CPLDebug("L1B", "Unknown data format \"%.2s\".",
                     abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF);
#endif
            return CE_Failure;
        }

        // Now read the dataset header record
        GByte abyRecHeader[L1B_NOAA9_HDR_REC_SIZE];
        if (VSIFSeekL(fp, L1B_NOAA9_HEADER_SIZE, SEEK_SET) < 0 ||
            VSIFReadL(abyRecHeader, 1, L1B_NOAA9_HDR_REC_SIZE, fp) <
                L1B_NOAA9_HDR_REC_SIZE)
        {
            CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
            return CE_Failure;
        }

        // Determine the spacecraft name
        switch (abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF])
        {
            case 4:
                eSpacecraftID = NOAA7;
                break;
            case 6:
                eSpacecraftID = NOAA8;
                break;
            case 7:
                eSpacecraftID = NOAA9;
                break;
            case 8:
                eSpacecraftID = NOAA10;
                break;
            case 1:
            {
                /* We could also use the time code to determine TIROS-N */
                if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
                    STARTS_WITH(pszFilename + 8, ".TN."))
                    eSpacecraftID = TIROSN;
                else
                    eSpacecraftID = NOAA11;
                break;
            }
            case 5:
                eSpacecraftID = NOAA12;
                break;
            case 2:
            {
                /* We could also use the time code to determine NOAA6 */
                if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
                    STARTS_WITH(pszFilename + 8, ".NA."))
                    eSpacecraftID = NOAA6;
                else
                    eSpacecraftID = NOAA13;
                break;
            }
            case 3:
                eSpacecraftID = NOAA14;
                break;
            default:
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Unknown spacecraft ID \"%d\".",
                         abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF]);

                eSpacecraftID = NOAA9_UNKNOWN;
                break;
        }

        // Determine the product data type
        int iWord = abyRecHeader[L1B_NOAA9_HDR_REC_PROD_OFF] >> 4;
        switch (iWord)
        {
            case 1:
                eProductType = LAC;
                break;
            case 2:
                eProductType = GAC;
                break;
            case 3:
                eProductType = HRPT;
                break;
            default:
#ifdef DEBUG
                CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
#endif
                return CE_Failure;
        }

        // Determine receiving station name
        iWord = (abyRecHeader[L1B_NOAA9_HDR_REC_DSTAT_OFF] & 0x60) >> 5;
        switch (iWord)
        {
            case 1:
                eSource = GC;
                break;
            case 2:
                eSource = WI;
                break;
            case 3:
                eSource = SO;
                break;
            default:
                eSource = UNKNOWN_STATION;
                break;
        }
    }

    else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
    {
        if (eL1BFormat == L1B_NOAA15)
        {
            GByte abyARSHeader[L1B_NOAA15_HEADER_SIZE];

            if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
                VSIFReadL(abyARSHeader, 1, L1B_NOAA15_HEADER_SIZE, fp) <
                    L1B_NOAA15_HEADER_SIZE)
            {
                CPLDebug("L1B", "Can't read NOAA-15 ARS header.");
                return CE_Failure;
            }

            // Determine number of bands
            int i;
            for (i = 0; i < L1B_NOAA15_HDR_CHAN_SIZE; i++)
            {
                if (abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 1 ||
                    abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 'Y')
                {
                    nBands++;
                    iChannelsMask |= (1 << i);
                }
            }
            if (nBands == 0 || nBands > 5)
            {
                nBands = 5;
                iChannelsMask = 0x1F;
            }

            // Determine data format (10-bit packed or 8/16-bit unpacked)
            if (STARTS_WITH_CI(
                    (const char *)abyARSHeader + L1B_NOAA15_HDR_WORD_OFF, "10"))
                iDataFormat = PACKED10BIT;
            else if (STARTS_WITH_CI((const char *)abyARSHeader +
                                        L1B_NOAA15_HDR_WORD_OFF,
                                    "16"))
                iDataFormat = UNPACKED16BIT;
            else if (STARTS_WITH_CI((const char *)abyARSHeader +
                                        L1B_NOAA15_HDR_WORD_OFF,
                                    "08"))
                iDataFormat = UNPACKED8BIT;
            else
            {
#ifdef DEBUG
                CPLDebug("L1B", "Unknown data format \"%.2s\".",
                         abyARSHeader + L1B_NOAA9_HDR_WORD_OFF);
#endif
                return CE_Failure;
            }
        }
        else
        {
            nBands = 5;
            iChannelsMask = 0x1F;
            iDataFormat = PACKED10BIT;
        }

        // Now read the dataset header record
        GByte abyRecHeader[L1B_NOAA15_HDR_REC_SIZE];
        if (VSIFSeekL(fp,
                      (eL1BFormat == L1B_NOAA15) ? L1B_NOAA15_HEADER_SIZE : 0,
                      SEEK_SET) < 0 ||
            VSIFReadL(abyRecHeader, 1, L1B_NOAA15_HDR_REC_SIZE, fp) <
                L1B_NOAA15_HDR_REC_SIZE)
        {
            CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
            return CE_Failure;
        }

        // Fetch dataset name
        memcpy(szDatasetName, abyRecHeader + L1B_NOAA15_HDR_REC_NAME_OFF,
               L1B_DATASET_NAME_SIZE);
        szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';

        // Determine processing center where the dataset was created
        if (STARTS_WITH_CI((const char *)abyRecHeader +
                               L1B_NOAA15_HDR_REC_SITE_OFF,
                           "CMS"))
            eProcCenter = CMS;
        else if (STARTS_WITH_CI((const char *)abyRecHeader +
                                    L1B_NOAA15_HDR_REC_SITE_OFF,
                                "DSS"))
            eProcCenter = DSS;
        else if (STARTS_WITH_CI((const char *)abyRecHeader +
                                    L1B_NOAA15_HDR_REC_SITE_OFF,
                                "NSS"))
            eProcCenter = NSS;
        else if (STARTS_WITH_CI((const char *)abyRecHeader +
                                    L1B_NOAA15_HDR_REC_SITE_OFF,
                                "UKM"))
            eProcCenter = UKM;
        else
            eProcCenter = UNKNOWN_CENTER;

        int nFormatVersionYear, nFormatVersionDayOfYear, nHeaderRecCount;

        /* Some products from NOAA-18 and NOAA-19 coming from 'ess' processing
         * station */
        /* have little-endian ordering. Try to detect it with some consistency
         * checks */
        int i = 0;
        do
        {
            nFormatVersionYear = GetUInt16(
                abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_YEAR_OFF);
            nFormatVersionDayOfYear = GetUInt16(
                abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_DAY_OFF);
            nHeaderRecCount =
                GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_HDR_REC_COUNT_OFF);
            if (i == 2)
                break;
            if (!(nFormatVersionYear >= 1980 && nFormatVersionYear <= 2100) &&
                !(nFormatVersionDayOfYear <= 366) && !(nHeaderRecCount == 1))
            {
                if (i == 0)
                    CPLDebug("L1B", "Trying little-endian ordering");
                else
                    CPLDebug("L1B", "Not completely convincing... Returning to "
                                    "big-endian order");
                bByteSwap = !bByteSwap;
            }
            else
                break;
            i++;
        } while (i <= 2);
        nRecordSizeFromHeader =
            GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_LOGICAL_REC_LENGTH_OFF);
        int nFormatVersion =
            GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_OFF);
        CPLDebug("L1B", "NOAA Level 1b Format Version Number = %d",
                 nFormatVersion);
        CPLDebug("L1B", "Level 1b Format Version Year = %d",
                 nFormatVersionYear);
        CPLDebug("L1B", "Level 1b Format Version Day of Year = %d",
                 nFormatVersionDayOfYear);
        CPLDebug("L1B",
                 "Logical Record Length of source Level 1b data set prior to "
                 "processing = %d",
                 nRecordSizeFromHeader);
        int nBlockSize =
            GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_BLOCK_SIZE_OFF);
        CPLDebug(
            "L1B",
            "Block Size of source Level 1b data set prior to processing = %d",
            nBlockSize);
        CPLDebug("L1B", "Count of Header Records in this Data Set = %d",
                 nHeaderRecCount);

        int nDataRecordCount =
            GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_DATA_RECORD_COUNT_OFF);
        CPLDebug("L1B", "Count of Data Records = %d", nDataRecordCount);

        int nCalibratedScanlineCount = GetUInt16(
            abyRecHeader + L1B_NOAA15_HDR_REC_CALIBRATED_SCANLINE_COUNT_OFF);
        CPLDebug("L1B", "Count of Calibrated, Earth Located Scan Lines = %d",
                 nCalibratedScanlineCount);

        int nMissingScanlineCount = GetUInt16(
            abyRecHeader + L1B_NOAA15_HDR_REC_MISSING_SCANLINE_COUNT_OFF);
        CPLDebug("L1B", "Count of Missing Scan Lines = %d",
                 nMissingScanlineCount);
        if (nMissingScanlineCount != 0)
            bExposeMaskBand = TRUE;

        char szEllipsoid[8 + 1];
        memcpy(szEllipsoid, abyRecHeader + L1B_NOAA15_HDR_REC_ELLIPSOID_OFF, 8);
        szEllipsoid[8] = '\0';
        CPLDebug("L1B", "Reference Ellipsoid Model ID = '%s'", szEllipsoid);
        if (EQUAL(szEllipsoid, "WGS-84  "))
        {
            m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
        }
        else if (EQUAL(szEllipsoid, "  GRS 80"))
        {
            m_oGCPSRS.importFromWkt(
                "GEOGCS[\"GRS 1980(IUGG, "
                "1980)\",DATUM[\"unknown\",SPHEROID[\"GRS80\",6378137,298."
                "257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],"
                "UNIT[\"degree\",0.0174532925199433]]");
        }

        // Determine the spacecraft name
        // See
        // http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83132-2.htm
        int iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_ID_OFF);
        switch (iWord)
        {
            case 2:
                eSpacecraftID = NOAA16;
                break;
            case 4:
                eSpacecraftID = NOAA15;
                break;
            case 6:
                eSpacecraftID = NOAA17;
                break;
            case 7:
                eSpacecraftID = NOAA18;
                break;
            case 8:
                eSpacecraftID = NOAA19;
                break;
            case 11:
                eSpacecraftID = METOP1;
                break;
            case 12:
                eSpacecraftID = METOP2;
                break;
            // METOP3 is not documented yet
            case 13:
                eSpacecraftID = METOP3;
                break;
            case 14:
                eSpacecraftID = METOP3;
                break;
            default:
#ifdef DEBUG
                CPLDebug("L1B", "Unknown spacecraft ID \"%d\".", iWord);
#endif
                return CE_Failure;
        }

        // Determine the product data type
        iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_PROD_OFF);
        switch (iWord)
        {
            case 1:
                eProductType = LAC;
                break;
            case 2:
                eProductType = GAC;
                break;
            case 3:
                eProductType = HRPT;
                break;
            case 4:   // XXX: documentation specifies the code '4'
            case 13:  // for FRAC but real datasets contain '13 here.'
                eProductType = FRAC;
                break;
            default:
#ifdef DEBUG
                CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
#endif
                return CE_Failure;
        }

        // Fetch instrument status. Helps to determine whether we have
        // 3A or 3B channel in the dataset.
        iInstrumentStatus =
            GetUInt32(abyRecHeader + L1B_NOAA15_HDR_REC_STAT_OFF);

        // Determine receiving station name
        iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_SRC_OFF);
        switch (iWord)
        {
            case 1:
                eSource = GC;
                break;
            case 2:
                eSource = WI;
                break;
            case 3:
                eSource = SO;
                break;
            case 4:
                eSource = SV;
                break;
            case 5:
                eSource = MO;
                break;
            default:
                eSource = UNKNOWN_STATION;
                break;
        }
    }
    else
        return CE_Failure;

    /* -------------------------------------------------------------------- */
    /*      Set fetched information as metadata records                     */
    /* -------------------------------------------------------------------- */

    SetMetadataItem("DATASET_NAME", szDatasetName);

    const char *pszText = nullptr;
    switch (eSpacecraftID)
    {
        case TIROSN:
            pszText = "TIROS-N";
            break;
        case NOAA6:
            pszText = "NOAA-6(A)";
            break;
        case NOAAB:
            pszText = "NOAA-B";
            break;
        case NOAA7:
            pszText = "NOAA-7(C)";
            break;
        case NOAA8:
            pszText = "NOAA-8(E)";
            break;
        case NOAA9_UNKNOWN:
            pszText = "UNKNOWN";
            break;
        case NOAA9:
            pszText = "NOAA-9(F)";
            break;
        case NOAA10:
            pszText = "NOAA-10(G)";
            break;
        case NOAA11:
            pszText = "NOAA-11(H)";
            break;
        case NOAA12:
            pszText = "NOAA-12(D)";
            break;
        case NOAA13:
            pszText = "NOAA-13(I)";
            break;
        case NOAA14:
            pszText = "NOAA-14(J)";
            break;
        case NOAA15:
            pszText = "NOAA-15(K)";
            break;
        case NOAA16:
            pszText = "NOAA-16(L)";
            break;
        case NOAA17:
            pszText = "NOAA-17(M)";
            break;
        case NOAA18:
            pszText = "NOAA-18(N)";
            break;
        case NOAA19:
            pszText = "NOAA-19(N')";
            break;
        case METOP2:
            pszText = "METOP-A(2)";
            break;
        case METOP1:
            pszText = "METOP-B(1)";
            break;
        case METOP3:
            pszText = "METOP-C(3)";
            break;
        default:
            pszText = "Unknown";
            break;
    }
    SetMetadataItem("SATELLITE", pszText);

    switch (eProductType)
    {
        case LAC:
            pszText = "AVHRR LAC";
            break;
        case HRPT:
            pszText = "AVHRR HRPT";
            break;
        case GAC:
            pszText = "AVHRR GAC";
            break;
        case FRAC:
            pszText = "AVHRR FRAC";
            break;
        default:
            pszText = "Unknown";
            break;
    }
    SetMetadataItem("DATA_TYPE", pszText);

    // Get revolution number as string, we don't need this value for processing
    char szRevolution[6];
    memcpy(szRevolution, szDatasetName + 32, 5);
    szRevolution[5] = '\0';
    SetMetadataItem("REVOLUTION", szRevolution);

    switch (eSource)
    {
        case DU:
            pszText = "Dundee, Scotland, UK";
            break;
        case GC:
            pszText = "Fairbanks, Alaska, USA (formerly Gilmore Creek)";
            break;
        case HO:
            pszText = "Honolulu, Hawaii, USA";
            break;
        case MO:
            pszText = "Monterey, California, USA";
            break;
        case WE:
            pszText = "Western Europe CDA, Lannion, France";
            break;
        case SO:
            pszText = "SOCC (Satellite Operations Control Center), Suitland, "
                      "Maryland, USA";
            break;
        case WI:
            pszText = "Wallops Island, Virginia, USA";
            break;
        default:
            pszText = "Unknown receiving station";
            break;
    }
    SetMetadataItem("SOURCE", pszText);

    switch (eProcCenter)
    {
        case CMS:
            pszText = "Centre de Meteorologie Spatiale - Lannion, France";
            break;
        case DSS:
            pszText =
                "Dundee Satellite Receiving Station - Dundee, Scotland, UK";
            break;
        case NSS:
            pszText = "NOAA/NESDIS - Suitland, Maryland, USA";
            break;
        case UKM:
            pszText =
                "United Kingdom Meteorological Office - Bracknell, England, UK";
            break;
        default:
            pszText = "Unknown processing center";
            break;
    }
    SetMetadataItem("PROCESSING_CENTER", pszText);

    return CE_None;
}

/************************************************************************/
/*                        ComputeFileOffsets()                          */
/************************************************************************/

int L1BDataset::ComputeFileOffsets()
{
    CPLDebug("L1B", "Data format = %s",
             (iDataFormat == PACKED10BIT)     ? "Packed 10 bit"
             : (iDataFormat == UNPACKED16BIT) ? "Unpacked 16 bit"
                                              : "Unpacked 8 bit");

    switch (eProductType)
    {
        case HRPT:
        case LAC:
        case FRAC:
            nRasterXSize = 2048;
            nBufferSize = 20484;
            iGCPStart =
                25 -
                1; /* http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
                    */
            iGCPStep = 40;
            nGCPsPerLine = 51;
            if (eL1BFormat == L1B_NOAA9)
            {
                if (iDataFormat == PACKED10BIT)
                {
                    nRecordSize = 14800;
                    nRecordDataEnd = 14104;
                }
                else if (iDataFormat == UNPACKED16BIT)
                {
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 4544;
                            nRecordDataEnd = 4544;
                            break;
                        case 2:
                            nRecordSize = 8640;
                            nRecordDataEnd = 8640;
                            break;
                        case 3:
                            nRecordSize = 12736;
                            nRecordDataEnd = 12736;
                            break;
                        case 4:
                            nRecordSize = 16832;
                            nRecordDataEnd = 16832;
                            break;
                        case 5:
                            nRecordSize = 20928;
                            nRecordDataEnd = 20928;
                            break;
                    }
                }
                else  // UNPACKED8BIT
                {
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 2496;
                            nRecordDataEnd = 2496;
                            break;
                        case 2:
                            nRecordSize = 4544;
                            nRecordDataEnd = 4544;
                            break;
                        case 3:
                            nRecordSize = 6592;
                            nRecordDataEnd = 6592;
                            break;
                        case 4:
                            nRecordSize = 8640;
                            nRecordDataEnd = 8640;
                            break;
                        case 5:
                            nRecordSize = 10688;
                            nRecordDataEnd = 10688;
                            break;
                    }
                }
                nDataStartOffset = nRecordSize + L1B_NOAA9_HEADER_SIZE;
                nRecordDataStart = 448;
                iGCPCodeOffset = 52;
                iGCPOffset = 104;
            }

            else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
            {
                if (iDataFormat == PACKED10BIT)
                {
                    nRecordSize = 15872;
                    nRecordDataEnd = 14920;
                    iCLAVRStart = 14984;
                }
                else if (iDataFormat == UNPACKED16BIT)
                { /* Table 8.3.1.3.3.1-3 */
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 6144;
                            nRecordDataEnd = 5360;
                            iCLAVRStart =
                                5368 + 56; /* guessed but not verified */
                            break;
                        case 2:
                            nRecordSize = 10240;
                            nRecordDataEnd = 9456;
                            iCLAVRStart =
                                9464 + 56; /* guessed but not verified */
                            break;
                        case 3:
                            nRecordSize = 14336;
                            nRecordDataEnd = 13552;
                            iCLAVRStart =
                                13560 + 56; /* guessed but not verified */
                            break;
                        case 4:
                            nRecordSize = 18432;
                            nRecordDataEnd = 17648;
                            iCLAVRStart =
                                17656 + 56; /* guessed but not verified */
                            break;
                        case 5:
                            nRecordSize = 22528;
                            nRecordDataEnd = 21744;
                            iCLAVRStart = 21752 + 56;
                            break;
                    }
                }
                else  // UNPACKED8BIT
                {     /* Table 8.3.1.3.3.1-2 */
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 4096;
                            nRecordDataEnd = 3312;
                            iCLAVRStart =
                                3320 + 56; /* guessed but not verified */
                            break;
                        case 2:
                            nRecordSize = 6144;
                            nRecordDataEnd = 5360;
                            iCLAVRStart =
                                5368 + 56; /* guessed but not verified */
                            break;
                        case 3:
                            nRecordSize = 8192;
                            nRecordDataEnd = 7408;
                            iCLAVRStart =
                                7416 + 56; /* guessed but not verified */
                            break;
                        case 4:
                            nRecordSize = 10240;
                            nRecordDataEnd = 9456;
                            iCLAVRStart =
                                9464 + 56; /* guessed but not verified */
                            break;
                        case 5:
                            nRecordSize = 12288;
                            nRecordDataEnd = 11504;
                            iCLAVRStart =
                                11512 + 56; /* guessed but not verified */
                            break;
                    }
                }
                nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
                                       ? nRecordDataEnd
                                       : nRecordSize + L1B_NOAA15_HEADER_SIZE;
                nRecordDataStart = 1264;
                iGCPCodeOffset = 0;  // XXX: not exist for NOAA15?
                iGCPOffset = 640;
            }
            else
                return 0;
            break;

        case GAC:
            nRasterXSize = 409;
            nBufferSize = 4092;
            iGCPStart = 5 - 1;  // FIXME: depends of scan direction
            iGCPStep = 8;
            nGCPsPerLine = 51;
            if (eL1BFormat == L1B_NOAA9)
            {
                if (iDataFormat == PACKED10BIT)
                {
                    nRecordSize = 3220;
                    nRecordDataEnd = 3176;
                }
                else if (iDataFormat == UNPACKED16BIT)
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 1268;
                            nRecordDataEnd = 1266;
                            break;
                        case 2:
                            nRecordSize = 2084;
                            nRecordDataEnd = 2084;
                            break;
                        case 3:
                            nRecordSize = 2904;
                            nRecordDataEnd = 2902;
                            break;
                        case 4:
                            nRecordSize = 3720;
                            nRecordDataEnd = 3720;
                            break;
                        case 5:
                            nRecordSize = 4540;
                            nRecordDataEnd = 4538;
                            break;
                    }
                else  // UNPACKED8BIT
                {
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 860;
                            nRecordDataEnd = 858;
                            break;
                        case 2:
                            nRecordSize = 1268;
                            nRecordDataEnd = 1266;
                            break;
                        case 3:
                            nRecordSize = 1676;
                            nRecordDataEnd = 1676;
                            break;
                        case 4:
                            nRecordSize = 2084;
                            nRecordDataEnd = 2084;
                            break;
                        case 5:
                            nRecordSize = 2496;
                            nRecordDataEnd = 2494;
                            break;
                    }
                }
                nDataStartOffset = nRecordSize * 2 + L1B_NOAA9_HEADER_SIZE;
                nRecordDataStart = 448;
                iGCPCodeOffset = 52;
                iGCPOffset = 104;
            }

            else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
            {
                if (iDataFormat == PACKED10BIT)
                {
                    nRecordSize = 4608;
                    nRecordDataEnd = 3992;
                    iCLAVRStart = 4056;
                }
                else if (iDataFormat == UNPACKED16BIT)
                { /* Table 8.3.1.4.3.1-3 */
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 2360;
                            nRecordDataEnd = 2082;
                            iCLAVRStart =
                                2088 + 56; /* guessed but not verified */
                            break;
                        case 2:
                            nRecordSize = 3176;
                            nRecordDataEnd = 2900;
                            iCLAVRStart =
                                2904 + 56; /* guessed but not verified */
                            break;
                        case 3:
                            nRecordSize = 3992;
                            nRecordDataEnd = 3718;
                            iCLAVRStart =
                                3720 + 56; /* guessed but not verified */
                            break;
                        case 4:
                            nRecordSize = 4816;
                            nRecordDataEnd = 4536;
                            iCLAVRStart =
                                4544 + 56; /* guessed but not verified */
                            break;
                        case 5:
                            nRecordSize = 5632;
                            nRecordDataEnd = 5354;
                            iCLAVRStart = 5360 + 56;
                            break;
                    }
                }
                else  // UNPACKED8BIT
                { /* Table 8.3.1.4.3.1-2 but record length is wrong in the table
                     ! */
                    switch (nBands)
                    {
                        case 1:
                            nRecordSize = 1952;
                            nRecordDataEnd = 1673;
                            iCLAVRStart =
                                1680 + 56; /* guessed but not verified */
                            break;
                        case 2:
                            nRecordSize = 2360;
                            nRecordDataEnd = 2082;
                            iCLAVRStart =
                                2088 + 56; /* guessed but not verified */
                            break;
                        case 3:
                            nRecordSize = 2768;
                            nRecordDataEnd = 2491;
                            iCLAVRStart =
                                2496 + 56; /* guessed but not verified */
                            break;
                        case 4:
                            nRecordSize = 3176;
                            nRecordDataEnd = 2900;
                            iCLAVRStart =
                                2904 + 56; /* guessed but not verified */
                            break;
                        case 5:
                            nRecordSize = 3584;
                            nRecordDataEnd = 3309;
                            iCLAVRStart =
                                3312 + 56; /* guessed but not verified */
                            break;
                    }
                }
                nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
                                       ? nRecordDataEnd
                                       : nRecordSize + L1B_NOAA15_HEADER_SIZE;
                nRecordDataStart = 1264;
                iGCPCodeOffset = 0;  // XXX: not exist for NOAA15?
                iGCPOffset = 640;
            }
            else
                return 0;
            break;
        default:
            return 0;
    }

    return 1;
}

/************************************************************************/
/*                       L1BGeolocDataset                               */
/************************************************************************/

class L1BGeolocDataset final : public GDALDataset
{
    friend class L1BGeolocRasterBand;

    L1BDataset *poL1BDS;
    int bInterpolGeolocationDS;

  public:
    L1BGeolocDataset(L1BDataset *poMainDS, int bInterpolGeolocationDS);
    virtual ~L1BGeolocDataset();

    static GDALDataset *CreateGeolocationDS(L1BDataset *poL1BDS,
                                            int bInterpolGeolocationDS);
};

/************************************************************************/
/*                       L1BGeolocRasterBand                            */
/************************************************************************/

class L1BGeolocRasterBand final : public GDALRasterBand
{
  public:
    L1BGeolocRasterBand(L1BGeolocDataset *poDS, int nBand);

    virtual CPLErr IReadBlock(int, int, void *) override;
    virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
};

/************************************************************************/
/*                        L1BGeolocDataset()                            */
/************************************************************************/

L1BGeolocDataset::L1BGeolocDataset(L1BDataset *poL1BDSIn,
                                   int bInterpolGeolocationDSIn)
    : poL1BDS(poL1BDSIn), bInterpolGeolocationDS(bInterpolGeolocationDSIn)
{
    if (bInterpolGeolocationDS)
        nRasterXSize = poL1BDS->nRasterXSize;
    else
        nRasterXSize = poL1BDS->nGCPsPerLine;
    nRasterYSize = poL1BDS->nRasterYSize;
}

/************************************************************************/
/*                       ~L1BGeolocDataset()                            */
/************************************************************************/

L1BGeolocDataset::~L1BGeolocDataset()
{
    delete poL1BDS;
}

/************************************************************************/
/*                        L1BGeolocRasterBand()                         */
/************************************************************************/

L1BGeolocRasterBand::L1BGeolocRasterBand(L1BGeolocDataset *poDSIn, int nBandIn)
{
    poDS = poDSIn;
    nBand = nBandIn;
    nRasterXSize = poDSIn->nRasterXSize;
    nRasterYSize = poDSIn->nRasterYSize;
    eDataType = GDT_Float64;
    nBlockXSize = nRasterXSize;
    nBlockYSize = 1;
    if (nBand == 1)
        SetDescription("GEOLOC X");
    else
        SetDescription("GEOLOC Y");
}

/************************************************************************/
/*                         LagrangeInterpol()                           */
/************************************************************************/

/* ----------------------------------------------------------------------------
 * Perform a Lagrangian interpolation through the given x,y coordinates
 * and return the interpolated y value for the given x value.
 * The array size and thus the polynomial order is defined by numpt.
 * Input: x[] and y[] are of size numpt,
 *  x0 is the x value for which we calculate the corresponding y
 * Returns: y value calculated for given x0.
 */
static double LagrangeInterpol(const double x[], const double y[], double x0,
                               int numpt)
{
    int i, j;
    double L;
    double y0 = 0;

    for (i = 0; i < numpt; i++)
    {
        L = 1.0;
        for (j = 0; j < numpt; j++)
        {
            if (i == j)
                continue;
            L = L * (x0 - x[j]) / (x[i] - x[j]);
        }
        y0 = y0 + L * y[i];
    }
    return y0;
}

/************************************************************************/
/*                         L1BInterpol()                                */
/************************************************************************/

/* ----------------------------------------------------------------------------
 * Interpolate an array of size numPoints where the only values set on input are
 * at knownFirst, and intervals of knownStep thereafter.
 * On return all the rest from 0..numPoints-1 will be filled in.
 * Uses the LagrangeInterpol() function to do the interpolation; 5-point for the
 * beginning and end of the array and 4-point for the rest.
 * To use this function for NOAA level 1B data extract the 51 latitude values
 * into their appropriate places in the vals array then call L1BInterpol to
 * calculate the rest of the values.  Do similarly for longitudes, solar zenith
 * angles, and any others which are present in the file.
 * Reference:
 *  http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
 */

#define MIDDLE_INTERP_ORDER 4
#define END_INTERP_ORDER 5 /* Ensure this is an odd number, 5 is suitable.*/

/* Convert number of known point to its index in full array */
#define IDX(N) ((N)*knownStep + knownFirst)

static void L1BInterpol(
    double vals[], int numKnown, /* Number of known points (typically 51) */
    int knownFirst, /* Index in full array of first known point (24) */
    int knownStep,  /* Interval to next and subsequent known points (40) */
    int numPoints /* Number of points in whole array (2048) */)
{
    int i, j;
    double x[END_INTERP_ORDER];
    double y[END_INTERP_ORDER];

    /* First extrapolate first 24 points */
    for (i = 0; i < END_INTERP_ORDER; i++)
    {
        x[i] = IDX(i);
        y[i] = vals[IDX(i)];
    }
    for (i = 0; i < knownFirst; i++)
    {
        vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
    }

    /* Next extrapolate last 23 points */
    for (i = 0; i < END_INTERP_ORDER; i++)
    {
        x[i] = IDX(numKnown - END_INTERP_ORDER + i);
        y[i] = vals[IDX(numKnown - END_INTERP_ORDER + i)];
    }
    for (i = IDX(numKnown - 1); i < numPoints; i++)
    {
        vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
    }

    /* Interpolate all intermediate points using two before and two after */
    for (i = knownFirst; i < IDX(numKnown - 1); i++)
    {
        double x2[MIDDLE_INTERP_ORDER];
        double y2[MIDDLE_INTERP_ORDER];
        int startpt;

        /* Find a suitable set of two known points before and two after */
        startpt = (i / knownStep) - MIDDLE_INTERP_ORDER / 2;
        if (startpt < 0)
            startpt = 0;
        if (startpt + MIDDLE_INTERP_ORDER - 1 >= numKnown)
            startpt = numKnown - MIDDLE_INTERP_ORDER;
        for (j = 0; j < MIDDLE_INTERP_ORDER; j++)
        {
            x2[j] = IDX(startpt + j);
            y2[j] = vals[IDX(startpt + j)];
        }
        vals[i] = LagrangeInterpol(x2, y2, i, MIDDLE_INTERP_ORDER);
    }
}

/************************************************************************/
/*                         IReadBlock()                                 */
/************************************************************************/

CPLErr L1BGeolocRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
                                       int nBlockYOff, void *pData)
{
    L1BGeolocDataset *poGDS = (L1BGeolocDataset *)poDS;
    L1BDataset *poL1BDS = poGDS->poL1BDS;
    GDAL_GCP *pasGCPList =
        (GDAL_GCP *)CPLCalloc(poL1BDS->nGCPsPerLine, sizeof(GDAL_GCP));
    GDALInitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);

    /* -------------------------------------------------------------------- */
    /*      Seek to data.                                                   */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));

    CPL_IGNORE_RET_VAL(
        VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordDataStart, poL1BDS->fp));

    /* Fetch the GCPs for the row */
    int nGotGCPs = poL1BDS->FetchGCPs(pasGCPList, pabyRecordHeader, nBlockYOff);
    double *padfData = (double *)pData;
    int i;
    if (poGDS->bInterpolGeolocationDS)
    {
        /* Fill the known position */
        for (i = 0; i < nGotGCPs; i++)
        {
            double dfVal =
                (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
            padfData[poL1BDS->iGCPStart + i * poL1BDS->iGCPStep] = dfVal;
        }

        if (nGotGCPs == poL1BDS->nGCPsPerLine)
        {
            /* And do Lagangian interpolation to fill the holes */
            L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
                        poL1BDS->iGCPStep, nRasterXSize);
        }
        else
        {
            int iFirstNonValid = 0;
            if (nGotGCPs > 5)
                iFirstNonValid = poL1BDS->iGCPStart +
                                 nGotGCPs * poL1BDS->iGCPStep +
                                 poL1BDS->iGCPStep / 2;
            for (i = iFirstNonValid; i < nRasterXSize; i++)
            {
                padfData[i] = GetNoDataValue(nullptr);
            }
            if (iFirstNonValid > 0)
            {
                L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
                            poL1BDS->iGCPStep, iFirstNonValid);
            }
        }
    }
    else
    {
        for (i = 0; i < nGotGCPs; i++)
        {
            padfData[i] =
                (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
        }
        for (i = nGotGCPs; i < nRasterXSize; i++)
            padfData[i] = GetNoDataValue(nullptr);
    }

    if (poL1BDS->eLocationIndicator == ASCEND)
    {
        for (i = 0; i < nRasterXSize / 2; i++)
        {
            double dfTmp = padfData[i];
            padfData[i] = padfData[nRasterXSize - 1 - i];
            padfData[nRasterXSize - 1 - i] = dfTmp;
        }
    }

    CPLFree(pabyRecordHeader);
    GDALDeinitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);
    CPLFree(pasGCPList);

    return CE_None;
}

/************************************************************************/
/*                        GetNoDataValue()                              */
/************************************************************************/

double L1BGeolocRasterBand::GetNoDataValue(int *pbSuccess)
{
    if (pbSuccess)
        *pbSuccess = TRUE;
    return -200.0;
}

/************************************************************************/
/*                      CreateGeolocationDS()                           */
/************************************************************************/

GDALDataset *L1BGeolocDataset::CreateGeolocationDS(L1BDataset *poL1BDS,
                                                   int bInterpolGeolocationDS)
{
    L1BGeolocDataset *poGeolocDS =
        new L1BGeolocDataset(poL1BDS, bInterpolGeolocationDS);
    for (int i = 1; i <= 2; i++)
    {
        poGeolocDS->SetBand(i, new L1BGeolocRasterBand(poGeolocDS, i));
    }
    return poGeolocDS;
}

/************************************************************************/
/*                    L1BSolarZenithAnglesDataset                       */
/************************************************************************/

class L1BSolarZenithAnglesDataset final : public GDALDataset
{
    friend class L1BSolarZenithAnglesRasterBand;

    L1BDataset *poL1BDS;

  public:
    explicit L1BSolarZenithAnglesDataset(L1BDataset *poMainDS);
    virtual ~L1BSolarZenithAnglesDataset();

    static GDALDataset *CreateSolarZenithAnglesDS(L1BDataset *poL1BDS);
};

/************************************************************************/
/*                  L1BSolarZenithAnglesRasterBand                      */
/************************************************************************/

class L1BSolarZenithAnglesRasterBand final : public GDALRasterBand
{
  public:
    L1BSolarZenithAnglesRasterBand(L1BSolarZenithAnglesDataset *poDS,
                                   int nBand);

    virtual CPLErr IReadBlock(int, int, void *) override;
    virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
};

/************************************************************************/
/*                  L1BSolarZenithAnglesDataset()                       */
/************************************************************************/

L1BSolarZenithAnglesDataset::L1BSolarZenithAnglesDataset(L1BDataset *poL1BDSIn)
{
    poL1BDS = poL1BDSIn;
    nRasterXSize = 51;
    nRasterYSize = poL1BDSIn->nRasterYSize;
}

/************************************************************************/
/*                  ~L1BSolarZenithAnglesDataset()                      */
/************************************************************************/

L1BSolarZenithAnglesDataset::~L1BSolarZenithAnglesDataset()
{
    delete poL1BDS;
}

/************************************************************************/
/*                  L1BSolarZenithAnglesRasterBand()                    */
/************************************************************************/

L1BSolarZenithAnglesRasterBand::L1BSolarZenithAnglesRasterBand(
    L1BSolarZenithAnglesDataset *poDSIn, int nBandIn)
{
    poDS = poDSIn;
    nBand = nBandIn;
    nRasterXSize = poDSIn->nRasterXSize;
    nRasterYSize = poDSIn->nRasterYSize;
    eDataType = GDT_Float32;
    nBlockXSize = nRasterXSize;
    nBlockYSize = 1;
}

/************************************************************************/
/*                         IReadBlock()                                 */
/************************************************************************/

CPLErr L1BSolarZenithAnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
                                                  int nBlockYOff, void *pData)
{
    L1BSolarZenithAnglesDataset *poGDS = (L1BSolarZenithAnglesDataset *)poDS;
    L1BDataset *poL1BDS = poGDS->poL1BDS;
    int i;

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);

    /* -------------------------------------------------------------------- */
    /*      Seek to data.                                                   */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));

    CPL_IGNORE_RET_VAL(
        VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));

    const int nValidValues =
        std::min(nRasterXSize,
                 static_cast<int>(pabyRecordHeader[poL1BDS->iGCPCodeOffset]));
    float *pafData = (float *)pData;

    int bHasFractional = (poL1BDS->nRecordDataEnd + 20 <= poL1BDS->nRecordSize);

#ifdef notdef
    if (bHasFractional)
    {
        for (i = 0; i < 20; i++)
        {
            GByte val = pabyRecordHeader[poL1BDS->nRecordDataEnd + i];
            for (int j = 0; j < 8; j++)
                fprintf(stderr, "%c", /*ok*/
                        ((val >> (7 - j)) & 1) ? '1' : '0');
            fprintf(stderr, " "); /*ok*/
        }
        fprintf(stderr, "\n"); /*ok*/
    }
#endif

    for (i = 0; i < nValidValues; i++)
    {
        pafData[i] = pabyRecordHeader[poL1BDS->iGCPCodeOffset + 1 + i] / 2.0f;

        if (bHasFractional)
        {
            /* Cf
             * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/l/app-l.htm#notl-2
             */
            /* This is not very clear on if bits must be counted from MSB or LSB
             */
            /* but when testing on n12gac10bit.l1b, it appears that the 3 bits
             * for i=0 are the 3 MSB bits */
            int nAddBitStart = i * 3;
            int nFractional;

#if 1
            if ((nAddBitStart % 8) + 3 <= 8)
            {
                nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
                                                (nAddBitStart / 8)] >>
                               (8 - ((nAddBitStart % 8) + 3))) &
                              0x7;
            }
            else
            {
                nFractional = (((pabyRecordHeader[poL1BDS->nRecordDataEnd +
                                                  (nAddBitStart / 8)]
                                 << 8) |
                                pabyRecordHeader[poL1BDS->nRecordDataEnd +
                                                 (nAddBitStart / 8) + 1]) >>
                               (16 - ((nAddBitStart % 8) + 3))) &
                              0x7;
            }
#else
            nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
                                            (nAddBitStart / 8)] >>
                           (nAddBitStart % 8)) &
                          0x7;
            if ((nAddBitStart % 8) + 3 > 8)
                nFractional |= (pabyRecordHeader[poL1BDS->nRecordDataEnd +
                                                 (nAddBitStart / 8) + 1] &
                                ((1 << (((nAddBitStart % 8) + 3 - 8))) - 1))
                               << (3 - ((((nAddBitStart % 8) + 3 - 8))));
            * /
#endif
            if (nFractional > 4)
            {
                CPLDebug("L1B",
                         "For nBlockYOff=%d, i=%d, wrong fractional value : %d",
                         nBlockYOff, i, nFractional);
            }

            pafData[i] += nFractional / 10.0f;
        }
    }

    for (; i < nRasterXSize; i++)
        pafData[i] = static_cast<float>(GetNoDataValue(nullptr));

    if (poL1BDS->eLocationIndicator == ASCEND)
    {
        for (i = 0; i < nRasterXSize / 2; i++)
        {
            float fTmp = pafData[i];
            pafData[i] = pafData[nRasterXSize - 1 - i];
            pafData[nRasterXSize - 1 - i] = fTmp;
        }
    }

    CPLFree(pabyRecordHeader);

    return CE_None;
}

/************************************************************************/
/*                        GetNoDataValue()                              */
/************************************************************************/

double L1BSolarZenithAnglesRasterBand::GetNoDataValue(int *pbSuccess)
{
    if (pbSuccess)
        *pbSuccess = TRUE;
    return -200.0;
}

/************************************************************************/
/*                      CreateSolarZenithAnglesDS()                     */
/************************************************************************/

GDALDataset *
L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(L1BDataset *poL1BDS)
{
    L1BSolarZenithAnglesDataset *poGeolocDS =
        new L1BSolarZenithAnglesDataset(poL1BDS);
    for (int i = 1; i <= 1; i++)
    {
        poGeolocDS->SetBand(i,
                            new L1BSolarZenithAnglesRasterBand(poGeolocDS, i));
    }
    return poGeolocDS;
}

/************************************************************************/
/*                     L1BNOAA15AnglesDataset                           */
/************************************************************************/

class L1BNOAA15AnglesDataset final : public GDALDataset
{
    friend class L1BNOAA15AnglesRasterBand;

    L1BDataset *poL1BDS;

  public:
    explicit L1BNOAA15AnglesDataset(L1BDataset *poMainDS);
    virtual ~L1BNOAA15AnglesDataset();

    static GDALDataset *CreateAnglesDS(L1BDataset *poL1BDS);
};

/************************************************************************/
/*                     L1BNOAA15AnglesRasterBand                        */
/************************************************************************/

class L1BNOAA15AnglesRasterBand final : public GDALRasterBand
{
  public:
    L1BNOAA15AnglesRasterBand(L1BNOAA15AnglesDataset *poDS, int nBand);

    virtual CPLErr IReadBlock(int, int, void *) override;
};

/************************************************************************/
/*                       L1BNOAA15AnglesDataset()                       */
/************************************************************************/

L1BNOAA15AnglesDataset::L1BNOAA15AnglesDataset(L1BDataset *poL1BDSIn)
    : poL1BDS(poL1BDSIn)
{
    nRasterXSize = 51;
    nRasterYSize = poL1BDS->nRasterYSize;
}

/************************************************************************/
/*                     ~L1BNOAA15AnglesDataset()                        */
/************************************************************************/

L1BNOAA15AnglesDataset::~L1BNOAA15AnglesDataset()
{
    delete poL1BDS;
}

/************************************************************************/
/*                      L1BNOAA15AnglesRasterBand()                     */
/************************************************************************/

L1BNOAA15AnglesRasterBand::L1BNOAA15AnglesRasterBand(
    L1BNOAA15AnglesDataset *poDSIn, int nBandIn)
{
    poDS = poDSIn;
    nBand = nBandIn;
    nRasterXSize = poDSIn->nRasterXSize;
    nRasterYSize = poDSIn->nRasterYSize;
    eDataType = GDT_Float32;
    nBlockXSize = nRasterXSize;
    nBlockYSize = 1;
    if (nBand == 1)
        SetDescription("Solar zenith angles");
    else if (nBand == 2)
        SetDescription("Satellite zenith angles");
    else
        SetDescription("Relative azimuth angles");
}

/************************************************************************/
/*                         IReadBlock()                                 */
/************************************************************************/

CPLErr L1BNOAA15AnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
                                             int nBlockYOff, void *pData)
{
    L1BNOAA15AnglesDataset *poGDS = (L1BNOAA15AnglesDataset *)poDS;
    L1BDataset *poL1BDS = poGDS->poL1BDS;
    int i;

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);

    /* -------------------------------------------------------------------- */
    /*      Seek to data.                                                   */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));

    CPL_IGNORE_RET_VAL(
        VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));

    float *pafData = (float *)pData;

    for (i = 0; i < nRasterXSize; i++)
    {
        GInt16 i16 =
            poL1BDS->GetInt16(pabyRecordHeader + 328 + 6 * i + 2 * (nBand - 1));
        pafData[i] = i16 / 100.0f;
    }

    if (poL1BDS->eLocationIndicator == ASCEND)
    {
        for (i = 0; i < nRasterXSize / 2; i++)
        {
            float fTmp = pafData[i];
            pafData[i] = pafData[nRasterXSize - 1 - i];
            pafData[nRasterXSize - 1 - i] = fTmp;
        }
    }

    CPLFree(pabyRecordHeader);

    return CE_None;
}

/************************************************************************/
/*                            CreateAnglesDS()                          */
/************************************************************************/

GDALDataset *L1BNOAA15AnglesDataset::CreateAnglesDS(L1BDataset *poL1BDS)
{
    L1BNOAA15AnglesDataset *poGeolocDS = new L1BNOAA15AnglesDataset(poL1BDS);
    for (int i = 1; i <= 3; i++)
    {
        poGeolocDS->SetBand(i, new L1BNOAA15AnglesRasterBand(poGeolocDS, i));
    }
    return poGeolocDS;
}

/************************************************************************/
/*                          L1BCloudsDataset                            */
/************************************************************************/

class L1BCloudsDataset final : public GDALDataset
{
    friend class L1BCloudsRasterBand;

    L1BDataset *poL1BDS;

  public:
    explicit L1BCloudsDataset(L1BDataset *poMainDS);
    virtual ~L1BCloudsDataset();

    static GDALDataset *CreateCloudsDS(L1BDataset *poL1BDS);
};

/************************************************************************/
/*                        L1BCloudsRasterBand                           */
/************************************************************************/

class L1BCloudsRasterBand final : public GDALRasterBand
{
  public:
    L1BCloudsRasterBand(L1BCloudsDataset *poDS, int nBand);

    virtual CPLErr IReadBlock(int, int, void *) override;
};

/************************************************************************/
/*                         L1BCloudsDataset()                           */
/************************************************************************/

L1BCloudsDataset::L1BCloudsDataset(L1BDataset *poL1BDSIn) : poL1BDS(poL1BDSIn)
{
    nRasterXSize = poL1BDSIn->nRasterXSize;
    nRasterYSize = poL1BDSIn->nRasterYSize;
}

/************************************************************************/
/*                        ~L1BCloudsDataset()                           */
/************************************************************************/

L1BCloudsDataset::~L1BCloudsDataset()
{
    delete poL1BDS;
}

/************************************************************************/
/*                         L1BCloudsRasterBand()                        */
/************************************************************************/

L1BCloudsRasterBand::L1BCloudsRasterBand(L1BCloudsDataset *poDSIn, int nBandIn)
{
    poDS = poDSIn;
    nBand = nBandIn;
    nRasterXSize = poDSIn->nRasterXSize;
    nRasterYSize = poDSIn->nRasterYSize;
    eDataType = GDT_Byte;
    nBlockXSize = nRasterXSize;
    nBlockYSize = 1;
}

/************************************************************************/
/*                         IReadBlock()                                 */
/************************************************************************/

CPLErr L1BCloudsRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
                                       int nBlockYOff, void *pData)
{
    L1BCloudsDataset *poGDS = (L1BCloudsDataset *)poDS;
    L1BDataset *poL1BDS = poGDS->poL1BDS;
    int i;

    GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);

    /* -------------------------------------------------------------------- */
    /*      Seek to data.                                                   */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(
        VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));

    CPL_IGNORE_RET_VAL(
        VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));

    GByte *pabyData = (GByte *)pData;

    for (i = 0; i < nRasterXSize; i++)
    {
        pabyData[i] = ((pabyRecordHeader[poL1BDS->iCLAVRStart + (i / 4)] >>
                        (8 - ((i % 4) * 2 + 2))) &
                       0x3);
    }

    if (poL1BDS->eLocationIndicator == ASCEND)
    {
        for (i = 0; i < nRasterXSize / 2; i++)
        {
            GByte byTmp = pabyData[i];
            pabyData[i] = pabyData[nRasterXSize - 1 - i];
            pabyData[nRasterXSize - 1 - i] = byTmp;
        }
    }

    CPLFree(pabyRecordHeader);

    return CE_None;
}

/************************************************************************/
/*                      CreateCloudsDS()                     */
/************************************************************************/

GDALDataset *L1BCloudsDataset::CreateCloudsDS(L1BDataset *poL1BDS)
{
    L1BCloudsDataset *poGeolocDS = new L1BCloudsDataset(poL1BDS);
    for (int i = 1; i <= 1; i++)
    {
        poGeolocDS->SetBand(i, new L1BCloudsRasterBand(poGeolocDS, i));
    }
    return poGeolocDS;
}

/************************************************************************/
/*                           DetectFormat()                             */
/************************************************************************/

L1BFileFormat L1BDataset::DetectFormat(const char *pszFilename,
                                       const GByte *pabyHeader,
                                       int nHeaderBytes)

{
    if (pabyHeader == nullptr || nHeaderBytes < L1B_NOAA9_HEADER_SIZE)
        return L1B_NONE;

    // try NOAA-18 formats
    if (*(pabyHeader + 0) == '\0' && *(pabyHeader + 1) == '\0' &&
        *(pabyHeader + 2) == '\0' && *(pabyHeader + 3) == '\0' &&
        *(pabyHeader + 4) == '\0' && *(pabyHeader + 5) == '\0' &&
        EQUALN((const char *)(pabyHeader + 22), "/N1BD/N18/", 10))
        return L1B_NOAA15_NOHDR;

    // We will try the NOAA-15 and later formats first
    if (nHeaderBytes > L1B_NOAA15_HEADER_SIZE + 61 &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 25) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 30) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 33) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 40) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 46) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 52) == '.' &&
        *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 61) == '.')
        return L1B_NOAA15;

    // Next try the NOAA-9/14 formats
    if (*(pabyHeader + 8 + 25) == '.' && *(pabyHeader + 8 + 30) == '.' &&
        *(pabyHeader + 8 + 33) == '.' && *(pabyHeader + 8 + 40) == '.' &&
        *(pabyHeader + 8 + 46) == '.' && *(pabyHeader + 8 + 52) == '.' &&
        *(pabyHeader + 8 + 61) == '.')
        return L1B_NOAA9;

    // Next try the NOAA-9/14 formats with dataset name in EBCDIC
    if (*(pabyHeader + 8 + 25) == 'K' && *(pabyHeader + 8 + 30) == 'K' &&
        *(pabyHeader + 8 + 33) == 'K' && *(pabyHeader + 8 + 40) == 'K' &&
        *(pabyHeader + 8 + 46) == 'K' && *(pabyHeader + 8 + 52) == 'K' &&
        *(pabyHeader + 8 + 61) == 'K')
        return L1B_NOAA9;

    // Finally try the AAPP formats
    if (*(pabyHeader + 25) == '.' && *(pabyHeader + 30) == '.' &&
        *(pabyHeader + 33) == '.' && *(pabyHeader + 40) == '.' &&
        *(pabyHeader + 46) == '.' && *(pabyHeader + 52) == '.' &&
        *(pabyHeader + 61) == '.')
        return L1B_NOAA15_NOHDR;

    // A few NOAA <= 9 datasets with no dataset name in TBM header
    if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE && pszFilename[3] == '.' &&
        pszFilename[8] == '.' && pszFilename[11] == '.' &&
        pszFilename[18] == '.' && pszFilename[24] == '.' &&
        pszFilename[30] == '.' && pszFilename[39] == '.' &&
        memcmp(pabyHeader + 30,
               "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0"
               "\0\0\0\0\0\0\0\0\0\0\0",
               L1B_DATASET_NAME_SIZE) == 0 &&
        (pabyHeader[75] == '+' || pabyHeader[75] == '-') &&
        (pabyHeader[78] == '+' || pabyHeader[78] == '-') &&
        (pabyHeader[81] == '+' || pabyHeader[81] == '-') &&
        (pabyHeader[85] == '+' || pabyHeader[85] == '-'))
        return L1B_NOAA9;

    return L1B_NONE;
}

/************************************************************************/
/*                              Identify()                              */
/************************************************************************/

int L1BDataset::Identify(GDALOpenInfo *poOpenInfo)

{
    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
        return TRUE;
    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
        return TRUE;
    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:"))
        return TRUE;
    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
        return TRUE;
    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
        return TRUE;

    if (DetectFormat(CPLGetFilename(poOpenInfo->pszFilename),
                     poOpenInfo->pabyHeader,
                     poOpenInfo->nHeaderBytes) == L1B_NONE)
        return FALSE;

    return TRUE;
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *L1BDataset::Open(GDALOpenInfo *poOpenInfo)

{
    GDALDataset *poOutDS = nullptr;
    VSILFILE *fp = nullptr;
    CPLString osFilename = poOpenInfo->pszFilename;
    int bAskGeolocationDS = FALSE;
    int bInterpolGeolocationDS = FALSE;
    int bAskSolarZenithAnglesDS = FALSE;
    int bAskAnglesDS = FALSE;
    int bAskCloudsDS = FALSE;
    L1BFileFormat eL1BFormat;

    if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:") ||
        STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:") ||
        STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:") ||
        STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:") ||
        STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
    {
        GByte abyHeader[1024];
        const char *pszFilename = nullptr;
        if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
        {
            bAskGeolocationDS = TRUE;
            bInterpolGeolocationDS = TRUE;
            pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS_INTERPOL:");
        }
        else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
        {
            bAskGeolocationDS = TRUE;
            pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS:");
        }
        else if (STARTS_WITH_CI(poOpenInfo->pszFilename,
                                "L1B_SOLAR_ZENITH_ANGLES:"))
        {
            bAskSolarZenithAnglesDS = TRUE;
            pszFilename =
                poOpenInfo->pszFilename + strlen("L1B_SOLAR_ZENITH_ANGLES:");
        }
        else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
        {
            bAskAnglesDS = TRUE;
            pszFilename = poOpenInfo->pszFilename + strlen("L1B_ANGLES:");
        }
        else
        {
            bAskCloudsDS = TRUE;
            pszFilename = poOpenInfo->pszFilename + strlen("L1B_CLOUDS:");
        }
        if (pszFilename[0] == '"')
            pszFilename++;
        osFilename = pszFilename;
        if (!osFilename.empty() && osFilename.back() == '"')
            osFilename.pop_back();
        fp = VSIFOpenL(osFilename, "rb");
        if (!fp)
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Can't open file \"%s\".",
                     osFilename.c_str());
            return nullptr;
        }
        CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 1, sizeof(abyHeader) - 1, fp));
        abyHeader[sizeof(abyHeader) - 1] = '\0';
        eL1BFormat = DetectFormat(CPLGetFilename(osFilename), abyHeader,
                                  sizeof(abyHeader));
    }
    else
        eL1BFormat =
            DetectFormat(CPLGetFilename(osFilename), poOpenInfo->pabyHeader,
                         poOpenInfo->nHeaderBytes);

    if (eL1BFormat == L1B_NONE)
    {
        if (fp != nullptr)
            CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Confirm the requested access is supported.                      */
    /* -------------------------------------------------------------------- */
    if (poOpenInfo->eAccess == GA_Update)
    {
        ReportUpdateNotSupportedByDriver("L1B");
        if (fp != nullptr)
            CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */
    L1BDataset *poDS;
    VSIStatBufL sStat;

    poDS = new L1BDataset(eL1BFormat);

    if (fp == nullptr)
        fp = VSIFOpenL(osFilename, "rb");
    poDS->fp = fp;
    if (!poDS->fp || VSIStatL(osFilename, &sStat) != 0)
    {
        CPLDebug("L1B", "Can't open file \"%s\".", osFilename.c_str());
        goto bad;
    }

    /* -------------------------------------------------------------------- */
    /*      Read the header.                                                */
    /* -------------------------------------------------------------------- */
    if (poDS->ProcessDatasetHeader(CPLGetFilename(osFilename)) != CE_None)
    {
        CPLDebug("L1B", "Error reading L1B record header.");
        goto bad;
    }

    if (poDS->eL1BFormat == L1B_NOAA15_NOHDR &&
        poDS->nRecordSizeFromHeader == 22016 &&
        (sStat.st_size % 22016 /* poDS->nRecordSizeFromHeader*/) == 0)
    {
        poDS->iDataFormat = UNPACKED16BIT;
        poDS->ComputeFileOffsets();
        poDS->nDataStartOffset = poDS->nRecordSizeFromHeader;
        poDS->nRecordSize = poDS->nRecordSizeFromHeader;
        poDS->iCLAVRStart = 0;
    }
    else if (poDS->bGuessDataFormat)
    {
        int nTempYSize;
        GUInt16 nScanlineNumber;
        int j;

        /* If the data format is unspecified, try each one of the 3 known data
         * formats */
        /* It is considered valid when the spacing between the first 5 scanline
         * numbers */
        /* is a constant */

        for (j = 0; j < 3; j++)
        {
            poDS->iDataFormat = (L1BDataFormat)(PACKED10BIT + j);
            if (!poDS->ComputeFileOffsets())
                goto bad;

            nTempYSize = static_cast<int>(
                (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);
            if (nTempYSize < 5)
                continue;

            int nLastScanlineNumber = 0;
            int nDiffLine = 0;
            int i;
            for (i = 0; i < 5; i++)
            {
                nScanlineNumber = 0;

                CPL_IGNORE_RET_VAL(VSIFSeekL(
                    poDS->fp, poDS->nDataStartOffset + i * poDS->nRecordSize,
                    SEEK_SET));
                CPL_IGNORE_RET_VAL(VSIFReadL(&nScanlineNumber, 1, 2, poDS->fp));
                nScanlineNumber = poDS->GetUInt16(&nScanlineNumber);

                if (i == 1)
                {
                    nDiffLine = nScanlineNumber - nLastScanlineNumber;
                    if (nDiffLine == 0)
                        break;
                }
                else if (i > 1)
                {
                    if (nDiffLine != nScanlineNumber - nLastScanlineNumber)
                        break;
                }

                nLastScanlineNumber = nScanlineNumber;
            }

            if (i == 5)
            {
                CPLDebug("L1B", "Guessed data format : %s",
                         (poDS->iDataFormat == PACKED10BIT)    ? "10"
                         : (poDS->iDataFormat == UNPACKED8BIT) ? "08"
                                                               : "16");
                break;
            }
        }

        if (j == 3)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Could not guess data format of L1B product");
            goto bad;
        }
    }
    else
    {
        if (!poDS->ComputeFileOffsets())
            goto bad;
    }

    CPLDebug("L1B", "nRecordDataStart = %d", poDS->nRecordDataStart);
    CPLDebug("L1B", "nRecordDataEnd = %d", poDS->nRecordDataEnd);
    CPLDebug("L1B", "nDataStartOffset = %d", poDS->nDataStartOffset);
    CPLDebug("L1B", "iCLAVRStart = %d", poDS->iCLAVRStart);
    CPLDebug("L1B", "nRecordSize = %d", poDS->nRecordSize);

    // Compute number of lines dynamically, so we can read partially
    // downloaded files.
    if (poDS->nDataStartOffset > sStat.st_size)
        goto bad;
    poDS->nRasterYSize = static_cast<int>(
        (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);

    /* -------------------------------------------------------------------- */
    /*      Deal with GCPs                                                  */
    /* -------------------------------------------------------------------- */
    poDS->ProcessRecordHeaders();

    if (bAskGeolocationDS)
    {
        return L1BGeolocDataset::CreateGeolocationDS(poDS,
                                                     bInterpolGeolocationDS);
    }
    else if (bAskSolarZenithAnglesDS)
    {
        if (eL1BFormat == L1B_NOAA9)
            return L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(poDS);
        else
        {
            delete poDS;
            return nullptr;
        }
    }
    else if (bAskAnglesDS)
    {
        if (eL1BFormat != L1B_NOAA9)
            return L1BNOAA15AnglesDataset::CreateAnglesDS(poDS);
        else
        {
            delete poDS;
            return nullptr;
        }
    }
    else if (bAskCloudsDS)
    {
        if (poDS->iCLAVRStart > 0)
            poOutDS = L1BCloudsDataset::CreateCloudsDS(poDS);
        else
        {
            delete poDS;
            return nullptr;
        }
    }
    else
    {
        poOutDS = poDS;
    }

    {
        CPLString osTMP;
        int bInterpol =
            CPLTestBool(CPLGetConfigOption("L1B_INTERPOL_GCPS", "TRUE"));

        char *pszWKT = nullptr;
        poDS->m_oGCPSRS.exportToWkt(&pszWKT);
        poOutDS->SetMetadataItem("SRS", pszWKT,
                                 "GEOLOCATION"); /* unused by gdalgeoloc.cpp */
        CPLFree(pszWKT);

        if (bInterpol)
            osTMP.Printf("L1BGCPS_INTERPOL:\"%s\"", osFilename.c_str());
        else
            osTMP.Printf("L1BGCPS:\"%s\"", osFilename.c_str());
        poOutDS->SetMetadataItem("X_DATASET", osTMP, "GEOLOCATION");
        poOutDS->SetMetadataItem("X_BAND", "1", "GEOLOCATION");
        poOutDS->SetMetadataItem("Y_DATASET", osTMP, "GEOLOCATION");
        poOutDS->SetMetadataItem("Y_BAND", "2", "GEOLOCATION");

        if (bInterpol)
        {
            poOutDS->SetMetadataItem("PIXEL_OFFSET", "0", "GEOLOCATION");
            poOutDS->SetMetadataItem("PIXEL_STEP", "1", "GEOLOCATION");
        }
        else
        {
            osTMP.Printf("%d", poDS->iGCPStart);
            poOutDS->SetMetadataItem("PIXEL_OFFSET", osTMP, "GEOLOCATION");
            osTMP.Printf("%d", poDS->iGCPStep);
            poOutDS->SetMetadataItem("PIXEL_STEP", osTMP, "GEOLOCATION");
        }

        poOutDS->SetMetadataItem("LINE_OFFSET", "0", "GEOLOCATION");
        poOutDS->SetMetadataItem("LINE_STEP", "1", "GEOLOCATION");
    }

    if (poOutDS != poDS)
        return poOutDS;

    if (eL1BFormat == L1B_NOAA9)
    {
        char **papszSubdatasets = nullptr;
        papszSubdatasets = CSLSetNameValue(
            papszSubdatasets, "SUBDATASET_1_NAME",
            CPLSPrintf("L1B_SOLAR_ZENITH_ANGLES:\"%s\"", osFilename.c_str()));
        papszSubdatasets = CSLSetNameValue(
            papszSubdatasets, "SUBDATASET_1_DESC", "Solar zenith angles");
        poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
        CSLDestroy(papszSubdatasets);
    }
    else
    {
        char **papszSubdatasets = nullptr;
        papszSubdatasets = CSLSetNameValue(
            papszSubdatasets, "SUBDATASET_1_NAME",
            CPLSPrintf("L1B_ANGLES:\"%s\"", osFilename.c_str()));
        papszSubdatasets =
            CSLSetNameValue(papszSubdatasets, "SUBDATASET_1_DESC",
                            "Solar zenith angles, satellite zenith angles and "
                            "relative azimuth angles");

        if (poDS->iCLAVRStart > 0)
        {
            papszSubdatasets = CSLSetNameValue(
                papszSubdatasets, "SUBDATASET_2_NAME",
                CPLSPrintf("L1B_CLOUDS:\"%s\"", osFilename.c_str()));
            papszSubdatasets =
                CSLSetNameValue(papszSubdatasets, "SUBDATASET_2_DESC",
                                "Clouds from AVHRR (CLAVR)");
        }

        poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
        CSLDestroy(papszSubdatasets);
    }

    /* -------------------------------------------------------------------- */
    /*      Create band information objects.                                */
    /* -------------------------------------------------------------------- */
    int iBand, i;

    for (iBand = 1, i = 0; iBand <= poDS->nBands; iBand++)
    {
        poDS->SetBand(iBand, new L1BRasterBand(poDS, iBand));

        // Channels descriptions
        if (poDS->eSpacecraftID >= NOAA6 && poDS->eSpacecraftID <= METOP3)
        {
            if (!(i & 0x01) && poDS->iChannelsMask & 0x01)
            {
                poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[0]);
                i |= 0x01;
                continue;
            }
            if (!(i & 0x02) && poDS->iChannelsMask & 0x02)
            {
                poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[1]);
                i |= 0x02;
                continue;
            }
            if (!(i & 0x04) && poDS->iChannelsMask & 0x04)
            {
                if (poDS->eSpacecraftID >= NOAA15 &&
                    poDS->eSpacecraftID <= METOP3)
                    if (poDS->iInstrumentStatus & 0x0400)
                        poDS->GetRasterBand(iBand)->SetDescription(
                            apszBandDesc[7]);
                    else
                        poDS->GetRasterBand(iBand)->SetDescription(
                            apszBandDesc[6]);
                else
                    poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[2]);
                i |= 0x04;
                continue;
            }
            if (!(i & 0x08) && poDS->iChannelsMask & 0x08)
            {
                poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
                i |= 0x08;
                continue;
            }
            if (!(i & 0x10) && poDS->iChannelsMask & 0x10)
            {
                if (poDS->eSpacecraftID == NOAA13)  // 5 NOAA-13
                    poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[5]);
                else if (poDS->eSpacecraftID == NOAA6 ||
                         poDS->eSpacecraftID == NOAA8 ||
                         poDS->eSpacecraftID == NOAA10)  // 4 NOAA-6,-8,-10
                    poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
                else
                    poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[4]);
                i |= 0x10;
                continue;
            }
        }
    }

    if (poDS->bExposeMaskBand)
        poDS->poMaskBand = new L1BMaskBand(poDS);

    /* -------------------------------------------------------------------- */
    /*      Initialize any PAM information.                                 */
    /* -------------------------------------------------------------------- */
    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->TryLoadXML();

    /* -------------------------------------------------------------------- */
    /*      Check for external overviews.                                   */
    /* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
                                poOpenInfo->GetSiblingFiles());

    /* -------------------------------------------------------------------- */
    /*      Fetch metadata in CSV file                                      */
    /* -------------------------------------------------------------------- */
    if (CPLTestBool(CPLGetConfigOption("L1B_FETCH_METADATA", "NO")))
    {
        poDS->FetchMetadata();
    }

    return poDS;

bad:
    delete poDS;
    return nullptr;
}

/************************************************************************/
/*                        GDALRegister_L1B()                            */
/************************************************************************/

void GDALRegister_L1B()

{
    if (GDALGetDriverByName("L1B") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("L1B");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
                              "NOAA Polar Orbiter Level 1b Data Set");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/l1b.html");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");

    poDriver->pfnOpen = L1BDataset::Open;
    poDriver->pfnIdentify = L1BDataset::Identify;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}