Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: HF2 driver
4 : * Purpose: GDALDataset driver for HF2/HFZ dataset.
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2010-2012, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "gdal_frmts.h"
15 : #include "gdal_pam.h"
16 : #include "ogr_spatialref.h"
17 :
18 : #include <cstdlib>
19 : #include <cmath>
20 :
21 : #include <algorithm>
22 : #include <limits>
23 :
24 : /************************************************************************/
25 : /* ==================================================================== */
26 : /* HF2Dataset */
27 : /* ==================================================================== */
28 : /************************************************************************/
29 :
30 : class HF2RasterBand;
31 :
32 : class HF2Dataset final : public GDALPamDataset
33 : {
34 : friend class HF2RasterBand;
35 :
36 : VSILFILE *fp;
37 : double adfGeoTransform[6];
38 : OGRSpatialReference m_oSRS{};
39 : vsi_l_offset *panBlockOffset; // tile 0 is a the bottom left
40 :
41 : int nTileSize;
42 : int bHasLoaderBlockMap;
43 : int LoadBlockMap();
44 :
45 : public:
46 : HF2Dataset();
47 : virtual ~HF2Dataset();
48 :
49 : virtual CPLErr GetGeoTransform(double *) override;
50 : const OGRSpatialReference *GetSpatialRef() const override;
51 :
52 : static GDALDataset *Open(GDALOpenInfo *);
53 : static int Identify(GDALOpenInfo *);
54 : static GDALDataset *CreateCopy(const char *pszFilename,
55 : GDALDataset *poSrcDS, int bStrict,
56 : char **papszOptions,
57 : GDALProgressFunc pfnProgress,
58 : void *pProgressData);
59 : };
60 :
61 : /************************************************************************/
62 : /* ==================================================================== */
63 : /* HF2RasterBand */
64 : /* ==================================================================== */
65 : /************************************************************************/
66 :
67 : class HF2RasterBand final : public GDALPamRasterBand
68 : {
69 : friend class HF2Dataset;
70 :
71 : float *pafBlockData;
72 : int nLastBlockYOffFromBottom;
73 :
74 : public:
75 : HF2RasterBand(HF2Dataset *, int, GDALDataType);
76 : virtual ~HF2RasterBand();
77 :
78 : virtual CPLErr IReadBlock(int, int, void *) override;
79 : };
80 :
81 : /************************************************************************/
82 : /* HF2RasterBand() */
83 : /************************************************************************/
84 :
85 23 : HF2RasterBand::HF2RasterBand(HF2Dataset *poDSIn, int nBandIn, GDALDataType eDT)
86 23 : : pafBlockData(nullptr), nLastBlockYOffFromBottom(-1)
87 : {
88 23 : poDS = poDSIn;
89 23 : nBand = nBandIn;
90 :
91 23 : eDataType = eDT;
92 :
93 23 : nBlockXSize = poDSIn->nTileSize;
94 23 : nBlockYSize = 1;
95 23 : }
96 :
97 : /************************************************************************/
98 : /* ~HF2RasterBand() */
99 : /************************************************************************/
100 :
101 46 : HF2RasterBand::~HF2RasterBand()
102 : {
103 23 : CPLFree(pafBlockData);
104 46 : }
105 :
106 : /************************************************************************/
107 : /* IReadBlock() */
108 : /************************************************************************/
109 :
110 764 : CPLErr HF2RasterBand::IReadBlock(int nBlockXOff, int nLineYOff, void *pImage)
111 :
112 : {
113 764 : HF2Dataset *poGDS = (HF2Dataset *)poDS;
114 :
115 764 : const int nXBlocks = DIV_ROUND_UP(nRasterXSize, poGDS->nTileSize);
116 :
117 764 : if (!poGDS->LoadBlockMap())
118 0 : return CE_Failure;
119 :
120 764 : const int nMaxTileHeight = std::min(poGDS->nTileSize, nRasterYSize);
121 764 : if (pafBlockData == nullptr)
122 : {
123 8 : if (nMaxTileHeight > 10 * 1024 * 1024 / nRasterXSize)
124 : {
125 0 : VSIFSeekL(poGDS->fp, 0, SEEK_END);
126 0 : vsi_l_offset nSize = VSIFTellL(poGDS->fp);
127 0 : if (nSize <
128 0 : static_cast<vsi_l_offset>(nMaxTileHeight) * nRasterXSize)
129 : {
130 0 : CPLError(CE_Failure, CPLE_FileIO, "File too short");
131 0 : return CE_Failure;
132 : }
133 : }
134 8 : pafBlockData =
135 8 : (float *)VSIMalloc3(sizeof(float), nRasterXSize, nMaxTileHeight);
136 8 : if (pafBlockData == nullptr)
137 0 : return CE_Failure;
138 : }
139 :
140 764 : const int nLineYOffFromBottom = nRasterYSize - 1 - nLineYOff;
141 :
142 764 : const int nBlockYOffFromBottom = nLineYOffFromBottom / nBlockXSize;
143 764 : const int nYOffInTile = nLineYOffFromBottom % nBlockXSize;
144 :
145 764 : if (nBlockYOffFromBottom != nLastBlockYOffFromBottom)
146 : {
147 10 : nLastBlockYOffFromBottom = nBlockYOffFromBottom;
148 :
149 10 : memset(pafBlockData, 0, sizeof(float) * nRasterXSize * nMaxTileHeight);
150 :
151 : /* 4 * nBlockXSize is the upper bound */
152 10 : void *pabyData = CPLMalloc(4 * nBlockXSize);
153 :
154 24 : for (int nxoff = 0; nxoff < nXBlocks; nxoff++)
155 : {
156 14 : VSIFSeekL(
157 : poGDS->fp,
158 14 : poGDS->panBlockOffset[nBlockYOffFromBottom * nXBlocks + nxoff],
159 : SEEK_SET);
160 : float fScale, fOff;
161 14 : VSIFReadL(&fScale, 4, 1, poGDS->fp);
162 14 : VSIFReadL(&fOff, 4, 1, poGDS->fp);
163 14 : CPL_LSBPTR32(&fScale);
164 14 : CPL_LSBPTR32(&fOff);
165 :
166 : const int nTileWidth =
167 14 : std::min(nBlockXSize, nRasterXSize - nxoff * nBlockXSize);
168 : const int nTileHeight = std::min(
169 14 : nBlockXSize, nRasterYSize - nBlockYOffFromBottom * nBlockXSize);
170 :
171 778 : for (int j = 0; j < nTileHeight; j++)
172 : {
173 : GByte nWordSize;
174 764 : VSIFReadL(&nWordSize, 1, 1, poGDS->fp);
175 764 : if (nWordSize != 1 && nWordSize != 2 && nWordSize != 4)
176 : {
177 0 : CPLError(CE_Failure, CPLE_AppDefined,
178 : "Unexpected word size : %d", (int)nWordSize);
179 0 : break;
180 : }
181 :
182 : GInt32 nVal;
183 764 : VSIFReadL(&nVal, 4, 1, poGDS->fp);
184 764 : CPL_LSBPTR32(&nVal);
185 764 : const size_t nToRead =
186 764 : static_cast<size_t>(nWordSize * (nTileWidth - 1));
187 764 : const size_t nRead = VSIFReadL(pabyData, 1, nToRead, poGDS->fp);
188 764 : if (nRead != nToRead)
189 : {
190 0 : CPLError(CE_Failure, CPLE_FileIO,
191 : "File too short: got %d, expected %d",
192 : static_cast<int>(nRead),
193 : static_cast<int>(nToRead));
194 0 : CPLFree(pabyData);
195 0 : return CE_Failure;
196 : }
197 : #if defined(CPL_MSB)
198 : if (nWordSize > 1)
199 : GDALSwapWords(pabyData, nWordSize, nTileWidth - 1,
200 : nWordSize);
201 : #endif
202 :
203 764 : double dfVal = nVal * (double)fScale + fOff;
204 764 : if (dfVal > std::numeric_limits<float>::max())
205 0 : dfVal = std::numeric_limits<float>::max();
206 764 : else if (dfVal < std::numeric_limits<float>::min())
207 402 : dfVal = std::numeric_limits<float>::min();
208 764 : pafBlockData[nxoff * nBlockXSize + j * nRasterXSize + 0] =
209 764 : static_cast<float>(dfVal);
210 111684 : for (int i = 1; i < nTileWidth; i++)
211 : {
212 : int nInc;
213 110920 : if (nWordSize == 1)
214 38920 : nInc = ((signed char *)pabyData)[i - 1];
215 72000 : else if (nWordSize == 2)
216 72000 : nInc = ((GInt16 *)pabyData)[i - 1];
217 : else
218 0 : nInc = ((GInt32 *)pabyData)[i - 1];
219 110920 : if ((nInc >= 0 && nVal > INT_MAX - nInc) ||
220 110920 : (nInc == INT_MIN && nVal < 0) ||
221 16066 : (nInc < 0 && nVal < INT_MIN - nInc))
222 : {
223 0 : CPLError(CE_Failure, CPLE_FileIO, "int32 overflow");
224 0 : CPLFree(pabyData);
225 0 : return CE_Failure;
226 : }
227 110920 : nVal += nInc;
228 110920 : dfVal = nVal * (double)fScale + fOff;
229 110920 : if (dfVal > std::numeric_limits<float>::max())
230 0 : dfVal = std::numeric_limits<float>::max();
231 110920 : else if (dfVal < std::numeric_limits<float>::min())
232 23778 : dfVal = std::numeric_limits<float>::min();
233 110920 : pafBlockData[nxoff * nBlockXSize + j * nRasterXSize + i] =
234 110920 : static_cast<float>(dfVal);
235 : }
236 : }
237 : }
238 :
239 10 : CPLFree(pabyData);
240 : }
241 :
242 : const int nTileWidth =
243 764 : std::min(nBlockXSize, nRasterXSize - nBlockXOff * nBlockXSize);
244 764 : memcpy(pImage,
245 764 : pafBlockData + nBlockXOff * nBlockXSize + nYOffInTile * nRasterXSize,
246 764 : nTileWidth * sizeof(float));
247 :
248 764 : return CE_None;
249 : }
250 :
251 : /************************************************************************/
252 : /* ~HF2Dataset() */
253 : /************************************************************************/
254 :
255 23 : HF2Dataset::HF2Dataset()
256 : : fp(nullptr), panBlockOffset(nullptr), nTileSize(0),
257 23 : bHasLoaderBlockMap(FALSE)
258 : {
259 23 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
260 23 : adfGeoTransform[0] = 0;
261 23 : adfGeoTransform[1] = 1;
262 23 : adfGeoTransform[2] = 0;
263 23 : adfGeoTransform[3] = 0;
264 23 : adfGeoTransform[4] = 0;
265 23 : adfGeoTransform[5] = 1;
266 23 : }
267 :
268 : /************************************************************************/
269 : /* ~HF2Dataset() */
270 : /************************************************************************/
271 :
272 46 : HF2Dataset::~HF2Dataset()
273 :
274 : {
275 23 : FlushCache(true);
276 23 : CPLFree(panBlockOffset);
277 23 : if (fp)
278 23 : VSIFCloseL(fp);
279 46 : }
280 :
281 : /************************************************************************/
282 : /* LoadBlockMap() */
283 : /************************************************************************/
284 :
285 764 : int HF2Dataset::LoadBlockMap()
286 : {
287 764 : if (bHasLoaderBlockMap)
288 756 : return panBlockOffset != nullptr;
289 :
290 8 : bHasLoaderBlockMap = TRUE;
291 :
292 8 : const int nXBlocks = (nRasterXSize + nTileSize - 1) / nTileSize;
293 8 : const int nYBlocks = (nRasterYSize + nTileSize - 1) / nTileSize;
294 8 : if (nXBlocks * nYBlocks > 1000000)
295 : {
296 0 : vsi_l_offset nCurOff = VSIFTellL(fp);
297 0 : VSIFSeekL(fp, 0, SEEK_END);
298 0 : vsi_l_offset nSize = VSIFTellL(fp);
299 0 : VSIFSeekL(fp, nCurOff, SEEK_SET);
300 : // Check that the file is big enough to have 8 bytes per block
301 0 : if (static_cast<vsi_l_offset>(nXBlocks) * nYBlocks > nSize / 8)
302 : {
303 0 : return FALSE;
304 : }
305 : }
306 8 : panBlockOffset =
307 8 : (vsi_l_offset *)VSIMalloc3(sizeof(vsi_l_offset), nXBlocks, nYBlocks);
308 8 : if (panBlockOffset == nullptr)
309 : {
310 0 : return FALSE;
311 : }
312 18 : for (int j = 0; j < nYBlocks; j++)
313 : {
314 24 : for (int i = 0; i < nXBlocks; i++)
315 : {
316 14 : vsi_l_offset nOff = VSIFTellL(fp);
317 14 : panBlockOffset[j * nXBlocks + i] = nOff;
318 : // VSIFSeekL(fp, 4 + 4, SEEK_CUR);
319 : float fScale, fOff;
320 14 : VSIFReadL(&fScale, 4, 1, fp);
321 14 : VSIFReadL(&fOff, 4, 1, fp);
322 14 : CPL_LSBPTR32(&fScale);
323 14 : CPL_LSBPTR32(&fOff);
324 : // printf("fScale = %f, fOff = %f\n", fScale, fOff);
325 14 : const int nCols = std::min(nTileSize, nRasterXSize - nTileSize * i);
326 : const int nLines =
327 14 : std::min(nTileSize, nRasterYSize - nTileSize * j);
328 778 : for (int k = 0; k < nLines; k++)
329 : {
330 : GByte nWordSize;
331 764 : if (VSIFReadL(&nWordSize, 1, 1, fp) != 1)
332 : {
333 0 : CPLError(CE_Failure, CPLE_FileIO, "File too short");
334 0 : VSIFree(panBlockOffset);
335 0 : panBlockOffset = nullptr;
336 0 : return FALSE;
337 : }
338 : // printf("nWordSize=%d\n", nWordSize);
339 764 : if (nWordSize == 1 || nWordSize == 2 || nWordSize == 4)
340 764 : VSIFSeekL(
341 : fp,
342 764 : static_cast<vsi_l_offset>(4 + nWordSize * (nCols - 1)),
343 : SEEK_CUR);
344 : else
345 : {
346 0 : CPLError(CE_Failure, CPLE_AppDefined,
347 : "Got unexpected byte depth (%d) for block (%d, "
348 : "%d) line %d",
349 : (int)nWordSize, i, j, k);
350 0 : VSIFree(panBlockOffset);
351 0 : panBlockOffset = nullptr;
352 0 : return FALSE;
353 : }
354 : }
355 : }
356 : }
357 :
358 8 : return TRUE;
359 : }
360 :
361 : /************************************************************************/
362 : /* GetSpatialRef() */
363 : /************************************************************************/
364 :
365 15 : const OGRSpatialReference *HF2Dataset::GetSpatialRef() const
366 :
367 : {
368 15 : if (!m_oSRS.IsEmpty())
369 15 : return &m_oSRS;
370 0 : return GDALPamDataset::GetSpatialRef();
371 : }
372 :
373 : /************************************************************************/
374 : /* Identify() */
375 : /************************************************************************/
376 :
377 50703 : int HF2Dataset::Identify(GDALOpenInfo *poOpenInfo)
378 : {
379 :
380 50703 : GDALOpenInfo *poOpenInfoToDelete = nullptr;
381 : /* GZipped .hf2 files are common, so automagically open them */
382 : /* if the /vsigzip/ has not been explicitly passed */
383 101402 : CPLString osFilename; // keep in that scope
384 50697 : if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
385 50695 : (strlen(poOpenInfo->pszFilename) > 6 &&
386 49838 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6,
387 101406 : "hf2.gz"))) &&
388 8 : !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
389 : {
390 7 : osFilename = "/vsigzip/";
391 7 : osFilename += poOpenInfo->pszFilename;
392 7 : poOpenInfo = poOpenInfoToDelete = new GDALOpenInfo(
393 7 : osFilename.c_str(), GA_ReadOnly, poOpenInfo->GetSiblingFiles());
394 : }
395 :
396 50703 : if (poOpenInfo->nHeaderBytes < 28)
397 : {
398 47847 : delete poOpenInfoToDelete;
399 47843 : return FALSE;
400 : }
401 :
402 2856 : if (memcmp(poOpenInfo->pabyHeader, "HF2\0\0\0\0", 6) != 0)
403 : {
404 2825 : delete poOpenInfoToDelete;
405 2825 : return FALSE;
406 : }
407 :
408 31 : delete poOpenInfoToDelete;
409 31 : return TRUE;
410 : }
411 :
412 : /************************************************************************/
413 : /* Open() */
414 : /************************************************************************/
415 :
416 23 : GDALDataset *HF2Dataset::Open(GDALOpenInfo *poOpenInfo)
417 :
418 : {
419 46 : CPLString osOriginalFilename(poOpenInfo->pszFilename);
420 :
421 23 : if (!Identify(poOpenInfo))
422 0 : return nullptr;
423 :
424 23 : GDALOpenInfo *poOpenInfoToDelete = nullptr;
425 : /* GZipped .hf2 files are common, so automagically open them */
426 : /* if the /vsigzip/ has not been explicitly passed */
427 46 : CPLString osFilename(poOpenInfo->pszFilename);
428 23 : if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
429 20 : (strlen(poOpenInfo->pszFilename) > 6 &&
430 20 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6,
431 46 : "hf2.gz"))) &&
432 3 : !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
433 : {
434 2 : osFilename = "/vsigzip/";
435 2 : osFilename += poOpenInfo->pszFilename;
436 2 : poOpenInfo = poOpenInfoToDelete = new GDALOpenInfo(
437 2 : osFilename.c_str(), GA_ReadOnly, poOpenInfo->GetSiblingFiles());
438 : }
439 :
440 : /* -------------------------------------------------------------------- */
441 : /* Parse header */
442 : /* -------------------------------------------------------------------- */
443 :
444 : int nXSize, nYSize;
445 23 : memcpy(&nXSize, poOpenInfo->pabyHeader + 6, 4);
446 23 : CPL_LSBPTR32(&nXSize);
447 23 : memcpy(&nYSize, poOpenInfo->pabyHeader + 10, 4);
448 23 : CPL_LSBPTR32(&nYSize);
449 :
450 : GUInt16 nTileSize;
451 23 : memcpy(&nTileSize, poOpenInfo->pabyHeader + 14, 2);
452 23 : CPL_LSBPTR16(&nTileSize);
453 :
454 : float fVertPres, fHorizScale;
455 23 : memcpy(&fVertPres, poOpenInfo->pabyHeader + 16, 4);
456 23 : CPL_LSBPTR32(&fVertPres);
457 23 : memcpy(&fHorizScale, poOpenInfo->pabyHeader + 20, 4);
458 23 : CPL_LSBPTR32(&fHorizScale);
459 :
460 : GUInt32 nExtendedHeaderLen;
461 23 : memcpy(&nExtendedHeaderLen, poOpenInfo->pabyHeader + 24, 4);
462 23 : CPL_LSBPTR32(&nExtendedHeaderLen);
463 :
464 23 : delete poOpenInfoToDelete;
465 23 : poOpenInfoToDelete = nullptr;
466 :
467 23 : if (nTileSize < 8)
468 0 : return nullptr;
469 23 : if (nXSize <= 0 || nXSize > INT_MAX - nTileSize || nYSize <= 0 ||
470 23 : nYSize > INT_MAX - nTileSize)
471 0 : return nullptr;
472 : /* To avoid later potential int overflows */
473 23 : if (nExtendedHeaderLen > 1024 * 65536)
474 0 : return nullptr;
475 :
476 23 : if (!GDALCheckDatasetDimensions(nXSize, nYSize))
477 : {
478 0 : return nullptr;
479 : }
480 23 : const int nXBlocks = (nXSize + nTileSize - 1) / nTileSize;
481 23 : const int nYBlocks = (nYSize + nTileSize - 1) / nTileSize;
482 23 : if (nXBlocks > INT_MAX / nYBlocks)
483 : {
484 0 : return nullptr;
485 : }
486 :
487 : /* -------------------------------------------------------------------- */
488 : /* Parse extended blocks */
489 : /* -------------------------------------------------------------------- */
490 :
491 23 : VSILFILE *fp = VSIFOpenL(osFilename.c_str(), "rb");
492 23 : if (fp == nullptr)
493 0 : return nullptr;
494 :
495 23 : VSIFSeekL(fp, 28, SEEK_SET);
496 :
497 23 : int bHasExtent = FALSE;
498 23 : double dfMinX = 0.0;
499 23 : double dfMaxX = 0.0;
500 23 : double dfMinY = 0.0;
501 23 : double dfMaxY = 0.0;
502 23 : int bHasUTMZone = FALSE;
503 23 : GInt16 nUTMZone = 0;
504 23 : int bHasEPSGDatumCode = FALSE;
505 23 : GInt16 nEPSGDatumCode = 0;
506 23 : int bHasEPSGCode = FALSE;
507 23 : GInt16 nEPSGCode = 0;
508 23 : int bHasRelativePrecision = FALSE;
509 23 : float fRelativePrecision = 0.0f;
510 23 : char szApplicationName[256] = {0};
511 :
512 23 : GUInt32 nExtendedHeaderOff = 0;
513 84 : while (nExtendedHeaderOff < nExtendedHeaderLen)
514 : {
515 : char pabyBlockHeader[24];
516 61 : VSIFReadL(pabyBlockHeader, 24, 1, fp);
517 :
518 : char szBlockName[16 + 1];
519 61 : memcpy(szBlockName, pabyBlockHeader + 4, 16);
520 61 : szBlockName[16] = 0;
521 : GUInt32 nBlockSize;
522 61 : memcpy(&nBlockSize, pabyBlockHeader + 20, 4);
523 61 : CPL_LSBPTR32(&nBlockSize);
524 61 : if (nBlockSize > 65536)
525 0 : break;
526 :
527 61 : nExtendedHeaderOff += 24 + nBlockSize;
528 :
529 61 : if (strcmp(szBlockName, "georef-extents") == 0 && nBlockSize == 34)
530 : {
531 : char pabyBlockData[34];
532 23 : VSIFReadL(pabyBlockData, 34, 1, fp);
533 :
534 23 : memcpy(&dfMinX, pabyBlockData + 2, 8);
535 23 : CPL_LSBPTR64(&dfMinX);
536 23 : memcpy(&dfMaxX, pabyBlockData + 2 + 8, 8);
537 23 : CPL_LSBPTR64(&dfMaxX);
538 23 : memcpy(&dfMinY, pabyBlockData + 2 + 8 + 8, 8);
539 23 : CPL_LSBPTR64(&dfMinY);
540 23 : memcpy(&dfMaxY, pabyBlockData + 2 + 8 + 8 + 8, 8);
541 23 : CPL_LSBPTR64(&dfMaxY);
542 :
543 23 : bHasExtent = TRUE;
544 : }
545 38 : else if (strcmp(szBlockName, "georef-utm") == 0 && nBlockSize == 2)
546 : {
547 9 : VSIFReadL(&nUTMZone, 2, 1, fp);
548 9 : CPL_LSBPTR16(&nUTMZone);
549 9 : CPLDebug("HF2", "UTM Zone = %d", nUTMZone);
550 :
551 9 : bHasUTMZone = TRUE;
552 : }
553 29 : else if (strcmp(szBlockName, "georef-datum") == 0 && nBlockSize == 2)
554 : {
555 23 : VSIFReadL(&nEPSGDatumCode, 2, 1, fp);
556 23 : CPL_LSBPTR16(&nEPSGDatumCode);
557 23 : CPLDebug("HF2", "EPSG Datum Code = %d", nEPSGDatumCode);
558 :
559 23 : bHasEPSGDatumCode = TRUE;
560 : }
561 6 : else if (strcmp(szBlockName, "georef-epsg-prj") == 0 && nBlockSize == 2)
562 : {
563 6 : VSIFReadL(&nEPSGCode, 2, 1, fp);
564 6 : CPL_LSBPTR16(&nEPSGCode);
565 6 : CPLDebug("HF2", "EPSG Code = %d", nEPSGCode);
566 :
567 6 : bHasEPSGCode = TRUE;
568 : }
569 0 : else if (strcmp(szBlockName, "precis-rel") == 0 && nBlockSize == 4)
570 : {
571 0 : VSIFReadL(&fRelativePrecision, 4, 1, fp);
572 0 : CPL_LSBPTR32(&fRelativePrecision);
573 :
574 0 : bHasRelativePrecision = TRUE;
575 : }
576 0 : else if (strcmp(szBlockName, "app-name") == 0 &&
577 0 : nBlockSize < sizeof(szApplicationName))
578 : {
579 0 : VSIFReadL(szApplicationName, nBlockSize, 1, fp);
580 0 : szApplicationName[nBlockSize] = 0;
581 : }
582 : else
583 : {
584 0 : CPLDebug("HF2", "Skipping block %s", szBlockName);
585 0 : VSIFSeekL(fp, nBlockSize, SEEK_CUR);
586 : }
587 : }
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Create a corresponding GDALDataset. */
591 : /* -------------------------------------------------------------------- */
592 23 : HF2Dataset *poDS = new HF2Dataset();
593 23 : poDS->fp = fp;
594 23 : poDS->nRasterXSize = nXSize;
595 23 : poDS->nRasterYSize = nYSize;
596 23 : poDS->nTileSize = nTileSize;
597 23 : CPLDebug("HF2", "nXSize = %d, nYSize = %d, nTileSize = %d", nXSize, nYSize,
598 : nTileSize);
599 23 : if (bHasExtent)
600 : {
601 23 : poDS->adfGeoTransform[0] = dfMinX;
602 23 : poDS->adfGeoTransform[3] = dfMaxY;
603 23 : poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nXSize;
604 23 : poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nYSize;
605 : }
606 : else
607 : {
608 0 : poDS->adfGeoTransform[1] = fHorizScale;
609 0 : poDS->adfGeoTransform[5] = fHorizScale;
610 : }
611 :
612 23 : if (bHasEPSGCode)
613 : {
614 6 : poDS->m_oSRS.importFromEPSG(nEPSGCode);
615 : }
616 : else
617 : {
618 17 : bool bHasSRS = false;
619 34 : OGRSpatialReference oSRS;
620 17 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
621 17 : oSRS.SetGeogCS("unknown", "unknown", "unknown", SRS_WGS84_SEMIMAJOR,
622 : SRS_WGS84_INVFLATTENING);
623 17 : if (bHasEPSGDatumCode)
624 : {
625 17 : if (nEPSGDatumCode == 23 || nEPSGDatumCode == 6326)
626 : {
627 14 : bHasSRS = true;
628 14 : oSRS.SetWellKnownGeogCS("WGS84");
629 : }
630 3 : else if (nEPSGDatumCode >= 6000)
631 : {
632 : char szName[32];
633 3 : snprintf(szName, sizeof(szName), "EPSG:%d",
634 3 : nEPSGDatumCode - 2000);
635 3 : oSRS.SetWellKnownGeogCS(szName);
636 3 : bHasSRS = true;
637 : }
638 : }
639 :
640 17 : if (bHasUTMZone && std::abs(nUTMZone) >= 1 && std::abs(nUTMZone) <= 60)
641 : {
642 3 : bHasSRS = true;
643 3 : oSRS.SetUTM(std::abs(static_cast<int>(nUTMZone)), nUTMZone > 0);
644 : }
645 17 : if (bHasSRS)
646 17 : poDS->m_oSRS = std::move(oSRS);
647 : }
648 :
649 : /* -------------------------------------------------------------------- */
650 : /* Create band information objects. */
651 : /* -------------------------------------------------------------------- */
652 23 : poDS->nBands = 1;
653 46 : for (int i = 0; i < poDS->nBands; i++)
654 : {
655 23 : poDS->SetBand(i + 1, new HF2RasterBand(poDS, i + 1, GDT_Float32));
656 23 : poDS->GetRasterBand(i + 1)->SetUnitType("m");
657 : }
658 :
659 23 : if (szApplicationName[0] != '\0')
660 0 : poDS->SetMetadataItem("APPLICATION_NAME", szApplicationName);
661 23 : poDS->SetMetadataItem("VERTICAL_PRECISION",
662 46 : CPLString().Printf("%f", fVertPres));
663 23 : if (bHasRelativePrecision)
664 0 : poDS->SetMetadataItem("RELATIVE_VERTICAL_PRECISION",
665 0 : CPLString().Printf("%f", fRelativePrecision));
666 :
667 : /* -------------------------------------------------------------------- */
668 : /* Initialize any PAM information. */
669 : /* -------------------------------------------------------------------- */
670 23 : poDS->SetDescription(osOriginalFilename.c_str());
671 23 : poDS->TryLoadXML();
672 :
673 : /* -------------------------------------------------------------------- */
674 : /* Support overviews. */
675 : /* -------------------------------------------------------------------- */
676 23 : poDS->oOvManager.Initialize(poDS, osOriginalFilename.c_str());
677 23 : return poDS;
678 : }
679 :
680 : /************************************************************************/
681 : /* GetGeoTransform() */
682 : /************************************************************************/
683 :
684 16 : CPLErr HF2Dataset::GetGeoTransform(double *padfTransform)
685 :
686 : {
687 16 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
688 :
689 16 : return CE_None;
690 : }
691 :
692 36065 : static void WriteShort(VSILFILE *fp, GInt16 val)
693 : {
694 36065 : CPL_LSBPTR16(&val);
695 36065 : VSIFWriteL(&val, 2, 1, fp);
696 36065 : }
697 :
698 572 : static void WriteInt(VSILFILE *fp, GInt32 val)
699 : {
700 572 : CPL_LSBPTR32(&val);
701 572 : VSIFWriteL(&val, 4, 1, fp);
702 572 : }
703 :
704 66 : static void WriteFloat(VSILFILE *fp, float val)
705 : {
706 66 : CPL_LSBPTR32(&val);
707 66 : VSIFWriteL(&val, 4, 1, fp);
708 66 : }
709 :
710 60 : static void WriteDouble(VSILFILE *fp, double val)
711 : {
712 60 : CPL_LSBPTR64(&val);
713 60 : VSIFWriteL(&val, 8, 1, fp);
714 60 : }
715 :
716 : /************************************************************************/
717 : /* CreateCopy() */
718 : /************************************************************************/
719 :
720 23 : GDALDataset *HF2Dataset::CreateCopy(const char *pszFilename,
721 : GDALDataset *poSrcDS, int bStrict,
722 : char **papszOptions,
723 : GDALProgressFunc pfnProgress,
724 : void *pProgressData)
725 : {
726 : /* -------------------------------------------------------------------- */
727 : /* Some some rudimentary checks */
728 : /* -------------------------------------------------------------------- */
729 23 : int nBands = poSrcDS->GetRasterCount();
730 23 : if (nBands == 0)
731 : {
732 1 : CPLError(
733 : CE_Failure, CPLE_NotSupported,
734 : "HF2 driver does not support source dataset with zero band.\n");
735 1 : return nullptr;
736 : }
737 :
738 22 : if (nBands != 1)
739 : {
740 4 : CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
741 : "HF2 driver only uses the first band of the dataset.\n");
742 4 : if (bStrict)
743 4 : return nullptr;
744 : }
745 :
746 18 : if (pfnProgress && !pfnProgress(0.0, nullptr, pProgressData))
747 0 : return nullptr;
748 :
749 : /* -------------------------------------------------------------------- */
750 : /* Get source dataset info */
751 : /* -------------------------------------------------------------------- */
752 :
753 18 : const int nXSize = poSrcDS->GetRasterXSize();
754 18 : const int nYSize = poSrcDS->GetRasterYSize();
755 : double adfGeoTransform[6];
756 18 : poSrcDS->GetGeoTransform(adfGeoTransform);
757 18 : int bHasGeoTransform =
758 18 : !(adfGeoTransform[0] == 0 && adfGeoTransform[1] == 1 &&
759 0 : adfGeoTransform[2] == 0 && adfGeoTransform[3] == 0 &&
760 0 : adfGeoTransform[4] == 0 && adfGeoTransform[5] == 1);
761 18 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
762 : {
763 0 : CPLError(CE_Failure, CPLE_NotSupported,
764 : "HF2 driver does not support CreateCopy() from skewed or "
765 : "rotated dataset.\n");
766 0 : return nullptr;
767 : }
768 :
769 18 : GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
770 : GDALDataType eReqDT;
771 18 : float fVertPres = (float)0.01;
772 18 : if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16)
773 : {
774 8 : fVertPres = 1;
775 8 : eReqDT = GDT_Int16;
776 : }
777 : else
778 10 : eReqDT = GDT_Float32;
779 :
780 : /* -------------------------------------------------------------------- */
781 : /* Read creation options */
782 : /* -------------------------------------------------------------------- */
783 18 : const char *pszCompressed = CSLFetchNameValue(papszOptions, "COMPRESS");
784 18 : bool bCompress = false;
785 18 : if (pszCompressed)
786 1 : bCompress = CPLTestBool(pszCompressed);
787 :
788 : const char *pszVerticalPrecision =
789 18 : CSLFetchNameValue(papszOptions, "VERTICAL_PRECISION");
790 18 : if (pszVerticalPrecision)
791 : {
792 0 : fVertPres = (float)CPLAtofM(pszVerticalPrecision);
793 0 : if (fVertPres <= 0)
794 : {
795 0 : CPLError(
796 : CE_Warning, CPLE_AppDefined,
797 : "Unsupported value for VERTICAL_PRECISION. Defaulting to 0.01");
798 0 : fVertPres = (float)0.01;
799 : }
800 0 : if (eReqDT == GDT_Int16 && fVertPres > 1)
801 0 : eReqDT = GDT_Float32;
802 : }
803 :
804 18 : const char *pszBlockSize = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
805 18 : int nTileSize = 256;
806 18 : if (pszBlockSize)
807 : {
808 1 : nTileSize = atoi(pszBlockSize);
809 1 : if (nTileSize < 8 || nTileSize > 4096)
810 : {
811 0 : CPLError(CE_Warning, CPLE_AppDefined,
812 : "Unsupported value for BLOCKSIZE. Defaulting to 256");
813 0 : nTileSize = 256;
814 : }
815 : }
816 :
817 : /* -------------------------------------------------------------------- */
818 : /* Parse source dataset georeferencing info */
819 : /* -------------------------------------------------------------------- */
820 :
821 18 : int nExtendedHeaderLen = 0;
822 18 : if (bHasGeoTransform)
823 18 : nExtendedHeaderLen += 58;
824 18 : const char *pszProjectionRef = poSrcDS->GetProjectionRef();
825 18 : int nDatumCode = -2;
826 18 : int nUTMZone = 0;
827 18 : int bNorth = FALSE;
828 18 : int nEPSGCode = 0;
829 18 : int nExtentUnits = 1;
830 18 : if (pszProjectionRef != nullptr && pszProjectionRef[0] != '\0')
831 : {
832 36 : OGRSpatialReference oSRS;
833 18 : if (oSRS.importFromWkt(pszProjectionRef) == OGRERR_NONE)
834 : {
835 18 : const char *pszValue = nullptr;
836 22 : if (oSRS.GetAuthorityName("GEOGCS|DATUM") != nullptr &&
837 4 : EQUAL(oSRS.GetAuthorityName("GEOGCS|DATUM"), "EPSG"))
838 4 : nDatumCode = atoi(oSRS.GetAuthorityCode("GEOGCS|DATUM"));
839 14 : else if ((pszValue = oSRS.GetAttrValue("GEOGCS|DATUM")) != nullptr)
840 : {
841 14 : if (strstr(pszValue, "WGS") && strstr(pszValue, "84"))
842 14 : nDatumCode = 6326;
843 : }
844 :
845 18 : nUTMZone = oSRS.GetUTMZone(&bNorth);
846 : }
847 20 : if (oSRS.GetAuthorityName("PROJCS") != nullptr &&
848 2 : EQUAL(oSRS.GetAuthorityName("PROJCS"), "EPSG"))
849 2 : nEPSGCode = atoi(oSRS.GetAuthorityCode("PROJCS"));
850 :
851 18 : if (oSRS.IsGeographic())
852 : {
853 15 : nExtentUnits = 0;
854 : }
855 : else
856 : {
857 3 : const double dfLinear = oSRS.GetLinearUnits();
858 :
859 3 : if (std::abs(dfLinear - 0.3048) < 0.0000001)
860 0 : nExtentUnits = 2;
861 3 : else if (std::abs(dfLinear - CPLAtof(SRS_UL_US_FOOT_CONV)) <
862 : 0.00000001)
863 0 : nExtentUnits = 3;
864 : else
865 3 : nExtentUnits = 1;
866 : }
867 : }
868 18 : if (nDatumCode != -2)
869 18 : nExtendedHeaderLen += 26;
870 18 : if (nUTMZone != 0)
871 3 : nExtendedHeaderLen += 26;
872 18 : if (nEPSGCode)
873 2 : nExtendedHeaderLen += 26;
874 :
875 : /* -------------------------------------------------------------------- */
876 : /* Create target file */
877 : /* -------------------------------------------------------------------- */
878 :
879 36 : CPLString osFilename;
880 18 : if (bCompress)
881 : {
882 1 : osFilename = "/vsigzip/";
883 1 : osFilename += pszFilename;
884 : }
885 : else
886 17 : osFilename = pszFilename;
887 18 : VSILFILE *fp = VSIFOpenL(osFilename.c_str(), "wb");
888 18 : if (fp == nullptr)
889 : {
890 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszFilename);
891 3 : return nullptr;
892 : }
893 :
894 : /* -------------------------------------------------------------------- */
895 : /* Write header */
896 : /* -------------------------------------------------------------------- */
897 :
898 15 : VSIFWriteL("HF2\0", 4, 1, fp);
899 15 : WriteShort(fp, 0);
900 15 : WriteInt(fp, nXSize);
901 15 : WriteInt(fp, nYSize);
902 15 : WriteShort(fp, (GInt16)nTileSize);
903 15 : WriteFloat(fp, fVertPres);
904 15 : const float fHorizScale =
905 15 : (float)((fabs(adfGeoTransform[1]) + fabs(adfGeoTransform[5])) / 2);
906 15 : WriteFloat(fp, fHorizScale);
907 15 : WriteInt(fp, nExtendedHeaderLen);
908 :
909 : /* -------------------------------------------------------------------- */
910 : /* Write extended header */
911 : /* -------------------------------------------------------------------- */
912 :
913 : char szBlockName[16 + 1];
914 15 : if (bHasGeoTransform)
915 : {
916 15 : VSIFWriteL("bin\0", 4, 1, fp);
917 15 : memset(szBlockName, 0, 16 + 1);
918 15 : strcpy(szBlockName, "georef-extents");
919 15 : VSIFWriteL(szBlockName, 16, 1, fp);
920 15 : WriteInt(fp, 34);
921 15 : WriteShort(fp, (GInt16)nExtentUnits);
922 15 : WriteDouble(fp, adfGeoTransform[0]);
923 15 : WriteDouble(fp, adfGeoTransform[0] + nXSize * adfGeoTransform[1]);
924 15 : WriteDouble(fp, adfGeoTransform[3] + nYSize * adfGeoTransform[5]);
925 15 : WriteDouble(fp, adfGeoTransform[3]);
926 : }
927 15 : if (nUTMZone != 0)
928 : {
929 3 : VSIFWriteL("bin\0", 4, 1, fp);
930 3 : memset(szBlockName, 0, 16 + 1);
931 3 : strcpy(szBlockName, "georef-utm");
932 3 : VSIFWriteL(szBlockName, 16, 1, fp);
933 3 : WriteInt(fp, 2);
934 3 : WriteShort(fp, (GInt16)((bNorth) ? nUTMZone : -nUTMZone));
935 : }
936 15 : if (nDatumCode != -2)
937 : {
938 15 : VSIFWriteL("bin\0", 4, 1, fp);
939 15 : memset(szBlockName, 0, 16 + 1);
940 15 : strcpy(szBlockName, "georef-datum");
941 15 : VSIFWriteL(szBlockName, 16, 1, fp);
942 15 : WriteInt(fp, 2);
943 15 : WriteShort(fp, (GInt16)nDatumCode);
944 : }
945 15 : if (nEPSGCode != 0)
946 : {
947 2 : VSIFWriteL("bin\0", 4, 1, fp);
948 2 : memset(szBlockName, 0, 16 + 1);
949 2 : strcpy(szBlockName, "georef-epsg-prj");
950 2 : VSIFWriteL(szBlockName, 16, 1, fp);
951 2 : WriteInt(fp, 2);
952 2 : WriteShort(fp, (GInt16)nEPSGCode);
953 : }
954 :
955 : /* -------------------------------------------------------------------- */
956 : /* Copy imagery */
957 : /* -------------------------------------------------------------------- */
958 15 : const int nXBlocks = (nXSize + nTileSize - 1) / nTileSize;
959 15 : const int nYBlocks = (nYSize + nTileSize - 1) / nTileSize;
960 :
961 15 : void *pTileBuffer = VSI_MALLOC3_VERBOSE(nTileSize, nTileSize,
962 : GDALGetDataTypeSizeBytes(eReqDT));
963 15 : if (pTileBuffer == nullptr)
964 : {
965 0 : VSIFCloseL(fp);
966 0 : return nullptr;
967 : }
968 :
969 15 : CPLErr eErr = CE_None;
970 31 : for (int j = 0; j < nYBlocks && eErr == CE_None; j++)
971 : {
972 34 : for (int i = 0; i < nXBlocks && eErr == CE_None; i++)
973 : {
974 18 : const int nReqXSize = std::min(nTileSize, nXSize - i * nTileSize);
975 18 : const int nReqYSize = std::min(nTileSize, nYSize - j * nTileSize);
976 36 : eErr = poSrcDS->GetRasterBand(1)->RasterIO(
977 : GF_Read, i * nTileSize,
978 18 : std::max(0, nYSize - (j + 1) * nTileSize), nReqXSize, nReqYSize,
979 : pTileBuffer, nReqXSize, nReqYSize, eReqDT, 0, 0, nullptr);
980 18 : if (eErr != CE_None)
981 0 : break;
982 :
983 18 : if (eReqDT == GDT_Int16)
984 : {
985 8 : WriteFloat(fp, 1); /* scale */
986 8 : WriteFloat(fp, 0); /* offset */
987 209 : for (int k = 0; k < nReqYSize; k++)
988 : {
989 201 : int nLastVal =
990 : ((short *)
991 201 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
992 201 : GByte nWordSize = 1;
993 15641 : for (int l = 1; l < nReqXSize; l++)
994 : {
995 15440 : const int nVal =
996 : ((short *)
997 15440 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize +
998 15440 : l];
999 15440 : const int nDiff = nVal - nLastVal;
1000 15440 : if (nDiff < -32768 || nDiff > 32767)
1001 : {
1002 0 : nWordSize = 4;
1003 0 : break;
1004 : }
1005 15440 : if (nDiff < -128 || nDiff > 127)
1006 0 : nWordSize = 2;
1007 15440 : nLastVal = nVal;
1008 : }
1009 :
1010 201 : VSIFWriteL(&nWordSize, 1, 1, fp);
1011 201 : nLastVal =
1012 : ((short *)
1013 201 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
1014 201 : WriteInt(fp, nLastVal);
1015 15641 : for (int l = 1; l < nReqXSize; l++)
1016 : {
1017 15440 : const int nVal =
1018 : ((short *)
1019 15440 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize +
1020 15440 : l];
1021 15440 : const int nDiff = nVal - nLastVal;
1022 15440 : if (nWordSize == 1)
1023 : {
1024 15440 : CPLAssert(nDiff >= -128 && nDiff <= 127);
1025 15440 : signed char chDiff = (signed char)nDiff;
1026 15440 : VSIFWriteL(&chDiff, 1, 1, fp);
1027 : }
1028 0 : else if (nWordSize == 2)
1029 : {
1030 0 : CPLAssert(nDiff >= -32768 && nDiff <= 32767);
1031 0 : WriteShort(fp, (short)nDiff);
1032 : }
1033 : else
1034 : {
1035 0 : WriteInt(fp, nDiff);
1036 : }
1037 15440 : nLastVal = nVal;
1038 : }
1039 : }
1040 : }
1041 : else
1042 : {
1043 10 : float fMinVal = ((float *)pTileBuffer)[0];
1044 10 : float fMaxVal = fMinVal;
1045 41301 : for (int k = 1; k < nReqYSize * nReqXSize; k++)
1046 : {
1047 41291 : float fVal = ((float *)pTileBuffer)[k];
1048 41291 : if (std::isnan(fVal))
1049 : {
1050 0 : CPLError(CE_Failure, CPLE_NotSupported,
1051 : "NaN value found");
1052 0 : eErr = CE_Failure;
1053 0 : break;
1054 : }
1055 41291 : if (fVal < fMinVal)
1056 0 : fMinVal = fVal;
1057 41291 : if (fVal > fMaxVal)
1058 180 : fMaxVal = fVal;
1059 : }
1060 10 : if (eErr == CE_Failure)
1061 0 : break;
1062 :
1063 10 : float fIntRange = (fMaxVal - fMinVal) / fVertPres;
1064 10 : if (fIntRange >
1065 10 : static_cast<float>(std::numeric_limits<int>::max()))
1066 : {
1067 0 : CPLError(CE_Failure, CPLE_NotSupported,
1068 : "VERTICAL_PRECISION too small regarding actual "
1069 : "range of values");
1070 0 : eErr = CE_Failure;
1071 0 : break;
1072 : }
1073 10 : float fScale =
1074 10 : (fMinVal == fMaxVal) ? 1 : (fMaxVal - fMinVal) / fIntRange;
1075 10 : if (fScale == 0.0f)
1076 : {
1077 0 : CPLError(CE_Failure, CPLE_NotSupported, "Scale == 0.0f");
1078 0 : eErr = CE_Failure;
1079 0 : break;
1080 : }
1081 10 : float fOffset = fMinVal;
1082 10 : WriteFloat(fp, fScale); /* scale */
1083 10 : WriteFloat(fp, fOffset); /* offset */
1084 301 : for (int k = 0; k < nReqYSize; k++)
1085 : {
1086 291 : float fLastVal =
1087 : ((float *)
1088 291 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
1089 291 : float fIntLastVal = (fLastVal - fOffset) / fScale;
1090 291 : CPLAssert(
1091 : fIntLastVal <=
1092 : static_cast<float>(std::numeric_limits<int>::max()));
1093 291 : int nLastVal = (int)fIntLastVal;
1094 291 : GByte nWordSize = 1;
1095 41301 : for (int l = 1; l < nReqXSize; l++)
1096 : {
1097 41010 : float fVal =
1098 : ((float *)
1099 41010 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize +
1100 41010 : l];
1101 41010 : float fIntVal = (fVal - fOffset) / fScale;
1102 41010 : CPLAssert(fIntVal <=
1103 : static_cast<float>(
1104 : std::numeric_limits<int>::max()));
1105 41010 : const int nVal = (int)fIntVal;
1106 41010 : const int nDiff = nVal - nLastVal;
1107 41010 : CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
1108 41010 : if (nDiff < -32768 || nDiff > 32767)
1109 : {
1110 0 : nWordSize = 4;
1111 0 : break;
1112 : }
1113 41010 : if (nDiff < -128 || nDiff > 127)
1114 201 : nWordSize = 2;
1115 41010 : nLastVal = nVal;
1116 : }
1117 :
1118 291 : VSIFWriteL(&nWordSize, 1, 1, fp);
1119 291 : fLastVal =
1120 : ((float *)
1121 291 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
1122 291 : fIntLastVal = (fLastVal - fOffset) / fScale;
1123 291 : nLastVal = (int)fIntLastVal;
1124 291 : WriteInt(fp, nLastVal);
1125 41301 : for (int l = 1; l < nReqXSize; l++)
1126 : {
1127 41010 : float fVal =
1128 : ((float *)
1129 41010 : pTileBuffer)[(nReqYSize - k - 1) * nReqXSize +
1130 41010 : l];
1131 41010 : float fIntVal = (fVal - fOffset) / fScale;
1132 41010 : int nVal = (int)fIntVal;
1133 41010 : int nDiff = nVal - nLastVal;
1134 41010 : CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
1135 41010 : if (nWordSize == 1)
1136 : {
1137 5010 : CPLAssert(nDiff >= -128 && nDiff <= 127);
1138 5010 : signed char chDiff = (signed char)nDiff;
1139 5010 : VSIFWriteL(&chDiff, 1, 1, fp);
1140 : }
1141 36000 : else if (nWordSize == 2)
1142 : {
1143 36000 : CPLAssert(nDiff >= -32768 && nDiff <= 32767);
1144 36000 : WriteShort(fp, (short)nDiff);
1145 : }
1146 : else
1147 : {
1148 0 : WriteInt(fp, nDiff);
1149 : }
1150 41010 : nLastVal = nVal;
1151 : }
1152 : }
1153 : }
1154 :
1155 36 : if (pfnProgress && !pfnProgress((j * nXBlocks + i + 1) * 1.0 /
1156 18 : (nXBlocks * nYBlocks),
1157 : nullptr, pProgressData))
1158 : {
1159 0 : eErr = CE_Failure;
1160 0 : break;
1161 : }
1162 : }
1163 : }
1164 :
1165 15 : CPLFree(pTileBuffer);
1166 :
1167 15 : VSIFCloseL(fp);
1168 :
1169 15 : if (eErr != CE_None)
1170 0 : return nullptr;
1171 :
1172 15 : GDALOpenInfo oOpenInfo(osFilename.c_str(), GA_ReadOnly);
1173 15 : HF2Dataset *poDS = reinterpret_cast<HF2Dataset *>(Open(&oOpenInfo));
1174 :
1175 15 : if (poDS)
1176 15 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1177 :
1178 15 : return poDS;
1179 : }
1180 :
1181 : /************************************************************************/
1182 : /* GDALRegister_HF2() */
1183 : /************************************************************************/
1184 :
1185 1595 : void GDALRegister_HF2()
1186 :
1187 : {
1188 1595 : if (GDALGetDriverByName("HF2") != nullptr)
1189 302 : return;
1190 :
1191 1293 : GDALDriver *poDriver = new GDALDriver();
1192 :
1193 1293 : poDriver->SetDescription("HF2");
1194 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1195 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "HF2/HFZ heightfield raster");
1196 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/hf2.html");
1197 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hf2");
1198 1293 : poDriver->SetMetadataItem(
1199 : GDAL_DMD_CREATIONOPTIONLIST,
1200 : "<CreationOptionList>"
1201 : " <Option name='VERTICAL_PRECISION' type='float' default='0.01' "
1202 : "description='Vertical precision.'/>"
1203 : " <Option name='COMPRESS' type='boolean' default='false' "
1204 : "description='Set to true to produce a GZip compressed file.'/>"
1205 : " <Option name='BLOCKSIZE' type='int' default='256' "
1206 : "description='Tile size.'/>"
1207 1293 : "</CreationOptionList>");
1208 :
1209 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1210 :
1211 1293 : poDriver->pfnOpen = HF2Dataset::Open;
1212 1293 : poDriver->pfnIdentify = HF2Dataset::Identify;
1213 1293 : poDriver->pfnCreateCopy = HF2Dataset::CreateCopy;
1214 :
1215 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1216 : }
|
