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