EVOLUTION-MANAGER

Edit File: msgndataset.cpp

/******************************************************************************
 *
 * Project:  MSG Native Reader
 * Purpose:  All code for EUMETSAT Archive format reader
 * Author:   Frans van den Bergh, fvdbergh@csir.co.za
 *
 ******************************************************************************
 * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
 * Copyright (c) 2008-2009, Even Rouault <even dot rouault at mines-paris dot org>
 *
 * 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 "gdal_frmts.h"
#include "gdal_priv.h"
#include "ogr_spatialref.h"

#include "msg_reader_core.h"

#include <algorithm>

using namespace msg_native_format;

CPL_CVSID("$Id: msgndataset.cpp 36501 2016-11-25 14:09:24Z rouault $");

typedef enum {
    MODE_VISIR,     // Visible and Infrared bands (1 through 11) in 10-bit raw mode
    MODE_HRV,       // Pan band (band 11) only, in 10-bit raw mode
    MODE_RAD     // Black-body temperature (K) for thermal bands only (4-10), 64-bit float
    } open_mode_type;

class MSGNRasterBand;

/************************************************************************/
/* ==================================================================== */
/*                            MSGNDataset                               */
/* ==================================================================== */
/************************************************************************/

class MSGNDataset : public GDALDataset
{
    friend class MSGNRasterBand;

FILE       *fp;

Msg_reader_core*    msg_reader_core;
    double adfGeoTransform[6];
    char   *pszProjection;

public:
        MSGNDataset();
               ~MSGNDataset();

static GDALDataset *Open( GDALOpenInfo * );

CPLErr     GetGeoTransform( double * padfTransform ) override;
    const char *GetProjectionRef() override;
};

/************************************************************************/
/* ==================================================================== */
/*                            MSGNRasterBand                            */
/* ==================================================================== */
/************************************************************************/

class MSGNRasterBand : public GDALRasterBand
{
    friend class MSGNDataset;

unsigned int packet_size;
    unsigned int bytes_per_line;
    unsigned int interline_spacing;
    unsigned int orig_band_no;      // The name of the band
    unsigned int band_in_file;      // The effective index of the band in the file
    open_mode_type open_mode;

double  GetNoDataValue (int *pbSuccess=NULL) override {
        if (pbSuccess) {
            *pbSuccess = 1;
        }
        return MSGN_NODATA_VALUE;
    }

double MSGN_NODATA_VALUE;

char band_description[30];

public:

MSGNRasterBand( MSGNDataset *, int , open_mode_type mode, int orig_band_no, int band_in_file);

virtual CPLErr IReadBlock( int, int, void * ) override;
    virtual double GetMinimum( int *pbSuccess = NULL ) override;
    virtual double GetMaximum(int *pbSuccess = NULL ) override;
    virtual const char* GetDescription() const override { return band_description; }
};

/************************************************************************/
/*                           MSGNRasterBand()                            */
/************************************************************************/

MSGNRasterBand::MSGNRasterBand( MSGNDataset *poDSIn, int nBandIn,
                                open_mode_type mode, int orig_band_noIn,
                                int band_in_fileIn ) :
    packet_size(0),
    bytes_per_line(0),
    interline_spacing(poDSIn->msg_reader_core->get_interline_spacing()),
    orig_band_no(orig_band_noIn),
    band_in_file(band_in_fileIn),
    open_mode(mode)
{
    poDS = poDSIn;
    nBand = nBandIn;  // GDAL's band number, i.e. always starts at 1.

snprintf(band_description, sizeof(band_description),
             "band %02u", orig_band_no);

if( mode != MODE_RAD )
    {
        eDataType = GDT_UInt16;
        MSGN_NODATA_VALUE = 0;
    }
    else
    {
        eDataType = GDT_Float64;
        MSGN_NODATA_VALUE = -1000;
    }

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

if( mode != MODE_HRV )
    {
        packet_size = poDSIn->msg_reader_core->get_visir_packet_size();
        bytes_per_line = poDSIn->msg_reader_core->get_visir_bytes_per_line();
    }
    else
    {
        packet_size = poDSIn->msg_reader_core->get_hrv_packet_size();
        bytes_per_line = poDSIn->msg_reader_core->get_hrv_bytes_per_line();
    }
}

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

CPLErr MSGNRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
                                   int nBlockYOff,
                                   void * pImage )

{
    MSGNDataset *poGDS = (MSGNDataset *) poDS;

// invert y position
    int i_nBlockYOff = poDS->GetRasterYSize() - 1 - nBlockYOff;

unsigned int data_length =  bytes_per_line + (unsigned int)sizeof(SUB_VISIRLINE);
    unsigned int data_offset = 0;

if (open_mode != MODE_HRV) {
        data_offset = poGDS->msg_reader_core->get_f_data_offset() +
            interline_spacing*i_nBlockYOff  + (band_in_file-1)*packet_size +
            (packet_size - data_length);
    } else {
        data_offset = poGDS->msg_reader_core->get_f_data_offset() +
            interline_spacing*(int(i_nBlockYOff/3) + 1) -
            packet_size*(3 - (i_nBlockYOff % 3)) + (packet_size - data_length);
    }

if( VSIFSeek( poGDS->fp, data_offset, SEEK_SET ) != 0 )
        return CE_Failure;

char *pszRecord = (char *) CPLMalloc(data_length);
    size_t nread = VSIFRead( pszRecord, 1, data_length, poGDS->fp );

SUB_VISIRLINE* p = (SUB_VISIRLINE*) pszRecord;
    to_native(*p);

if (p->lineValidity != 1) {
        for (int c=0; c < nBlockXSize; c++) {
            if (open_mode != MODE_RAD) {
                ((GUInt16 *)pImage)[c] = (GUInt16)MSGN_NODATA_VALUE;
            } else {
                ((double *)pImage)[c] = MSGN_NODATA_VALUE;
            }
        }
    }

if ( nread != data_length ||
        ( open_mode != MODE_HRV && (p->lineNumberInVisirGrid - poGDS->msg_reader_core->get_line_start()) != (unsigned int)i_nBlockYOff )
       ) { // no sophisticated checking for HRV at the moment
        CPLFree( pszRecord );

CPLError( CE_Failure, CPLE_AppDefined, "MSGN Scanline corrupt." );

return CE_Failure;
    }

// unpack the 10-bit values into 16-bit unsigned short ints
    unsigned char *cptr = (unsigned char*)pszRecord +
        (data_length - bytes_per_line);
    int bitsLeft = 8;
    unsigned short value = 0;

if (open_mode != MODE_RAD) {
        for (int c=0; c < nBlockXSize; c++) {
            value = 0;
            for (int bit=0; bit < 10; bit++) {
                value <<= 1;
                if (*cptr & 128) {
                    value |= 1;
                }
                *cptr <<= 1;
                bitsLeft--;
                if (bitsLeft == 0) {
                    cptr++;
                bitsLeft = 8;
                }
            }
            ((GUInt16 *)pImage)[nBlockXSize-1 - c] = value;
        }
    } else {
        // radiance mode
        for (int c=0; c < nBlockXSize; c++) {
            value = 0;
            for (int bit=0; bit < 10; bit++) {
                value <<= 1;
                if (*cptr & 128) {
                    value |= 1;
                }
                *cptr <<= 1;
                bitsLeft--;
                if (bitsLeft == 0) {
                    cptr++;
                bitsLeft = 8;
                }
            }
            double dvalue = double(value);
            double bbvalue = dvalue * poGDS->msg_reader_core->get_calibration_parameters()[orig_band_no-1].cal_slope +
                poGDS->msg_reader_core->get_calibration_parameters()[orig_band_no-1].cal_offset;

((double *)pImage)[nBlockXSize-1 -c] = bbvalue;
        }
    }
    CPLFree( pszRecord );
    return CE_None;
}

/************************************************************************/
/*                             GetMinimum()                             */
/************************************************************************/
double MSGNRasterBand::GetMinimum( int *pbSuccess ) {
    if (pbSuccess) {
        *pbSuccess = 1;
    }
    return open_mode != MODE_RAD ? 1 : GDALRasterBand::GetMinimum(pbSuccess);
}

/************************************************************************/
/*                             GetMaximum()                             */
/************************************************************************/
double MSGNRasterBand::GetMaximum(int *pbSuccess ) {
    if (pbSuccess) {
        *pbSuccess = 1;
    }
    return open_mode != MODE_RAD ? 1023 : GDALRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/* ==================================================================== */
/*                             MSGNDataset                             */
/* ==================================================================== */
/************************************************************************/

MSGNDataset::MSGNDataset() :
    fp(NULL),
    msg_reader_core(NULL),
    pszProjection(CPLStrdup(""))
{
    std::fill_n(adfGeoTransform, CPL_ARRAYSIZE(adfGeoTransform), 0);
}

/************************************************************************/
/*                            ~MSGNDataset()                             */
/************************************************************************/

MSGNDataset::~MSGNDataset()

{
    if( fp != NULL )
        VSIFClose( fp );

if (msg_reader_core) {
        delete msg_reader_core;
    }

CPLFree(pszProjection);
}

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

CPLErr MSGNDataset::GetGeoTransform( double * padfTransform )

{

for (int i=0; i < 6; i++) {
        padfTransform[i] = adfGeoTransform[i];
    }

return CE_None;
}

/************************************************************************/
/*                          GetProjectionRef()                          */
/************************************************************************/

const char *MSGNDataset::GetProjectionRef()

{
    return pszProjection;
}

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

GDALDataset *MSGNDataset::Open( GDALOpenInfo * poOpenInfo )

{
    open_mode_type open_mode = MODE_VISIR;
    GDALOpenInfo* open_info = poOpenInfo;

if (!poOpenInfo->bStatOK) {
        if ( STARTS_WITH_CI(poOpenInfo->pszFilename, "HRV:") ) {
            open_info = new GDALOpenInfo(&poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
            open_mode = MODE_HRV;
        } else
        if ( STARTS_WITH_CI(poOpenInfo->pszFilename, "RAD:") ) {
            open_info = new GDALOpenInfo(&poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
            open_mode = MODE_RAD;
        }
    }

/* -------------------------------------------------------------------- */
/*      Before trying MSGNOpen() we first verify that there is at        */
/*      least one "\n#keyword" type signature in the first chunk of     */
/*      the file.                                                       */
/* -------------------------------------------------------------------- */
    if( open_info->fpL == NULL || open_info->nHeaderBytes < 50 ) {
        if (open_info != poOpenInfo) {
            delete open_info;
        }
        return NULL;
    }

/* check if this is a "NATIVE" MSG format image */
    if( !STARTS_WITH_CI((char *)open_info->pabyHeader, "FormatName                  : NATIVE") )
    {
        if (open_info != poOpenInfo) {
            delete open_info;
        }
        return NULL;
    }

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

/* -------------------------------------------------------------------- */
/*      Create a corresponding GDALDataset.                             */
/* -------------------------------------------------------------------- */
    FILE* fp = VSIFOpen( open_info->pszFilename, "rb" );
    if( fp == NULL ) {
        if (open_info != poOpenInfo) {
            delete open_info;
        }
        return NULL;
    }

MSGNDataset *poDS = new MSGNDataset();

poDS->fp = fp;

/* -------------------------------------------------------------------- */
/*      Read the header.                                                */
/* -------------------------------------------------------------------- */
    // first reset the file pointer, then hand over to the msg_reader_core
    CPL_IGNORE_RET_VAL(VSIFSeek( poDS->fp, 0, SEEK_SET ));

poDS->msg_reader_core = new Msg_reader_core(poDS->fp);

if (!poDS->msg_reader_core->get_open_success()) {
        if (open_info != poOpenInfo) {
            delete open_info;
        }
        delete poDS;
        return NULL;
    }

poDS->nRasterXSize = poDS->msg_reader_core->get_columns();
    poDS->nRasterYSize = poDS->msg_reader_core->get_lines();

if (open_mode == MODE_HRV) {
        poDS->nRasterXSize *= 3;
        poDS->nRasterYSize *= 3;
    }

/* -------------------------------------------------------------------- */
/*      Create band information objects.                                */
/* -------------------------------------------------------------------- */
    unsigned int i;
    unsigned int band_count = 1;
    unsigned int missing_band_count = 0;
    unsigned char* bands = poDS->msg_reader_core->get_band_map();
    unsigned char band_map[MSG_NUM_CHANNELS+1];   // map GDAL band numbers to MSG channels
    for (i=0; i < MSG_NUM_CHANNELS; i++) {
        if (bands[i]) {
            bool ok_to_add = false;
            switch (open_mode) {
                case MODE_VISIR:
                    ok_to_add = i < MSG_NUM_CHANNELS - 1;
                    break;
                case MODE_RAD:
                    ok_to_add = (i <= 2) || (Msg_reader_core::Blackbody_LUT[i+1].B != 0);
                    break;
                case MODE_HRV:
                    ok_to_add = i == MSG_NUM_CHANNELS - 1;
                    break;
            }
            if (ok_to_add) {
                poDS->SetBand( band_count, new MSGNRasterBand( poDS, band_count, open_mode, i+1, i+1 - missing_band_count));
                band_map[band_count] = (unsigned char) (i+1);
                band_count++;
            }
        } else {
            missing_band_count++;
        }
    }

double pixel_gsd_x;
    double pixel_gsd_y;
    double origin_x;
    double origin_y;

if (open_mode != MODE_HRV) {
        pixel_gsd_x = 1000 * poDS->msg_reader_core->get_col_dir_step();  // convert from km to m
        pixel_gsd_y = 1000 * poDS->msg_reader_core->get_line_dir_step(); // convert from km to m
        origin_x = -pixel_gsd_x * (-(Conversions::nlines / 2.0) + poDS->msg_reader_core->get_col_start());
        origin_y = -pixel_gsd_y * ((Conversions::nlines / 2.0) - poDS->msg_reader_core->get_line_start());
    } else {
        pixel_gsd_x = 1000 * poDS->msg_reader_core->get_col_dir_step() / 3.0;  // convert from km to m, approximate for HRV
        pixel_gsd_y = 1000 * poDS->msg_reader_core->get_line_dir_step() / 3.0; // convert from km to m, approximate for HRV
        origin_x = -pixel_gsd_x * (-(3*Conversions::nlines / 2.0) + 3*poDS->msg_reader_core->get_col_start());
        origin_y = -pixel_gsd_y * ((3*Conversions::nlines / 2.0) - 3*poDS->msg_reader_core->get_line_start());
    }

poDS->adfGeoTransform[0] = origin_x;
    poDS->adfGeoTransform[1] = pixel_gsd_x;
    poDS->adfGeoTransform[2] = 0.0;

poDS->adfGeoTransform[3] = origin_y;
    poDS->adfGeoTransform[4] = 0.0;
    poDS->adfGeoTransform[5] = -pixel_gsd_y;

OGRSpatialReference oSRS;

oSRS.SetProjCS("Geostationary projection (MSG)");
    oSRS.SetGEOS(  0, 35785831, 0, 0 );
    oSRS.SetGeogCS(
        "MSG Ellipsoid",
        "MSG_DATUM",
        "MSG_SPHEROID",
        Conversions::rpol * 1000.0,
        1 / ( 1 - Conversions::rpol/Conversions::req)
    );

oSRS.exportToWkt( &(poDS->pszProjection) );

CALIBRATION* cal = poDS->msg_reader_core->get_calibration_parameters();
    char tagname[30];
    char field[300];

poDS->SetMetadataItem("Radiometric parameters format", "offset slope");
    for (i=1; i < band_count; i++) {
        snprintf(tagname, sizeof(tagname), "ch%02u_cal", band_map[i]);
        CPLsnprintf(field, sizeof(field), "%.12e %.12e", cal[band_map[i]-1].cal_offset, cal[band_map[i]-1].cal_slope);
        poDS->SetMetadataItem(tagname, field);
    }

snprintf(field, sizeof(field), "%04u%02u%02u/%02u:%02u",
        poDS->msg_reader_core->get_year(),
        poDS->msg_reader_core->get_month(),
        poDS->msg_reader_core->get_day(),
        poDS->msg_reader_core->get_hour(),
        poDS->msg_reader_core->get_minute()
    );
    poDS->SetMetadataItem("Date/Time", field);

snprintf(field, sizeof(field), "%u %u",
         poDS->msg_reader_core->get_line_start(),
         poDS->msg_reader_core->get_col_start()
    );
    poDS->SetMetadataItem("Origin", field);

if (open_info != poOpenInfo) {
        delete open_info;
    }

return poDS;
}

/************************************************************************/
/*                          GDALRegister_MSGN()                         */
/************************************************************************/

void GDALRegister_MSGN()

{
    if( GDALGetDriverByName( "MSGN" ) != NULL )
        return;

GDALDriver *poDriver = new GDALDriver();

poDriver->SetDescription( "MSGN" );
    poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
    poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
                               "EUMETSAT Archive native (.nat)" );
    poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_msgn.html" );
    poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nat" );

poDriver->pfnOpen = MSGNDataset::Open;

GetGDALDriverManager()->RegisterDriver( poDriver );
}