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