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 : 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 2 : ~HKVRasterBand() override
44 1 : {
45 2 : }
46 : };
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 : double adfGeoTransform[6];
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(double *) 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 : adfGeoTransform[0] = 0.0;
309 1 : adfGeoTransform[1] = 1.0;
310 1 : adfGeoTransform[2] = 0.0;
311 1 : adfGeoTransform[3] = 0.0;
312 1 : adfGeoTransform[4] = 0.0;
313 1 : adfGeoTransform[5] = 1.0;
314 1 : }
315 :
316 : /************************************************************************/
317 : /* ~HKVDataset() */
318 : /************************************************************************/
319 :
320 2 : HKVDataset::~HKVDataset()
321 :
322 : {
323 1 : HKVDataset::Close();
324 2 : }
325 :
326 : /************************************************************************/
327 : /* Close() */
328 : /************************************************************************/
329 :
330 2 : CPLErr HKVDataset::Close()
331 : {
332 2 : CPLErr eErr = CE_None;
333 2 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
334 : {
335 1 : if (HKVDataset::FlushCache(true) != CE_None)
336 0 : eErr = CE_Failure;
337 :
338 1 : if (fpBlob)
339 : {
340 1 : if (VSIFCloseL(fpBlob) != 0)
341 : {
342 0 : eErr = CE_Failure;
343 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
344 : }
345 : }
346 :
347 1 : if (nGCPCount > 0)
348 : {
349 1 : GDALDeinitGCPs(nGCPCount, pasGCPList);
350 1 : CPLFree(pasGCPList);
351 : }
352 :
353 1 : CPLFree(pszPath);
354 1 : CSLDestroy(papszGeoref);
355 1 : CSLDestroy(papszAttrib);
356 :
357 1 : if (GDALPamDataset::Close() != CE_None)
358 0 : eErr = CE_Failure;
359 : }
360 2 : return eErr;
361 : }
362 :
363 : /************************************************************************/
364 : /* GetGeoTransform() */
365 : /************************************************************************/
366 :
367 0 : CPLErr HKVDataset::GetGeoTransform(double *padfTransform)
368 :
369 : {
370 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
371 0 : return CE_None;
372 : }
373 :
374 : /************************************************************************/
375 : /* GetGCP() */
376 : /************************************************************************/
377 :
378 0 : const GDAL_GCP *HKVDataset::GetGCPs()
379 :
380 : {
381 0 : return pasGCPList;
382 : }
383 :
384 : /************************************************************************/
385 : /* ProcessGeorefGCP() */
386 : /************************************************************************/
387 :
388 5 : void HKVDataset::ProcessGeorefGCP(char **papszGeorefIn, const char *pszBase,
389 : double dfRasterX, double dfRasterY)
390 :
391 : {
392 : /* -------------------------------------------------------------------- */
393 : /* Fetch the GCP from the string list. */
394 : /* -------------------------------------------------------------------- */
395 5 : char szFieldName[128] = {'\0'};
396 5 : snprintf(szFieldName, sizeof(szFieldName), "%s.latitude", pszBase);
397 5 : double dfLat = 0.0;
398 5 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
399 0 : return;
400 : else
401 5 : dfLat = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
402 :
403 5 : snprintf(szFieldName, sizeof(szFieldName), "%s.longitude", pszBase);
404 5 : double dfLong = 0.0;
405 5 : if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
406 0 : return;
407 : else
408 5 : dfLong = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
409 :
410 : /* -------------------------------------------------------------------- */
411 : /* Add the gcp to the internal list. */
412 : /* -------------------------------------------------------------------- */
413 5 : GDALInitGCPs(1, pasGCPList + nGCPCount);
414 :
415 5 : CPLFree(pasGCPList[nGCPCount].pszId);
416 :
417 5 : pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
418 :
419 5 : pasGCPList[nGCPCount].dfGCPX = dfLong;
420 5 : pasGCPList[nGCPCount].dfGCPY = dfLat;
421 5 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
422 :
423 5 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
424 5 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
425 :
426 5 : nGCPCount++;
427 : }
428 :
429 : /************************************************************************/
430 : /* ProcessGeoref() */
431 : /************************************************************************/
432 :
433 1 : void HKVDataset::ProcessGeoref(const char *pszFilename)
434 :
435 : {
436 : /* -------------------------------------------------------------------- */
437 : /* Load the georef file, and boil white space away from around */
438 : /* the equal sign. */
439 : /* -------------------------------------------------------------------- */
440 1 : CSLDestroy(papszGeoref);
441 1 : papszGeoref = CSLLoad(pszFilename);
442 1 : if (papszGeoref == nullptr)
443 0 : return;
444 :
445 1 : HKVSpheroidList *hkvEllipsoids = new HKVSpheroidList;
446 :
447 14 : for (int i = 0; papszGeoref[i] != nullptr; i++)
448 : {
449 13 : int iDst = 0;
450 13 : char *pszLine = papszGeoref[i];
451 :
452 433 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
453 : {
454 420 : if (pszLine[iSrc] != ' ')
455 : {
456 420 : pszLine[iDst++] = pszLine[iSrc];
457 : }
458 : }
459 13 : pszLine[iDst] = '\0';
460 : }
461 :
462 : /* -------------------------------------------------------------------- */
463 : /* Try to get GCPs, in lat/longs . */
464 : /* -------------------------------------------------------------------- */
465 1 : nGCPCount = 0;
466 1 : pasGCPList = reinterpret_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), 5));
467 :
468 1 : if (MFF2version > 1.0)
469 : {
470 1 : ProcessGeorefGCP(papszGeoref, "top_left", 0, 0);
471 1 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize(), 0);
472 1 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0, GetRasterYSize());
473 1 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize(),
474 1 : GetRasterYSize());
475 1 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
476 1 : GetRasterYSize() / 2.0);
477 : }
478 : else
479 : {
480 0 : ProcessGeorefGCP(papszGeoref, "top_left", 0.5, 0.5);
481 0 : ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize() - 0.5, 0.5);
482 0 : ProcessGeorefGCP(papszGeoref, "bottom_left", 0.5,
483 0 : GetRasterYSize() - 0.5);
484 0 : ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize() - 0.5,
485 0 : GetRasterYSize() - 0.5);
486 0 : ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
487 0 : GetRasterYSize() / 2.0);
488 : }
489 :
490 1 : if (nGCPCount == 0)
491 : {
492 0 : CPLFree(pasGCPList);
493 0 : pasGCPList = nullptr;
494 : }
495 :
496 : /* -------------------------------------------------------------------- */
497 : /* Do we have a recognised projection? */
498 : /* -------------------------------------------------------------------- */
499 1 : const char *pszProjName = CSLFetchNameValue(papszGeoref, "projection.name");
500 : const char *pszOriginLong =
501 1 : CSLFetchNameValue(papszGeoref, "projection.origin_longitude");
502 : const char *pszSpheroidName =
503 1 : CSLFetchNameValue(papszGeoref, "spheroid.name");
504 :
505 2 : if (pszSpheroidName != nullptr &&
506 1 : hkvEllipsoids->SpheroidInList(pszSpheroidName))
507 : {
508 : #if 0
509 : // TODO(schwehr): Enable in trunk after 2.1 branch and fix.
510 : // Breaks tests on some platforms.
511 : CPLError( CE_Failure, CPLE_AppDefined,
512 : "Unrecognized ellipsoid. Not handled. "
513 : "Spheroid name not in spheroid list: '%s'",
514 : pszSpheroidName );
515 : #endif
516 : // Why were eq_radius and inv_flattening never used?
517 : // eq_radius = hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
518 : // inv_flattening =
519 : // hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
520 : }
521 0 : else if (pszProjName != nullptr)
522 : {
523 0 : CPLError(CE_Warning, CPLE_AppDefined,
524 : "Unrecognized ellipsoid. Not handled.");
525 : // TODO(schwehr): This error is was never what was happening.
526 : // CPLError( CE_Warning, CPLE_AppDefined,
527 : // "Unrecognized ellipsoid. Using wgs-84 parameters.");
528 : // eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84"); */
529 : // inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
530 : }
531 :
532 1 : if (pszProjName != nullptr && EQUAL(pszProjName, "utm") && nGCPCount == 5)
533 : {
534 : // int nZone = (int)((CPLAtof(pszOriginLong)+184.5) / 6.0);
535 1 : int nZone = 31; // TODO(schwehr): Where does 31 come from?
536 :
537 1 : if (pszOriginLong == nullptr)
538 : {
539 : // If origin not specified, assume 0.0.
540 0 : CPLError(
541 : CE_Warning, CPLE_AppDefined,
542 : "No projection origin longitude specified. Assuming 0.0.");
543 : }
544 : else
545 : {
546 1 : nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
547 : }
548 :
549 2 : OGRSpatialReference oUTM;
550 :
551 1 : if (pasGCPList[4].dfGCPY < 0)
552 0 : oUTM.SetUTM(nZone, 0);
553 : else
554 1 : oUTM.SetUTM(nZone, 1);
555 :
556 2 : OGRSpatialReference oLL;
557 1 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
558 1 : if (pszOriginLong != nullptr)
559 : {
560 1 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
561 1 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
562 : }
563 :
564 1 : if ((pszSpheroidName == nullptr) ||
565 1 : (EQUAL(pszSpheroidName, "wgs-84")) ||
566 1 : (EQUAL(pszSpheroidName, "wgs_84")))
567 : {
568 0 : oUTM.SetWellKnownGeogCS("WGS84");
569 0 : oLL.SetWellKnownGeogCS("WGS84");
570 : }
571 : else
572 : {
573 1 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
574 : {
575 1 : oUTM.SetGeogCS(
576 : "unknown", "unknown", pszSpheroidName,
577 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
578 : hkvEllipsoids->GetSpheroidInverseFlattening(
579 : pszSpheroidName));
580 1 : oLL.SetGeogCS(
581 : "unknown", "unknown", pszSpheroidName,
582 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
583 : hkvEllipsoids->GetSpheroidInverseFlattening(
584 : pszSpheroidName));
585 : }
586 : else
587 : {
588 0 : CPLError(CE_Warning, CPLE_AppDefined,
589 : "Unrecognized ellipsoid. Using wgs-84 parameters.");
590 0 : oUTM.SetWellKnownGeogCS("WGS84");
591 0 : oLL.SetWellKnownGeogCS("WGS84");
592 : }
593 : }
594 :
595 : OGRCoordinateTransformation *poTransform =
596 1 : OGRCreateCoordinateTransformation(&oLL, &oUTM);
597 :
598 1 : bool bSuccess = true;
599 1 : if (poTransform == nullptr)
600 : {
601 0 : CPLErrorReset();
602 0 : bSuccess = false;
603 : }
604 :
605 1 : double dfUtmX[5] = {0.0};
606 1 : double dfUtmY[5] = {0.0};
607 :
608 1 : if (poTransform != nullptr)
609 : {
610 6 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
611 : {
612 5 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
613 5 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
614 :
615 5 : if (bSuccess && !poTransform->Transform(1, &(dfUtmX[gcp_index]),
616 : &(dfUtmY[gcp_index])))
617 0 : bSuccess = false;
618 : }
619 : }
620 :
621 1 : if (bSuccess)
622 : {
623 : // Update GCPS to proper projection.
624 6 : for (int gcp_index = 0; gcp_index < 5; gcp_index++)
625 : {
626 5 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
627 5 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
628 : }
629 :
630 1 : m_oGCPSRS = oUTM;
631 :
632 1 : bool transform_ok = CPL_TO_BOOL(
633 1 : GDALGCPsToGeoTransform(5, pasGCPList, adfGeoTransform, 0));
634 :
635 1 : if (!transform_ok)
636 : {
637 : // Transform may not be sufficient in all cases (slant range
638 : // projection).
639 0 : adfGeoTransform[0] = 0.0;
640 0 : adfGeoTransform[1] = 1.0;
641 0 : adfGeoTransform[2] = 0.0;
642 0 : adfGeoTransform[3] = 0.0;
643 0 : adfGeoTransform[4] = 0.0;
644 0 : adfGeoTransform[5] = 1.0;
645 0 : m_oGCPSRS.Clear();
646 : }
647 : else
648 : {
649 1 : m_oSRS = std::move(oUTM);
650 : }
651 : }
652 :
653 1 : if (poTransform != nullptr)
654 2 : delete poTransform;
655 : }
656 0 : else if (pszProjName != nullptr && nGCPCount == 5)
657 : {
658 0 : OGRSpatialReference oLL;
659 0 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
660 :
661 0 : if (pszOriginLong != nullptr)
662 : {
663 0 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
664 : }
665 :
666 0 : if (pszSpheroidName == nullptr ||
667 0 : EQUAL(pszSpheroidName, "wgs-84") || // Dash.
668 0 : EQUAL(pszSpheroidName, "wgs_84")) // Underscore.
669 : {
670 0 : oLL.SetWellKnownGeogCS("WGS84");
671 : }
672 : else
673 : {
674 0 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
675 : {
676 0 : oLL.SetGeogCS(
677 : "", "", pszSpheroidName,
678 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
679 : hkvEllipsoids->GetSpheroidInverseFlattening(
680 : pszSpheroidName));
681 : }
682 : else
683 : {
684 0 : CPLError(CE_Warning, CPLE_AppDefined,
685 : "Unrecognized ellipsoid. "
686 : "Using wgs-84 parameters.");
687 0 : oLL.SetWellKnownGeogCS("WGS84");
688 : }
689 : }
690 :
691 0 : const bool transform_ok = CPL_TO_BOOL(
692 0 : GDALGCPsToGeoTransform(5, pasGCPList, adfGeoTransform, 0));
693 :
694 0 : m_oSRS.Clear();
695 :
696 0 : if (!transform_ok)
697 : {
698 0 : adfGeoTransform[0] = 0.0;
699 0 : adfGeoTransform[1] = 1.0;
700 0 : adfGeoTransform[2] = 0.0;
701 0 : adfGeoTransform[3] = 0.0;
702 0 : adfGeoTransform[4] = 0.0;
703 0 : adfGeoTransform[5] = 1.0;
704 : }
705 : else
706 : {
707 0 : m_oSRS = oLL;
708 : }
709 :
710 0 : m_oGCPSRS = std::move(oLL);
711 : }
712 :
713 1 : delete hkvEllipsoids;
714 : }
715 :
716 : /************************************************************************/
717 : /* Open() */
718 : /************************************************************************/
719 :
720 29531 : GDALDataset *HKVDataset::Open(GDALOpenInfo *poOpenInfo)
721 :
722 : {
723 : /* -------------------------------------------------------------------- */
724 : /* We assume the dataset is passed as a directory. Check for */
725 : /* an attrib and blob file as a minimum. */
726 : /* -------------------------------------------------------------------- */
727 29531 : if (!poOpenInfo->bIsDirectory)
728 29297 : return nullptr;
729 :
730 : std::string osFilename =
731 468 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "image_data", nullptr);
732 : VSIStatBuf sStat;
733 234 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
734 : osFilename =
735 233 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "blob", nullptr);
736 234 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
737 233 : return nullptr;
738 :
739 : osFilename =
740 1 : CPLFormFilenameSafe(poOpenInfo->pszFilename, "attrib", nullptr);
741 1 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
742 0 : return nullptr;
743 :
744 : /* -------------------------------------------------------------------- */
745 : /* Load the attrib file, and boil white space away from around */
746 : /* the equal sign. */
747 : /* -------------------------------------------------------------------- */
748 1 : char **papszAttrib = CSLLoad(osFilename.c_str());
749 1 : if (papszAttrib == nullptr)
750 0 : return nullptr;
751 :
752 10 : for (int i = 0; papszAttrib[i] != nullptr; i++)
753 : {
754 9 : int iDst = 0;
755 9 : char *pszLine = papszAttrib[i];
756 :
757 252 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
758 : {
759 243 : if (pszLine[iSrc] != ' ')
760 : {
761 211 : pszLine[iDst++] = pszLine[iSrc];
762 : }
763 : }
764 9 : pszLine[iDst] = '\0';
765 : }
766 :
767 : /* -------------------------------------------------------------------- */
768 : /* Create a corresponding GDALDataset. */
769 : /* -------------------------------------------------------------------- */
770 2 : auto poDS = std::make_unique<HKVDataset>();
771 :
772 1 : poDS->pszPath = CPLStrdup(poOpenInfo->pszFilename);
773 1 : poDS->papszAttrib = papszAttrib;
774 :
775 1 : poDS->eAccess = poOpenInfo->eAccess;
776 :
777 : /* -------------------------------------------------------------------- */
778 : /* Set some dataset wide information. */
779 : /* -------------------------------------------------------------------- */
780 1 : bool bNative = false;
781 1 : bool bComplex = false;
782 1 : int nRawBands = 0;
783 :
784 2 : if (CSLFetchNameValue(papszAttrib, "extent.cols") == nullptr ||
785 1 : CSLFetchNameValue(papszAttrib, "extent.rows") == nullptr)
786 : {
787 0 : return nullptr;
788 : }
789 :
790 1 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib, "extent.cols"));
791 1 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib, "extent.rows"));
792 :
793 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
794 : {
795 0 : return nullptr;
796 : }
797 :
798 1 : const char *pszValue = CSLFetchNameValue(papszAttrib, "pixel.order");
799 1 : if (pszValue == nullptr)
800 0 : bNative = true;
801 : else
802 : {
803 : #ifdef CPL_MSB
804 : bNative = strstr(pszValue, "*msbf") != NULL;
805 : #else
806 1 : bNative = strstr(pszValue, "*lsbf") != nullptr;
807 : #endif
808 : }
809 :
810 1 : bool bNoDataSet = false;
811 1 : double dfNoDataValue = 0.0;
812 1 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.no_data");
813 1 : if (pszValue != nullptr)
814 : {
815 0 : bNoDataSet = true;
816 0 : dfNoDataValue = CPLAtof(pszValue);
817 : }
818 :
819 1 : pszValue = CSLFetchNameValue(papszAttrib, "channel.enumeration");
820 1 : if (pszValue != nullptr)
821 1 : nRawBands = atoi(pszValue);
822 : else
823 0 : nRawBands = 1;
824 :
825 1 : if (!GDALCheckBandCount(nRawBands, TRUE))
826 : {
827 0 : return nullptr;
828 : }
829 :
830 1 : pszValue = CSLFetchNameValue(papszAttrib, "pixel.field");
831 1 : if (pszValue != nullptr && strstr(pszValue, "*complex") != nullptr)
832 0 : bComplex = true;
833 : else
834 1 : bComplex = false;
835 :
836 : /* Get the version number, if present (if not, assume old version. */
837 : /* Versions differ in their interpretation of corner coordinates. */
838 :
839 1 : if (CSLFetchNameValue(papszAttrib, "version") != nullptr)
840 1 : poDS->SetVersion(static_cast<float>(
841 1 : CPLAtof(CSLFetchNameValue(papszAttrib, "version"))));
842 : else
843 0 : poDS->SetVersion(1.0);
844 :
845 : /* -------------------------------------------------------------------- */
846 : /* Figure out the datatype */
847 : /* -------------------------------------------------------------------- */
848 1 : const char *pszEncoding = CSLFetchNameValue(papszAttrib, "pixel.encoding");
849 1 : if (pszEncoding == nullptr)
850 0 : pszEncoding = "{ *unsigned }";
851 :
852 1 : int nSize = 1;
853 1 : if (CSLFetchNameValue(papszAttrib, "pixel.size") != nullptr)
854 1 : nSize = atoi(CSLFetchNameValue(papszAttrib, "pixel.size")) / 8;
855 : #if 0
856 : int nPseudoBands;
857 : if( bComplex )
858 : nPseudoBands = 2;
859 : else
860 : nPseudoBands = 1;
861 : #endif
862 :
863 : GDALDataType eType;
864 1 : if (nSize == 1)
865 1 : eType = GDT_Byte;
866 0 : else if (nSize == 2 && strstr(pszEncoding, "*unsigned") != nullptr)
867 0 : eType = GDT_UInt16;
868 0 : else if (nSize == 4 && bComplex)
869 0 : eType = GDT_CInt16;
870 0 : else if (nSize == 2)
871 0 : eType = GDT_Int16;
872 0 : else if (nSize == 4 && strstr(pszEncoding, "*unsigned") != nullptr)
873 0 : eType = GDT_UInt32;
874 0 : else if (nSize == 8 && strstr(pszEncoding, "*two") != nullptr && bComplex)
875 0 : eType = GDT_CInt32;
876 0 : else if (nSize == 4 && strstr(pszEncoding, "*two") != nullptr)
877 0 : eType = GDT_Int32;
878 0 : else if (nSize == 8 && bComplex)
879 0 : eType = GDT_CFloat32;
880 0 : else if (nSize == 4)
881 0 : eType = GDT_Float32;
882 0 : else if (nSize == 16 && bComplex)
883 0 : eType = GDT_CFloat64;
884 0 : else if (nSize == 8)
885 0 : eType = GDT_Float64;
886 : else
887 : {
888 0 : CPLError(CE_Failure, CPLE_AppDefined,
889 : "Unsupported pixel data type in %s.\n"
890 : "pixel.size=%d pixel.encoding=%s",
891 0 : poDS->pszPath, nSize, pszEncoding);
892 0 : return nullptr;
893 : }
894 :
895 : /* -------------------------------------------------------------------- */
896 : /* Open the blob file. */
897 : /* -------------------------------------------------------------------- */
898 1 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "image_data", nullptr);
899 1 : if (VSIStat(osFilename.c_str(), &sStat) != 0)
900 0 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "blob", nullptr);
901 1 : if (poOpenInfo->eAccess == GA_ReadOnly)
902 : {
903 1 : poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb");
904 1 : if (poDS->fpBlob == nullptr)
905 : {
906 0 : CPLError(CE_Failure, CPLE_OpenFailed,
907 : "Unable to open file %s for read access.",
908 : osFilename.c_str());
909 0 : return nullptr;
910 : }
911 : }
912 : else
913 : {
914 0 : poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb+");
915 0 : if (poDS->fpBlob == nullptr)
916 : {
917 0 : CPLError(CE_Failure, CPLE_OpenFailed,
918 : "Unable to open file %s for update access.",
919 : osFilename.c_str());
920 0 : return nullptr;
921 : }
922 : }
923 :
924 : /* -------------------------------------------------------------------- */
925 : /* Build the overview filename, as blob file = "_ovr". */
926 : /* -------------------------------------------------------------------- */
927 2 : std::string osOvrFilename(osFilename);
928 1 : osOvrFilename += "_ovr";
929 :
930 : /* -------------------------------------------------------------------- */
931 : /* Define the bands. */
932 : /* -------------------------------------------------------------------- */
933 1 : const int nPixelOffset = nRawBands * nSize;
934 1 : const int nLineOffset = nPixelOffset * poDS->GetRasterXSize();
935 1 : int nOffset = 0;
936 :
937 2 : for (int iRawBand = 0; iRawBand < nRawBands; iRawBand++)
938 : {
939 : auto poBand = std::make_unique<HKVRasterBand>(
940 1 : poDS.get(), poDS->GetRasterCount() + 1, poDS->fpBlob, nOffset,
941 1 : nPixelOffset, nLineOffset, eType, bNative);
942 1 : if (!poBand->IsValid())
943 0 : return nullptr;
944 :
945 1 : if (bNoDataSet)
946 0 : poBand->SetNoDataValue(dfNoDataValue);
947 1 : poDS->SetBand(poDS->GetRasterCount() + 1, std::move(poBand));
948 1 : nOffset += GDALGetDataTypeSizeBytes(eType);
949 : }
950 :
951 1 : poDS->eRasterType = eType;
952 :
953 : /* -------------------------------------------------------------------- */
954 : /* Process the georef file if there is one. */
955 : /* -------------------------------------------------------------------- */
956 1 : osFilename = CPLFormFilenameSafe(poDS->pszPath, "georef", nullptr);
957 1 : if (VSIStat(osFilename.c_str(), &sStat) == 0)
958 1 : poDS->ProcessGeoref(osFilename.c_str());
959 :
960 : /* -------------------------------------------------------------------- */
961 : /* Initialize any PAM information. */
962 : /* -------------------------------------------------------------------- */
963 1 : poDS->SetDescription(osOvrFilename.c_str());
964 1 : poDS->TryLoadXML();
965 :
966 : /* -------------------------------------------------------------------- */
967 : /* Handle overviews. */
968 : /* -------------------------------------------------------------------- */
969 1 : poDS->oOvManager.Initialize(poDS.get(), osOvrFilename.c_str(), nullptr,
970 : TRUE);
971 :
972 1 : return poDS.release();
973 : }
974 :
975 : /************************************************************************/
976 : /* GDALRegister_HKV() */
977 : /************************************************************************/
978 :
979 1667 : void GDALRegister_HKV()
980 :
981 : {
982 1667 : if (GDALGetDriverByName("MFF2") != nullptr)
983 282 : return;
984 :
985 1385 : GDALDriver *poDriver = new GDALDriver();
986 :
987 1385 : poDriver->SetDescription("MFF2");
988 1385 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
989 1385 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF2 (HKV) Raster");
990 1385 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff2.html");
991 :
992 1385 : poDriver->pfnOpen = HKVDataset::Open;
993 :
994 1385 : GetGDALDriverManager()->RegisterDriver(poDriver);
995 : }
|
