/****************************************************************************** * * 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 /************************************************************************/ /* 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(), static_cast(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; }