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