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