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 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include <ctype.h>
15 :
16 : #include "atlsci_spheroid.h"
17 : #include "cpl_string.h"
18 : #include "gdal_frmts.h"
19 : #include "gdal_priv.h"
20 : #include "ogr_spatialref.h"
21 : #include "rawdataset.h"
22 :
23 : #include <cmath>
24 :
25 : #include <algorithm>
26 :
27 : /************************************************************************/
28 : /* ==================================================================== */
29 : /* HKVRasterBand */
30 : /* ==================================================================== */
31 : /************************************************************************/
32 :
33 : class HKVDataset;
34 :
35 2 : class HKVRasterBand final : public RawRasterBand
36 : {
37 : friend class HKVDataset;
38 :
39 : public:
40 : HKVRasterBand(HKVDataset *poDS, int nBand, VSILFILE *fpRaw,
41 : unsigned int nImgOffset, int nPixelOffset, int nLineOffset,
42 : GDALDataType eDataType, int bNativeOrder);
43 :
44 : ~HKVRasterBand() override;
45 : };
46 :
47 : HKVRasterBand::~HKVRasterBand() = default;
48 :
49 : /************************************************************************/
50 : /* HKV Spheroids */
51 : /************************************************************************/
52 :
53 : class HKVSpheroidList final : public SpheroidList
54 : {
55 : public:
56 : HKVSpheroidList();
57 : };
58 :
59 1 : HKVSpheroidList ::HKVSpheroidList()
60 : {
61 1 : num_spheroids = 58;
62 1 : epsilonR = 0.1;
63 1 : epsilonI = 0.000001;
64 :
65 1 : spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830", 6377563.396,
66 : 299.3249646);
67 1 : spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",
68 : 6377340.189, 299.3249646);
69 1 : spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",
70 : 6378160, 298.25);
71 1 : spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",
72 : 6377483.865, 299.1528128);
73 1 : spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841", 6377397.155,
74 : 299.1528128);
75 1 : spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858", 6378294.0,
76 : 294.297);
77 1 : spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866", 6378206.4,
78 : 294.9786982);
79 1 : spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880", 6378249.145,
80 : 293.465);
81 1 : spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",
82 : 6377276.345, 300.8017);
83 1 : spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",
84 : 6377298.556, 300.8017);
85 1 : spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",
86 : 6377301.243, 300.8017);
87 1 : spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",
88 : 6377295.664, 300.8017);
89 1 : spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",
90 : 6377304.063, 300.8017);
91 1 : spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",
92 : 6377309.613, 300.8017);
93 1 : spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",
94 : 6378155, 298.3);
95 1 : spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906", 6378200,
96 : 298.3);
97 1 : spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960", 6378270,
98 : 297);
99 1 : spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes", 6378273.0,
100 : 298.279);
101 1 : spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",
102 : 6378160, 298.247);
103 1 : spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",
104 : 6378388, 297);
105 1 : spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67", 6378160.0,
106 : 298.254);
107 1 : spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75", 6378140.0,
108 : 298.25298);
109 1 : spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",
110 : 6378245, 298.3);
111 1 : spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula", 6378165.0,
112 : 292.308);
113 1 : spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80", 6378137,
114 : 298.257222101);
115 1 : spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",
116 : 6378160, 298.25);
117 1 : spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72", 6378135,
118 : 298.26);
119 1 : spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84", 6378137,
120 : 298.257223563);
121 1 : spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84", 6378137.0,
122 : 298.252841);
123 1 : spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel", 6377397.0,
124 : 299.1976073);
125 :
126 1 : spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830", 6377563.396,
127 : 299.3249646);
128 1 : spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",
129 : 6377340.189, 299.3249646);
130 1 : spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",
131 : 6378160, 298.25);
132 1 : spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",
133 : 6377483.865, 299.1528128);
134 1 : spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",
135 : 6377397.155, 299.1528128);
136 1 : spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858", 6378294.0,
137 : 294.297);
138 1 : spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866", 6378206.4,
139 : 294.9786982);
140 1 : spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",
141 : 6378249.145, 293.465);
142 1 : spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",
143 : 6377276.345, 300.8017);
144 1 : spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",
145 : 6377298.556, 300.8017);
146 1 : spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",
147 : 6377301.243, 300.8017);
148 1 : spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",
149 : 6377295.664, 300.8017);
150 1 : spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",
151 : 6377304.063, 300.8017);
152 1 : spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",
153 : 6377309.613, 300.8017);
154 1 : spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",
155 : 6378155, 298.3);
156 1 : spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906", 6378200,
157 : 298.3);
158 1 : spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960", 6378270,
159 : 297);
160 1 : spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",
161 : 6378160, 298.247);
162 1 : spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",
163 : 6378388, 297);
164 1 : spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67", 6378160.0,
165 : 298.254);
166 1 : spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75", 6378140.0,
167 : 298.25298);
168 1 : spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",
169 : 6378245, 298.3);
170 1 : spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80", 6378137,
171 : 298.257222101);
172 1 : spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",
173 : 6378160, 298.25);
174 1 : spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72", 6378135,
175 : 298.26);
176 1 : spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84", 6378137,
177 : 298.257223563);
178 1 : spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84", 6378137.0,
179 : 298.252841);
180 1 : spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel", 6377397.0,
181 : 299.1976073);
182 1 : }
183 :
184 : /************************************************************************/
185 : /* ==================================================================== */
186 : /* HKVDataset */
187 : /* ==================================================================== */
188 : /************************************************************************/
189 :
190 : class HKVDataset final : public RawDataset
191 : {
192 : friend class HKVRasterBand;
193 :
194 : char *pszPath;
195 : VSILFILE *fpBlob;
196 :
197 : int nGCPCount;
198 : GDAL_GCP *pasGCPList;
199 :
200 : void ProcessGeoref(const char *);
201 : void ProcessGeorefGCP(char **, const char *, double, double);
202 :
203 1 : void SetVersion(float version_number)
204 : {
205 : // Update stored info.
206 1 : MFF2version = version_number;
207 1 : }
208 :
209 : float MFF2version;
210 :
211 : GDALDataType eRasterType;
212 :
213 : void SetNoDataValue(double);
214 :
215 : OGRSpatialReference m_oSRS{};
216 : OGRSpatialReference m_oGCPSRS{};
217 : GDALGeoTransform m_gt{};
218 :
219 : char **papszAttrib;
220 :
221 : char **papszGeoref;
222 :
223 : // NOTE: The MFF2 format goes against GDAL's API in that nodata values are
224 : // set per-dataset rather than per-band. To compromise, for writing out,
225 : // the dataset's nodata value will be set to the last value set on any of
226 : // the raster bands.
227 :
228 : bool bNoDataSet;
229 : double dfNoDataValue;
230 :
231 : CPL_DISALLOW_COPY_ASSIGN(HKVDataset)
232 :
233 : CPLErr Close() override;
234 :
235 : public:
236 : HKVDataset();
237 : ~HKVDataset() override;
238 :
239 0 : int GetGCPCount() override /* const */
240 : {
241 0 : return nGCPCount;
242 : }
243 :
244 0 : const OGRSpatialReference *GetGCPSpatialRef() const override
245 : {
246 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
247 : }
248 :
249 : const GDAL_GCP *GetGCPs() override;
250 :
251 0 : const OGRSpatialReference *GetSpatialRef() const override
252 : {
253 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
254 : }
255 :
256 : CPLErr GetGeoTransform(GDALGeoTransform &) const override;
257 :
258 : static GDALDataset *Open(GDALOpenInfo *);
259 : };
260 :
261 : /************************************************************************/
262 : /* ==================================================================== */
263 : /* HKVRasterBand */
264 : /* ==================================================================== */
265 : /************************************************************************/
266 :
267 : /************************************************************************/
268 : /* HKVRasterBand() */
269 : /************************************************************************/
270 :
271 1 : HKVRasterBand::HKVRasterBand(HKVDataset *poDSIn, int nBandIn, VSILFILE *fpRawIn,
272 : unsigned int nImgOffsetIn, int nPixelOffsetIn,
273 : int nLineOffsetIn, GDALDataType eDataTypeIn,
274 1 : int bNativeOrderIn)
275 : : RawRasterBand(GDALDataset::FromHandle(poDSIn), nBandIn, fpRawIn,
276 : nImgOffsetIn, nPixelOffsetIn, nLineOffsetIn, eDataTypeIn,
277 1 : bNativeOrderIn, RawRasterBand::OwnFP::NO)
278 :
279 : {
280 1 : poDS = poDSIn;
281 1 : nBand = nBandIn;
282 :
283 1 : nBlockXSize = poDS->GetRasterXSize();
284 1 : nBlockYSize = 1;
285 1 : }
286 :
287 : /************************************************************************/
288 : /* ==================================================================== */
289 : /* HKVDataset */
290 : /* ==================================================================== */
291 : /************************************************************************/
292 :
293 : /************************************************************************/
294 : /* HKVDataset() */
295 : /************************************************************************/
296 :
297 1 : HKVDataset::HKVDataset()
298 : : pszPath(nullptr), fpBlob(nullptr), nGCPCount(0), pasGCPList(nullptr),
299 : // Initialize datasets to new version; change if necessary.
300 : MFF2version(1.1f), eRasterType(GDT_Unknown), papszAttrib(nullptr),
301 1 : papszGeoref(nullptr), bNoDataSet(false), dfNoDataValue(0.0)
302 : {
303 1 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
304 1 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
305 1 : }
306 :
307 : /************************************************************************/
308 : /* ~HKVDataset() */
309 : /************************************************************************/
310 :
311 2 : HKVDataset::~HKVDataset()
312 :
313 : {
314 1 : HKVDataset::Close();
315 2 : }
316 :
317 : /************************************************************************/
318 : /* Close() */
319 : /************************************************************************/
320 :
321 2 : CPLErr HKVDataset::Close()
322 : {
323 2 : CPLErr eErr = CE_None;
324 2 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
325 : {
326 1 : if (HKVDataset::FlushCache(true) != CE_None)
327 0 : eErr = CE_Failure;
328 :
329 1 : if (fpBlob)
330 : {
331 1 : if (VSIFCloseL(fpBlob) != 0)
332 : {
333 0 : eErr = CE_Failure;
334 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
335 : }
336 : }
337 :
338 1 : if (nGCPCount > 0)
339 : {
340 1 : GDALDeinitGCPs(nGCPCount, pasGCPList);
341 1 : CPLFree(pasGCPList);
342 : }
343 :
344 1 : CPLFree(pszPath);
345 1 : CSLDestroy(papszGeoref);
346 1 : CSLDestroy(papszAttrib);
347 :
348 1 : if (GDALPamDataset::Close() != CE_None)
349 0 : eErr = CE_Failure;
350 : }
351 2 : return eErr;
352 : }
353 :
354 : /************************************************************************/
355 : /* GetGeoTransform() */
356 : /************************************************************************/
357 :
358 0 : CPLErr HKVDataset::GetGeoTransform(GDALGeoTransform >) const
359 :
360 : {
361 0 : gt = m_gt;
362 0 : return CE_None;
363 : }
364 :
365 : /************************************************************************/
366 : /* GetGCP() */
367 : /************************************************************************/
368 :
369 0 : const GDAL_GCP *HKVDataset::GetGCPs()
370 :
371 : {
372 0 : return pasGCPList;
373 : }
374 :
375 : /************************************************************************/
376 : /* ProcessGeorefGCP() */
377 : /************************************************************************/
378 :
379 5 : void HKVDataset::ProcessGeorefGCP(char **papszGeorefIn, const char *pszBase,
380 : double dfRasterX, double dfRasterY)
381 :
382 : {
383 : /* -------------------------------------------------------------------- */
384 : /* Fetch the GCP from the string list. */
385 : /* -------------------------------------------------------------------- */
386 5 : char szFieldName[128] = {'\0'};
387 5 : snprintf(szFieldName, sizeof(szFieldName), "%s.latitude", pszBase);
388 5 : double dfLat = 0.0;
389 5 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
390 0 : return;
391 : else
392 5 : dfLat = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
393 :
394 5 : snprintf(szFieldName, sizeof(szFieldName), "%s.longitude", pszBase);
395 5 : double dfLong = 0.0;
396 5 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
397 0 : return;
398 : else
399 5 : dfLong = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
400 :
401 : /* -------------------------------------------------------------------- */
402 : /* Add the gcp to the internal list. */
403 : /* -------------------------------------------------------------------- */
404 5 : GDALInitGCPs(1, pasGCPList + nGCPCount);
405 :
406 5 : CPLFree(pasGCPList[nGCPCount].pszId);
407 :
408 5 : pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
409 :
410 5 : pasGCPList[nGCPCount].dfGCPX = dfLong;
411 5 : pasGCPList[nGCPCount].dfGCPY = dfLat;
412 5 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
413 :
414 5 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
415 5 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
416 :
417 5 : nGCPCount++;
418 : }
419 :
420 : /************************************************************************/
421 : /* ProcessGeoref() */
422 : /************************************************************************/
423 :
424 1 : void HKVDataset::ProcessGeoref(const char *pszFilename)
425 :
426 : {
427 : /* -------------------------------------------------------------------- */
428 : /* Load the georef file, and boil white space away from around */
429 : /* the equal sign. */
430 : /* -------------------------------------------------------------------- */
431 1 : CSLDestroy(papszGeoref);
432 1 : papszGeoref = CSLLoad(pszFilename);
433 1 : if (papszGeoref == nullptr)
434 0 : return;
435 :
436 1 : HKVSpheroidList *hkvEllipsoids = new HKVSpheroidList;
437 :
438 14 : for (int i = 0; papszGeoref[i] != nullptr; i++)
439 : {
440 13 : int iDst = 0;
441 13 : char *pszLine = papszGeoref[i];
442 :
443 433 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
444 : {
445 420 : if (pszLine[iSrc] != ' ')
446 : {
447 420 : pszLine[iDst++] = pszLine[iSrc];
448 : }
449 : }
450 13 : pszLine[iDst] = '\0';
451 : }
452 :
453 : /* -------------------------------------------------------------------- */
454 : /* Try to get GCPs, in lat/longs . */
455 : /* -------------------------------------------------------------------- */
456 1 : nGCPCount = 0;
457 1 : pasGCPList = static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), 5));
458 :
459 1 : if (MFF2version > 1.0)
460 : {
461 1 : ProcessGeorefGCP(papszGeoref, "top_left", 0, 0);
462 1 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize(), 0);
463 1 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0, GetRasterYSize());
464 1 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize(),
465 1 : GetRasterYSize());
466 1 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
467 1 : GetRasterYSize() / 2.0);
468 : }
469 : else
470 : {
471 0 : ProcessGeorefGCP(papszGeoref, "top_left", 0.5, 0.5);
472 0 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize() - 0.5, 0.5);
473 0 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0.5,
474 0 : GetRasterYSize() - 0.5);
475 0 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize() - 0.5,
476 0 : GetRasterYSize() - 0.5);
477 0 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
478 0 : GetRasterYSize() / 2.0);
479 : }
480 :
481 1 : if (nGCPCount == 0)
482 : {
483 0 : CPLFree(pasGCPList);
484 0 : pasGCPList = nullptr;
485 : }
486 :
487 : /* -------------------------------------------------------------------- */
488 : /* Do we have a recognised projection? */
489 : /* -------------------------------------------------------------------- */
490 1 : const char *pszProjName = CSLFetchNameValue(papszGeoref, "projection.name");
491 : const char *pszOriginLong =
492 1 : CSLFetchNameValue(papszGeoref, "projection.origin_longitude");
493 : const char *pszSpheroidName =
494 1 : CSLFetchNameValue(papszGeoref, "spheroid.name");
495 :
496 2 : if (pszSpheroidName != nullptr &&
497 1 : hkvEllipsoids->SpheroidInList(pszSpheroidName))
498 : {
499 : #if 0
500 : // TODO(schwehr): Enable in trunk after 2.1 branch and fix.
501 : // Breaks tests on some platforms.
502 : CPLError( CE_Failure, CPLE_AppDefined,
503 : "Unrecognized ellipsoid. Not handled. "
504 : "Spheroid name not in spheroid list: '%s'",
505 : pszSpheroidName );
506 : #endif
507 : // Why were eq_radius and inv_flattening never used?
508 : // eq_radius = hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
509 : // inv_flattening =
510 : // hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
511 : }
512 0 : else if (pszProjName != nullptr)
513 : {
514 0 : CPLError(CE_Warning, CPLE_AppDefined,
515 : "Unrecognized ellipsoid. Not handled.");
516 : // TODO(schwehr): This error is was never what was happening.
517 : // CPLError( CE_Warning, CPLE_AppDefined,
518 : // "Unrecognized ellipsoid. Using wgs-84 parameters.");
519 : // eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84"); */
520 : // inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
521 : }
522 :
523 1 : if (pszProjName != nullptr && EQUAL(pszProjName, "utm") && nGCPCount == 5)
524 : {
525 : // int nZone = (int)((CPLAtof(pszOriginLong)+184.5) / 6.0);
526 1 : int nZone = 31; // TODO(schwehr): Where does 31 come from?
527 :
528 1 : if (pszOriginLong == nullptr)
529 : {
530 : // If origin not specified, assume 0.0.
531 0 : CPLError(
532 : CE_Warning, CPLE_AppDefined,
533 : "No projection origin longitude specified. Assuming 0.0.");
534 : }
535 : else
536 : {
537 1 : nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
538 : }
539 :
540 2 : OGRSpatialReference oUTM;
541 :
542 1 : if (pasGCPList[4].dfGCPY < 0)
543 0 : oUTM.SetUTM(nZone, 0);
544 : else
545 1 : oUTM.SetUTM(nZone, 1);
546 :
547 2 : OGRSpatialReference oLL;
548 1 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
549 1 : if (pszOriginLong != nullptr)
550 : {
551 1 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
552 1 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
553 : }
554 :
555 1 : if ((pszSpheroidName == nullptr) ||
556 1 : (EQUAL(pszSpheroidName, "wgs-84")) ||
557 1 : (EQUAL(pszSpheroidName, "wgs_84")))
558 : {
559 0 : oUTM.SetWellKnownGeogCS("WGS84");
560 0 : oLL.SetWellKnownGeogCS("WGS84");
561 : }
562 : else
563 : {
564 1 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
565 : {
566 1 : oUTM.SetGeogCS(
567 : "unknown", "unknown", pszSpheroidName,
568 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
569 : hkvEllipsoids->GetSpheroidInverseFlattening(
570 : pszSpheroidName));
571 1 : oLL.SetGeogCS(
572 : "unknown", "unknown", pszSpheroidName,
573 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
574 : hkvEllipsoids->GetSpheroidInverseFlattening(
575 : pszSpheroidName));
576 : }
577 : else
578 : {
579 0 : CPLError(CE_Warning, CPLE_AppDefined,
580 : "Unrecognized ellipsoid. Using wgs-84 parameters.");
581 0 : oUTM.SetWellKnownGeogCS("WGS84");
582 0 : oLL.SetWellKnownGeogCS("WGS84");
583 : }
584 : }
585 :
586 : OGRCoordinateTransformation *poTransform =
587 1 : OGRCreateCoordinateTransformation(&oLL, &oUTM);
588 :
589 1 : bool bSuccess = true;
590 1 : if (poTransform == nullptr)
591 : {
592 0 : CPLErrorReset();
593 0 : bSuccess = false;
594 : }
595 :
596 1 : double dfUtmX[5] = {0.0};
597 1 : double dfUtmY[5] = {0.0};
598 :
599 1 : if (poTransform != nullptr)
600 : {
601 6 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
602 : {
603 5 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
604 5 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
605 :
606 5 : if (bSuccess && !poTransform->Transform(1, &(dfUtmX[gcp_index]),
607 : &(dfUtmY[gcp_index])))
608 0 : bSuccess = false;
609 : }
610 : }
611 :
612 1 : if (bSuccess)
613 : {
614 : // Update GCPS to proper projection.
615 6 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
616 : {
617 5 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
618 5 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
619 : }
620 :
621 1 : m_oGCPSRS = oUTM;
622 :
623 2 : bool transform_ok = CPL_TO_BOOL(
624 1 : GDALGCPsToGeoTransform(5, pasGCPList, m_gt.data(), 0));
625 :
626 1 : if (!transform_ok)
627 : {
628 : // Transform may not be sufficient in all cases (slant range
629 : // projection).
630 0 : m_gt = GDALGeoTransform();
631 0 : m_oGCPSRS.Clear();
632 : }
633 : else
634 : {
635 1 : m_oSRS = std::move(oUTM);
636 : }
637 : }
638 :
639 1 : if (poTransform != nullptr)
640 2 : delete poTransform;
641 : }
642 0 : else if (pszProjName != nullptr && nGCPCount == 5)
643 : {
644 0 : OGRSpatialReference oLL;
645 0 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
646 :
647 0 : if (pszOriginLong != nullptr)
648 : {
649 0 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
650 : }
651 :
652 0 : if (pszSpheroidName == nullptr ||
653 0 : EQUAL(pszSpheroidName, "wgs-84") || // Dash.
654 0 : EQUAL(pszSpheroidName, "wgs_84")) // Underscore.
655 : {
656 0 : oLL.SetWellKnownGeogCS("WGS84");
657 : }
658 : else
659 : {
660 0 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
661 : {
662 0 : oLL.SetGeogCS(
663 : "", "", pszSpheroidName,
664 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
665 : hkvEllipsoids->GetSpheroidInverseFlattening(
666 : pszSpheroidName));
667 : }
668 : else
669 : {
670 0 : CPLError(CE_Warning, CPLE_AppDefined,
671 : "Unrecognized ellipsoid. "
672 : "Using wgs-84 parameters.");
673 0 : oLL.SetWellKnownGeogCS("WGS84");
674 : }
675 : }
676 :
677 : const bool transform_ok =
678 0 : CPL_TO_BOOL(GDALGCPsToGeoTransform(5, pasGCPList, m_gt.data(), 0));
679 :
680 0 : m_oSRS.Clear();
681 :
682 0 : if (!transform_ok)
683 : {
684 0 : m_gt = GDALGeoTransform();
685 : }
686 : else
687 : {
688 0 : m_oSRS = oLL;
689 : }
690 :
691 0 : m_oGCPSRS = std::move(oLL);
692 : }
693 :
694 1 : delete hkvEllipsoids;
695 : }
696 :
697 : /************************************************************************/
698 : /* Open() */
699 : /************************************************************************/
700 :
701 34417 : GDALDataset *HKVDataset::Open(GDALOpenInfo *poOpenInfo)
702 :
703 : {
704 : /* -------------------------------------------------------------------- */
705 : /* We assume the dataset is passed as a directory. Check for */
706 : /* an attrib and blob file as a minimum. */
707 : /* -------------------------------------------------------------------- */
708 34417 : if (!poOpenInfo->bIsDirectory)
709 34060 : return nullptr;
710 :
711 : std::string osFilename =
712 714 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "image_data", nullptr);
713 : VSIStatBuf sStat;
714 357 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
715 : osFilename =
716 356 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "blob", nullptr);
717 357 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
718 356 : return nullptr;
719 :
720 : osFilename =
721 1 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "attrib", nullptr);
722 1 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
723 0 : return nullptr;
724 :
725 : /* -------------------------------------------------------------------- */
726 : /* Load the attrib file, and boil white space away from around */
727 : /* the equal sign. */
728 : /* -------------------------------------------------------------------- */
729 1 : char **papszAttrib = CSLLoad(osFilename.c_str());
730 1 : if (papszAttrib == nullptr)
731 0 : return nullptr;
732 :
733 10 : for (int i = 0; papszAttrib[i] != nullptr; i++)
734 : {
735 9 : int iDst = 0;
736 9 : char *pszLine = papszAttrib[i];
737 :
738 252 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
739 : {
740 243 : if (pszLine[iSrc] != ' ')
741 : {
742 211 : pszLine[iDst++] = pszLine[iSrc];
743 : }
744 : }
745 9 : pszLine[iDst] = '\0';
746 : }
747 :
748 : /* -------------------------------------------------------------------- */
749 : /* Create a corresponding GDALDataset. */
750 : /* -------------------------------------------------------------------- */
751 2 : auto poDS = std::make_unique<HKVDataset>();
752 :
753 1 : poDS->pszPath = CPLStrdup(poOpenInfo->pszFilename);
754 1 : poDS->papszAttrib = papszAttrib;
755 :
756 1 : poDS->eAccess = poOpenInfo->eAccess;
757 :
758 : /* -------------------------------------------------------------------- */
759 : /* Set some dataset wide information. */
760 : /* -------------------------------------------------------------------- */
761 1 : bool bNative = false;
762 1 : bool bComplex = false;
763 1 : int nRawBands = 0;
764 :
765 2 : if (CSLFetchNameValue(papszAttrib, "extent.cols") == nullptr ||
766 1 : CSLFetchNameValue(papszAttrib, "extent.rows") == nullptr)
767 : {
768 0 : return nullptr;
769 : }
770 :
771 1 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib, "extent.cols"));
772 1 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib, "extent.rows"));
773 :
774 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
775 : {
776 0 : return nullptr;
777 : }
778 :
779 1 : const char *pszValue = CSLFetchNameValue(papszAttrib, "pixel.order");
780 1 : if (pszValue == nullptr)
781 0 : bNative = true;
782 : else
783 : {
784 : #ifdef CPL_MSB
785 : bNative = strstr(pszValue, "*msbf") != NULL;
786 : #else
787 1 : bNative = strstr(pszValue, "*lsbf") != nullptr;
788 : #endif
789 : }
790 :
791 1 : bool bNoDataSet = false;
792 1 : double dfNoDataValue = 0.0;
793 1 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.no_data");
794 1 : if (pszValue != nullptr)
795 : {
796 0 : bNoDataSet = true;
797 0 : dfNoDataValue = CPLAtof(pszValue);
798 : }
799 :
800 1 : pszValue = CSLFetchNameValue(papszAttrib, "channel.enumeration");
801 1 : if (pszValue != nullptr)
802 1 : nRawBands = atoi(pszValue);
803 : else
804 0 : nRawBands = 1;
805 :
806 1 : if (!GDALCheckBandCount(nRawBands, TRUE))
807 : {
808 0 : return nullptr;
809 : }
810 :
811 1 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.field");
812 1 : if (pszValue != nullptr && strstr(pszValue, "*complex") != nullptr)
813 0 : bComplex = true;
814 : else
815 1 : bComplex = false;
816 :
817 : /* Get the version number, if present (if not, assume old version. */
818 : /* Versions differ in their interpretation of corner coordinates. */
819 :
820 1 : if (CSLFetchNameValue(papszAttrib, "version") != nullptr)
821 1 : poDS->SetVersion(static_cast<float>(
822 1 : CPLAtof(CSLFetchNameValue(papszAttrib, "version"))));
823 : else
824 0 : poDS->SetVersion(1.0);
825 :
826 : /* -------------------------------------------------------------------- */
827 : /* Figure out the datatype */
828 : /* -------------------------------------------------------------------- */
829 1 : const char *pszEncoding = CSLFetchNameValue(papszAttrib, "pixel.encoding");
830 1 : if (pszEncoding == nullptr)
831 0 : pszEncoding = "{ *unsigned }";
832 :
833 1 : int nSize = 1;
834 1 : if (CSLFetchNameValue(papszAttrib, "pixel.size") != nullptr)
835 1 : nSize = atoi(CSLFetchNameValue(papszAttrib, "pixel.size")) / 8;
836 : #if 0
837 : int nPseudoBands;
838 : if( bComplex )
839 : nPseudoBands = 2;
840 : else
841 : nPseudoBands = 1;
842 : #endif
843 :
844 : GDALDataType eType;
845 1 : if (nSize == 1)
846 1 : eType = GDT_Byte;
847 0 : else if (nSize == 2 && strstr(pszEncoding, "*unsigned") != nullptr)
848 0 : eType = GDT_UInt16;
849 0 : else if (nSize == 4 && bComplex)
850 0 : eType = GDT_CInt16;
851 0 : else if (nSize == 2)
852 0 : eType = GDT_Int16;
853 0 : else if (nSize == 4 && strstr(pszEncoding, "*unsigned") != nullptr)
854 0 : eType = GDT_UInt32;
855 0 : else if (nSize == 8 && strstr(pszEncoding, "*two") != nullptr && bComplex)
856 0 : eType = GDT_CInt32;
857 0 : else if (nSize == 4 && strstr(pszEncoding, "*two") != nullptr)
858 0 : eType = GDT_Int32;
859 0 : else if (nSize == 8 && bComplex)
860 0 : eType = GDT_CFloat32;
861 0 : else if (nSize == 4)
862 0 : eType = GDT_Float32;
863 0 : else if (nSize == 16 && bComplex)
864 0 : eType = GDT_CFloat64;
865 0 : else if (nSize == 8)
866 0 : eType = GDT_Float64;
867 : else
868 : {
869 0 : CPLError(CE_Failure, CPLE_AppDefined,
870 : "Unsupported pixel data type in %s.\n"
871 : "pixel.size=%d pixel.encoding=%s",
872 0 : poDS->pszPath, nSize, pszEncoding);
873 0 : return nullptr;
874 : }
875 :
876 : /* -------------------------------------------------------------------- */
877 : /* Open the blob file. */
878 : /* -------------------------------------------------------------------- */
879 1 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "image_data", nullptr);
880 1 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
881 0 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "blob", nullptr);
882 1 : if (poOpenInfo->eAccess == GA_ReadOnly)
883 : {
884 1 : poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb");
885 1 : if (poDS->fpBlob == nullptr)
886 : {
887 0 : CPLError(CE_Failure, CPLE_OpenFailed,
888 : "Unable to open file %s for read access.",
889 : osFilename.c_str());
890 0 : return nullptr;
891 : }
892 : }
893 : else
894 : {
895 0 : poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb+");
896 0 : if (poDS->fpBlob == nullptr)
897 : {
898 0 : CPLError(CE_Failure, CPLE_OpenFailed,
899 : "Unable to open file %s for update access.",
900 : osFilename.c_str());
901 0 : return nullptr;
902 : }
903 : }
904 :
905 : /* -------------------------------------------------------------------- */
906 : /* Build the overview filename, as blob file = "_ovr". */
907 : /* -------------------------------------------------------------------- */
908 2 : std::string osOvrFilename(osFilename);
909 1 : osOvrFilename += "_ovr";
910 :
911 : /* -------------------------------------------------------------------- */
912 : /* Define the bands. */
913 : /* -------------------------------------------------------------------- */
914 1 : const int nPixelOffset = nRawBands * nSize;
915 1 : const int nLineOffset = nPixelOffset * poDS->GetRasterXSize();
916 1 : int nOffset = 0;
917 :
918 2 : for (int iRawBand = 0; iRawBand < nRawBands; iRawBand++)
919 : {
920 : auto poBand = std::make_unique<HKVRasterBand>(
921 1 : poDS.get(), poDS->GetRasterCount() + 1, poDS->fpBlob, nOffset,
922 1 : nPixelOffset, nLineOffset, eType, bNative);
923 1 : if (!poBand->IsValid())
924 0 : return nullptr;
925 :
926 1 : if (bNoDataSet)
927 0 : poBand->SetNoDataValue(dfNoDataValue);
928 1 : poDS->SetBand(poDS->GetRasterCount() + 1, std::move(poBand));
929 1 : nOffset += GDALGetDataTypeSizeBytes(eType);
930 : }
931 :
932 1 : poDS->eRasterType = eType;
933 :
934 : /* -------------------------------------------------------------------- */
935 : /* Process the georef file if there is one. */
936 : /* -------------------------------------------------------------------- */
937 1 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "georef", nullptr);
938 1 : if (VSIStat(osFilename.c_str(), &sStat) == 0)
939 1 : poDS->ProcessGeoref(osFilename.c_str());
940 :
941 : /* -------------------------------------------------------------------- */
942 : /* Initialize any PAM information. */
943 : /* -------------------------------------------------------------------- */
944 1 : poDS->SetDescription(osOvrFilename.c_str());
945 1 : poDS->TryLoadXML();
946 :
947 : /* -------------------------------------------------------------------- */
948 : /* Handle overviews. */
949 : /* -------------------------------------------------------------------- */
950 1 : poDS->oOvManager.Initialize(poDS.get(), osOvrFilename.c_str(), nullptr,
951 : TRUE);
952 :
953 1 : return poDS.release();
954 : }
955 :
956 : /************************************************************************/
957 : /* GDALRegister_HKV() */
958 : /************************************************************************/
959 :
960 2038 : void GDALRegister_HKV()
961 :
962 : {
963 2038 : if (GDALGetDriverByName("MFF2") != nullptr)
964 283 : return;
965 :
966 1755 : GDALDriver *poDriver = new GDALDriver();
967 :
968 1755 : poDriver->SetDescription("MFF2");
969 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
970 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF2 (HKV) Raster");
971 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff2.html");
972 :
973 1755 : poDriver->pfnOpen = HKVDataset::Open;
974 :
975 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
976 : }
|
