Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GView
4 : * Purpose: Implementation of Atlantis HKV labelled blob support
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Frank Warmerdam
9 : * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include <ctype.h>
31 :
32 : #include "atlsci_spheroid.h"
33 : #include "cpl_string.h"
34 : #include "gdal_frmts.h"
35 : #include "ogr_spatialref.h"
36 : #include "rawdataset.h"
37 :
38 : #include <cmath>
39 :
40 : #include <algorithm>
41 :
42 : /************************************************************************/
43 : /* ==================================================================== */
44 : /* HKVRasterBand */
45 : /* ==================================================================== */
46 : /************************************************************************/
47 :
48 : class HKVDataset;
49 :
50 : class HKVRasterBand final : public RawRasterBand
51 : {
52 : friend class HKVDataset;
53 :
54 : public:
55 : HKVRasterBand(HKVDataset *poDS, int nBand, VSILFILE *fpRaw,
56 : unsigned int nImgOffset, int nPixelOffset, int nLineOffset,
57 : GDALDataType eDataType, int bNativeOrder);
58 :
59 188 : ~HKVRasterBand() override
60 94 : {
61 188 : }
62 :
63 : CPLErr SetNoDataValue(double) override;
64 : };
65 :
66 : /************************************************************************/
67 : /* HKV Spheroids */
68 : /************************************************************************/
69 :
70 : class HKVSpheroidList : public SpheroidList
71 : {
72 : public:
73 : HKVSpheroidList();
74 :
75 44 : ~HKVSpheroidList()
76 44 : {
77 44 : }
78 : };
79 :
80 44 : HKVSpheroidList ::HKVSpheroidList()
81 : {
82 44 : num_spheroids = 58;
83 44 : epsilonR = 0.1;
84 44 : epsilonI = 0.000001;
85 :
86 44 : spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830", 6377563.396,
87 : 299.3249646);
88 44 : spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",
89 : 6377340.189, 299.3249646);
90 44 : spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",
91 : 6378160, 298.25);
92 44 : spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",
93 : 6377483.865, 299.1528128);
94 44 : spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841", 6377397.155,
95 : 299.1528128);
96 44 : spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858", 6378294.0,
97 : 294.297);
98 44 : spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866", 6378206.4,
99 : 294.9786982);
100 44 : spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880", 6378249.145,
101 : 293.465);
102 44 : spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",
103 : 6377276.345, 300.8017);
104 44 : spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",
105 : 6377298.556, 300.8017);
106 44 : spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",
107 : 6377301.243, 300.8017);
108 44 : spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",
109 : 6377295.664, 300.8017);
110 44 : spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",
111 : 6377304.063, 300.8017);
112 44 : spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",
113 : 6377309.613, 300.8017);
114 44 : spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",
115 : 6378155, 298.3);
116 44 : spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906", 6378200,
117 : 298.3);
118 44 : spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960", 6378270,
119 : 297);
120 44 : spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes", 6378273.0,
121 : 298.279);
122 44 : spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",
123 : 6378160, 298.247);
124 44 : spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",
125 : 6378388, 297);
126 44 : spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67", 6378160.0,
127 : 298.254);
128 44 : spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75", 6378140.0,
129 : 298.25298);
130 44 : spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",
131 : 6378245, 298.3);
132 44 : spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula", 6378165.0,
133 : 292.308);
134 44 : spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80", 6378137,
135 : 298.257222101);
136 44 : spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",
137 : 6378160, 298.25);
138 44 : spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72", 6378135,
139 : 298.26);
140 44 : spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84", 6378137,
141 : 298.257223563);
142 44 : spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84", 6378137.0,
143 : 298.252841);
144 44 : spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel", 6377397.0,
145 : 299.1976073);
146 :
147 44 : spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830", 6377563.396,
148 : 299.3249646);
149 44 : spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",
150 : 6377340.189, 299.3249646);
151 44 : spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",
152 : 6378160, 298.25);
153 44 : spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",
154 : 6377483.865, 299.1528128);
155 44 : spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",
156 : 6377397.155, 299.1528128);
157 44 : spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858", 6378294.0,
158 : 294.297);
159 44 : spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866", 6378206.4,
160 : 294.9786982);
161 44 : spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",
162 : 6378249.145, 293.465);
163 44 : spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",
164 : 6377276.345, 300.8017);
165 44 : spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",
166 : 6377298.556, 300.8017);
167 44 : spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",
168 : 6377301.243, 300.8017);
169 44 : spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",
170 : 6377295.664, 300.8017);
171 44 : spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",
172 : 6377304.063, 300.8017);
173 44 : spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",
174 : 6377309.613, 300.8017);
175 44 : spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",
176 : 6378155, 298.3);
177 44 : spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906", 6378200,
178 : 298.3);
179 44 : spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960", 6378270,
180 : 297);
181 44 : spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",
182 : 6378160, 298.247);
183 44 : spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",
184 : 6378388, 297);
185 44 : spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67", 6378160.0,
186 : 298.254);
187 44 : spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75", 6378140.0,
188 : 298.25298);
189 44 : spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",
190 : 6378245, 298.3);
191 44 : spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80", 6378137,
192 : 298.257222101);
193 44 : spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",
194 : 6378160, 298.25);
195 44 : spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72", 6378135,
196 : 298.26);
197 44 : spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84", 6378137,
198 : 298.257223563);
199 44 : spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84", 6378137.0,
200 : 298.252841);
201 44 : spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel", 6377397.0,
202 : 299.1976073);
203 44 : }
204 :
205 : CPLErr SaveHKVAttribFile(const char *pszFilenameIn, int nXSize, int nYSize,
206 : int nBands, GDALDataType eType, int bNoDataSet,
207 : double dfNoDataValue);
208 :
209 : /************************************************************************/
210 : /* ==================================================================== */
211 : /* HKVDataset */
212 : /* ==================================================================== */
213 : /************************************************************************/
214 :
215 : class HKVDataset final : public RawDataset
216 : {
217 : friend class HKVRasterBand;
218 :
219 : char *pszPath;
220 : VSILFILE *fpBlob;
221 :
222 : int nGCPCount;
223 : GDAL_GCP *pasGCPList;
224 :
225 : void ProcessGeoref(const char *);
226 : void ProcessGeorefGCP(char **, const char *, double, double);
227 :
228 44 : void SetVersion(float version_number)
229 : {
230 : // Update stored info.
231 44 : MFF2version = version_number;
232 44 : }
233 :
234 : float MFF2version;
235 :
236 : GDALDataType eRasterType;
237 :
238 : void SetNoDataValue(double);
239 :
240 : OGRSpatialReference m_oSRS{};
241 : OGRSpatialReference m_oGCPSRS{};
242 : double adfGeoTransform[6];
243 :
244 : char **papszAttrib;
245 :
246 : bool bGeorefChanged;
247 : char **papszGeoref;
248 :
249 : // NOTE: The MFF2 format goes against GDAL's API in that nodata values are
250 : // set per-dataset rather than per-band. To compromise, for writing out,
251 : // the dataset's nodata value will be set to the last value set on any of
252 : // the raster bands.
253 :
254 : bool bNoDataSet;
255 : bool bNoDataChanged;
256 : double dfNoDataValue;
257 :
258 : CPL_DISALLOW_COPY_ASSIGN(HKVDataset)
259 :
260 : CPLErr Close() override;
261 :
262 : public:
263 : HKVDataset();
264 : ~HKVDataset() override;
265 :
266 2 : int GetGCPCount() override /* const */
267 : {
268 2 : return nGCPCount;
269 : }
270 :
271 0 : const OGRSpatialReference *GetGCPSpatialRef() const override
272 : {
273 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
274 : }
275 :
276 : const GDAL_GCP *GetGCPs() override;
277 :
278 15 : const OGRSpatialReference *GetSpatialRef() const override
279 : {
280 15 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
281 : }
282 :
283 : CPLErr GetGeoTransform(double *) override;
284 :
285 : CPLErr SetGeoTransform(double *) override;
286 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
287 :
288 : static GDALDataset *Open(GDALOpenInfo *);
289 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
290 : int nBands, GDALDataType eType,
291 : char **papszParamList);
292 : static GDALDataset *CreateCopy(const char *pszFilename,
293 : GDALDataset *poSrcDS, int bStrict,
294 : char **papszOptions,
295 : GDALProgressFunc pfnProgress,
296 : void *pProgressData);
297 :
298 : static CPLErr Delete(const char *pszName);
299 : };
300 :
301 : /************************************************************************/
302 : /* ==================================================================== */
303 : /* HKVRasterBand */
304 : /* ==================================================================== */
305 : /************************************************************************/
306 :
307 : /************************************************************************/
308 : /* HKVRasterBand() */
309 : /************************************************************************/
310 :
311 94 : HKVRasterBand::HKVRasterBand(HKVDataset *poDSIn, int nBandIn, VSILFILE *fpRawIn,
312 : unsigned int nImgOffsetIn, int nPixelOffsetIn,
313 : int nLineOffsetIn, GDALDataType eDataTypeIn,
314 94 : int bNativeOrderIn)
315 : : RawRasterBand(GDALDataset::FromHandle(poDSIn), nBandIn, fpRawIn,
316 : nImgOffsetIn, nPixelOffsetIn, nLineOffsetIn, eDataTypeIn,
317 94 : bNativeOrderIn, RawRasterBand::OwnFP::NO)
318 :
319 : {
320 94 : poDS = poDSIn;
321 94 : nBand = nBandIn;
322 :
323 94 : nBlockXSize = poDS->GetRasterXSize();
324 94 : nBlockYSize = 1;
325 94 : }
326 :
327 : /************************************************************************/
328 : /* SetNoDataValue() */
329 : /************************************************************************/
330 :
331 0 : CPLErr HKVRasterBand::SetNoDataValue(double dfNewValue)
332 :
333 : {
334 0 : HKVDataset *poHKVDS = reinterpret_cast<HKVDataset *>(poDS);
335 0 : RawRasterBand::SetNoDataValue(dfNewValue);
336 0 : poHKVDS->SetNoDataValue(dfNewValue);
337 :
338 0 : return CE_None;
339 : }
340 :
341 : /************************************************************************/
342 : /* ==================================================================== */
343 : /* HKVDataset */
344 : /* ==================================================================== */
345 : /************************************************************************/
346 :
347 : /************************************************************************/
348 : /* HKVDataset() */
349 : /************************************************************************/
350 :
351 44 : HKVDataset::HKVDataset()
352 : : pszPath(nullptr), fpBlob(nullptr), nGCPCount(0), pasGCPList(nullptr),
353 : // Initialize datasets to new version; change if necessary.
354 : MFF2version(1.1f), eRasterType(GDT_Unknown), papszAttrib(nullptr),
355 : bGeorefChanged(false), papszGeoref(nullptr), bNoDataSet(false),
356 44 : bNoDataChanged(false), dfNoDataValue(0.0)
357 : {
358 44 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
359 44 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
360 44 : adfGeoTransform[0] = 0.0;
361 44 : adfGeoTransform[1] = 1.0;
362 44 : adfGeoTransform[2] = 0.0;
363 44 : adfGeoTransform[3] = 0.0;
364 44 : adfGeoTransform[4] = 0.0;
365 44 : adfGeoTransform[5] = 1.0;
366 44 : }
367 :
368 : /************************************************************************/
369 : /* ~HKVDataset() */
370 : /************************************************************************/
371 :
372 88 : HKVDataset::~HKVDataset()
373 :
374 : {
375 44 : HKVDataset::Close();
376 88 : }
377 :
378 : /************************************************************************/
379 : /* Close() */
380 : /************************************************************************/
381 :
382 88 : CPLErr HKVDataset::Close()
383 : {
384 88 : CPLErr eErr = CE_None;
385 88 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
386 : {
387 44 : if (HKVDataset::FlushCache(true) != CE_None)
388 0 : eErr = CE_Failure;
389 :
390 44 : if (bGeorefChanged)
391 : {
392 : const char *pszFilename =
393 26 : CPLFormFilename(pszPath, "georef", nullptr);
394 26 : CSLSave(papszGeoref, pszFilename);
395 : }
396 :
397 44 : if (bNoDataChanged)
398 : {
399 0 : SaveHKVAttribFile(pszPath, nRasterXSize, nRasterYSize, nBands,
400 0 : eRasterType, bNoDataSet, dfNoDataValue);
401 : }
402 :
403 44 : if (fpBlob)
404 : {
405 44 : if (VSIFCloseL(fpBlob) != 0)
406 : {
407 0 : eErr = CE_Failure;
408 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
409 : }
410 : }
411 :
412 44 : if (nGCPCount > 0)
413 : {
414 14 : GDALDeinitGCPs(nGCPCount, pasGCPList);
415 14 : CPLFree(pasGCPList);
416 : }
417 :
418 44 : CPLFree(pszPath);
419 44 : CSLDestroy(papszGeoref);
420 44 : CSLDestroy(papszAttrib);
421 :
422 44 : if (GDALPamDataset::Close() != CE_None)
423 0 : eErr = CE_Failure;
424 : }
425 88 : return eErr;
426 : }
427 :
428 : /************************************************************************/
429 : /* SetNoDataValue() */
430 : /************************************************************************/
431 :
432 0 : void HKVDataset::SetNoDataValue(double dfNewValue)
433 :
434 : {
435 0 : bNoDataSet = true;
436 0 : bNoDataChanged = true;
437 0 : dfNoDataValue = dfNewValue;
438 0 : }
439 :
440 : /************************************************************************/
441 : /* SaveHKVAttribFile() */
442 : /************************************************************************/
443 :
444 26 : CPLErr SaveHKVAttribFile(const char *pszFilenameIn, int nXSize, int nYSize,
445 : int nBands, GDALDataType eType, int bNoDataSet,
446 : double dfNoDataValue)
447 :
448 : {
449 26 : const char *pszFilename = CPLFormFilename(pszFilenameIn, "attrib", nullptr);
450 :
451 26 : FILE *fp = VSIFOpen(pszFilename, "wt");
452 26 : if (fp == nullptr)
453 : {
454 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.",
455 : pszFilename);
456 0 : return CE_Failure;
457 : }
458 :
459 26 : fprintf(fp, "channel.enumeration = %d\n", nBands);
460 26 : fprintf(fp, "channel.interleave = { *pixel tile sequential }\n");
461 26 : fprintf(fp, "extent.cols = %d\n", nXSize);
462 26 : fprintf(fp, "extent.rows = %d\n", nYSize);
463 :
464 26 : switch (eType)
465 : {
466 11 : case GDT_Byte:
467 11 : fprintf(fp, "pixel.encoding = "
468 : "{ *unsigned twos-complement ieee-754 }\n");
469 11 : break;
470 :
471 3 : case GDT_UInt16:
472 3 : fprintf(fp, "pixel.encoding = "
473 : "{ *unsigned twos-complement ieee-754 }\n");
474 3 : break;
475 :
476 6 : case GDT_CInt16:
477 : case GDT_Int16:
478 6 : fprintf(fp, "pixel.encoding = "
479 : "{ unsigned *twos-complement ieee-754 }\n");
480 6 : break;
481 :
482 6 : case GDT_CFloat32:
483 : case GDT_Float32:
484 6 : fprintf(fp, "pixel.encoding = "
485 : "{ unsigned twos-complement *ieee-754 }\n");
486 6 : break;
487 :
488 0 : default:
489 0 : CPLAssert(false);
490 : }
491 :
492 26 : fprintf(fp, "pixel.size = %d\n", GDALGetDataTypeSizeBits(eType));
493 26 : if (GDALDataTypeIsComplex(eType))
494 6 : fprintf(fp, "pixel.field = { real *complex }\n");
495 : else
496 20 : fprintf(fp, "pixel.field = { *real complex }\n");
497 :
498 : #ifdef CPL_MSB
499 : fprintf(fp, "pixel.order = { lsbf *msbf }\n");
500 : #else
501 26 : fprintf(fp, "pixel.order = { *lsbf msbf }\n");
502 : #endif
503 :
504 26 : if (bNoDataSet)
505 0 : fprintf(fp, "pixel.no_data = %s\n", CPLSPrintf("%f", dfNoDataValue));
506 :
507 : // Version information- only create the new style.
508 26 : fprintf(fp, "version = 1.1");
509 :
510 26 : if (VSIFClose(fp) != 0)
511 0 : return CE_Failure;
512 26 : return CE_None;
513 : }
514 :
515 : /************************************************************************/
516 : /* GetGeoTransform() */
517 : /************************************************************************/
518 :
519 29 : CPLErr HKVDataset::GetGeoTransform(double *padfTransform)
520 :
521 : {
522 29 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
523 29 : return CE_None;
524 : }
525 :
526 : /************************************************************************/
527 : /* SetGeoTransform() */
528 : /************************************************************************/
529 :
530 26 : CPLErr HKVDataset::SetGeoTransform(double *padfTransform)
531 :
532 : {
533 : // NOTE: Geotransform coordinates must match the current projection
534 : // of the dataset being changed (not the geotransform source).
535 : // i.e. be in lat/longs for LL projected; UTM for UTM projected.
536 : // SET PROJECTION BEFORE SETTING GEOTRANSFORM TO AVOID SYNCHRONIZATION
537 : // PROBLEMS.
538 :
539 : // Update the geotransform itself.
540 26 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
541 :
542 : // Clear previous gcps.
543 26 : if (nGCPCount > 0)
544 : {
545 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
546 0 : CPLFree(pasGCPList);
547 : }
548 26 : nGCPCount = 0;
549 26 : pasGCPList = nullptr;
550 :
551 : // Return if the identity transform is set.
552 26 : if (adfGeoTransform[0] == 0.0 && adfGeoTransform[1] == 1.0 &&
553 0 : adfGeoTransform[2] == 0.0 && adfGeoTransform[3] == 0.0 &&
554 0 : adfGeoTransform[4] == 0.0 && adfGeoTransform[5] == 1.0)
555 0 : return CE_None;
556 :
557 : // Update georef text info for saving later, and update GCPs to match
558 : // geotransform.
559 :
560 26 : OGRCoordinateTransformation *poTransform = nullptr;
561 26 : bool bSuccess = true;
562 :
563 : // Projection parameter checking will have been done in SetProjection.
564 37 : if ((CSLFetchNameValue(papszGeoref, "projection.name") != nullptr) &&
565 11 : (EQUAL(CSLFetchNameValue(papszGeoref, "projection.name"), "UTM")))
566 : {
567 1 : auto poLLSRS = m_oSRS.CloneGeogCS();
568 1 : if (poLLSRS)
569 : {
570 1 : poLLSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
571 1 : poTransform = OGRCreateCoordinateTransformation(&m_oSRS, poLLSRS);
572 1 : delete poLLSRS;
573 1 : if (poTransform == nullptr)
574 : {
575 0 : bSuccess = false;
576 0 : CPLErrorReset();
577 : }
578 : }
579 : else
580 : {
581 0 : bSuccess = false;
582 : }
583 : }
584 25 : else if (((CSLFetchNameValue(papszGeoref, "projection.name") != nullptr) &&
585 10 : (!EQUAL(CSLFetchNameValue(papszGeoref, "projection.name"),
586 50 : "LL"))) ||
587 25 : (CSLFetchNameValue(papszGeoref, "projection.name") == nullptr))
588 : {
589 15 : return CE_Failure;
590 : }
591 :
592 11 : nGCPCount = 0;
593 11 : pasGCPList = static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), 5));
594 :
595 : /* -------------------------------------------------------------------- */
596 : /* top left */
597 : /* -------------------------------------------------------------------- */
598 11 : GDALInitGCPs(1, pasGCPList + nGCPCount);
599 11 : CPLFree(pasGCPList[nGCPCount].pszId);
600 11 : pasGCPList[nGCPCount].pszId = CPLStrdup("top_left");
601 :
602 11 : double temp_lat = 0.0;
603 11 : double temp_long = 0.0;
604 11 : if (MFF2version > 1.0)
605 : {
606 11 : temp_lat = padfTransform[3];
607 11 : temp_long = padfTransform[0];
608 11 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
609 11 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
610 : }
611 : else
612 : {
613 0 : temp_lat =
614 0 : padfTransform[3] + 0.5 * padfTransform[4] + 0.5 * padfTransform[5];
615 0 : temp_long =
616 0 : padfTransform[0] + 0.5 * padfTransform[1] + 0.5 * padfTransform[2];
617 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
618 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
619 : }
620 11 : pasGCPList[nGCPCount].dfGCPX = temp_long;
621 11 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
622 11 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
623 11 : nGCPCount++;
624 :
625 11 : if (poTransform != nullptr)
626 : {
627 1 : if (!bSuccess || !poTransform->Transform(1, &temp_long, &temp_lat))
628 0 : bSuccess = false;
629 : }
630 :
631 11 : if (bSuccess)
632 : {
633 11 : char szValue[128] = {'\0'};
634 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_lat);
635 11 : papszGeoref =
636 11 : CSLSetNameValue(papszGeoref, "top_left.latitude", szValue);
637 :
638 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_long);
639 11 : papszGeoref =
640 11 : CSLSetNameValue(papszGeoref, "top_left.longitude", szValue);
641 : }
642 :
643 : /* -------------------------------------------------------------------- */
644 : /* top_right */
645 : /* -------------------------------------------------------------------- */
646 11 : GDALInitGCPs(1, pasGCPList + nGCPCount);
647 11 : CPLFree(pasGCPList[nGCPCount].pszId);
648 11 : pasGCPList[nGCPCount].pszId = CPLStrdup("top_right");
649 :
650 11 : if (MFF2version > 1.0)
651 : {
652 11 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4];
653 11 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1];
654 11 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
655 11 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
656 : }
657 : else
658 : {
659 0 : temp_lat = padfTransform[3] +
660 0 : (GetRasterXSize() - 0.5) * padfTransform[4] +
661 0 : 0.5 * padfTransform[5];
662 0 : temp_long = padfTransform[0] +
663 0 : (GetRasterXSize() - 0.5) * padfTransform[1] +
664 0 : 0.5 * padfTransform[2];
665 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize() - 0.5;
666 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
667 : }
668 11 : pasGCPList[nGCPCount].dfGCPX = temp_long;
669 11 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
670 11 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
671 11 : nGCPCount++;
672 :
673 11 : if (poTransform != nullptr)
674 : {
675 1 : if (!bSuccess || !poTransform->Transform(1, &temp_long, &temp_lat))
676 0 : bSuccess = false;
677 : }
678 :
679 11 : if (bSuccess)
680 : {
681 11 : char szValue[128] = {'\0'};
682 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_lat);
683 11 : papszGeoref =
684 11 : CSLSetNameValue(papszGeoref, "top_right.latitude", szValue);
685 :
686 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_long);
687 11 : papszGeoref =
688 11 : CSLSetNameValue(papszGeoref, "top_right.longitude", szValue);
689 : }
690 :
691 : /* -------------------------------------------------------------------- */
692 : /* bottom_left */
693 : /* -------------------------------------------------------------------- */
694 11 : GDALInitGCPs(1, pasGCPList + nGCPCount);
695 11 : CPLFree(pasGCPList[nGCPCount].pszId);
696 11 : pasGCPList[nGCPCount].pszId = CPLStrdup("bottom_left");
697 :
698 11 : if (MFF2version > 1.0)
699 : {
700 11 : temp_lat = padfTransform[3] + GetRasterYSize() * padfTransform[5];
701 11 : temp_long = padfTransform[0] + GetRasterYSize() * padfTransform[2];
702 11 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
703 11 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
704 : }
705 : else
706 : {
707 0 : temp_lat = padfTransform[3] + 0.5 * padfTransform[4] +
708 0 : (GetRasterYSize() - 0.5) * padfTransform[5];
709 0 : temp_long = padfTransform[0] + 0.5 * padfTransform[1] +
710 0 : (GetRasterYSize() - 0.5) * padfTransform[2];
711 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
712 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize() - 0.5;
713 : }
714 11 : pasGCPList[nGCPCount].dfGCPX = temp_long;
715 11 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
716 11 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
717 11 : nGCPCount++;
718 :
719 11 : if (poTransform != nullptr)
720 : {
721 1 : if (!bSuccess || !poTransform->Transform(1, &temp_long, &temp_lat))
722 0 : bSuccess = false;
723 : }
724 :
725 11 : if (bSuccess)
726 : {
727 11 : char szValue[128] = {'\0'};
728 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_lat);
729 11 : papszGeoref =
730 11 : CSLSetNameValue(papszGeoref, "bottom_left.latitude", szValue);
731 :
732 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_long);
733 11 : papszGeoref =
734 11 : CSLSetNameValue(papszGeoref, "bottom_left.longitude", szValue);
735 : }
736 :
737 : /* -------------------------------------------------------------------- */
738 : /* bottom_right */
739 : /* -------------------------------------------------------------------- */
740 11 : GDALInitGCPs(1, pasGCPList + nGCPCount);
741 11 : CPLFree(pasGCPList[nGCPCount].pszId);
742 11 : pasGCPList[nGCPCount].pszId = CPLStrdup("bottom_right");
743 :
744 11 : if (MFF2version > 1.0)
745 : {
746 11 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] +
747 11 : GetRasterYSize() * padfTransform[5];
748 11 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] +
749 11 : GetRasterYSize() * padfTransform[2];
750 11 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
751 11 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
752 : }
753 : else
754 : {
755 0 : temp_lat = padfTransform[3] +
756 0 : (GetRasterXSize() - 0.5) * padfTransform[4] +
757 0 : (GetRasterYSize() - 0.5) * padfTransform[5];
758 0 : temp_long = padfTransform[0] +
759 0 : (GetRasterXSize() - 0.5) * padfTransform[1] +
760 0 : (GetRasterYSize() - 0.5) * padfTransform[2];
761 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize() - 0.5;
762 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize() - 0.5;
763 : }
764 11 : pasGCPList[nGCPCount].dfGCPX = temp_long;
765 11 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
766 11 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
767 11 : nGCPCount++;
768 :
769 11 : if (poTransform != nullptr)
770 : {
771 1 : if (!bSuccess || !poTransform->Transform(1, &temp_long, &temp_lat))
772 0 : bSuccess = false;
773 : }
774 :
775 11 : if (bSuccess)
776 : {
777 11 : char szValue[128] = {'\0'};
778 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_lat);
779 11 : papszGeoref =
780 11 : CSLSetNameValue(papszGeoref, "bottom_right.latitude", szValue);
781 :
782 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_long);
783 11 : papszGeoref =
784 11 : CSLSetNameValue(papszGeoref, "bottom_right.longitude", szValue);
785 : }
786 :
787 : /* -------------------------------------------------------------------- */
788 : /* Center */
789 : /* -------------------------------------------------------------------- */
790 11 : GDALInitGCPs(1, pasGCPList + nGCPCount);
791 11 : CPLFree(pasGCPList[nGCPCount].pszId);
792 11 : pasGCPList[nGCPCount].pszId = CPLStrdup("centre");
793 :
794 11 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
795 11 : GetRasterYSize() * padfTransform[5] * 0.5;
796 11 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
797 11 : GetRasterYSize() * padfTransform[2] * 0.5;
798 11 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize() / 2.0;
799 11 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize() / 2.0;
800 :
801 11 : pasGCPList[nGCPCount].dfGCPX = temp_long;
802 11 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
803 11 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
804 11 : nGCPCount++;
805 :
806 11 : if (poTransform != nullptr)
807 : {
808 1 : if (!bSuccess || !poTransform->Transform(1, &temp_long, &temp_lat))
809 0 : bSuccess = false;
810 : }
811 :
812 11 : if (bSuccess)
813 : {
814 11 : char szValue[128] = {'\0'};
815 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_lat);
816 11 : papszGeoref = CSLSetNameValue(papszGeoref, "centre.latitude", szValue);
817 :
818 11 : CPLsnprintf(szValue, sizeof(szValue), "%.10f", temp_long);
819 11 : papszGeoref = CSLSetNameValue(papszGeoref, "centre.longitude", szValue);
820 : }
821 :
822 11 : if (!bSuccess)
823 : {
824 0 : CPLError(CE_Warning, CPLE_AppDefined,
825 : "Error setting header info in SetGeoTransform. "
826 : "Changes may not be saved properly.");
827 : }
828 :
829 11 : if (poTransform != nullptr)
830 1 : delete poTransform;
831 :
832 11 : bGeorefChanged = true;
833 :
834 11 : return CE_None;
835 : }
836 :
837 : /************************************************************************/
838 : /* SetSpatialRef() */
839 : /* */
840 : /* We provide very limited support for setting the projection. */
841 : /************************************************************************/
842 :
843 26 : CPLErr HKVDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
844 :
845 : {
846 : // Update a georef file.
847 26 : if (poSRS == nullptr)
848 : {
849 0 : m_oSRS.Clear();
850 0 : return CE_None;
851 : }
852 26 : m_oSRS = *poSRS;
853 :
854 27 : if ((m_oSRS.GetAttrValue("PROJECTION") != nullptr) &&
855 1 : (EQUAL(m_oSRS.GetAttrValue("PROJECTION"), SRS_PT_TRANSVERSE_MERCATOR)))
856 : {
857 1 : papszGeoref = CSLSetNameValue(papszGeoref, "projection.name", "utm");
858 1 : OGRErr ogrerrorOl = OGRERR_NONE;
859 1 : papszGeoref = CSLSetNameValue(
860 : papszGeoref, "projection.origin_longitude",
861 : CPLSPrintf("%f", m_oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0,
862 : &ogrerrorOl)));
863 : }
864 50 : else if (m_oSRS.GetAttrValue("PROJECTION") == nullptr &&
865 25 : m_oSRS.IsGeographic())
866 : {
867 25 : papszGeoref = CSLSetNameValue(papszGeoref, "projection.name", "LL");
868 : }
869 : else
870 : {
871 0 : CPLError(CE_Warning, CPLE_AppDefined, "Unrecognized projection.");
872 0 : return CE_Failure;
873 : }
874 :
875 26 : OGRErr ogrerrorEq = OGRERR_NONE;
876 26 : const double eq_radius = m_oSRS.GetSemiMajor(&ogrerrorEq);
877 :
878 26 : OGRErr ogrerrorInvf = OGRERR_NONE;
879 26 : const double inv_flattening = m_oSRS.GetInvFlattening(&ogrerrorInvf);
880 :
881 26 : if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE))
882 : {
883 26 : HKVSpheroidList *hkvEllipsoids = new HKVSpheroidList;
884 : char *spheroid_name =
885 26 : hkvEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(
886 : eq_radius, inv_flattening);
887 26 : if (spheroid_name != nullptr)
888 : {
889 26 : papszGeoref =
890 26 : CSLSetNameValue(papszGeoref, "spheroid.name", spheroid_name);
891 : }
892 26 : CPLFree(spheroid_name);
893 26 : delete hkvEllipsoids;
894 : }
895 : else
896 : {
897 : // Default to previous behavior if spheroid not found by radius and
898 : // inverse flattening.
899 0 : char *pszProjection = nullptr;
900 0 : m_oSRS.exportToWkt(&pszProjection);
901 0 : if (pszProjection && strstr(pszProjection, "Bessel") != nullptr)
902 : {
903 0 : papszGeoref =
904 0 : CSLSetNameValue(papszGeoref, "spheroid.name", "ev-bessel");
905 : }
906 : else
907 : {
908 0 : papszGeoref =
909 0 : CSLSetNameValue(papszGeoref, "spheroid.name", "ev-wgs-84");
910 : }
911 0 : CPLFree(pszProjection);
912 : }
913 26 : bGeorefChanged = true;
914 26 : return CE_None;
915 : }
916 :
917 : /************************************************************************/
918 : /* GetGCP() */
919 : /************************************************************************/
920 :
921 0 : const GDAL_GCP *HKVDataset::GetGCPs()
922 :
923 : {
924 0 : return pasGCPList;
925 : }
926 :
927 : /************************************************************************/
928 : /* ProcessGeorefGCP() */
929 : /************************************************************************/
930 :
931 90 : void HKVDataset::ProcessGeorefGCP(char **papszGeorefIn, const char *pszBase,
932 : double dfRasterX, double dfRasterY)
933 :
934 : {
935 : /* -------------------------------------------------------------------- */
936 : /* Fetch the GCP from the string list. */
937 : /* -------------------------------------------------------------------- */
938 90 : char szFieldName[128] = {'\0'};
939 90 : snprintf(szFieldName, sizeof(szFieldName), "%s.latitude", pszBase);
940 90 : double dfLat = 0.0;
941 90 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
942 75 : return;
943 : else
944 15 : dfLat = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
945 :
946 15 : snprintf(szFieldName, sizeof(szFieldName), "%s.longitude", pszBase);
947 15 : double dfLong = 0.0;
948 15 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
949 0 : return;
950 : else
951 15 : dfLong = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
952 :
953 : /* -------------------------------------------------------------------- */
954 : /* Add the gcp to the internal list. */
955 : /* -------------------------------------------------------------------- */
956 15 : GDALInitGCPs(1, pasGCPList + nGCPCount);
957 :
958 15 : CPLFree(pasGCPList[nGCPCount].pszId);
959 :
960 15 : pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
961 :
962 15 : pasGCPList[nGCPCount].dfGCPX = dfLong;
963 15 : pasGCPList[nGCPCount].dfGCPY = dfLat;
964 15 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
965 :
966 15 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
967 15 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
968 :
969 15 : nGCPCount++;
970 : }
971 :
972 : /************************************************************************/
973 : /* ProcessGeoref() */
974 : /************************************************************************/
975 :
976 18 : void HKVDataset::ProcessGeoref(const char *pszFilename)
977 :
978 : {
979 : /* -------------------------------------------------------------------- */
980 : /* Load the georef file, and boil white space away from around */
981 : /* the equal sign. */
982 : /* -------------------------------------------------------------------- */
983 18 : CSLDestroy(papszGeoref);
984 18 : papszGeoref = CSLLoad(pszFilename);
985 18 : if (papszGeoref == nullptr)
986 0 : return;
987 :
988 18 : HKVSpheroidList *hkvEllipsoids = new HKVSpheroidList;
989 :
990 87 : for (int i = 0; papszGeoref[i] != nullptr; i++)
991 : {
992 69 : int iDst = 0;
993 69 : char *pszLine = papszGeoref[i];
994 :
995 1899 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
996 : {
997 1830 : if (pszLine[iSrc] != ' ')
998 : {
999 1830 : pszLine[iDst++] = pszLine[iSrc];
1000 : }
1001 : }
1002 69 : pszLine[iDst] = '\0';
1003 : }
1004 :
1005 : /* -------------------------------------------------------------------- */
1006 : /* Try to get GCPs, in lat/longs . */
1007 : /* -------------------------------------------------------------------- */
1008 18 : nGCPCount = 0;
1009 18 : pasGCPList = reinterpret_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), 5));
1010 :
1011 18 : if (MFF2version > 1.0)
1012 : {
1013 18 : ProcessGeorefGCP(papszGeoref, "top_left", 0, 0);
1014 18 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize(), 0);
1015 18 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0, GetRasterYSize());
1016 18 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize(),
1017 18 : GetRasterYSize());
1018 18 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
1019 18 : GetRasterYSize() / 2.0);
1020 : }
1021 : else
1022 : {
1023 0 : ProcessGeorefGCP(papszGeoref, "top_left", 0.5, 0.5);
1024 0 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize() - 0.5, 0.5);
1025 0 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0.5,
1026 0 : GetRasterYSize() - 0.5);
1027 0 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize() - 0.5,
1028 0 : GetRasterYSize() - 0.5);
1029 0 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
1030 0 : GetRasterYSize() / 2.0);
1031 : }
1032 :
1033 18 : if (nGCPCount == 0)
1034 : {
1035 15 : CPLFree(pasGCPList);
1036 15 : pasGCPList = nullptr;
1037 : }
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Do we have a recognised projection? */
1041 : /* -------------------------------------------------------------------- */
1042 18 : const char *pszProjName = CSLFetchNameValue(papszGeoref, "projection.name");
1043 : const char *pszOriginLong =
1044 18 : CSLFetchNameValue(papszGeoref, "projection.origin_longitude");
1045 : const char *pszSpheroidName =
1046 18 : CSLFetchNameValue(papszGeoref, "spheroid.name");
1047 :
1048 36 : if (pszSpheroidName != nullptr &&
1049 18 : hkvEllipsoids->SpheroidInList(pszSpheroidName))
1050 : {
1051 : #if 0
1052 : // TODO(schwehr): Enable in trunk after 2.1 branch and fix.
1053 : // Breaks tests on some platforms.
1054 : CPLError( CE_Failure, CPLE_AppDefined,
1055 : "Unrecognized ellipsoid. Not handled. "
1056 : "Spheroid name not in spheroid list: '%s'",
1057 : pszSpheroidName );
1058 : #endif
1059 : // Why were eq_radius and inv_flattening never used?
1060 : // eq_radius = hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
1061 : // inv_flattening =
1062 : // hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
1063 : }
1064 0 : else if (pszProjName != nullptr)
1065 : {
1066 0 : CPLError(CE_Warning, CPLE_AppDefined,
1067 : "Unrecognized ellipsoid. Not handled.");
1068 : // TODO(schwehr): This error is was never what was happening.
1069 : // CPLError( CE_Warning, CPLE_AppDefined,
1070 : // "Unrecognized ellipsoid. Using wgs-84 parameters.");
1071 : // eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84"); */
1072 : // inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
1073 : }
1074 :
1075 18 : if (pszProjName != nullptr && EQUAL(pszProjName, "utm") && nGCPCount == 5)
1076 : {
1077 : // int nZone = (int)((CPLAtof(pszOriginLong)+184.5) / 6.0);
1078 3 : int nZone = 31; // TODO(schwehr): Where does 31 come from?
1079 :
1080 3 : if (pszOriginLong == nullptr)
1081 : {
1082 : // If origin not specified, assume 0.0.
1083 0 : CPLError(
1084 : CE_Warning, CPLE_AppDefined,
1085 : "No projection origin longitude specified. Assuming 0.0.");
1086 : }
1087 : else
1088 : {
1089 3 : nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
1090 : }
1091 :
1092 6 : OGRSpatialReference oUTM;
1093 :
1094 3 : if (pasGCPList[4].dfGCPY < 0)
1095 0 : oUTM.SetUTM(nZone, 0);
1096 : else
1097 3 : oUTM.SetUTM(nZone, 1);
1098 :
1099 6 : OGRSpatialReference oLL;
1100 3 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1101 3 : if (pszOriginLong != nullptr)
1102 : {
1103 3 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
1104 3 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
1105 : }
1106 :
1107 3 : if ((pszSpheroidName == nullptr) ||
1108 3 : (EQUAL(pszSpheroidName, "wgs-84")) ||
1109 3 : (EQUAL(pszSpheroidName, "wgs_84")))
1110 : {
1111 0 : oUTM.SetWellKnownGeogCS("WGS84");
1112 0 : oLL.SetWellKnownGeogCS("WGS84");
1113 : }
1114 : else
1115 : {
1116 3 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1117 : {
1118 3 : oUTM.SetGeogCS(
1119 : "unknown", "unknown", pszSpheroidName,
1120 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1121 : hkvEllipsoids->GetSpheroidInverseFlattening(
1122 : pszSpheroidName));
1123 3 : oLL.SetGeogCS(
1124 : "unknown", "unknown", pszSpheroidName,
1125 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1126 : hkvEllipsoids->GetSpheroidInverseFlattening(
1127 : pszSpheroidName));
1128 : }
1129 : else
1130 : {
1131 0 : CPLError(CE_Warning, CPLE_AppDefined,
1132 : "Unrecognized ellipsoid. Using wgs-84 parameters.");
1133 0 : oUTM.SetWellKnownGeogCS("WGS84");
1134 0 : oLL.SetWellKnownGeogCS("WGS84");
1135 : }
1136 : }
1137 :
1138 : OGRCoordinateTransformation *poTransform =
1139 3 : OGRCreateCoordinateTransformation(&oLL, &oUTM);
1140 :
1141 3 : bool bSuccess = true;
1142 3 : if (poTransform == nullptr)
1143 : {
1144 0 : CPLErrorReset();
1145 0 : bSuccess = false;
1146 : }
1147 :
1148 3 : double dfUtmX[5] = {0.0};
1149 3 : double dfUtmY[5] = {0.0};
1150 :
1151 3 : if (poTransform != nullptr)
1152 : {
1153 18 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
1154 : {
1155 15 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
1156 15 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
1157 :
1158 15 : if (bSuccess && !poTransform->Transform(1, &(dfUtmX[gcp_index]),
1159 : &(dfUtmY[gcp_index])))
1160 0 : bSuccess = false;
1161 : }
1162 : }
1163 :
1164 3 : if (bSuccess)
1165 : {
1166 : // Update GCPS to proper projection.
1167 18 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
1168 : {
1169 15 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
1170 15 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
1171 : }
1172 :
1173 3 : m_oGCPSRS = oUTM;
1174 :
1175 3 : bool transform_ok = CPL_TO_BOOL(
1176 3 : GDALGCPsToGeoTransform(5, pasGCPList, adfGeoTransform, 0));
1177 :
1178 3 : if (!transform_ok)
1179 : {
1180 : // Transform may not be sufficient in all cases (slant range
1181 : // projection).
1182 0 : adfGeoTransform[0] = 0.0;
1183 0 : adfGeoTransform[1] = 1.0;
1184 0 : adfGeoTransform[2] = 0.0;
1185 0 : adfGeoTransform[3] = 0.0;
1186 0 : adfGeoTransform[4] = 0.0;
1187 0 : adfGeoTransform[5] = 1.0;
1188 0 : m_oGCPSRS.Clear();
1189 : }
1190 : else
1191 : {
1192 3 : m_oSRS = std::move(oUTM);
1193 : }
1194 : }
1195 :
1196 3 : if (poTransform != nullptr)
1197 6 : delete poTransform;
1198 : }
1199 15 : else if (pszProjName != nullptr && nGCPCount == 5)
1200 : {
1201 0 : OGRSpatialReference oLL;
1202 0 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1203 :
1204 0 : if (pszOriginLong != nullptr)
1205 : {
1206 0 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
1207 : }
1208 :
1209 0 : if (pszSpheroidName == nullptr ||
1210 0 : EQUAL(pszSpheroidName, "wgs-84") || // Dash.
1211 0 : EQUAL(pszSpheroidName, "wgs_84")) // Underscore.
1212 : {
1213 0 : oLL.SetWellKnownGeogCS("WGS84");
1214 : }
1215 : else
1216 : {
1217 0 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1218 : {
1219 0 : oLL.SetGeogCS(
1220 : "", "", pszSpheroidName,
1221 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1222 : hkvEllipsoids->GetSpheroidInverseFlattening(
1223 : pszSpheroidName));
1224 : }
1225 : else
1226 : {
1227 0 : CPLError(CE_Warning, CPLE_AppDefined,
1228 : "Unrecognized ellipsoid. "
1229 : "Using wgs-84 parameters.");
1230 0 : oLL.SetWellKnownGeogCS("WGS84");
1231 : }
1232 : }
1233 :
1234 0 : const bool transform_ok = CPL_TO_BOOL(
1235 0 : GDALGCPsToGeoTransform(5, pasGCPList, adfGeoTransform, 0));
1236 :
1237 0 : m_oSRS.Clear();
1238 :
1239 0 : if (!transform_ok)
1240 : {
1241 0 : adfGeoTransform[0] = 0.0;
1242 0 : adfGeoTransform[1] = 1.0;
1243 0 : adfGeoTransform[2] = 0.0;
1244 0 : adfGeoTransform[3] = 0.0;
1245 0 : adfGeoTransform[4] = 0.0;
1246 0 : adfGeoTransform[5] = 1.0;
1247 : }
1248 : else
1249 : {
1250 0 : m_oSRS = oLL;
1251 : }
1252 :
1253 0 : m_oGCPSRS = std::move(oLL);
1254 : }
1255 :
1256 18 : delete hkvEllipsoids;
1257 : }
1258 :
1259 : /************************************************************************/
1260 : /* Open() */
1261 : /************************************************************************/
1262 :
1263 29263 : GDALDataset *HKVDataset::Open(GDALOpenInfo *poOpenInfo)
1264 :
1265 : {
1266 : /* -------------------------------------------------------------------- */
1267 : /* We assume the dataset is passed as a directory. Check for */
1268 : /* an attrib and blob file as a minimum. */
1269 : /* -------------------------------------------------------------------- */
1270 29263 : if (!poOpenInfo->bIsDirectory)
1271 29002 : return nullptr;
1272 :
1273 : const char *pszFilename =
1274 261 : CPLFormFilename(poOpenInfo->pszFilename, "image_data", nullptr);
1275 : VSIStatBuf sStat;
1276 261 : if (VSIStat(pszFilename, &sStat) != 0)
1277 217 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "blob", nullptr);
1278 261 : if (VSIStat(pszFilename, &sStat) != 0)
1279 217 : return nullptr;
1280 :
1281 44 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "attrib", nullptr);
1282 44 : if (VSIStat(pszFilename, &sStat) != 0)
1283 0 : return nullptr;
1284 :
1285 : /* -------------------------------------------------------------------- */
1286 : /* Load the attrib file, and boil white space away from around */
1287 : /* the equal sign. */
1288 : /* -------------------------------------------------------------------- */
1289 44 : char **papszAttrib = CSLLoad(pszFilename);
1290 44 : if (papszAttrib == nullptr)
1291 0 : return nullptr;
1292 :
1293 440 : for (int i = 0; papszAttrib[i] != nullptr; i++)
1294 : {
1295 396 : int iDst = 0;
1296 396 : char *pszLine = papszAttrib[i];
1297 :
1298 11173 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
1299 : {
1300 10777 : if (pszLine[iSrc] != ' ')
1301 : {
1302 9369 : pszLine[iDst++] = pszLine[iSrc];
1303 : }
1304 : }
1305 396 : pszLine[iDst] = '\0';
1306 : }
1307 :
1308 : /* -------------------------------------------------------------------- */
1309 : /* Create a corresponding GDALDataset. */
1310 : /* -------------------------------------------------------------------- */
1311 88 : auto poDS = std::make_unique<HKVDataset>();
1312 :
1313 44 : poDS->pszPath = CPLStrdup(poOpenInfo->pszFilename);
1314 44 : poDS->papszAttrib = papszAttrib;
1315 :
1316 44 : poDS->eAccess = poOpenInfo->eAccess;
1317 :
1318 : /* -------------------------------------------------------------------- */
1319 : /* Set some dataset wide information. */
1320 : /* -------------------------------------------------------------------- */
1321 44 : bool bNative = false;
1322 44 : bool bComplex = false;
1323 44 : int nRawBands = 0;
1324 :
1325 88 : if (CSLFetchNameValue(papszAttrib, "extent.cols") == nullptr ||
1326 44 : CSLFetchNameValue(papszAttrib, "extent.rows") == nullptr)
1327 : {
1328 0 : return nullptr;
1329 : }
1330 :
1331 44 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib, "extent.cols"));
1332 44 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib, "extent.rows"));
1333 :
1334 44 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
1335 : {
1336 0 : return nullptr;
1337 : }
1338 :
1339 44 : const char *pszValue = CSLFetchNameValue(papszAttrib, "pixel.order");
1340 44 : if (pszValue == nullptr)
1341 0 : bNative = true;
1342 : else
1343 : {
1344 : #ifdef CPL_MSB
1345 : bNative = strstr(pszValue, "*msbf") != NULL;
1346 : #else
1347 44 : bNative = strstr(pszValue, "*lsbf") != nullptr;
1348 : #endif
1349 : }
1350 :
1351 44 : bool bNoDataSet = false;
1352 44 : double dfNoDataValue = 0.0;
1353 44 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.no_data");
1354 44 : if (pszValue != nullptr)
1355 : {
1356 0 : bNoDataSet = true;
1357 0 : dfNoDataValue = CPLAtof(pszValue);
1358 : }
1359 :
1360 44 : pszValue = CSLFetchNameValue(papszAttrib, "channel.enumeration");
1361 44 : if (pszValue != nullptr)
1362 44 : nRawBands = atoi(pszValue);
1363 : else
1364 0 : nRawBands = 1;
1365 :
1366 44 : if (!GDALCheckBandCount(nRawBands, TRUE))
1367 : {
1368 0 : return nullptr;
1369 : }
1370 :
1371 44 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.field");
1372 44 : if (pszValue != nullptr && strstr(pszValue, "*complex") != nullptr)
1373 10 : bComplex = true;
1374 : else
1375 34 : bComplex = false;
1376 :
1377 : /* Get the version number, if present (if not, assume old version. */
1378 : /* Versions differ in their interpretation of corner coordinates. */
1379 :
1380 44 : if (CSLFetchNameValue(papszAttrib, "version") != nullptr)
1381 44 : poDS->SetVersion(static_cast<float>(
1382 44 : CPLAtof(CSLFetchNameValue(papszAttrib, "version"))));
1383 : else
1384 0 : poDS->SetVersion(1.0);
1385 :
1386 : /* -------------------------------------------------------------------- */
1387 : /* Figure out the datatype */
1388 : /* -------------------------------------------------------------------- */
1389 44 : const char *pszEncoding = CSLFetchNameValue(papszAttrib, "pixel.encoding");
1390 44 : if (pszEncoding == nullptr)
1391 0 : pszEncoding = "{ *unsigned }";
1392 :
1393 44 : int nSize = 1;
1394 44 : if (CSLFetchNameValue(papszAttrib, "pixel.size") != nullptr)
1395 44 : nSize = atoi(CSLFetchNameValue(papszAttrib, "pixel.size")) / 8;
1396 : #if 0
1397 : int nPseudoBands;
1398 : if( bComplex )
1399 : nPseudoBands = 2;
1400 : else
1401 : nPseudoBands = 1;
1402 : #endif
1403 :
1404 : GDALDataType eType;
1405 44 : if (nSize == 1)
1406 19 : eType = GDT_Byte;
1407 25 : else if (nSize == 2 && strstr(pszEncoding, "*unsigned") != nullptr)
1408 5 : eType = GDT_UInt16;
1409 20 : else if (nSize == 4 && bComplex)
1410 5 : eType = GDT_CInt16;
1411 15 : else if (nSize == 2)
1412 5 : eType = GDT_Int16;
1413 10 : else if (nSize == 4 && strstr(pszEncoding, "*unsigned") != nullptr)
1414 0 : eType = GDT_UInt32;
1415 10 : else if (nSize == 8 && strstr(pszEncoding, "*two") != nullptr && bComplex)
1416 0 : eType = GDT_CInt32;
1417 10 : else if (nSize == 4 && strstr(pszEncoding, "*two") != nullptr)
1418 0 : eType = GDT_Int32;
1419 10 : else if (nSize == 8 && bComplex)
1420 5 : eType = GDT_CFloat32;
1421 5 : else if (nSize == 4)
1422 5 : eType = GDT_Float32;
1423 0 : else if (nSize == 16 && bComplex)
1424 0 : eType = GDT_CFloat64;
1425 0 : else if (nSize == 8)
1426 0 : eType = GDT_Float64;
1427 : else
1428 : {
1429 0 : CPLError(CE_Failure, CPLE_AppDefined,
1430 : "Unsupported pixel data type in %s.\n"
1431 : "pixel.size=%d pixel.encoding=%s",
1432 0 : poDS->pszPath, nSize, pszEncoding);
1433 0 : return nullptr;
1434 : }
1435 :
1436 : /* -------------------------------------------------------------------- */
1437 : /* Open the blob file. */
1438 : /* -------------------------------------------------------------------- */
1439 44 : pszFilename = CPLFormFilename(poDS->pszPath, "image_data", nullptr);
1440 44 : if (VSIStat(pszFilename, &sStat) != 0)
1441 0 : pszFilename = CPLFormFilename(poDS->pszPath, "blob", nullptr);
1442 44 : if (poOpenInfo->eAccess == GA_ReadOnly)
1443 : {
1444 18 : poDS->fpBlob = VSIFOpenL(pszFilename, "rb");
1445 18 : if (poDS->fpBlob == nullptr)
1446 : {
1447 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1448 : "Unable to open file %s for read access.", pszFilename);
1449 0 : return nullptr;
1450 : }
1451 : }
1452 : else
1453 : {
1454 26 : poDS->fpBlob = VSIFOpenL(pszFilename, "rb+");
1455 26 : if (poDS->fpBlob == nullptr)
1456 : {
1457 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1458 : "Unable to open file %s for update access.", pszFilename);
1459 0 : return nullptr;
1460 : }
1461 : }
1462 :
1463 : /* -------------------------------------------------------------------- */
1464 : /* Build the overview filename, as blob file = "_ovr". */
1465 : /* -------------------------------------------------------------------- */
1466 88 : std::string osOvrFilename(pszFilename);
1467 44 : osOvrFilename += "_ovr";
1468 :
1469 : /* -------------------------------------------------------------------- */
1470 : /* Define the bands. */
1471 : /* -------------------------------------------------------------------- */
1472 44 : const int nPixelOffset = nRawBands * nSize;
1473 44 : const int nLineOffset = nPixelOffset * poDS->GetRasterXSize();
1474 44 : int nOffset = 0;
1475 :
1476 138 : for (int iRawBand = 0; iRawBand < nRawBands; iRawBand++)
1477 : {
1478 : auto poBand = std::make_unique<HKVRasterBand>(
1479 94 : poDS.get(), poDS->GetRasterCount() + 1, poDS->fpBlob, nOffset,
1480 94 : nPixelOffset, nLineOffset, eType, bNative);
1481 94 : if (!poBand->IsValid())
1482 0 : return nullptr;
1483 :
1484 94 : if (bNoDataSet)
1485 0 : poBand->SetNoDataValue(dfNoDataValue);
1486 94 : poDS->SetBand(poDS->GetRasterCount() + 1, std::move(poBand));
1487 94 : nOffset += GDALGetDataTypeSizeBytes(eType);
1488 : }
1489 :
1490 44 : poDS->eRasterType = eType;
1491 :
1492 : /* -------------------------------------------------------------------- */
1493 : /* Process the georef file if there is one. */
1494 : /* -------------------------------------------------------------------- */
1495 44 : pszFilename = CPLFormFilename(poDS->pszPath, "georef", nullptr);
1496 44 : if (VSIStat(pszFilename, &sStat) == 0)
1497 18 : poDS->ProcessGeoref(pszFilename);
1498 :
1499 : /* -------------------------------------------------------------------- */
1500 : /* Initialize any PAM information. */
1501 : /* -------------------------------------------------------------------- */
1502 44 : poDS->SetDescription(osOvrFilename.c_str());
1503 44 : poDS->TryLoadXML();
1504 :
1505 : /* -------------------------------------------------------------------- */
1506 : /* Handle overviews. */
1507 : /* -------------------------------------------------------------------- */
1508 44 : poDS->oOvManager.Initialize(poDS.get(), osOvrFilename.c_str(), nullptr,
1509 : TRUE);
1510 :
1511 44 : return poDS.release();
1512 : }
1513 :
1514 : /************************************************************************/
1515 : /* Create() */
1516 : /************************************************************************/
1517 :
1518 51 : GDALDataset *HKVDataset::Create(const char *pszFilenameIn, int nXSize,
1519 : int nYSize, int nBandsIn, GDALDataType eType,
1520 : char ** /* papszParamList */)
1521 :
1522 : {
1523 : /* -------------------------------------------------------------------- */
1524 : /* Verify input options. */
1525 : /* -------------------------------------------------------------------- */
1526 51 : if (nBandsIn <= 0)
1527 : {
1528 1 : CPLError(CE_Failure, CPLE_NotSupported,
1529 : "HKV driver does not support %d bands.", nBandsIn);
1530 1 : return nullptr;
1531 : }
1532 :
1533 50 : if (eType != GDT_Byte && eType != GDT_UInt16 && eType != GDT_Int16 &&
1534 27 : eType != GDT_CInt16 && eType != GDT_Float32 && eType != GDT_CFloat32)
1535 : {
1536 21 : CPLError(CE_Failure, CPLE_AppDefined,
1537 : "Attempt to create HKV file with currently unsupported\n"
1538 : "data type (%s).",
1539 : GDALGetDataTypeName(eType));
1540 :
1541 21 : return nullptr;
1542 : }
1543 :
1544 : /* -------------------------------------------------------------------- */
1545 : /* Establish the name of the directory we will be creating the */
1546 : /* new HKV directory in. Verify that this is a directory. */
1547 : /* -------------------------------------------------------------------- */
1548 29 : char *pszBaseDir = nullptr;
1549 :
1550 29 : if (strlen(CPLGetPath(pszFilenameIn)) == 0)
1551 0 : pszBaseDir = CPLStrdup(".");
1552 : else
1553 29 : pszBaseDir = CPLStrdup(CPLGetPath(pszFilenameIn));
1554 :
1555 : VSIStatBuf sStat;
1556 29 : if (CPLStat(pszBaseDir, &sStat) != 0 || !VSI_ISDIR(sStat.st_mode))
1557 : {
1558 3 : CPLError(CE_Failure, CPLE_AppDefined,
1559 : "Attempt to create HKV dataset under %s,\n"
1560 : "but this is not a valid directory.",
1561 : pszBaseDir);
1562 3 : CPLFree(pszBaseDir);
1563 3 : return nullptr;
1564 : }
1565 :
1566 26 : CPLFree(pszBaseDir);
1567 26 : pszBaseDir = nullptr;
1568 :
1569 26 : if (VSIMkdir(pszFilenameIn, 0755) != 0)
1570 : {
1571 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unable to create directory %s.",
1572 : pszFilenameIn);
1573 0 : return nullptr;
1574 : }
1575 :
1576 : /* -------------------------------------------------------------------- */
1577 : /* Create the header file. */
1578 : /* -------------------------------------------------------------------- */
1579 26 : CPLErr CEHeaderCreated = SaveHKVAttribFile(pszFilenameIn, nXSize, nYSize,
1580 : nBandsIn, eType, FALSE, 0.0);
1581 :
1582 26 : if (CEHeaderCreated != CE_None)
1583 0 : return nullptr;
1584 :
1585 : /* -------------------------------------------------------------------- */
1586 : /* Create the blob file. */
1587 : /* -------------------------------------------------------------------- */
1588 :
1589 : const char *pszFilename =
1590 26 : CPLFormFilename(pszFilenameIn, "image_data", nullptr);
1591 26 : FILE *fp = VSIFOpen(pszFilename, "wb");
1592 26 : if (fp == nullptr)
1593 : {
1594 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.\n",
1595 : pszFilename);
1596 0 : return nullptr;
1597 : }
1598 :
1599 26 : bool bOK = VSIFWrite(reinterpret_cast<void *>(const_cast<char *>("")), 1, 1,
1600 26 : fp) == 1;
1601 26 : if (VSIFClose(fp) != 0)
1602 0 : bOK &= false;
1603 :
1604 26 : if (!bOK)
1605 0 : return nullptr;
1606 : /* -------------------------------------------------------------------- */
1607 : /* Open the dataset normally. */
1608 : /* -------------------------------------------------------------------- */
1609 26 : return GDALDataset::FromHandle(GDALOpen(pszFilenameIn, GA_Update));
1610 : }
1611 :
1612 : /************************************************************************/
1613 : /* Delete() */
1614 : /* */
1615 : /* An HKV Blob dataset consists of a bunch of files in a */
1616 : /* directory. Try to delete all the files, then the */
1617 : /* directory. */
1618 : /************************************************************************/
1619 :
1620 1 : CPLErr HKVDataset::Delete(const char *pszName)
1621 :
1622 : {
1623 : VSIStatBuf sStat;
1624 1 : if (CPLStat(pszName, &sStat) != 0 || !VSI_ISDIR(sStat.st_mode))
1625 : {
1626 0 : CPLError(CE_Failure, CPLE_AppDefined,
1627 : "%s does not appear to be an HKV Dataset, as it is not "
1628 : "a path to a directory.",
1629 : pszName);
1630 0 : return CE_Failure;
1631 : }
1632 :
1633 1 : char **papszFiles = VSIReadDir(pszName);
1634 6 : for (int i = 0; i < CSLCount(papszFiles); i++)
1635 : {
1636 5 : if (EQUAL(papszFiles[i], ".") || EQUAL(papszFiles[i], ".."))
1637 2 : continue;
1638 :
1639 : const char *pszTarget =
1640 3 : CPLFormFilename(pszName, papszFiles[i], nullptr);
1641 3 : if (VSIUnlink(pszTarget) != 0)
1642 : {
1643 0 : CPLError(CE_Failure, CPLE_AppDefined,
1644 : "Unable to delete file %s,"
1645 : "HKVDataset Delete(%s) failed.",
1646 : pszTarget, pszName);
1647 0 : CSLDestroy(papszFiles);
1648 0 : return CE_Failure;
1649 : }
1650 : }
1651 :
1652 1 : CSLDestroy(papszFiles);
1653 :
1654 1 : if (VSIRmdir(pszName) != 0)
1655 : {
1656 0 : CPLError(CE_Failure, CPLE_AppDefined,
1657 : "Unable to delete directory %s,"
1658 : "HKVDataset Delete() failed.",
1659 : pszName);
1660 0 : return CE_Failure;
1661 : }
1662 :
1663 1 : return CE_None;
1664 : }
1665 :
1666 : /************************************************************************/
1667 : /* CreateCopy() */
1668 : /************************************************************************/
1669 :
1670 20 : GDALDataset *HKVDataset::CreateCopy(const char *pszFilename,
1671 : GDALDataset *poSrcDS,
1672 : CPL_UNUSED int bStrict, char **papszOptions,
1673 : GDALProgressFunc pfnProgress,
1674 : void *pProgressData)
1675 : {
1676 20 : int nBands = poSrcDS->GetRasterCount();
1677 20 : if (nBands == 0)
1678 : {
1679 1 : CPLError(CE_Failure, CPLE_NotSupported,
1680 : "HKV driver does not support source dataset with zero band.");
1681 1 : return nullptr;
1682 : }
1683 :
1684 19 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1685 :
1686 19 : if (!pfnProgress(0.0, nullptr, pProgressData))
1687 0 : return nullptr;
1688 :
1689 : /* check that other bands match type- sets type */
1690 : /* to unknown if they differ. */
1691 29 : for (int iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++)
1692 : {
1693 10 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
1694 10 : eType = GDALDataTypeUnion(eType, poBand->GetRasterDataType());
1695 : }
1696 :
1697 19 : HKVDataset *poDS = reinterpret_cast<HKVDataset *>(Create(
1698 : pszFilename, poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
1699 : poSrcDS->GetRasterCount(), eType, papszOptions));
1700 :
1701 : /* Check that Create worked- return Null if it didn't */
1702 19 : if (poDS == nullptr)
1703 8 : return nullptr;
1704 :
1705 : /* -------------------------------------------------------------------- */
1706 : /* Copy the image data. */
1707 : /* -------------------------------------------------------------------- */
1708 11 : const int nXSize = poDS->GetRasterXSize();
1709 11 : const int nYSize = poDS->GetRasterYSize();
1710 :
1711 : int nBlockXSize, nBlockYSize;
1712 11 : poDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
1713 :
1714 11 : const int nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize) *
1715 11 : ((nYSize + nBlockYSize - 1) / nBlockYSize) *
1716 11 : poSrcDS->GetRasterCount();
1717 :
1718 11 : int nBlocksDone = 0;
1719 32 : for (int iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++)
1720 : {
1721 21 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand + 1);
1722 21 : GDALRasterBand *poDstBand = poDS->GetRasterBand(iBand + 1);
1723 :
1724 : /* Get nodata value, if relevant */
1725 21 : int pbSuccess = FALSE;
1726 21 : double dfSrcNoDataValue = poSrcBand->GetNoDataValue(&pbSuccess);
1727 21 : if (pbSuccess)
1728 0 : poDS->SetNoDataValue(dfSrcNoDataValue);
1729 :
1730 63 : void *pData = CPLMalloc(nBlockXSize * nBlockYSize *
1731 21 : GDALGetDataTypeSize(eType) / 8);
1732 :
1733 241 : for (int iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize)
1734 : {
1735 440 : for (int iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize)
1736 : {
1737 220 : if (!pfnProgress((nBlocksDone++) /
1738 220 : static_cast<float>(nBlockTotal),
1739 : nullptr, pProgressData))
1740 : {
1741 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1742 0 : delete poDS;
1743 0 : CPLFree(pData);
1744 :
1745 : GDALDriver *poHKVDriver = reinterpret_cast<GDALDriver *>(
1746 0 : GDALGetDriverByName("MFF2"));
1747 0 : poHKVDriver->Delete(pszFilename);
1748 0 : return nullptr;
1749 : }
1750 :
1751 220 : const int nTBXSize = std::min(nBlockXSize, nXSize - iXOffset);
1752 220 : const int nTBYSize = std::min(nBlockYSize, nYSize - iYOffset);
1753 :
1754 220 : CPLErr eErr = poSrcBand->RasterIO(
1755 : GF_Read, iXOffset, iYOffset, nTBXSize, nTBYSize, pData,
1756 : nTBXSize, nTBYSize, eType, 0, 0, nullptr);
1757 220 : if (eErr != CE_None)
1758 : {
1759 0 : delete poDS;
1760 0 : CPLFree(pData);
1761 0 : return nullptr;
1762 : }
1763 :
1764 220 : eErr = poDstBand->RasterIO(GF_Write, iXOffset, iYOffset,
1765 : nTBXSize, nTBYSize, pData, nTBXSize,
1766 : nTBYSize, eType, 0, 0, nullptr);
1767 :
1768 220 : if (eErr != CE_None)
1769 : {
1770 0 : delete poDS;
1771 0 : CPLFree(pData);
1772 0 : return nullptr;
1773 : }
1774 : }
1775 : }
1776 :
1777 21 : CPLFree(pData);
1778 : }
1779 :
1780 : /* -------------------------------------------------------------------- */
1781 : /* Copy georeferencing information, if enough is available. */
1782 : /* Only copy geotransform-style info (won't work for slant range). */
1783 : /* -------------------------------------------------------------------- */
1784 :
1785 : double *tempGeoTransform =
1786 11 : static_cast<double *>(CPLMalloc(6 * sizeof(double)));
1787 :
1788 22 : if ((poSrcDS->GetGeoTransform(tempGeoTransform) == CE_None) &&
1789 11 : (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0 ||
1790 0 : tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0 ||
1791 0 : tempGeoTransform[4] != 0.0 || std::abs(tempGeoTransform[5]) != 1.0))
1792 : {
1793 11 : auto poSrcSRS = poSrcDS->GetSpatialRef();
1794 11 : if (poSrcSRS)
1795 : {
1796 11 : poDS->SetSpatialRef(poSrcSRS);
1797 11 : poDS->m_oGCPSRS = *poSrcSRS;
1798 : }
1799 11 : poDS->SetGeoTransform(tempGeoTransform);
1800 :
1801 11 : CPLFree(tempGeoTransform);
1802 :
1803 : // georef file will be saved automatically when dataset is deleted
1804 : // because SetProjection sets a flag to indicate it is necessary.
1805 : }
1806 : else
1807 : {
1808 0 : CPLFree(tempGeoTransform);
1809 : }
1810 :
1811 : // Make sure image data gets flushed.
1812 32 : for (int iBand = 0; iBand < poDS->GetRasterCount(); iBand++)
1813 : {
1814 : RawRasterBand *poDstBand =
1815 21 : reinterpret_cast<RawRasterBand *>(poDS->GetRasterBand(iBand + 1));
1816 21 : poDstBand->FlushCache(false);
1817 : }
1818 :
1819 11 : if (!pfnProgress(1.0, nullptr, pProgressData))
1820 : {
1821 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1822 0 : delete poDS;
1823 :
1824 : GDALDriver *poHKVDriver =
1825 0 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("MFF2"));
1826 0 : poHKVDriver->Delete(pszFilename);
1827 0 : return nullptr;
1828 : }
1829 :
1830 11 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1831 :
1832 11 : return poDS;
1833 : }
1834 :
1835 : /************************************************************************/
1836 : /* GDALRegister_HKV() */
1837 : /************************************************************************/
1838 :
1839 1520 : void GDALRegister_HKV()
1840 :
1841 : {
1842 1520 : if (GDALGetDriverByName("MFF2") != nullptr)
1843 301 : return;
1844 :
1845 1219 : GDALDriver *poDriver = new GDALDriver();
1846 :
1847 1219 : poDriver->SetDescription("MFF2");
1848 1219 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1849 1219 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF2 (HKV) Raster");
1850 1219 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff2.html");
1851 1219 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1852 : "Byte Int16 UInt16 Int32 UInt32 CInt16 "
1853 1219 : "CInt32 Float32 Float64 CFloat32 CFloat64");
1854 :
1855 1219 : poDriver->pfnOpen = HKVDataset::Open;
1856 1219 : poDriver->pfnCreate = HKVDataset::Create;
1857 1219 : poDriver->pfnDelete = HKVDataset::Delete;
1858 1219 : poDriver->pfnCreateCopy = HKVDataset::CreateCopy;
1859 :
1860 1219 : GetGDALDriverManager()->RegisterDriver(poDriver);
1861 : }
|
