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