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