frmts/sdts/sdtsrasterreader.cpp

/******************************************************************************
 *
 * Project:  SDTS Translator
 * Purpose:  Implementation of SDTSRasterReader class.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 1999, Frank Warmerdam
 * Copyright (c) 2008-2013, 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 "sdts_al.h"

#include <algorithm>

CPL_CVSID("$Id: sdtsrasterreader.cpp b1c9c12ad373e40b955162b45d704070d4ebf7b0 2019-06-19 16:50:15 +0200 Even Rouault $")

/************************************************************************/
/*                          SDTSRasterReader()                          */
/************************************************************************/

SDTSRasterReader::SDTSRasterReader() :
    nXSize(0),
    nYSize(0),
    nXBlockSize(0),
    nYBlockSize(0),
    nXStart(0),
    nYStart(0)
{
    strcpy( szINTR, "CE" );
    memset( szModule, 0, sizeof(szModule) );
    memset( adfTransform, 0, sizeof(adfTransform) );
    memset( szFMT, 0, sizeof(szFMT) );
    memset( szUNITS, 0, sizeof(szUNITS) );
    memset( szLabel, 0, sizeof(szLabel) );
}

/************************************************************************/
/*                             ~SDTSRasterReader()                     */
/************************************************************************/

SDTSRasterReader::~SDTSRasterReader() {}

/************************************************************************/
/*                               Close()                                */
/************************************************************************/

void SDTSRasterReader::Close()

{
    oDDFModule.Close();
}

/************************************************************************/
/*                                Open()                                */
/*                                                                      */
/*      Open the requested cell file, and collect required              */
/*      information.                                                    */
/************************************************************************/

int SDTSRasterReader::Open( SDTS_CATD * poCATD, SDTS_IREF * poIREF,
                            const char * pszModule )

{
    snprintf( szModule, sizeof(szModule), "%s", pszModule );

/* ==================================================================== */
/*      Search the LDEF module for the requested cell module.           */
/* ==================================================================== */

/* -------------------------------------------------------------------- */
/*      Open the LDEF module, and report failure if it is missing.      */
/* -------------------------------------------------------------------- */
    if( poCATD->GetModuleFilePath("LDEF") == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find LDEF entry in CATD module ... "
                  "can't treat as raster.\n" );
        return FALSE;
    }

    DDFModule oLDEF;
    if( !oLDEF.Open( poCATD->GetModuleFilePath("LDEF") ) )
        return FALSE;

/* -------------------------------------------------------------------- */
/*      Read each record, till we find what we want.                    */
/* -------------------------------------------------------------------- */
    DDFRecord *poRecord = nullptr;
    while( (poRecord = oLDEF.ReadRecord() ) != nullptr )
    {
        const char* pszCandidateModule = poRecord->GetStringSubfield("LDEF",0,"CMNM",0);
        if( pszCandidateModule == nullptr )
        {
            poRecord = nullptr;
            break;
        }
        if( EQUAL(pszCandidateModule, pszModule) )
            break;
    }

    if( poRecord == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find module `%s' in LDEF file.\n",
                  pszModule );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Extract raster dimensions, and origin offset (0/1).             */
/* -------------------------------------------------------------------- */
    nXSize = poRecord->GetIntSubfield( "LDEF", 0, "NCOL", 0 );
    nYSize = poRecord->GetIntSubfield( "LDEF", 0, "NROW", 0 );

    nXStart = poRecord->GetIntSubfield( "LDEF", 0, "SOCI", 0 );
    nYStart = poRecord->GetIntSubfield( "LDEF", 0, "SORI", 0 );

/* -------------------------------------------------------------------- */
/*      Get the point in the pixel that the origin defines.  We only    */
/*      support top left and center.                                    */
/* -------------------------------------------------------------------- */
    const char* pszINTR = poRecord->GetStringSubfield(  "LDEF", 0, "INTR", 0 );
    if( pszINTR == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined, "Can't find INTR subfield of LDEF field" );
        return FALSE;
    }
    snprintf( szINTR, sizeof(szINTR), "%s", pszINTR );
    if( EQUAL(szINTR,"") )
        snprintf( szINTR, sizeof(szINTR), "%s", "CE" );

    if( !EQUAL(szINTR,"CE") && !EQUAL(szINTR,"TL") )
    {
        CPLError( CE_Warning, CPLE_AppDefined,
                  "Unsupported INTR value of `%s', assume CE.\n"
                  "Positions may be off by one pixel.\n",
                  szINTR );
        snprintf( szINTR, sizeof(szINTR), "%s", "CE" );
    }

/* -------------------------------------------------------------------- */
/*      Record the LDEF record number we used so we can find the        */
/*      corresponding RSDF record.                                      */
/* -------------------------------------------------------------------- */
    int nLDEF_RCID = poRecord->GetIntSubfield( "LDEF", 0, "RCID", 0 );

    oLDEF.Close();

/* ==================================================================== */
/*      Search the RSDF module for the requested cell module.           */
/* ==================================================================== */

/* -------------------------------------------------------------------- */
/*      Open the RSDF module, and report failure if it is missing.      */
/* -------------------------------------------------------------------- */
    if( poCATD->GetModuleFilePath("RSDF") == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find RSDF entry in CATD module ... "
                  "can't treat as raster.\n" );
        return FALSE;
    }

    DDFModule oRSDF;
    if( !oRSDF.Open( poCATD->GetModuleFilePath("RSDF") ) )
        return FALSE;

/* -------------------------------------------------------------------- */
/*      Read each record, till we find what we want.                    */
/* -------------------------------------------------------------------- */
    while( (poRecord = oRSDF.ReadRecord() ) != nullptr )
    {
        if( poRecord->GetIntSubfield("LYID",0,"RCID",0) == nLDEF_RCID )
            break;
    }

    if( poRecord == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find LDEF:%d record in RSDF file.\n",
                  nLDEF_RCID );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Establish the raster pixel/line to georef transformation.       */
/* -------------------------------------------------------------------- */

    if( poRecord->FindField( "SADR" ) == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find SADR field in RSDF record.\n" );
        return FALSE;
    }

    double dfZ;
    poIREF->GetSADR( poRecord->FindField( "SADR" ), 1,
                     adfTransform + 0, adfTransform + 3, &dfZ );

    adfTransform[1] = poIREF->dfXRes;
    adfTransform[2] = 0.0;
    adfTransform[4] = 0.0;
    adfTransform[5] = -1 * poIREF->dfYRes;

/* -------------------------------------------------------------------- */
/*      If the origin is the center of the pixel, then shift it back    */
/*      half a pixel to the top left of the top left.                   */
/* -------------------------------------------------------------------- */
    if( EQUAL(szINTR,"CE") )
    {
        adfTransform[0] -= adfTransform[1] * 0.5;
        adfTransform[3] -= adfTransform[5] * 0.5;
    }

/* -------------------------------------------------------------------- */
/*      Verify some other assumptions.                                  */
/* -------------------------------------------------------------------- */
    const char *pszString = poRecord->GetStringSubfield( "RSDF", 0, "OBRP", 0);
    if( pszString == nullptr )
        pszString = "";
    if( !EQUAL(pszString, "G2") )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "OBRP value of `%s' not expected 2D raster code (G2).\n",
                  pszString );
        return FALSE;
    }

    pszString = poRecord->GetStringSubfield( "RSDF", 0, "SCOR", 0);
    if( pszString == nullptr )
        pszString = "";
    if( !EQUAL(pszString,"TL") )
    {
        CPLError( CE_Warning, CPLE_AppDefined,
                  "SCOR (origin) is `%s' instead of expected top left.\n"
                  "Georef coordinates will likely be incorrect.\n",
                  pszString );
    }

    oRSDF.Close();

/* -------------------------------------------------------------------- */
/*      For now we will assume that the block size is one scanline.     */
/*      We will blow a gasket later while reading the cell file if      */
/*      this isn't the case.                                            */
/*                                                                      */
/*      This isn't a very flexible raster implementation!               */
/* -------------------------------------------------------------------- */
    nXBlockSize = nXSize;
    nYBlockSize = 1;

/* ==================================================================== */
/*      Fetch the data type used for the raster, and the units from     */
/*      the data dictionary/schema record (DDSH).                       */
/* ==================================================================== */

/* -------------------------------------------------------------------- */
/*      Open the DDSH module, and report failure if it is missing.      */
/* -------------------------------------------------------------------- */
    if( poCATD->GetModuleFilePath("DDSH") == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find DDSH entry in CATD module ... "
                  "can't treat as raster.\n" );
        return FALSE;
    }

    DDFModule oDDSH;
    if( !oDDSH.Open( poCATD->GetModuleFilePath("DDSH") ) )
        return FALSE;

/* -------------------------------------------------------------------- */
/*      Read each record, till we find what we want.                    */
/* -------------------------------------------------------------------- */
    while( (poRecord = oDDSH.ReadRecord() ) != nullptr )
    {
        const char* pszName = poRecord->GetStringSubfield("DDSH",0,"NAME",0);
        if( pszName == nullptr )
        {
            poRecord = nullptr;
            break;
        }
        if( EQUAL(pszName,pszModule) )
            break;
    }

    if( poRecord == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Can't find DDSH record for %s.\n",
                  pszModule );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Get some values we are interested in.                           */
/* -------------------------------------------------------------------- */
    if( poRecord->GetStringSubfield("DDSH",0,"FMT",0) != nullptr )
        snprintf( szFMT, sizeof(szFMT), "%s", poRecord->GetStringSubfield("DDSH",0,"FMT",0) );
    else
        snprintf( szFMT, sizeof(szFMT), "%s", "BI16" );
    if( !EQUAL(szFMT, "BI16") && !EQUAL(szFMT, "BFP32") )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Unhandled FMT=%s", szFMT );
        return FALSE;
    }

    if( poRecord->GetStringSubfield("DDSH",0,"UNIT",0) != nullptr )
        snprintf( szUNITS, sizeof(szUNITS), "%s", poRecord->GetStringSubfield("DDSH",0,"UNIT",0) );
    else
        snprintf( szUNITS, sizeof(szUNITS), "%s", "METERS" );

    if( poRecord->GetStringSubfield("DDSH",0,"ATLB",0) != nullptr )
        snprintf( szLabel, sizeof(szLabel), "%s", poRecord->GetStringSubfield("DDSH",0,"ATLB",0) );
    else
        strcpy( szLabel, "" );

/* -------------------------------------------------------------------- */
/*      Open the cell file.                                             */
/* -------------------------------------------------------------------- */
    return oDDFModule.Open( poCATD->GetModuleFilePath(pszModule) );
}

/************************************************************************/
/*                              GetBlock()                              */
/*                                                                      */
/*      Read a requested block of raster data from the file.            */
/*                                                                      */
/*      Currently we will always use sequential access.  In the         */
/*      future we should modify the iso8211 library to support          */
/*      seeking, and modify this to seek directly to the right          */
/*      record once its location is known.                              */
/************************************************************************/

/**
  Read a block of raster data from the file.

  @param nXOffset X block offset into the file.  Normally zero for scanline
  organized raster files.

  @param nYOffset Y block offset into the file.  Normally the scanline offset
  from top of raster for scanline organized raster files.

  @param pData pointer to GInt16 (signed short) buffer of data into which to
  read the raster.

  @return TRUE on success and FALSE on error.

  */

int SDTSRasterReader::GetBlock( CPL_UNUSED int nXOffset,
                                int nYOffset,
                                void * pData )
{
    CPLAssert( nXOffset == 0 );

/* -------------------------------------------------------------------- */
/*      Analyse the datatype.                                           */
/* -------------------------------------------------------------------- */
    CPLAssert( EQUAL(szFMT,"BI16") || EQUAL(szFMT,"BFP32") );

    int nBytesPerValue;
    if( EQUAL(szFMT,"BI16") )
        nBytesPerValue = 2;
    else
        nBytesPerValue = 4;

    DDFRecord *poRecord = nullptr;

    for( int iTry = 0; iTry < 2; iTry++ )
    {
    /* -------------------------------------------------------------------- */
    /*      Read through till we find the desired record.                   */
    /* -------------------------------------------------------------------- */
        CPLErrorReset();
        while( (poRecord = oDDFModule.ReadRecord()) != nullptr )
        {
            if( poRecord->GetIntSubfield( "CELL", 0, "ROWI", 0 )
                == nYOffset + nYStart )
            {
                break;
            }
        }

        if( CPLGetLastErrorType() == CE_Failure )
            return FALSE;

    /* -------------------------------------------------------------------- */
    /*      If we didn't get what we needed just start over.                */
    /* -------------------------------------------------------------------- */
        if( poRecord == nullptr )
        {
            if (iTry == 0)
                oDDFModule.Rewind();
            else
            {
                CPLError( CE_Failure, CPLE_AppDefined,
                          "Cannot read scanline %d.  Raster access failed.\n",
                          nYOffset );
                return FALSE;
            }
        }
        else
        {
            break;
        }
    }

/* -------------------------------------------------------------------- */
/*      Validate the records size.  Does it represent exactly one       */
/*      scanline?                                                       */
/* -------------------------------------------------------------------- */
    DDFField *poCVLS = poRecord->FindField( "CVLS" );
    if( poCVLS == nullptr )
        return FALSE;

    if( poCVLS->GetRepeatCount() != nXSize )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Cell record is %d long, but we expected %d, the number\n"
                  "of pixels in a scanline.  Raster access failed.\n",
                  poCVLS->GetRepeatCount(), nXSize );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Does the CVLS field consist of exactly 1 B(16) field?           */
/* -------------------------------------------------------------------- */
    if( poCVLS->GetDataSize() < nBytesPerValue * nXSize
        || poCVLS->GetDataSize() > nBytesPerValue * nXSize + 1 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Cell record is not of expected format.  Raster access "
                  "failed.\n" );

        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Copy the data to the application buffer, and byte swap if       */
/*      required.                                                       */
/* -------------------------------------------------------------------- */
    memcpy( pData, poCVLS->GetData(), nXSize * nBytesPerValue );

#ifdef CPL_LSB
    if( nBytesPerValue == 2 )
    {
        for( int i = 0; i < nXSize; i++ )
        {
            reinterpret_cast<GInt16 *>( pData )[i] = CPL_MSBWORD16(
                reinterpret_cast<GInt16 *>( pData )[i] );
        }
    }
    else
    {
        for( int i = 0; i < nXSize; i++ )
        {
            CPL_MSBPTR32( reinterpret_cast<GByte *>( pData ) + i*4 );
        }
    }
#endif

    return TRUE;
}

/************************************************************************/
/*                            GetTransform()                            */
/************************************************************************/

/**
  Fetch the transformation between pixel/line coordinates and georeferenced
  coordinates.

  @param padfTransformOut pointer to an array of six doubles which will be
  filled with the georeferencing transform.

  @return TRUE is returned, indicating success.

  The padfTransformOut array consists of six values.  The pixel/line coordinate
  (Xp,Yp) can be related to a georeferenced coordinate (Xg,Yg) or (Easting,
  Northing).

  <pre>
  Xg = padfTransformOut[0] + Xp * padfTransform[1] + Yp * padfTransform[2]
  Yg = padfTransformOut[3] + Xp * padfTransform[4] + Yp * padfTransform[5]
  </pre>

  In other words, for a north up image the top left corner of the top left
  pixel is at georeferenced coordinate (padfTransform[0],padfTransform[3])
  the pixel width is padfTransform[1], the pixel height is padfTransform[5]
  and padfTransform[2] and padfTransform[4] will be zero.

  */

int SDTSRasterReader::GetTransform( double * padfTransformOut )

{
    memcpy( padfTransformOut, adfTransform, sizeof(double)*6 );

    return TRUE;
}

/************************************************************************/
/*                           GetRasterType()                            */
/************************************************************************/

/**
 * Fetch the pixel data type.
 *
 * Returns one of SDTS_RT_INT16 (1) or SDTS_RT_FLOAT32 (6) indicating the
 * type of buffer that should be passed to GetBlock().
 */

int SDTSRasterReader::GetRasterType()

{
    if( EQUAL(szFMT,"BFP32") )
        return SDTS_RT_FLOAT32;

    return SDTS_RT_INT16;
}

/************************************************************************/
/*                             GetMinMax()                              */
/************************************************************************/

/**
 * Fetch the minimum and maximum raster values that occur in the file.
 *
 * Note this operation current results in a scan of the entire file.
 *
 * @param pdfMin variable in which the minimum value encountered is returned.
 * @param pdfMax variable in which the maximum value encountered is returned.
 * @param dfNoData a value to ignore when computing min/max, defaults to
 * -32766.
 *
 * @return TRUE on success, or FALSE if an error occurs.
 */

int SDTSRasterReader::GetMinMax( double * pdfMin, double * pdfMax,
                                 double dfNoData )

{
    CPLAssert( GetBlockXSize() == GetXSize() && GetBlockYSize() == 1 );

    bool bFirst = true;
    const bool b32Bit = GetRasterType() == SDTS_RT_FLOAT32;
    void *pBuffer = CPLMalloc(sizeof(float) * GetXSize());

    for( int iLine = 0; iLine < GetYSize(); iLine++ )
    {
        if( !GetBlock( 0, iLine, pBuffer ) )
        {
            CPLFree( pBuffer );
            return FALSE;
        }

        for( int iPixel = 0; iPixel < GetXSize(); iPixel++ )
        {
            double dfValue;

            if( b32Bit )
                dfValue = reinterpret_cast<float *>( pBuffer )[iPixel];
            else
                dfValue = reinterpret_cast<short *>( pBuffer )[iPixel];

            if( dfValue != dfNoData )
            {
                if( bFirst )
                {
                    *pdfMin = dfValue;
                    *pdfMax = dfValue;
                    bFirst = false;
                }
                else
                {
                    *pdfMin = std::min( *pdfMin, dfValue );
                    *pdfMax = std::max( *pdfMax, dfValue );
                }
            }
        }
    }

    CPLFree( pBuffer );

    return !bFirst;
}