/******************************************************************************

 *

 * 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 "gdal_driver.h"

#include "gdal_drivermanager.h"

#include "gdal_openinfo.h"

#include "gdal_cpp_functions.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);

    ~L1BDataset() override;

    int GetGCPCount() override;

    const OGRSpatialReference *GetGCPSpatialRef() const override;

    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 );

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

    GDALRasterBand *GetMaskBand() override;

    int GetMaskFlags() override;

};

/************************************************************************/

/* ==================================================================== */

/*                            L1BMaskBand                               */

/* ==================================================================== */

/************************************************************************/

class L1BMaskBand final : public GDALPamRasterBand

{

    friend class L1BDataset;

  public:

    explicit L1BMaskBand(L1BDataset *);

    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_UInt8;

    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 = cpl::down_cast<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 = cpl::down_cast<L1BDataset *>(poDS);

    if (poGDS->poMaskBand)

        return poGDS->poMaskBand;

    return GDALPamRasterBand::GetMaskBand();

}

/************************************************************************/

/*                            GetMaskFlags()                            */

/************************************************************************/

int L1BRasterBand::GetMaskFlags()

{

    L1BDataset *poGDS = cpl::down_cast<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 = cpl::down_cast<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()

{

    if (nRasterYSize == 0)

        return;

    void *pRecordHeader = CPLCalloc(1, nRecordDataStart);

    CPL_IGNORE_RET_VAL(

        VSIFSeekL(fp, static_cast<vsi_l_offset>(nDataStartOffset), SEEK_SET));

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

    FetchTimeCode(&sStartTime, pRecordHeader, &eLocationIndicator);

    CPL_IGNORE_RET_VAL(

        VSIFSeekL(fp,

                  nDataStartOffset +

                      static_cast<vsi_l_offset>(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)