Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:    Scaled Integer Gridded DEM .sigdem Driver
       4             :  * Purpose:    Implementation of Scaled Integer Gridded DEM
       5             :  * Author:     Paul Austin, paul.austin@revolsys.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2018, Paul Austin <paul.austin@revolsys.com>
       9             :  *
      10             :  * Permission is hereby granted, free of charge, to any person obtaining a
      11             :  * copy of this software and associated documentation files (the "Software"),
      12             :  * to deal in the Software without restriction, including without limitation
      13             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14             :  * and/or sell copies of the Software, and to permit persons to whom the
      15             :  * Software is furnished to do so, subject to the following conditions:
      16             :  *
      17             :  * The above copyright notice and this permission notice shall be included
      18             :  * in all copies or substantial portions of the Software.
      19             :  *
      20             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26             :  * DEALINGS IN THE SOFTWARE.
      27             :  ****************************************************************************/
      28             : 
      29             : #include "sigdemdataset.h"
      30             : #include "rawdataset.h"
      31             : 
      32             : #include <algorithm>
      33             : #include <limits>
      34             : 
      35             : #ifdef CPL_IS_LSB
      36             : #define SWAP_SIGDEM_HEADER(abyHeader)                                          \
      37             :     {                                                                          \
      38             :         CPL_SWAP16PTR(abyHeader + 6);                                          \
      39             :         CPL_SWAP32PTR(abyHeader + 8);                                          \
      40             :         CPL_SWAP64PTR(abyHeader + 12);                                         \
      41             :         CPL_SWAP64PTR(abyHeader + 20);                                         \
      42             :         CPL_SWAP64PTR(abyHeader + 28);                                         \
      43             :         CPL_SWAP64PTR(abyHeader + 36);                                         \
      44             :         CPL_SWAP64PTR(abyHeader + 44);                                         \
      45             :         CPL_SWAP64PTR(abyHeader + 52);                                         \
      46             :         CPL_SWAP64PTR(abyHeader + 60);                                         \
      47             :         CPL_SWAP64PTR(abyHeader + 68);                                         \
      48             :         CPL_SWAP64PTR(abyHeader + 76);                                         \
      49             :         CPL_SWAP64PTR(abyHeader + 84);                                         \
      50             :         CPL_SWAP64PTR(abyHeader + 92);                                         \
      51             :         CPL_SWAP64PTR(abyHeader + 100);                                        \
      52             :         CPL_SWAP32PTR(abyHeader + 108);                                        \
      53             :         CPL_SWAP32PTR(abyHeader + 112);                                        \
      54             :         CPL_SWAP64PTR(abyHeader + 116);                                        \
      55             :         CPL_SWAP64PTR(abyHeader + 124);                                        \
      56             :     }
      57             : #else
      58             : #define SWAP_SIGDEM_HEADER(abyHeader)
      59             : #endif
      60             : 
      61             : constexpr int CELL_SIZE_FILE = 4;
      62             : 
      63             : constexpr int CELL_SIZE_MEM = 8;
      64             : 
      65             : constexpr vsi_l_offset HEADER_LENGTH = 132;
      66             : 
      67             : constexpr int32_t NO_DATA = 0x80000000;
      68             : 
      69             : constexpr char SIGDEM_FILE_TYPE[6] = {'S', 'I', 'G', 'D', 'E', 'M'};
      70             : 
      71          24 : static OGRSpatialReference *BuildSRS(const char *pszWKT)
      72             : {
      73          24 :     OGRSpatialReference *poSRS = new OGRSpatialReference();
      74          24 :     if (poSRS->importFromWkt(pszWKT) != OGRERR_NONE)
      75             :     {
      76           0 :         delete poSRS;
      77           0 :         return nullptr;
      78             :     }
      79             :     else
      80             :     {
      81          24 :         if (poSRS->AutoIdentifyEPSG() != OGRERR_NONE)
      82             :         {
      83          22 :             auto poSRSMatch = poSRS->FindBestMatch(100);
      84          22 :             if (poSRSMatch)
      85             :             {
      86           1 :                 poSRS->Release();
      87           1 :                 poSRS = poSRSMatch;
      88           1 :                 poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      89             :             }
      90             :         }
      91          24 :         return poSRS;
      92             :     }
      93             : }
      94             : 
      95        1520 : void GDALRegister_SIGDEM()
      96             : {
      97        1520 :     if (GDALGetDriverByName("SIGDEM") == nullptr)
      98             :     {
      99        1219 :         GDALDriver *poDriver = new GDALDriver();
     100             : 
     101        1219 :         poDriver->SetDescription("SIGDEM");
     102        1219 :         poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     103        1219 :         poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     104        1219 :                                   "Scaled Integer Gridded DEM .sigdem");
     105        1219 :         poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
     106        1219 :                                   "drivers/raster/sigdem.html");
     107        1219 :         poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "sigdem");
     108             : 
     109        1219 :         poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     110        1219 :         poDriver->pfnCreateCopy = SIGDEMDataset::CreateCopy;
     111        1219 :         poDriver->pfnIdentify = SIGDEMDataset::Identify;
     112        1219 :         poDriver->pfnOpen = SIGDEMDataset::Open;
     113             : 
     114        1219 :         GetGDALDriverManager()->RegisterDriver(poDriver);
     115             :     }
     116        1520 : }
     117             : 
     118          24 : static int32_t GetCoordinateSystemId(const char *pszProjection)
     119             : {
     120          24 :     int32_t coordinateSystemId = 0;
     121          24 :     OGRSpatialReference *poSRS = BuildSRS(pszProjection);
     122          24 :     if (poSRS != nullptr)
     123             :     {
     124          48 :         std::string pszRoot;
     125          24 :         if (poSRS->IsProjected())
     126             :         {
     127           3 :             pszRoot = "PROJCS";
     128             :         }
     129             :         else
     130             :         {
     131          21 :             pszRoot = "GEOCS";
     132             :         }
     133          24 :         const char *pszAuthName = poSRS->GetAuthorityName(pszRoot.c_str());
     134          24 :         const char *pszAuthCode = poSRS->GetAuthorityCode(pszRoot.c_str());
     135          24 :         if (pszAuthName != nullptr && EQUAL(pszAuthName, "EPSG") &&
     136             :             pszAuthCode != nullptr)
     137             :         {
     138           3 :             coordinateSystemId = atoi(pszAuthCode);
     139             :         }
     140             :     }
     141          24 :     delete poSRS;
     142          24 :     return coordinateSystemId;
     143             : }
     144             : 
     145          20 : SIGDEMDataset::SIGDEMDataset(const SIGDEMHeader &sHeaderIn)
     146          20 :     : fpImage(nullptr), sHeader(sHeaderIn)
     147             : {
     148          20 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     149          20 :     this->nRasterXSize = sHeader.nCols;
     150          20 :     this->nRasterYSize = sHeader.nRows;
     151             : 
     152          20 :     this->adfGeoTransform[0] = sHeader.dfMinX;
     153          20 :     this->adfGeoTransform[1] = sHeader.dfXDim;
     154          20 :     this->adfGeoTransform[2] = 0.0;
     155          20 :     this->adfGeoTransform[3] = sHeader.dfMaxY;
     156          20 :     this->adfGeoTransform[4] = 0.0;
     157          20 :     this->adfGeoTransform[5] = -sHeader.dfYDim;
     158          20 : }
     159             : 
     160          40 : SIGDEMDataset::~SIGDEMDataset()
     161             : {
     162          20 :     FlushCache(true);
     163             : 
     164          20 :     if (fpImage != nullptr)
     165             :     {
     166          20 :         if (VSIFCloseL(fpImage) != 0)
     167             :         {
     168           0 :             CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     169             :         }
     170             :     }
     171          40 : }
     172             : 
     173          32 : GDALDataset *SIGDEMDataset::CreateCopy(const char *pszFilename,
     174             :                                        GDALDataset *poSrcDS, int /*bStrict*/,
     175             :                                        char ** /*papszOptions*/,
     176             :                                        GDALProgressFunc pfnProgress,
     177             :                                        void *pProgressData)
     178             : {
     179          32 :     const int nBands = poSrcDS->GetRasterCount();
     180          32 :     double adfGeoTransform[6] = {};
     181          32 :     if (poSrcDS->GetGeoTransform(adfGeoTransform) != CE_None)
     182             :     {
     183           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     184             :                  "SIGDEM driver requires a valid GeoTransform.");
     185           0 :         return nullptr;
     186             :     }
     187             : 
     188          32 :     if (nBands != 1)
     189             :     {
     190           5 :         CPLError(CE_Failure, CPLE_NotSupported,
     191             :                  "SIGDEM driver doesn't support %d bands.  Must be 1 band.",
     192             :                  nBands);
     193           5 :         return nullptr;
     194             :     }
     195             : 
     196          27 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
     197          27 :     if (fp == nullptr)
     198             :     {
     199           3 :         CPLError(CE_Failure, CPLE_OpenFailed,
     200             :                  "Attempt to create file `%s' failed.", pszFilename);
     201           3 :         return nullptr;
     202             :     }
     203             : 
     204          24 :     GDALRasterBand *band = poSrcDS->GetRasterBand(1);
     205          24 :     const char *pszProjection = poSrcDS->GetProjectionRef();
     206             : 
     207          24 :     int32_t nCols = poSrcDS->GetRasterXSize();
     208          24 :     int32_t nRows = poSrcDS->GetRasterYSize();
     209          24 :     int32_t nCoordinateSystemId = GetCoordinateSystemId(pszProjection);
     210             : 
     211          24 :     SIGDEMHeader sHeader;
     212          24 :     sHeader.nCoordinateSystemId = nCoordinateSystemId;
     213          24 :     sHeader.dfMinX = adfGeoTransform[0];
     214          24 :     const char *pszMin = band->GetMetadataItem("STATISTICS_MINIMUM");
     215          24 :     if (pszMin == nullptr)
     216             :     {
     217           3 :         sHeader.dfMinZ = -10000;
     218             :     }
     219             :     else
     220             :     {
     221          21 :         sHeader.dfMinZ = CPLAtof(pszMin);
     222             :     }
     223          24 :     sHeader.dfMaxY = adfGeoTransform[3];
     224          24 :     const char *pszMax = band->GetMetadataItem("STATISTICS_MAXIMUM");
     225          24 :     if (pszMax == nullptr)
     226             :     {
     227           3 :         sHeader.dfMaxZ = 10000;
     228             :     }
     229             :     else
     230             :     {
     231          21 :         sHeader.dfMaxZ = CPLAtof(pszMax);
     232             :     }
     233          24 :     sHeader.nCols = poSrcDS->GetRasterXSize();
     234          24 :     sHeader.nRows = poSrcDS->GetRasterYSize();
     235          24 :     sHeader.dfXDim = adfGeoTransform[1];
     236          24 :     sHeader.dfYDim = -adfGeoTransform[5];
     237          24 :     sHeader.dfMaxX = sHeader.dfMinX + sHeader.nCols * sHeader.dfXDim;
     238          24 :     sHeader.dfMinY = sHeader.dfMaxY - sHeader.nRows * sHeader.dfYDim;
     239          24 :     sHeader.dfOffsetX = sHeader.dfMinX;
     240          24 :     sHeader.dfOffsetY = sHeader.dfMinY;
     241             : 
     242          24 :     if (!sHeader.Write(fp))
     243             :     {
     244           3 :         VSIUnlink(pszFilename);
     245           3 :         VSIFCloseL(fp);
     246           3 :         return nullptr;
     247             :     }
     248             : 
     249             :     // Write fill with all NO_DATA values
     250             :     int32_t *row =
     251          21 :         static_cast<int32_t *>(VSI_MALLOC2_VERBOSE(nCols, sizeof(int32_t)));
     252          21 :     if (!row)
     253             :     {
     254           0 :         VSIUnlink(pszFilename);
     255           0 :         VSIFCloseL(fp);
     256           0 :         return nullptr;
     257             :     }
     258          21 :     std::fill(row, row + nCols, CPL_MSBWORD32(NO_DATA));
     259         236 :     for (int i = 0; i < nRows; i++)
     260             :     {
     261         222 :         if (VSIFWriteL(row, CELL_SIZE_FILE, nCols, fp) !=
     262         222 :             static_cast<size_t>(nCols))
     263             :         {
     264           7 :             VSIFree(row);
     265           7 :             VSIUnlink(pszFilename);
     266           7 :             VSIFCloseL(fp);
     267           7 :             return nullptr;
     268             :         }
     269             :     }
     270          14 :     VSIFree(row);
     271             : 
     272          14 :     if (VSIFCloseL(fp) != 0)
     273             :     {
     274           0 :         return nullptr;
     275             :     }
     276             : 
     277          14 :     if (nCoordinateSystemId <= 0)
     278             :     {
     279          11 :         if (!EQUAL(pszProjection, ""))
     280             :         {
     281          22 :             CPLString osPrjFilename = CPLResetExtension(pszFilename, "prj");
     282          11 :             VSILFILE *fpProj = VSIFOpenL(osPrjFilename, "wt");
     283          11 :             if (fpProj != nullptr)
     284             :             {
     285          22 :                 OGRSpatialReference oSRS;
     286          11 :                 oSRS.importFromWkt(pszProjection);
     287          11 :                 oSRS.morphToESRI();
     288          11 :                 char *pszESRIProjection = nullptr;
     289          11 :                 oSRS.exportToWkt(&pszESRIProjection);
     290          11 :                 CPL_IGNORE_RET_VAL(VSIFWriteL(
     291             :                     pszESRIProjection, 1, strlen(pszESRIProjection), fpProj));
     292             : 
     293          11 :                 CPL_IGNORE_RET_VAL(VSIFCloseL(fpProj));
     294          11 :                 CPLFree(pszESRIProjection);
     295             :             }
     296             :             else
     297             :             {
     298           0 :                 CPLError(CE_Failure, CPLE_FileIO, "Unable to create file %s.",
     299             :                          osPrjFilename.c_str());
     300             :             }
     301             :         }
     302             :     }
     303          28 :     GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
     304          14 :     GDALDataset *poDstDS = Open(&oOpenInfo);
     305          28 :     if (poDstDS != nullptr &&
     306          14 :         GDALDatasetCopyWholeRaster(poSrcDS, poDstDS, nullptr, pfnProgress,
     307             :                                    pProgressData) == OGRERR_NONE)
     308             :     {
     309          14 :         return poDstDS;
     310             :     }
     311             :     else
     312             :     {
     313           0 :         VSIUnlink(pszFilename);
     314           0 :         return nullptr;
     315             :     }
     316             : }
     317             : 
     318           2 : CPLErr SIGDEMDataset::GetGeoTransform(double *padfTransform)
     319             : {
     320           2 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     321           2 :     return CE_None;
     322             : }
     323             : 
     324           2 : const OGRSpatialReference *SIGDEMDataset::GetSpatialRef() const
     325             : {
     326           2 :     return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
     327             : }
     328             : 
     329       47981 : int SIGDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
     330             : {
     331       47981 :     if (poOpenInfo->nHeaderBytes < static_cast<int>(HEADER_LENGTH))
     332             :     {
     333       45486 :         return FALSE;
     334             :     }
     335        2495 :     return memcmp(poOpenInfo->pabyHeader, SIGDEM_FILE_TYPE,
     336        2495 :                   sizeof(SIGDEM_FILE_TYPE)) == 0;
     337             : }
     338             : 
     339          20 : GDALDataset *SIGDEMDataset::Open(GDALOpenInfo *poOpenInfo)
     340             : {
     341          20 :     VSILFILE *fp = poOpenInfo->fpL;
     342             : 
     343          20 :     SIGDEMHeader sHeader;
     344          20 :     if (SIGDEMDataset::Identify(poOpenInfo) != TRUE || fp == nullptr)
     345             :     {
     346           0 :         return nullptr;
     347             :     }
     348             : 
     349          20 :     sHeader.Read(poOpenInfo->pabyHeader);
     350             : 
     351          20 :     if (!GDALCheckDatasetDimensions(sHeader.nCols, sHeader.nRows))
     352             :     {
     353           0 :         return nullptr;
     354             :     }
     355             : 
     356          40 :     OGRSpatialReference oSRS;
     357          20 :     oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     358             : 
     359          20 :     if (sHeader.nCoordinateSystemId > 0)
     360             :     {
     361           9 :         if (oSRS.importFromEPSG(sHeader.nCoordinateSystemId) != OGRERR_NONE)
     362             :         {
     363           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     364             :                      "SIGDEM unable to find coordinateSystemId=%d.",
     365             :                      sHeader.nCoordinateSystemId);
     366           0 :             return nullptr;
     367             :         }
     368             :     }
     369             :     else
     370             :     {
     371             :         CPLString osPrjFilename =
     372          11 :             CPLResetExtension(poOpenInfo->pszFilename, "prj");
     373             :         VSIStatBufL sStatBuf;
     374          11 :         int nRet = VSIStatL(osPrjFilename, &sStatBuf);
     375          11 :         if (nRet != 0 && VSIIsCaseSensitiveFS(osPrjFilename))
     376             :         {
     377           0 :             osPrjFilename = CPLResetExtension(poOpenInfo->pszFilename, "PRJ");
     378           0 :             nRet = VSIStatL(osPrjFilename, &sStatBuf);
     379             :         }
     380             : 
     381          11 :         if (nRet == 0)
     382             :         {
     383          11 :             char **papszPrj = CSLLoad(osPrjFilename);
     384          11 :             if (oSRS.importFromESRI(papszPrj) != OGRERR_NONE)
     385             :             {
     386           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
     387             :                          "SIGDEM unable to read projection from %s.",
     388             :                          osPrjFilename.c_str());
     389           0 :                 CSLDestroy(papszPrj);
     390           0 :                 return nullptr;
     391             :             }
     392          11 :             CSLDestroy(papszPrj);
     393             :         }
     394             :         else
     395             :         {
     396           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     397             :                      "SIGDEM unable to find projection.");
     398           0 :             return nullptr;
     399             :         }
     400             :     }
     401             : 
     402          20 :     if (sHeader.nCols > std::numeric_limits<int>::max() / (CELL_SIZE_MEM))
     403             :     {
     404           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
     405           0 :         return nullptr;
     406             :     }
     407             : 
     408          20 :     if (!RAWDatasetCheckMemoryUsage(sHeader.nCols, sHeader.nRows, 1, 4, 4,
     409          20 :                                     4 * sHeader.nCols, 0, 0, poOpenInfo->fpL))
     410             :     {
     411           0 :         return nullptr;
     412             :     }
     413          20 :     SIGDEMDataset *poDS = new SIGDEMDataset(sHeader);
     414             : 
     415          20 :     poDS->m_oSRS = std::move(oSRS);
     416             : 
     417          20 :     poDS->fpImage = poOpenInfo->fpL;
     418          20 :     poOpenInfo->fpL = nullptr;
     419          20 :     poDS->eAccess = poOpenInfo->eAccess;
     420             : 
     421          20 :     poDS->SetDescription(poOpenInfo->pszFilename);
     422          20 :     poDS->PamInitialize();
     423             : 
     424          20 :     poDS->nBands = 1;
     425          20 :     CPLErrorReset();
     426             :     SIGDEMRasterBand *poBand = new SIGDEMRasterBand(
     427          20 :         poDS, poDS->fpImage, sHeader.dfMinZ, sHeader.dfMaxZ);
     428             : 
     429          20 :     poDS->SetBand(1, poBand);
     430          20 :     if (CPLGetLastErrorType() != CE_None)
     431             :     {
     432           0 :         poDS->nBands = 1;
     433           0 :         delete poDS;
     434           0 :         return nullptr;
     435             :     }
     436             : 
     437             :     // Initialize any PAM information.
     438          20 :     poDS->TryLoadXML();
     439             : 
     440             :     // Check for overviews.
     441          20 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
     442             : 
     443          20 :     return poDS;
     444             : }
     445             : 
     446          44 : SIGDEMHeader::SIGDEMHeader()
     447             : {
     448          44 : }
     449             : 
     450          20 : bool SIGDEMHeader::Read(const GByte *pabyHeader)
     451             : {
     452             :     GByte abyHeader[HEADER_LENGTH];
     453          20 :     memcpy(abyHeader, pabyHeader, HEADER_LENGTH);
     454             : 
     455          20 :     SWAP_SIGDEM_HEADER(abyHeader);
     456          20 :     memcpy(&(this->version), abyHeader + 6, 2);
     457          20 :     memcpy(&(this->nCoordinateSystemId), abyHeader + 8, 4);
     458          20 :     memcpy(&(this->dfOffsetX), abyHeader + 12, 8);
     459          20 :     memcpy(&(this->dfScaleFactorX), abyHeader + 20, 8);
     460          20 :     memcpy(&(this->dfOffsetY), abyHeader + 28, 8);
     461          20 :     memcpy(&(this->dfScaleFactorY), abyHeader + 36, 8);
     462          20 :     memcpy(&(this->dfOffsetZ), abyHeader + 44, 8);
     463          20 :     memcpy(&(this->dfScaleFactorZ), abyHeader + 52, 8);
     464          20 :     memcpy(&(this->dfMinX), abyHeader + 60, 8);
     465          20 :     memcpy(&(this->dfMinY), abyHeader + 68, 8);
     466          20 :     memcpy(&(this->dfMinZ), abyHeader + 76, 8);
     467          20 :     memcpy(&(this->dfMaxX), abyHeader + 84, 8);
     468          20 :     memcpy(&(this->dfMaxY), abyHeader + 92, 8);
     469          20 :     memcpy(&(this->dfMaxZ), abyHeader + 100, 8);
     470          20 :     memcpy(&(this->nCols), abyHeader + 108, 4);
     471          20 :     memcpy(&(this->nRows), abyHeader + 112, 4);
     472          20 :     memcpy(&(this->dfXDim), abyHeader + 116, 8);
     473          20 :     memcpy(&(this->dfYDim), abyHeader + 124, 8);
     474             : 
     475          20 :     return true;
     476             : }
     477             : 
     478          24 : bool SIGDEMHeader::Write(VSILFILE *fp)
     479             : {
     480             :     GByte abyHeader[HEADER_LENGTH];
     481             : 
     482          24 :     memcpy(abyHeader, &(SIGDEM_FILE_TYPE), 6);
     483          24 :     memcpy(abyHeader + 6, &(this->version), 2);
     484          24 :     memcpy(abyHeader + 8, &(this->nCoordinateSystemId), 4);
     485          24 :     memcpy(abyHeader + 12, &(this->dfOffsetX), 8);
     486          24 :     memcpy(abyHeader + 20, &(this->dfScaleFactorX), 8);
     487          24 :     memcpy(abyHeader + 28, &(this->dfOffsetY), 8);
     488          24 :     memcpy(abyHeader + 36, &(this->dfScaleFactorY), 8);
     489          24 :     memcpy(abyHeader + 44, &(this->dfOffsetZ), 8);
     490          24 :     memcpy(abyHeader + 52, &(this->dfScaleFactorZ), 8);
     491          24 :     memcpy(abyHeader + 60, &(this->dfMinX), 8);
     492          24 :     memcpy(abyHeader + 68, &(this->dfMinY), 8);
     493          24 :     memcpy(abyHeader + 76, &(this->dfMinZ), 8);
     494          24 :     memcpy(abyHeader + 84, &(this->dfMaxX), 8);
     495          24 :     memcpy(abyHeader + 92, &(this->dfMaxY), 8);
     496          24 :     memcpy(abyHeader + 100, &(this->dfMaxZ), 8);
     497          24 :     memcpy(abyHeader + 108, &(this->nCols), 4);
     498          24 :     memcpy(abyHeader + 112, &(this->nRows), 4);
     499          24 :     memcpy(abyHeader + 116, &(this->dfXDim), 8);
     500          24 :     memcpy(abyHeader + 124, &(this->dfYDim), 8);
     501          24 :     SWAP_SIGDEM_HEADER(abyHeader);
     502          24 :     return VSIFWriteL(&abyHeader, HEADER_LENGTH, 1, fp) == 1;
     503             : }
     504             : 
     505          20 : SIGDEMRasterBand::SIGDEMRasterBand(SIGDEMDataset *poDSIn, VSILFILE *fpRawIn,
     506          20 :                                    double dfMinZ, double dfMaxZ)
     507          20 :     : dfOffsetZ(poDSIn->sHeader.dfOffsetZ),
     508          20 :       dfScaleFactorZ(poDSIn->sHeader.dfScaleFactorZ), fpRawL(fpRawIn)
     509             : {
     510          20 :     this->poDS = poDSIn;
     511          20 :     this->nBand = 1;
     512          20 :     this->nRasterXSize = poDSIn->GetRasterXSize();
     513          20 :     this->nRasterYSize = poDSIn->GetRasterYSize();
     514          20 :     this->nBlockXSize = this->nRasterXSize;
     515          20 :     this->nBlockYSize = 1;
     516          20 :     this->eDataType = GDT_Float64;
     517             : 
     518          20 :     this->nBlockSizeBytes = nRasterXSize * CELL_SIZE_FILE;
     519             : 
     520          20 :     this->pBlockBuffer = static_cast<int32_t *>(
     521          20 :         VSI_MALLOC2_VERBOSE(nRasterXSize, sizeof(int32_t)));
     522          20 :     SetNoDataValue(-9999);
     523          40 :     CPLString osValue;
     524          20 :     SetMetadataItem("STATISTICS_MINIMUM", osValue.Printf("%.15g", dfMinZ));
     525          20 :     SetMetadataItem("STATISTICS_MAXIMUM", osValue.Printf("%.15g", dfMaxZ));
     526          20 : }
     527             : 
     528          75 : CPLErr SIGDEMRasterBand::IReadBlock(int /*nBlockXOff*/, int nBlockYOff,
     529             :                                     void *pImage)
     530             : {
     531             : 
     532          75 :     const int nBlockIndex = nRasterYSize - nBlockYOff - 1;
     533             : 
     534          75 :     if (nLoadedBlockIndex == nBlockIndex)
     535             :     {
     536           0 :         return CE_None;
     537             :     }
     538          75 :     const vsi_l_offset nReadStart =
     539             :         HEADER_LENGTH +
     540          75 :         static_cast<vsi_l_offset>(nBlockSizeBytes) * nBlockIndex;
     541             : 
     542             :     // Seek to the correct line.
     543          75 :     if (VSIFSeekL(fpRawL, nReadStart, SEEK_SET) == -1)
     544             :     {
     545           0 :         if (poDS != nullptr && poDS->GetAccess() == GA_ReadOnly)
     546             :         {
     547           0 :             CPLError(CE_Failure, CPLE_FileIO,
     548             :                      "Failed to seek to block %d @ " CPL_FRMT_GUIB ".",
     549             :                      nBlockIndex, nReadStart);
     550           0 :             return CE_Failure;
     551             :         }
     552             :         else
     553             :         {
     554           0 :             std::fill(pBlockBuffer, pBlockBuffer + nRasterXSize, 0);
     555           0 :             nLoadedBlockIndex = nBlockIndex;
     556           0 :             return CE_None;
     557             :         }
     558             :     }
     559             :     const size_t nCellReadCount =
     560          75 :         VSIFReadL(pBlockBuffer, CELL_SIZE_FILE, nRasterXSize, fpRawL);
     561          75 :     if (nCellReadCount < static_cast<size_t>(nRasterXSize))
     562             :     {
     563           0 :         if (poDS != nullptr && poDS->GetAccess() == GA_ReadOnly)
     564             :         {
     565           0 :             CPLError(CE_Failure, CPLE_FileIO, "Failed to read block %d.",
     566             :                      nBlockIndex);
     567           0 :             return CE_Failure;
     568             :         }
     569           0 :         std::fill(pBlockBuffer + nCellReadCount, pBlockBuffer + nRasterXSize,
     570             :                   NO_DATA);
     571             :     }
     572             : 
     573          75 :     nLoadedBlockIndex = nBlockIndex;
     574             : 
     575          75 :     const int32_t *pnSourceValues = pBlockBuffer;
     576          75 :     double *padfDestValues = static_cast<double *>(pImage);
     577          75 :     double dfOffset = this->dfOffsetZ;
     578          75 :     const double dfInvScaleFactor =
     579          75 :         dfScaleFactorZ != 0.0 ? 1.0 / dfScaleFactorZ : 0.0;
     580          75 :     int nCellCount = this->nRasterXSize;
     581        1960 :     for (int i = 0; i < nCellCount; i++)
     582             :     {
     583        1885 :         int32_t nValue = CPL_MSBWORD32(*pnSourceValues);
     584        1885 :         if (nValue == NO_DATA)
     585             :         {
     586           0 :             *padfDestValues = -9999;
     587             :         }
     588             :         else
     589             :         {
     590        1885 :             *padfDestValues = dfOffset + nValue * dfInvScaleFactor;
     591             :         }
     592             : 
     593        1885 :         pnSourceValues++;
     594        1885 :         padfDestValues++;
     595             :     }
     596             : 
     597          75 :     return CE_None;
     598             : }
     599             : 
     600         185 : CPLErr SIGDEMRasterBand::IWriteBlock(int /*nBlockXOff*/, int nBlockYOff,
     601             :                                      void *pImage)
     602             : {
     603             : 
     604         185 :     const int nBlockIndex = nRasterYSize - nBlockYOff - 1;
     605             : 
     606         185 :     const double *padfSourceValues = static_cast<double *>(pImage);
     607         185 :     int32_t *pnDestValues = pBlockBuffer;
     608         185 :     double dfOffset = this->dfOffsetZ;
     609         185 :     double dfScaleFactor = this->dfScaleFactorZ;
     610         185 :     int nCellCount = this->nRasterXSize;
     611        3170 :     for (int i = 0; i < nCellCount; i++)
     612             :     {
     613        2985 :         double dfValue = *padfSourceValues;
     614             :         int32_t nValue;
     615        2985 :         if (dfValue == -9999)
     616             :         {
     617           0 :             nValue = NO_DATA;
     618             :         }
     619             :         else
     620             :         {
     621        2985 :             nValue = static_cast<int32_t>(
     622        2985 :                 std::round((dfValue - dfOffset) * dfScaleFactor));
     623             :         }
     624        2985 :         *pnDestValues = CPL_MSBWORD32(nValue);
     625        2985 :         padfSourceValues++;
     626        2985 :         pnDestValues++;
     627             :     }
     628             : 
     629         185 :     const vsi_l_offset nWriteStart =
     630             :         HEADER_LENGTH +
     631         185 :         static_cast<vsi_l_offset>(nBlockSizeBytes) * nBlockIndex;
     632             : 
     633         370 :     if (VSIFSeekL(fpRawL, nWriteStart, SEEK_SET) == -1 ||
     634         185 :         VSIFWriteL(pBlockBuffer, CELL_SIZE_FILE, nRasterXSize, fpRawL) <
     635         185 :             static_cast<size_t>(nRasterXSize))
     636             :     {
     637           0 :         CPLError(CE_Failure, CPLE_FileIO, "Failed to write block %d to file.",
     638             :                  nBlockIndex);
     639             : 
     640           0 :         return CE_Failure;
     641             :     }
     642         185 :     return CE_None;
     643             : }
     644             : 
     645          40 : SIGDEMRasterBand::~SIGDEMRasterBand()
     646             : {
     647          20 :     SIGDEMRasterBand::FlushCache(true);
     648          20 :     VSIFree(pBlockBuffer);
     649          40 : }