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