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