/******************************************************************************
 *
 * Project:  MSG Driver
 * Purpose:  GDALDataset driver for MSG translator for read support.
 * Author:   Bas Retsios, retsios@itc.nl
 *
 ******************************************************************************
 * Copyright (c) 2004, ITC
 * Copyright (c) 2009, Even Rouault <even dot rouault at spatialys.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************/
#include "cpl_port.h"  // Must be first.

#include "gdal_frmts.h"
#include "msgdataset.h"
#include "prologue.h"
#include "xritheaderparser.h"
#include "reflectancecalculator.h"

#include "PublicDecompWT_headers.h"

#include <memory>
#include <vector>

#if _MSC_VER > 1000
#include <io.h>
#else
#include <stdio.h>
#endif

const double MSGDataset::rCentralWvl[12] = {0.635,  0.810,  1.640,  3.900,
                                            6.250,  7.350,  8.701,  9.660,
                                            10.800, 12.000, 13.400, 0.750};
const double MSGDataset::rVc[12] = {-1,       -1,       -1,       2569.094,
                                    1598.566, 1362.142, 1149.083, 1034.345,
                                    930.659,  839.661,  752.381,  -1};
const double MSGDataset::rA[12] = {-1,     -1,     -1,     0.9959,
                                   0.9963, 0.9991, 0.9996, 0.9999,
                                   0.9983, 0.9988, 0.9981, -1};
const double MSGDataset::rB[12] = {-1,    -1,    -1,    3.471, 2.219, 0.485,
                                   0.181, 0.060, 0.627, 0.397, 0.576, -1};
const int MSGDataset::iCentralPixelVIS_IR = 1856;  // center pixel VIS and IR
const int MSGDataset::iCentralPixelHRV = 5566;     // center pixel HRV
int MSGDataset::iCurrentSatelliteHint =
    1;  // satellite number hint 1,2,3,4 for MSG1, MSG2, MSG3 and MSG4
const char *MSGDataset::metadataDomain = "msg";  // the metadata domain

#define MAX_SATELLITES 4

/************************************************************************/
/*                    MSGDataset()                                     */
/************************************************************************/

MSGDataset::MSGDataset()

{
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    adfGeoTransform[0] = 0.0;
    adfGeoTransform[1] = 1.0;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = 0.0;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = 1.0;
}

/************************************************************************/
/*                    ~MSGDataset()                                     */
/************************************************************************/

MSGDataset::~MSGDataset()

{
    delete poTransform;
}

/************************************************************************/
/*                          GetGeoTransform()                           */
/************************************************************************/

CPLErr MSGDataset::GetGeoTransform(double *padfTransform)

{
    memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
    return CE_None;
}

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

GDALDataset *MSGDataset::Open(GDALOpenInfo *poOpenInfo)

{
    /* -------------------------------------------------------------------- */
    /*     Does this look like a MSG file                                   */
    /* -------------------------------------------------------------------- */
    // if( poOpenInfo->fp == nullptr)
    //   return nullptr;
    //  Do not touch the fp .. it will close by itself if not null after we
    //  return (whether it is recognized as HRIT or not)

    std::string command_line(poOpenInfo->pszFilename);

    MSGCommand command;
    std::string sErr = command.parse(command_line);
    if (sErr.length() > 0)
    {
        if (sErr.compare("-") !=
            0)  // this driver does not recognize this format .. be silent and
            // return false so that another driver can try
            CPLError(CE_Failure, CPLE_AppDefined, "%s", (sErr + "\n").c_str());
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Read the prologue.                                              */
    /* -------------------------------------------------------------------- */
    Prologue pp;

    int iCurrentSatellite =
        iCurrentSatelliteHint;  // Start with the hint. It is nice to have but
                                // don't rely on it...

    std::string sPrologueFileName =
        command.sPrologueFileName(iCurrentSatellite, 1);
    bool fPrologueExists = (access(sPrologueFileName.c_str(), 0) == 0);

    // Make sure we're testing for MSG1,2,3 or 4 exactly once, start with the
    // hint which is the most recently used.
    int iTries = 1;
    while (!fPrologueExists && (iTries < MAX_SATELLITES))
    {
        iCurrentSatellite = 1 + iCurrentSatellite % MAX_SATELLITES;
        sPrologueFileName = command.sPrologueFileName(iCurrentSatellite, 1);
        fPrologueExists = (access(sPrologueFileName.c_str(), 0) == 0);
        ++iTries;
    }

    if (fPrologueExists)
    {
        iCurrentSatelliteHint = iCurrentSatellite;  // set the hint
        std::ifstream p_file(sPrologueFileName.c_str(),
                             std::ios::in | std::ios::binary);
        XRITHeaderParser xhp(p_file);
        if (xhp.isValid() && xhp.isPrologue())
            pp.read(p_file);
        p_file.close();
    }
    else
    {
        std::string l_sErr = "The prologue of the data set could not be found "
                             "at the location specified:\n" +
                             sPrologueFileName + "\n";
        CPLError(CE_Failure, CPLE_AppDefined, "%s", l_sErr.c_str());
        return nullptr;
    }

    // We're confident the string is formatted as an MSG command_line

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */

    MSGDataset *poDS = new MSGDataset();
    poDS->command = command;  // copy it

    /* -------------------------------------------------------------------- */
    /*      Set the current satellite for this DS                           */
    /* -------------------------------------------------------------------- */

    poDS->iCurrentSatellite = iCurrentSatellite;  // copy it

    /* -------------------------------------------------------------------- */
    /*      Capture raster size from MSG prologue and submit it to GDAL     */
    /* -------------------------------------------------------------------- */

    if (command.channel[11] != 0)  // the HRV band
    {
        poDS->nRasterXSize = pp.idr()->ReferenceGridHRV->NumberOfColumns;
        poDS->nRasterYSize =
            abs(pp.idr()->PlannedCoverageHRV->UpperNorthLinePlanned -
                pp.idr()->PlannedCoverageHRV->LowerSouthLinePlanned) +
            1;
    }
    else
    {
        poDS->nRasterXSize =
            abs(pp.idr()->PlannedCoverageVIS_IR->WesternColumnPlanned -
                pp.idr()->PlannedCoverageVIS_IR->EasternColumnPlanned) +
            1;
        poDS->nRasterYSize =
            abs(pp.idr()->PlannedCoverageVIS_IR->NorthernLinePlanned -
                pp.idr()->PlannedCoverageVIS_IR->SouthernLinePlanned) +
            1;
    }

    /* -------------------------------------------------------------------- */
    /*      Set Georeference Information                                    */
    /* -------------------------------------------------------------------- */

    double rPixelSizeX;
    double rPixelSizeY;
    double rMinX;
    double rMaxY;

    if (command.channel[11] != 0)
    {
        rPixelSizeX = 1000 * pp.idr()->ReferenceGridHRV->ColumnDirGridStep;
        rPixelSizeY = 1000 * pp.idr()->ReferenceGridHRV->LineDirGridStep;
        // The MSG Level 1.5 Image Data Format Description (page 22f) defines
        // the center of pixel (5566,5566) from SE as sub satellite point for
        // the HRV channel (0°,0° for operational MSG).
        rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridHRV->NumberOfColumns -
                                iCentralPixelHRV + 0.5);
        rMaxY = rPixelSizeY *
                (pp.idr()->ReferenceGridHRV->NumberOfLines - iCentralPixelHRV +
                 0.5);  // scan direction south -> north
    }
    else
    {
        rPixelSizeX = 1000 * pp.idr()->ReferenceGridVIS_IR->ColumnDirGridStep;
        rPixelSizeY = 1000 * pp.idr()->ReferenceGridVIS_IR->LineDirGridStep;
        // The MSG Level 1.5 Image Data Format Description (page 22f) defines
        // the center of pixel (1856,1856) from SE as sub satellite point for
        // the VIS_IR channels (0°,0° for operational MSG).
        rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridVIS_IR->NumberOfColumns -
                                iCentralPixelVIS_IR + 0.5);
        rMaxY = rPixelSizeY *
                (pp.idr()->ReferenceGridVIS_IR->NumberOfLines -
                 iCentralPixelVIS_IR +
                 0.5);  // The y scan direction is always south -> north
    }
    poDS->adfGeoTransform[0] = rMinX;
    poDS->adfGeoTransform[3] = rMaxY;
    poDS->adfGeoTransform[1] = rPixelSizeX;
    poDS->adfGeoTransform[5] = -rPixelSizeY;
    poDS->adfGeoTransform[2] = 0.0;
    poDS->adfGeoTransform[4] = 0.0;

    /* -------------------------------------------------------------------- */
    /*      Set Projection Information                                      */
    /* -------------------------------------------------------------------- */

    poDS->m_oSRS.SetGEOS(0, 35785831, 0, 0);
    poDS->m_oSRS.SetWellKnownGeogCS(
        "WGS84");  // Temporary line to satisfy ERDAS (otherwise the ellips is
                   // "unnamed"). Eventually this should become the custom a and
                   // b ellips (CGMS).

    // The following are 3 different try-outs for also setting the ellips a and
    // b parameters. We leave them out for now however because this does not
    // work. In gdalwarp, when choosing some specific target SRS, the result is
    // an error message:
    //
    // ERROR 1: geocentric transformation missing z or ellps
    // ERROR 1: GDALWarperOperation::ComputeSourceWindow() failed because
    // the pfnTransformer failed.
    //
    // I can't explain the reason for the message at this time (could be a
    // problem in the way the SRS is set here, but also a bug in Proj.4 or GDAL.
    /*
    poDS->m_oSRS.SetGeogCS( NULL, NULL, NULL, 6378169, 295.488065897, NULL, 0,
    NULL, 0 );

    poDS->m_oSRS.SetGeogCS( "unnamed ellipse", "unknown", "unnamed", 6378169,
    295.488065897, "Greenwich", 0.0);

    poDS->m_oSRS.importFromProj4("+proj=geos +h=35785831 +a=6378169
    +b=6356583.8");
    */

    /* -------------------------------------------------------------------- */
    /*   Create a transformer to LatLon (only for Reflectance calculation)  */
    /* -------------------------------------------------------------------- */

    OGRSpatialReference *poSRSLongLat = poDS->m_oSRS.CloneGeogCS();
    if (poSRSLongLat)
    {
        poSRSLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
        poDS->poTransform =
            OGRCreateCoordinateTransformation(&(poDS->m_oSRS), poSRSLongLat);
        delete poSRSLongLat;
    }
    /* -------------------------------------------------------------------- */
    /*      Set the radiometric calibration parameters.                     */
    /* -------------------------------------------------------------------- */

    memcpy(poDS->rCalibrationOffset, pp.rpr()->Cal_Offset, sizeof(double) * 12);
    memcpy(poDS->rCalibrationSlope, pp.rpr()->Cal_Slope, sizeof(double) * 12);

    /* -------------------------------------------------------------------- */
    /*      Create band information objects.                                */
    /* -------------------------------------------------------------------- */
    poDS->nBands = command.iNrChannels() * command.iNrCycles;
    for (int iBand = 0; iBand < poDS->nBands; iBand++)
    {
        poDS->SetBand(iBand + 1, new MSGRasterBand(poDS, iBand + 1));
    }

    /* -------------------------------------------------------------------- */
    /*      Confirm the requested access is supported.                      */
    /* -------------------------------------------------------------------- */
    if (poOpenInfo->eAccess == GA_Update)
    {
        delete poDS;
        CPLError(CE_Failure, CPLE_NotSupported,
                 "The MSG driver does not support update access to existing"
                 " datasets.\n");
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*                 Set DataSet metadata information                    */
    /* -------------------------------------------------------------------- */
    CPLString metadataValue;
    metadataValue.Printf("%d", poDS->iCurrentSatellite);
    poDS->SetMetadataItem("satellite_number", metadataValue.c_str(),
                          metadataDomain);

    return poDS;
}

/************************************************************************/
/*                         MSGRasterBand()                              */
/************************************************************************/

const double MSGRasterBand::rRTOA[12] = {20.76, 23.24, 19.85, -1, -1, -1,
                                         -1,    -1,    -1,    -1, -1, 25.11};

MSGRasterBand::MSGRasterBand(MSGDataset *poDSIn, int nBandIn)
    : fScanNorth(false), iLowerShift(0), iSplitLine(0),
      iLowerWestColumnPlanned(0)

{
    this->poDS = poDSIn;
    this->nBand = nBandIn;

    // Find if we're dealing with MSG1, MSG2, MSG3 or MSG4
    // Doing this per band is the only way to guarantee time-series when the
    // satellite is changed

    int iCurrentSatellite =
        poDSIn->iCurrentSatellite;  // Start with the satellite selected by the
                                    // dataset

    std::string sPrologueFileName =
        poDSIn->command.sPrologueFileName(iCurrentSatellite, 1);
    bool fPrologueExists = (access(sPrologueFileName.c_str(), 0) == 0);

    // Make sure we're testing for MSG1,2,3 or 4 exactly once, start with the
    // one of the dataset
    int iTries = 1;
    while (!fPrologueExists && (iTries < MAX_SATELLITES))
    {
        iCurrentSatellite = 1 + iCurrentSatellite % MAX_SATELLITES;
        sPrologueFileName =
            poDSIn->command.sPrologueFileName(iCurrentSatellite, 1);
        fPrologueExists = (access(sPrologueFileName.c_str(), 0) == 0);
        ++iTries;
    }

    iSatellite =
        iCurrentSatellite;  // From here on, the satellite that corresponds to
                            // this band is settled to the current satellite

    nBlockXSize = poDSIn->GetRasterXSize();
    nBlockYSize = poDSIn->GetRasterYSize();

    /* -------------------------------------------------------------------- */
    /*      Open an input file and capture the header for the nr. of bits.  */
    /* -------------------------------------------------------------------- */
    int iStrip = 1;
    int iChannel = poDSIn->command.iChannel(
        1 + ((nBand - 1) % poDSIn->command.iNrChannels()));
    std::string input_file =
        poDSIn->command.sFileName(iSatellite, nBand, iStrip);
    while ((access(input_file.c_str(), 0) != 0) &&
           (iStrip <= poDSIn->command.iNrStrips(
                          iChannel)))  // compensate for missing strips
        input_file = poDSIn->command.sFileName(iSatellite, nBand, ++iStrip);

    if (iStrip <= poDSIn->command.iNrStrips(iChannel))
    {
        std::ifstream i_file(input_file.c_str(),
                             std::ios::in | std::ios::binary);

        if (i_file.good())
        {
            XRITHeaderParser xhp(i_file);

            if (xhp.isValid())
            {
                // Data type is either 8 or 16 bits .. we tell this to GDAL here
                eDataType = GDT_Byte;  // default .. always works
                if (xhp.nrBitsPerPixel() > 8)
                {
                    if (poDSIn->command.cDataConversion == 'N')
                        eDataType =
                            GDT_UInt16;  // normal case: MSG 10 bits data
                    else if (poDSIn->command.cDataConversion == 'B')
                        eDataType = GDT_Byte;  // output data type Byte
                    else
                        eDataType = GDT_Float32;  // Radiometric calibration
                }

                // make IReadBlock be called once per file
                nBlockYSize = xhp.nrRows();

                // remember the scan direction

                fScanNorth = xhp.isScannedNorth();
            }
        }

        i_file.close();
    }
    else if (nBand > 1)
    {
        // missing entire band .. take data from first band
        MSGRasterBand *pFirstRasterBand =
            (MSGRasterBand *)poDSIn->GetRasterBand(1);
        eDataType = pFirstRasterBand->eDataType;
        nBlockYSize = pFirstRasterBand->nBlockYSize;
        fScanNorth = pFirstRasterBand->fScanNorth;
    }
    else  // also first band is missing .. do something for fail-safety
    {
        if (poDSIn->command.cDataConversion == 'N')
            eDataType = GDT_UInt16;  // normal case: MSG 10 bits data
        else if (poDSIn->command.cDataConversion == 'B')
            eDataType = GDT_Byte;  // output data type Byte
        else
            eDataType = GDT_Float32;  // Radiometric calibration

        // nBlockYSize : default
        // fScanNorth : default
    }
    /* -------------------------------------------------------------------- */
    /*      For the HRV band, read the prologue for shift and splitline.    */
    /* -------------------------------------------------------------------- */

    if (iChannel == 12)
    {
        if (fPrologueExists)
        {
            std::ifstream p_file(sPrologueFileName.c_str(),
                                 std::ios::in | std::ios::binary);
            XRITHeaderParser xhp(p_file);
            Prologue pp;
            if (xhp.isValid() && xhp.isPrologue())
                pp.read(p_file);
            p_file.close();

            iLowerShift = pp.idr()->PlannedCoverageHRV->UpperWestColumnPlanned -
                          pp.idr()->PlannedCoverageHRV->LowerWestColumnPlanned;
            iSplitLine =
                abs(pp.idr()->PlannedCoverageHRV->UpperNorthLinePlanned -
                    pp.idr()->PlannedCoverageHRV->LowerNorthLinePlanned) +
                1;  // without the "+ 1" the image of 1-Jan-2005 splits
                    // incorrectly
            iLowerWestColumnPlanned =
                pp.idr()->PlannedCoverageHRV->LowerWestColumnPlanned;
        }
    }

    /* -------------------------------------------------------------------- */
    /*  Initialize the ReflectanceCalculator with the band-dependent info.  */
    /* -------------------------------------------------------------------- */

    int iCycle = 1 + (nBand - 1) / poDSIn->command.iNrChannels();
    std::string sTimeStamp = poDSIn->command.sCycle(iCycle);

    m_rc = new ReflectanceCalculator(sTimeStamp, rRTOA[iChannel - 1]);

    /* -------------------------------------------------------------------- */
    /*  Set DataSet metadata information                                   */
    /* -------------------------------------------------------------------- */
    CPLString metadataValue;
    metadataValue.Printf("%.10f", poDSIn->rCalibrationOffset[iChannel - 1]);
    SetMetadataItem("calibration_offset", metadataValue.c_str(),
                    poDSIn->metadataDomain);
    metadataValue.Printf("%.10f", poDSIn->rCalibrationSlope[iChannel - 1]);
    SetMetadataItem("calibration_slope", metadataValue.c_str(),
                    poDSIn->metadataDomain);
    metadataValue.Printf("%d", iChannel);
    SetMetadataItem("channel_number", metadataValue.c_str(),
                    poDSIn->metadataDomain);
}

/************************************************************************/
/*                           ~MSGRasterBand()                           */
/************************************************************************/
MSGRasterBand::~MSGRasterBand()
{
    delete m_rc;
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/
CPLErr MSGRasterBand::IReadBlock(int /*nBlockXOff*/, int nBlockYOff,
                                 void *pImage)

{

    MSGDataset *poGDS = (MSGDataset *)poDS;

    int iBytesPerPixel = 1;
    if (eDataType == GDT_UInt16)
        iBytesPerPixel = 2;
    else if (eDataType == GDT_Float32)
        iBytesPerPixel = 4;
    /* -------------------------------------------------------------------- */
    /*      Calculate the correct input file name based on nBlockYOff       */
    /* -------------------------------------------------------------------- */

    int strip_number;
    int iChannel = poGDS->command.iChannel(
        1 + ((nBand - 1) % poGDS->command.iNrChannels()));

    if (fScanNorth)
        strip_number = nBlockYOff + 1;
    else
        strip_number = poGDS->command.iNrStrips(iChannel) - nBlockYOff;

    std::string strip_input_file =
        poGDS->command.sFileName(iSatellite, nBand, strip_number);

    /* -------------------------------------------------------------------- */
    /*      Open the input file                                             */
    /* -------------------------------------------------------------------- */
    if (access(strip_input_file.c_str(), 0) == 0)  // does it exist?
    {
        std::ifstream i_file(strip_input_file.c_str(),
                             std::ios::in | std::ios::binary);

        if (i_file.good())
        {
            XRITHeaderParser xhp(i_file);

            if (xhp.isValid())
            {
                std::vector<short> QualityInfo;
                unsigned short chunk_height =
                    static_cast<unsigned short>(xhp.nrRows());
                unsigned short chunk_bpp =
                    static_cast<unsigned short>(xhp.nrBitsPerPixel());
                unsigned short chunk_width =
                    static_cast<unsigned short>(xhp.nrColumns());
                unsigned __int8 NR = (unsigned __int8)chunk_bpp;
                unsigned int nb_ibytes =
                    static_cast<unsigned int>(xhp.dataSize());
                int iShift = 0;
                bool fSplitStrip =
                    false;  // in the split strip the "shift" only happens
                            // before the split "row"
                int iSplitRow = 0;
                if (iChannel == 12)
                {
                    iSplitRow = iSplitLine % xhp.nrRows();
                    int iSplitBlock = iSplitLine / xhp.nrRows();
                    fSplitStrip =
                        (nBlockYOff ==
                         iSplitBlock);  // in the split strip the "shift" only
                                        // happens before the split "row"

                    // When iLowerShift > 0, the lower HRV image is shifted to
                    // the right When iLowerShift < 0, the lower HRV image is
                    // shifted to the left The available raster may be wider
                    // than needed, so that time series don't fall outside the
                    // raster.

                    if (nBlockYOff <= iSplitBlock)
                        iShift = -iLowerShift;
                    // iShift < 0 means upper image moves to the left
                    // iShift > 0 means upper image moves to the right
                }

                unsigned char *ibuf =
                    new (std::nothrow) unsigned char[nb_ibytes];
                if (ibuf == nullptr)
                {
                    CPLError(
                        CE_Failure, CPLE_AppDefined,
                        "Not enough memory to perform wavelet decompression\n");
                    return CE_Failure;
                }

                i_file.read((char *)ibuf, nb_ibytes);

                Util::CDataFieldCompressedImage img_compressed(
                    ibuf, nb_ibytes * 8, (unsigned char)chunk_bpp, chunk_width,
                    chunk_height);

                Util::CDataFieldUncompressedImage img_uncompressed;

                //****************************************************
                //*** Here comes the wavelets decompression routine
                COMP::DecompressWT(img_compressed, NR, img_uncompressed,
                                   QualityInfo);
                //****************************************************

                // convert:
                // Depth:
                // 8 bits -> 8 bits
                // 10 bits -> 16 bits (img_uncompressed contains the 10 bits
                // data in packed form) Geometry: chunk_width x chunk_height to
                // nBlockXSize x nBlockYSize

                // cases:
                // combination of the following:
                // - scan direction can be north or south
                // - eDataType can be GDT_Byte, GDT_UInt16 or GDT_Float32
                // - nBlockXSize == chunk_width or nBlockXSize > chunk_width
                // - when nBlockXSize > chunk_width, fSplitStrip can be true or
                // false we won't distinguish the following cases:
                // - NR can be == 8 or != 8
                // - when nBlockXSize > chunk_width, iShift iMinCOff-iMaxCOff <=
                // iShift <= 0

                int nBlockSize = nBlockXSize * nBlockYSize;
                int y = chunk_width * chunk_height;
                int iStep = -1;
                if (fScanNorth)  // image is the other way around
                {
                    y = -1;  // See how y is used below: += happens first, the
                             // result is used in the []
                    iStep = 1;
                }

                COMP::CImage cimg(img_uncompressed);  // unpack
                if (eDataType == GDT_Byte)
                {
                    if (nBlockXSize == chunk_width)  // optimized version
                    {
                        if (poGDS->command.cDataConversion == 'B')
                        {
                            for (int i = 0; i < nBlockSize; ++i)
                                ((GByte *)pImage)[i] = static_cast<GByte>(
                                    cimg.Get()[y += iStep] / 4);
                        }
                        else
                        {
                            for (int i = 0; i < nBlockSize; ++i)
                                ((GByte *)pImage)[i] =
                                    static_cast<GByte>(cimg.Get()[y += iStep]);
                        }
                    }
                    else
                    {
                        // initialize to 0's (so that it does not have to be
                        // done in an 'else' statement <performance>)
                        memset(pImage, 0,
                               nBlockXSize * nBlockYSize * iBytesPerPixel);
                        if (poGDS->command.cDataConversion == 'B')
                        {
                            for (int j = 0; j < chunk_height;
                                 ++j)  // assumption: nBlockYSize ==
                                       // chunk_height
                            {
                                int iXOffset = j * nBlockXSize + iShift;
                                iXOffset += nBlockXSize -
                                            iLowerWestColumnPlanned -
                                            1;  // Position the HRV part in the
                                                // frame; -1 to compensate the
                                                // pre-increment in the for-loop
                                if (fSplitStrip &&
                                    (j >= iSplitRow))  // In splitstrip, below
                                                       // splitline, thus do not
                                                       // shift!!
                                    iXOffset -= iShift;
                                for (int i = 0; i < chunk_width; ++i)
                                    ((GByte *)pImage)[++iXOffset] =
                                        static_cast<GByte>(
                                            cimg.Get()[y += iStep] / 4);
                            }
                        }
                        else
                        {
                            for (int j = 0; j < chunk_height;
                                 ++j)  // assumption: nBlockYSize ==
                                       // chunk_height
                            {
                                int iXOffset = j * nBlockXSize + iShift;
                                iXOffset += nBlockXSize -
                                            iLowerWestColumnPlanned -
                                            1;  // Position the HRV part in the
                                                // frame; -1 to compensate the
                                                // pre-increment in the for-loop
                                if (fSplitStrip &&
                                    (j >= iSplitRow))  // In splitstrip, below
                                                       // splitline, thus do not
                                                       // shift!!
                                    iXOffset -= iShift;
                                for (int i = 0; i < chunk_width; ++i)
                                    ((GByte *)pImage)[++iXOffset] =
                                        static_cast<GByte>(
                                            cimg.Get()[y += iStep]);
                            }
                        }
                    }
                }
                else if (eDataType ==
                         GDT_UInt16)  // this is our "normal case" if scan
                                      // direction is South: 10 bit MSG data
                                      // became 2 bytes per pixel
                {
                    if (nBlockXSize == chunk_width)  // optimized version
                    {
                        for (int i = 0; i < nBlockSize; ++i)
                            ((GUInt16 *)pImage)[i] = cimg.Get()[y += iStep];
                    }
                    else
                    {
                        // initialize to 0's (so that it does not have to be
                        // done in an 'else' statement <performance>)
                        memset(pImage, 0,
                               nBlockXSize * nBlockYSize * iBytesPerPixel);
                        for (int j = 0; j < chunk_height;
                             ++j)  // assumption: nBlockYSize == chunk_height
                        {
                            int iXOffset = j * nBlockXSize + iShift;
                            iXOffset += nBlockXSize - iLowerWestColumnPlanned -
                                        1;  // Position the HRV part in the
                                            // frame; -1 to compensate the
                                            // pre-increment in the for-loop
                            if (fSplitStrip &&
                                (j >=
                                 iSplitRow))  // In splitstrip, below splitline,
                                              // thus do not shift!!
                                iXOffset -= iShift;
                            for (int i = 0; i < chunk_width; ++i)
                                ((GUInt16 *)pImage)[++iXOffset] =
                                    cimg.Get()[y += iStep];
                        }
                    }
                }
                else if (eDataType ==
                         GDT_Float32)  // radiometric calibration is requested
                {
                    if (nBlockXSize == chunk_width)  // optimized version
                    {
                        for (int i = 0; i < nBlockSize; ++i)
                            ((float *)pImage)[i] =
                                (float)rRadiometricCorrection(
                                    cimg.Get()[y += iStep], iChannel,
                                    nBlockYOff * nBlockYSize + i / nBlockXSize,
                                    i % nBlockXSize, poGDS);
                    }
                    else
                    {
                        // initialize to 0's (so that it does not have to be
                        // done in an 'else' statement <performance>)
                        memset(pImage, 0,
                               nBlockXSize * nBlockYSize * iBytesPerPixel);
                        for (int j = 0; j < chunk_height;
                             ++j)  // assumption: nBlockYSize == chunk_height
                        {
                            int iXOffset = j * nBlockXSize + iShift;
                            iXOffset += nBlockXSize - iLowerWestColumnPlanned -
                                        1;  // Position the HRV part in the
                                            // frame; -1 to compensate the
                                            // pre-increment in the for-loop
                            if (fSplitStrip &&
                                (j >=
                                 iSplitRow))  // In splitstrip, below splitline,
                                              // thus do not shift!!
                                iXOffset -= iShift;
                            int iXFrom =
                                nBlockXSize - iLowerWestColumnPlanned +
                                iShift;  // i is used as the iCol parameter in
                                         // rRadiometricCorrection
                            int iXTo = nBlockXSize - iLowerWestColumnPlanned +
                                       chunk_width + iShift;
                            for (int i = iXFrom; i < iXTo;
                                 ++i)  // range always equal to chunk_width ..
                                       // this is to utilize i to get iCol
                                ((float *)pImage)[++iXOffset] =
                                    (float)rRadiometricCorrection(
                                        cimg.Get()[y += iStep], iChannel,
                                        nBlockYOff * nBlockYSize + j,
                                        (fSplitStrip && (j >= iSplitRow))
                                            ? (i - iShift)
                                            : i,
                                        poGDS);
                        }
                    }
                }
            }
            else  // header could not be opened .. make sure block contains 0's
                memset(pImage, 0, nBlockXSize * nBlockYSize * iBytesPerPixel);
        }
        else  // file could not be opened .. make sure block contains 0's
            memset(pImage, 0, nBlockXSize * nBlockYSize * iBytesPerPixel);

        i_file.close();
    }
    else  // file does not exist .. make sure block contains 0's
        memset(pImage, 0, nBlockXSize * nBlockYSize * iBytesPerPixel);

    return CE_None;
}

double MSGRasterBand::rRadiometricCorrection(unsigned int iDN, int iChannel,
                                             int iRow, int iCol,
                                             MSGDataset *poGDS) const
{
    int iIndex = iChannel - 1;  // just for speed optimization

    double rSlope = poGDS->rCalibrationSlope[iIndex];
    double rOffset = poGDS->rCalibrationOffset[iIndex];

    // Reflectance for visual bands, temperature for IR bands.
    if (poGDS->command.cDataConversion == 'T')
    {
        double rRadiance = rOffset + (iDN * rSlope);

        if (iChannel >= 4 &&
            iChannel <= 11)  // Channels 4 to 11 (infrared): Temperature
        {
            const double rC1 = 1.19104e-5;
            const double rC2 = 1.43877e+0;

            double cc2 = rC2 * poGDS->rVc[iIndex];
            double cc1 = rC1 * pow(poGDS->rVc[iIndex], 3) / rRadiance;
            double rTemperature =
                ((cc2 / log1p(cc1)) - poGDS->rB[iIndex]) / poGDS->rA[iIndex];
            return rTemperature;
        }
        else  // Channels 1,2,3 and 12 (visual): Reflectance
        {
            double rLon =
                poGDS->adfGeoTransform[0] +
                iCol * poGDS->adfGeoTransform[1];  // X, in "geos" meters
            double rLat =
                poGDS->adfGeoTransform[3] +
                iRow * poGDS->adfGeoTransform[5];  // Y, in "geos" meters
            if ((poGDS->poTransform != nullptr) &&
                poGDS->poTransform->Transform(1, &rLon,
                                              &rLat))  // transform it to latlon
                return m_rc->rGetReflectance(rRadiance, rLat, rLon);
            else
                return 0;
        }
    }
    else  // radiometric
    {
        if (poGDS->command.cDataConversion == 'R')
            return rOffset + (iDN * rSlope);
        else
        {
            double rFactor = 10 / pow(poGDS->rCentralWvl[iIndex], 2);
            return rFactor * (rOffset + (iDN * rSlope));
        }
    }
}

/************************************************************************/
/*                      GDALRegister_MSG()                              */
/************************************************************************/

void GDALRegister_MSG()

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

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("MSG");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "MSG HRIT Data");

    poDriver->pfnOpen = MSGDataset::Open;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}