/****************************************************************************** * * 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 * * 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 CPL_CVSID("$Id: sdtsrasterreader.cpp 1d00777992650359787f9a0d393afbd596f5d2a5 2018-01-15 14:47:44Z 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( pData )[i] = CPL_MSBWORD16( reinterpret_cast( pData )[i] ); } } else { for( int i = 0; i < nXSize; i++ ) { CPL_MSBPTR32( reinterpret_cast( 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). 
  Xg = padfTransformOut[0] + Xp * padfTransform[1] + Yp * padfTransform[2]
  Yg = padfTransformOut[3] + Xp * padfTransform[4] + Yp * padfTransform[5]
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( pBuffer )[iPixel]; else dfValue = reinterpret_cast( 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; }