Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GView
4 : * Purpose: Implementation of Atlantis MFF Support
5 : * Author: Frank Warmerdam, warmerda@home.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 "atlsci_spheroid.h"
15 : #include "cpl_string.h"
16 : #include "gdal_frmts.h"
17 : #include "gdal_priv.h"
18 : #include "ogr_spatialref.h"
19 : #include "rawdataset.h"
20 :
21 : #include <cctype>
22 : #include <cmath>
23 : #include <algorithm>
24 :
25 : enum
26 : {
27 : MFFPRJ_NONE,
28 : MFFPRJ_LL,
29 : MFFPRJ_UTM,
30 : MFFPRJ_UNRECOGNIZED
31 : };
32 :
33 : /************************************************************************/
34 : /* ==================================================================== */
35 : /* MFFDataset */
36 : /* ==================================================================== */
37 : /************************************************************************/
38 :
39 : class MFFDataset final : public RawDataset
40 : {
41 : int nGCPCount;
42 : GDAL_GCP *pasGCPList;
43 :
44 : OGRSpatialReference m_oSRS{};
45 : OGRSpatialReference m_oGCPSRS{};
46 : GDALGeoTransform m_gt{};
47 : char **m_papszFileList;
48 :
49 : void ScanForGCPs();
50 : void ScanForProjectionInfo();
51 :
52 : CPL_DISALLOW_COPY_ASSIGN(MFFDataset)
53 :
54 : CPLErr Close() override;
55 :
56 : public:
57 : MFFDataset();
58 : ~MFFDataset() override;
59 :
60 : char **papszHdrLines;
61 :
62 : VSILFILE **pafpBandFiles;
63 :
64 : char **GetFileList() override;
65 :
66 : int GetGCPCount() override;
67 :
68 0 : const OGRSpatialReference *GetGCPSpatialRef() const override
69 : {
70 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
71 : }
72 :
73 : const GDAL_GCP *GetGCPs() override;
74 :
75 0 : const OGRSpatialReference *GetSpatialRef() const override
76 : {
77 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
78 : }
79 :
80 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
81 :
82 : static GDALDataset *Open(GDALOpenInfo *);
83 : };
84 :
85 : /************************************************************************/
86 : /* ==================================================================== */
87 : /* MFFTiledBand */
88 : /* ==================================================================== */
89 : /************************************************************************/
90 :
91 : class MFFTiledBand final : public GDALRasterBand
92 : {
93 : friend class MFFDataset;
94 :
95 : VSILFILE *fpRaw;
96 : RawRasterBand::ByteOrder eByteOrder;
97 :
98 : CPL_DISALLOW_COPY_ASSIGN(MFFTiledBand)
99 :
100 : public:
101 : MFFTiledBand(MFFDataset *, int, VSILFILE *, int, int, GDALDataType,
102 : RawRasterBand::ByteOrder);
103 : ~MFFTiledBand() override;
104 :
105 : CPLErr IReadBlock(int, int, void *) override;
106 : };
107 :
108 : /************************************************************************/
109 : /* MFFTiledBand() */
110 : /************************************************************************/
111 :
112 2 : MFFTiledBand::MFFTiledBand(MFFDataset *poDSIn, int nBandIn, VSILFILE *fp,
113 : int nTileXSize, int nTileYSize,
114 : GDALDataType eDataTypeIn,
115 2 : RawRasterBand::ByteOrder eByteOrderIn)
116 2 : : fpRaw(fp), eByteOrder(eByteOrderIn)
117 : {
118 2 : poDS = poDSIn;
119 2 : nBand = nBandIn;
120 :
121 2 : eDataType = eDataTypeIn;
122 :
123 2 : nBlockXSize = nTileXSize;
124 2 : nBlockYSize = nTileYSize;
125 2 : }
126 :
127 : /************************************************************************/
128 : /* ~MFFTiledBand() */
129 : /************************************************************************/
130 :
131 4 : MFFTiledBand::~MFFTiledBand()
132 :
133 : {
134 2 : if (VSIFCloseL(fpRaw) != 0)
135 : {
136 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
137 : }
138 4 : }
139 :
140 : /************************************************************************/
141 : /* IReadBlock() */
142 : /************************************************************************/
143 :
144 1 : CPLErr MFFTiledBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
145 :
146 : {
147 1 : const int nTilesPerRow = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
148 1 : const int nWordSize = GDALGetDataTypeSizeBytes(eDataType);
149 1 : const int nBlockSize = nWordSize * nBlockXSize * nBlockYSize;
150 :
151 1 : const vsi_l_offset nOffset =
152 1 : nBlockSize *
153 1 : (nBlockXOff + static_cast<vsi_l_offset>(nBlockYOff) * nTilesPerRow);
154 :
155 2 : if (VSIFSeekL(fpRaw, nOffset, SEEK_SET) == -1 ||
156 1 : VSIFReadL(pImage, 1, nBlockSize, fpRaw) < 1)
157 : {
158 0 : CPLError(CE_Failure, CPLE_FileIO,
159 : "Read of tile %d/%d failed with fseek or fread error.",
160 : nBlockXOff, nBlockYOff);
161 0 : return CE_Failure;
162 : }
163 :
164 1 : if (eByteOrder != RawRasterBand::NATIVE_BYTE_ORDER && nWordSize > 1)
165 : {
166 0 : if (GDALDataTypeIsComplex(eDataType))
167 : {
168 0 : GDALSwapWords(pImage, nWordSize / 2, nBlockXSize * nBlockYSize,
169 : nWordSize);
170 0 : GDALSwapWords(reinterpret_cast<GByte *>(pImage) + nWordSize / 2,
171 0 : nWordSize / 2, nBlockXSize * nBlockYSize, nWordSize);
172 : }
173 : else
174 0 : GDALSwapWords(pImage, nWordSize, nBlockXSize * nBlockYSize,
175 : nWordSize);
176 : }
177 :
178 1 : return CE_None;
179 : }
180 :
181 : /************************************************************************/
182 : /* MFF Spheroids */
183 : /************************************************************************/
184 :
185 : class MFFSpheroidList final : public SpheroidList
186 : {
187 : public:
188 : MFFSpheroidList();
189 : };
190 :
191 6 : MFFSpheroidList ::MFFSpheroidList()
192 : {
193 6 : num_spheroids = 18;
194 :
195 6 : epsilonR = 0.1;
196 6 : epsilonI = 0.000001;
197 :
198 6 : spheroids[0].SetValuesByRadii("SPHERE", 6371007.0, 6371007.0);
199 6 : spheroids[1].SetValuesByRadii("EVEREST", 6377304.0, 6356103.0);
200 6 : spheroids[2].SetValuesByRadii("BESSEL", 6377397.0, 6356082.0);
201 6 : spheroids[3].SetValuesByRadii("AIRY", 6377563.0, 6356300.0);
202 6 : spheroids[4].SetValuesByRadii("CLARKE_1858", 6378294.0, 6356621.0);
203 6 : spheroids[5].SetValuesByRadii("CLARKE_1866", 6378206.4, 6356583.8);
204 6 : spheroids[6].SetValuesByRadii("CLARKE_1880", 6378249.0, 6356517.0);
205 6 : spheroids[7].SetValuesByRadii("HAYFORD", 6378388.0, 6356915.0);
206 6 : spheroids[8].SetValuesByRadii("KRASOVSKI", 6378245.0, 6356863.0);
207 6 : spheroids[9].SetValuesByRadii("HOUGH", 6378270.0, 6356794.0);
208 6 : spheroids[10].SetValuesByRadii("FISHER_60", 6378166.0, 6356784.0);
209 6 : spheroids[11].SetValuesByRadii("KAULA", 6378165.0, 6356345.0);
210 6 : spheroids[12].SetValuesByRadii("IUGG_67", 6378160.0, 6356775.0);
211 6 : spheroids[13].SetValuesByRadii("FISHER_68", 6378150.0, 6356330.0);
212 6 : spheroids[14].SetValuesByRadii("WGS_72", 6378135.0, 6356751.0);
213 6 : spheroids[15].SetValuesByRadii("IUGG_75", 6378140.0, 6356755.0);
214 6 : spheroids[16].SetValuesByRadii("WGS_84", 6378137.0, 6356752.0);
215 6 : spheroids[17].SetValuesByRadii("HUGHES", 6378273.0, 6356889.4);
216 6 : }
217 :
218 : /************************************************************************/
219 : /* ==================================================================== */
220 : /* MFFDataset */
221 : /* ==================================================================== */
222 : /************************************************************************/
223 :
224 : /************************************************************************/
225 : /* MFFDataset() */
226 : /************************************************************************/
227 :
228 8 : MFFDataset::MFFDataset()
229 : : nGCPCount(0), pasGCPList(nullptr), m_papszFileList(nullptr),
230 8 : papszHdrLines(nullptr), pafpBandFiles(nullptr)
231 : {
232 8 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
233 8 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
234 8 : }
235 :
236 : /************************************************************************/
237 : /* ~MFFDataset() */
238 : /************************************************************************/
239 :
240 16 : MFFDataset::~MFFDataset()
241 :
242 : {
243 8 : MFFDataset::Close();
244 16 : }
245 :
246 : /************************************************************************/
247 : /* Close() */
248 : /************************************************************************/
249 :
250 15 : CPLErr MFFDataset::Close()
251 : {
252 15 : CPLErr eErr = CE_None;
253 15 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
254 : {
255 8 : if (MFFDataset::FlushCache(true) != CE_None)
256 0 : eErr = CE_Failure;
257 :
258 8 : CSLDestroy(papszHdrLines);
259 8 : if (pafpBandFiles)
260 : {
261 0 : for (int i = 0; i < GetRasterCount(); i++)
262 : {
263 0 : if (pafpBandFiles[i])
264 : {
265 0 : if (VSIFCloseL(pafpBandFiles[i]) != 0)
266 : {
267 0 : eErr = CE_Failure;
268 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
269 : }
270 : }
271 : }
272 0 : CPLFree(pafpBandFiles);
273 : }
274 :
275 8 : if (nGCPCount > 0)
276 : {
277 6 : GDALDeinitGCPs(nGCPCount, pasGCPList);
278 : }
279 8 : CPLFree(pasGCPList);
280 8 : CSLDestroy(m_papszFileList);
281 :
282 8 : if (GDALPamDataset::Close() != CE_None)
283 0 : eErr = CE_Failure;
284 : }
285 15 : return eErr;
286 : }
287 :
288 : /************************************************************************/
289 : /* GetFileList() */
290 : /************************************************************************/
291 :
292 3 : char **MFFDataset::GetFileList()
293 : {
294 3 : char **papszFileList = RawDataset::GetFileList();
295 3 : papszFileList = CSLInsertStrings(papszFileList, -1, m_papszFileList);
296 3 : return papszFileList;
297 : }
298 :
299 : /************************************************************************/
300 : /* GetGCPCount() */
301 : /************************************************************************/
302 :
303 0 : int MFFDataset::GetGCPCount()
304 :
305 : {
306 0 : return nGCPCount;
307 : }
308 :
309 : /************************************************************************/
310 : /* GetGeoTransform() */
311 : /************************************************************************/
312 :
313 0 : CPLErr MFFDataset::GetGeoTransform(GDALGeoTransform >) const
314 : {
315 0 : gt = m_gt;
316 0 : return CE_None;
317 : }
318 :
319 : /************************************************************************/
320 : /* GetGCP() */
321 : /************************************************************************/
322 :
323 0 : const GDAL_GCP *MFFDataset::GetGCPs()
324 :
325 : {
326 0 : return pasGCPList;
327 : }
328 :
329 : /************************************************************************/
330 : /* ScanForGCPs() */
331 : /************************************************************************/
332 :
333 7 : void MFFDataset::ScanForGCPs()
334 :
335 : {
336 7 : int NUM_GCPS = 0;
337 :
338 7 : if (CSLFetchNameValue(papszHdrLines, "NUM_GCPS") != nullptr)
339 4 : NUM_GCPS = atoi(CSLFetchNameValue(papszHdrLines, "NUM_GCPS"));
340 7 : if (NUM_GCPS < 0)
341 0 : return;
342 :
343 7 : nGCPCount = 0;
344 7 : pasGCPList =
345 7 : static_cast<GDAL_GCP *>(VSICalloc(sizeof(GDAL_GCP), 5 + NUM_GCPS));
346 7 : if (pasGCPList == nullptr)
347 0 : return;
348 :
349 42 : for (int nCorner = 0; nCorner < 5; nCorner++)
350 : {
351 35 : const char *pszBase = nullptr;
352 35 : double dfRasterX = 0.0;
353 35 : double dfRasterY = 0.0;
354 :
355 35 : if (nCorner == 0)
356 : {
357 7 : dfRasterX = 0.5;
358 7 : dfRasterY = 0.5;
359 7 : pszBase = "TOP_LEFT_CORNER";
360 : }
361 28 : else if (nCorner == 1)
362 : {
363 7 : dfRasterX = GetRasterXSize() - 0.5;
364 7 : dfRasterY = 0.5;
365 7 : pszBase = "TOP_RIGHT_CORNER";
366 : }
367 21 : else if (nCorner == 2)
368 : {
369 7 : dfRasterX = GetRasterXSize() - 0.5;
370 7 : dfRasterY = GetRasterYSize() - 0.5;
371 7 : pszBase = "BOTTOM_RIGHT_CORNER";
372 : }
373 14 : else if (nCorner == 3)
374 : {
375 7 : dfRasterX = 0.5;
376 7 : dfRasterY = GetRasterYSize() - 0.5;
377 7 : pszBase = "BOTTOM_LEFT_CORNER";
378 : }
379 : else /* if( nCorner == 4 ) */
380 : {
381 7 : dfRasterX = GetRasterXSize() / 2.0;
382 7 : dfRasterY = GetRasterYSize() / 2.0;
383 7 : pszBase = "CENTRE";
384 : }
385 :
386 35 : char szLatName[40] = {'\0'};
387 35 : char szLongName[40] = {'\0'};
388 35 : snprintf(szLatName, sizeof(szLatName), "%s_LATITUDE", pszBase);
389 35 : snprintf(szLongName, sizeof(szLongName), "%s_LONGITUDE", pszBase);
390 :
391 45 : if (CSLFetchNameValue(papszHdrLines, szLatName) != nullptr &&
392 10 : CSLFetchNameValue(papszHdrLines, szLongName) != nullptr)
393 : {
394 10 : GDALInitGCPs(1, pasGCPList + nGCPCount);
395 :
396 10 : CPLFree(pasGCPList[nGCPCount].pszId);
397 :
398 10 : pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
399 :
400 20 : pasGCPList[nGCPCount].dfGCPX =
401 10 : CPLAtof(CSLFetchNameValue(papszHdrLines, szLongName));
402 20 : pasGCPList[nGCPCount].dfGCPY =
403 10 : CPLAtof(CSLFetchNameValue(papszHdrLines, szLatName));
404 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
405 :
406 10 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
407 10 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
408 :
409 10 : nGCPCount++;
410 : }
411 : }
412 :
413 : /* -------------------------------------------------------------------- */
414 : /* Collect standalone GCPs. They look like: */
415 : /* */
416 : /* GCPn = row, col, lat, long */
417 : /* GCP1 = 1, 1, 45.0, -75.0 */
418 : /* -------------------------------------------------------------------- */
419 11 : for (int i = 0; i < NUM_GCPS; i++)
420 : {
421 4 : char szName[25] = {'\0'};
422 4 : snprintf(szName, sizeof(szName), "GCP%d", i + 1);
423 4 : if (CSLFetchNameValue(papszHdrLines, szName) == nullptr)
424 0 : continue;
425 :
426 4 : char **papszTokens = CSLTokenizeStringComplex(
427 4 : CSLFetchNameValue(papszHdrLines, szName), ",", FALSE, FALSE);
428 4 : if (CSLCount(papszTokens) == 4)
429 : {
430 4 : GDALInitGCPs(1, pasGCPList + nGCPCount);
431 :
432 4 : CPLFree(pasGCPList[nGCPCount].pszId);
433 4 : pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
434 :
435 4 : pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[3]);
436 4 : pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[2]);
437 4 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
438 4 : pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]) + 0.5;
439 4 : pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[0]) + 0.5;
440 :
441 4 : nGCPCount++;
442 : }
443 :
444 4 : CSLDestroy(papszTokens);
445 : }
446 : }
447 :
448 : /************************************************************************/
449 : /* ScanForProjectionInfo */
450 : /************************************************************************/
451 :
452 7 : void MFFDataset::ScanForProjectionInfo()
453 : {
454 : const char *pszProjName =
455 7 : CSLFetchNameValue(papszHdrLines, "PROJECTION_NAME");
456 : const char *pszOriginLong =
457 7 : CSLFetchNameValue(papszHdrLines, "PROJECTION_ORIGIN_LONGITUDE");
458 : const char *pszSpheroidName =
459 7 : CSLFetchNameValue(papszHdrLines, "SPHEROID_NAME");
460 :
461 7 : if (pszProjName == nullptr)
462 : {
463 1 : m_oSRS.Clear();
464 1 : m_oGCPSRS.Clear();
465 1 : return;
466 : }
467 6 : else if ((!EQUAL(pszProjName, "utm")) && (!EQUAL(pszProjName, "ll")))
468 : {
469 0 : CPLError(CE_Warning, CPLE_AppDefined,
470 : "Only utm and lat/long projections are currently supported.");
471 0 : m_oSRS.Clear();
472 0 : m_oGCPSRS.Clear();
473 0 : return;
474 : }
475 6 : MFFSpheroidList *mffEllipsoids = new MFFSpheroidList;
476 :
477 12 : OGRSpatialReference oProj;
478 6 : oProj.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
479 6 : if (EQUAL(pszProjName, "utm"))
480 : {
481 : int nZone;
482 :
483 6 : if (pszOriginLong == nullptr)
484 : {
485 : // If origin not specified, assume 0.0.
486 4 : CPLError(
487 : CE_Warning, CPLE_AppDefined,
488 : "No projection origin longitude specified. Assuming 0.0.");
489 4 : nZone = 31;
490 : }
491 : else
492 2 : nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
493 :
494 6 : if (nGCPCount >= 5 && pasGCPList[4].dfGCPY < 0)
495 0 : oProj.SetUTM(nZone, 0);
496 : else
497 6 : oProj.SetUTM(nZone, 1);
498 :
499 6 : if (pszOriginLong != nullptr)
500 2 : oProj.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
501 : }
502 :
503 12 : OGRSpatialReference oLL;
504 6 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
505 6 : if (pszOriginLong != nullptr)
506 2 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
507 :
508 6 : if (pszSpheroidName == nullptr)
509 : {
510 4 : CPLError(CE_Warning, CPLE_AppDefined,
511 : "Unspecified ellipsoid. Using wgs-84 parameters.\n");
512 :
513 4 : oProj.SetWellKnownGeogCS("WGS84");
514 4 : oLL.SetWellKnownGeogCS("WGS84");
515 : }
516 : else
517 : {
518 2 : if (mffEllipsoids->SpheroidInList(pszSpheroidName))
519 : {
520 2 : oProj.SetGeogCS(
521 : "unknown", "unknown", pszSpheroidName,
522 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
523 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName));
524 2 : oLL.SetGeogCS(
525 : "unknown", "unknown", pszSpheroidName,
526 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
527 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName));
528 : }
529 0 : else if (EQUAL(pszSpheroidName, "USER_DEFINED"))
530 : {
531 : const char *pszSpheroidEqRadius =
532 0 : CSLFetchNameValue(papszHdrLines, "SPHEROID_EQUATORIAL_RADIUS");
533 : const char *pszSpheroidPolarRadius =
534 0 : CSLFetchNameValue(papszHdrLines, "SPHEROID_POLAR_RADIUS");
535 0 : if ((pszSpheroidEqRadius != nullptr) &&
536 : (pszSpheroidPolarRadius != nullptr))
537 : {
538 0 : const double eq_radius = CPLAtof(pszSpheroidEqRadius);
539 0 : const double polar_radius = CPLAtof(pszSpheroidPolarRadius);
540 0 : oProj.SetGeogCS("unknown", "unknown", "unknown", eq_radius,
541 0 : eq_radius / (eq_radius - polar_radius));
542 0 : oLL.SetGeogCS("unknown", "unknown", "unknown", eq_radius,
543 0 : eq_radius / (eq_radius - polar_radius));
544 : }
545 : else
546 : {
547 0 : CPLError(CE_Warning, CPLE_AppDefined,
548 : "Radii not specified for user-defined ellipsoid. "
549 : "Using wgs-84 parameters.");
550 0 : oProj.SetWellKnownGeogCS("WGS84");
551 0 : oLL.SetWellKnownGeogCS("WGS84");
552 : }
553 : }
554 : else
555 : {
556 0 : CPLError(CE_Warning, CPLE_AppDefined,
557 : "Unrecognized ellipsoid. Using wgs-84 parameters.");
558 0 : oProj.SetWellKnownGeogCS("WGS84");
559 0 : oLL.SetWellKnownGeogCS("WGS84");
560 : }
561 : }
562 :
563 : /* If a geotransform is sufficient to represent the GCP's (i.e. each */
564 : /* estimated gcp is within 0.25*pixel size of the actual value- this */
565 : /* is the test applied by GDALGCPsToGeoTransform), store the */
566 : /* geotransform. */
567 6 : bool transform_ok = false;
568 :
569 6 : if (EQUAL(pszProjName, "LL"))
570 : {
571 0 : transform_ok = CPL_TO_BOOL(
572 0 : GDALGCPsToGeoTransform(nGCPCount, pasGCPList, m_gt.data(), 0));
573 : }
574 : else
575 : {
576 : OGRCoordinateTransformation *poTransform =
577 6 : OGRCreateCoordinateTransformation(&oLL, &oProj);
578 6 : bool bSuccess = true;
579 6 : if (poTransform == nullptr)
580 : {
581 0 : CPLErrorReset();
582 0 : bSuccess = FALSE;
583 : }
584 :
585 : double *dfPrjX =
586 6 : static_cast<double *>(CPLMalloc(nGCPCount * sizeof(double)));
587 : double *dfPrjY =
588 6 : static_cast<double *>(CPLMalloc(nGCPCount * sizeof(double)));
589 :
590 20 : for (int gcp_index = 0; gcp_index < nGCPCount; gcp_index++)
591 : {
592 14 : dfPrjX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
593 14 : dfPrjY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
594 :
595 28 : if (bSuccess && !poTransform->Transform(1, &(dfPrjX[gcp_index]),
596 14 : &(dfPrjY[gcp_index])))
597 0 : bSuccess = FALSE;
598 : }
599 :
600 6 : if (bSuccess)
601 : {
602 20 : for (int gcp_index = 0; gcp_index < nGCPCount; gcp_index++)
603 : {
604 14 : pasGCPList[gcp_index].dfGCPX = dfPrjX[gcp_index];
605 14 : pasGCPList[gcp_index].dfGCPY = dfPrjY[gcp_index];
606 : }
607 12 : transform_ok = CPL_TO_BOOL(
608 6 : GDALGCPsToGeoTransform(nGCPCount, pasGCPList, m_gt.data(), 0));
609 : }
610 :
611 6 : if (poTransform)
612 6 : delete poTransform;
613 :
614 6 : CPLFree(dfPrjX);
615 6 : CPLFree(dfPrjY);
616 : }
617 :
618 6 : m_oSRS = oProj;
619 6 : m_oGCPSRS = std::move(oProj);
620 :
621 6 : if (!transform_ok)
622 : {
623 : /* transform is sufficient in some cases (slant range, standalone gcps)
624 : */
625 4 : m_gt = GDALGeoTransform();
626 4 : m_oSRS.Clear();
627 : }
628 :
629 6 : delete mffEllipsoids;
630 : }
631 :
632 : /************************************************************************/
633 : /* Open() */
634 : /************************************************************************/
635 :
636 34420 : GDALDataset *MFFDataset::Open(GDALOpenInfo *poOpenInfo)
637 :
638 : {
639 : /* -------------------------------------------------------------------- */
640 : /* We assume the user is pointing to the header file. */
641 : /* -------------------------------------------------------------------- */
642 34420 : if (poOpenInfo->nHeaderBytes < 17 || poOpenInfo->fpL == nullptr)
643 30973 : return nullptr;
644 :
645 3447 : if (!poOpenInfo->IsExtensionEqualToCI("hdr"))
646 3429 : return nullptr;
647 :
648 : /* -------------------------------------------------------------------- */
649 : /* Load the .hdr file, and compress white space out around the */
650 : /* equal sign. */
651 : /* -------------------------------------------------------------------- */
652 18 : char **papszHdrLines = CSLLoad(poOpenInfo->pszFilename);
653 18 : if (papszHdrLines == nullptr)
654 0 : return nullptr;
655 :
656 : // Remove spaces. e.g.
657 : // SPHEROID_NAME = CLARKE_1866 -> SPHEROID_NAME=CLARKE_1866
658 577 : for (int i = 0; papszHdrLines[i] != nullptr; i++)
659 : {
660 559 : int iDst = 0;
661 559 : char *pszLine = papszHdrLines[i];
662 :
663 16291 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
664 : {
665 15732 : if (pszLine[iSrc] != ' ')
666 : {
667 13778 : pszLine[iDst++] = pszLine[iSrc];
668 : }
669 : }
670 559 : pszLine[iDst] = '\0';
671 : }
672 :
673 : /* -------------------------------------------------------------------- */
674 : /* Verify it is an MFF file. */
675 : /* -------------------------------------------------------------------- */
676 26 : if (CSLFetchNameValue(papszHdrLines, "IMAGE_FILE_FORMAT") != nullptr &&
677 8 : !EQUAL(CSLFetchNameValue(papszHdrLines, "IMAGE_FILE_FORMAT"), "MFF"))
678 : {
679 0 : CSLDestroy(papszHdrLines);
680 0 : return nullptr;
681 : }
682 :
683 18 : if ((CSLFetchNameValue(papszHdrLines, "IMAGE_LINES") == nullptr ||
684 28 : CSLFetchNameValue(papszHdrLines, "LINE_SAMPLES") == nullptr) &&
685 10 : (CSLFetchNameValue(papszHdrLines, "no_rows") == nullptr ||
686 0 : CSLFetchNameValue(papszHdrLines, "no_columns") == nullptr))
687 : {
688 10 : CSLDestroy(papszHdrLines);
689 10 : return nullptr;
690 : }
691 :
692 : /* -------------------------------------------------------------------- */
693 : /* Create a corresponding GDALDataset. */
694 : /* -------------------------------------------------------------------- */
695 16 : auto poDS = std::make_unique<MFFDataset>();
696 :
697 8 : poDS->papszHdrLines = papszHdrLines;
698 :
699 8 : poDS->eAccess = poOpenInfo->eAccess;
700 :
701 : /* -------------------------------------------------------------------- */
702 : /* Set some dataset wide information. */
703 : /* -------------------------------------------------------------------- */
704 11 : if (CSLFetchNameValue(papszHdrLines, "no_rows") != nullptr &&
705 3 : CSLFetchNameValue(papszHdrLines, "no_columns") != nullptr)
706 : {
707 0 : poDS->nRasterXSize =
708 0 : atoi(CSLFetchNameValue(papszHdrLines, "no_columns"));
709 0 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszHdrLines, "no_rows"));
710 : }
711 : else
712 : {
713 16 : poDS->nRasterXSize =
714 8 : atoi(CSLFetchNameValue(papszHdrLines, "LINE_SAMPLES"));
715 8 : poDS->nRasterYSize =
716 8 : atoi(CSLFetchNameValue(papszHdrLines, "IMAGE_LINES"));
717 : }
718 :
719 8 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
720 : {
721 0 : return nullptr;
722 : }
723 :
724 8 : RawRasterBand::ByteOrder eByteOrder = RawRasterBand::NATIVE_BYTE_ORDER;
725 :
726 8 : const char *pszByteOrder = CSLFetchNameValue(papszHdrLines, "BYTE_ORDER");
727 8 : if (pszByteOrder)
728 : {
729 3 : eByteOrder = EQUAL(pszByteOrder, "LSB")
730 3 : ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
731 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN;
732 : }
733 :
734 : /* -------------------------------------------------------------------- */
735 : /* Get some information specific to APP tiled files. */
736 : /* -------------------------------------------------------------------- */
737 8 : int nTileXSize = 0;
738 8 : int nTileYSize = 0;
739 8 : const char *pszRefinedType = CSLFetchNameValue(papszHdrLines, "type");
740 8 : const bool bTiled = CSLFetchNameValue(papszHdrLines, "no_rows") != nullptr;
741 :
742 8 : if (bTiled)
743 : {
744 3 : if (CSLFetchNameValue(papszHdrLines, "tile_size_rows"))
745 3 : nTileYSize =
746 3 : atoi(CSLFetchNameValue(papszHdrLines, "tile_size_rows"));
747 3 : if (CSLFetchNameValue(papszHdrLines, "tile_size_columns"))
748 3 : nTileXSize =
749 3 : atoi(CSLFetchNameValue(papszHdrLines, "tile_size_columns"));
750 :
751 3 : if (nTileXSize <= 0 || nTileYSize <= 0 ||
752 9 : poDS->nRasterXSize - 1 > INT_MAX - nTileXSize ||
753 3 : poDS->nRasterYSize - 1 > INT_MAX - nTileYSize)
754 : {
755 0 : return nullptr;
756 : }
757 : }
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Read the directory to find matching band files. */
761 : /* -------------------------------------------------------------------- */
762 16 : const std::string osTargetPath = CPLGetPathSafe(poOpenInfo->pszFilename);
763 : const std::string osTargetBase =
764 16 : CPLGetBasenameSafe(poOpenInfo->pszFilename);
765 : const CPLStringList aosDirFiles(
766 16 : VSIReadDir(CPLGetPathSafe(poOpenInfo->pszFilename).c_str()));
767 8 : if (aosDirFiles.empty())
768 : {
769 0 : return nullptr;
770 : }
771 :
772 8 : int nSkipped = 0;
773 8 : for (int nRawBand = 0; true; nRawBand++)
774 : {
775 : /* Find the next raw band file. */
776 :
777 15 : int i = 0; // Used after for.
778 89 : for (; i < aosDirFiles.size(); i++)
779 : {
780 82 : if (!EQUAL(CPLGetBasenameSafe(aosDirFiles[i]).c_str(),
781 : osTargetBase.c_str()))
782 57 : continue;
783 :
784 25 : const std::string osExtension = CPLGetExtensionSafe(aosDirFiles[i]);
785 25 : if (osExtension.size() >= 2 &&
786 25 : isdigit(static_cast<unsigned char>(osExtension[1])) &&
787 58 : atoi(osExtension.c_str() + 1) == nRawBand &&
788 8 : strchr("bBcCiIjJrRxXzZ", osExtension[0]) != nullptr)
789 8 : break;
790 : }
791 :
792 15 : if (i == aosDirFiles.size())
793 7 : break;
794 :
795 : /* open the file for required level of access */
796 : const std::string osRawFilename =
797 8 : CPLFormFilenameSafe(osTargetPath.c_str(), aosDirFiles[i], nullptr);
798 :
799 8 : VSILFILE *fpRaw = nullptr;
800 8 : if (poOpenInfo->eAccess == GA_Update)
801 0 : fpRaw = VSIFOpenL(osRawFilename.c_str(), "rb+");
802 : else
803 8 : fpRaw = VSIFOpenL(osRawFilename.c_str(), "rb");
804 :
805 8 : if (fpRaw == nullptr)
806 : {
807 0 : CPLError(CE_Warning, CPLE_OpenFailed,
808 : "Unable to open %s ... skipping.", osRawFilename.c_str());
809 0 : nSkipped++;
810 0 : continue;
811 : }
812 16 : poDS->m_papszFileList =
813 8 : CSLAddString(poDS->m_papszFileList, osRawFilename.c_str());
814 :
815 8 : GDALDataType eDataType = GDT_Unknown;
816 8 : const std::string osExt = CPLGetExtensionSafe(aosDirFiles[i]);
817 8 : if (pszRefinedType != nullptr)
818 : {
819 0 : if (EQUAL(pszRefinedType, "C*4"))
820 0 : eDataType = GDT_CFloat32;
821 0 : else if (EQUAL(pszRefinedType, "C*8"))
822 0 : eDataType = GDT_CFloat64;
823 0 : else if (EQUAL(pszRefinedType, "R*4"))
824 0 : eDataType = GDT_Float32;
825 0 : else if (EQUAL(pszRefinedType, "R*8"))
826 0 : eDataType = GDT_Float64;
827 0 : else if (EQUAL(pszRefinedType, "I*1"))
828 0 : eDataType = GDT_Byte;
829 0 : else if (EQUAL(pszRefinedType, "I*2"))
830 0 : eDataType = GDT_Int16;
831 0 : else if (EQUAL(pszRefinedType, "I*4"))
832 0 : eDataType = GDT_Int32;
833 0 : else if (EQUAL(pszRefinedType, "U*2"))
834 0 : eDataType = GDT_UInt16;
835 0 : else if (EQUAL(pszRefinedType, "U*4"))
836 0 : eDataType = GDT_UInt32;
837 0 : else if (EQUAL(pszRefinedType, "J*1"))
838 : {
839 0 : CPLError(
840 : CE_Warning, CPLE_OpenFailed,
841 : "Unable to open band %d because type J*1 is not handled. "
842 : "Skipping.",
843 : nRawBand + 1);
844 0 : nSkipped++;
845 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
846 0 : continue; // Does not support 1 byte complex.
847 : }
848 0 : else if (EQUAL(pszRefinedType, "J*2"))
849 0 : eDataType = GDT_CInt16;
850 0 : else if (EQUAL(pszRefinedType, "K*4"))
851 0 : eDataType = GDT_CInt32;
852 : else
853 : {
854 0 : CPLError(
855 : CE_Warning, CPLE_OpenFailed,
856 : "Unable to open band %d because type %s is not handled. "
857 : "Skipping.\n",
858 : nRawBand + 1, pszRefinedType);
859 0 : nSkipped++;
860 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
861 0 : continue;
862 : }
863 : }
864 8 : else if (STARTS_WITH_CI(osExt.c_str(), "b"))
865 : {
866 7 : eDataType = GDT_Byte;
867 : }
868 1 : else if (STARTS_WITH_CI(osExt.c_str(), "i"))
869 : {
870 0 : eDataType = GDT_UInt16;
871 : }
872 1 : else if (STARTS_WITH_CI(osExt.c_str(), "j"))
873 : {
874 0 : eDataType = GDT_CInt16;
875 : }
876 1 : else if (STARTS_WITH_CI(osExt.c_str(), "r"))
877 : {
878 1 : eDataType = GDT_Float32;
879 : }
880 0 : else if (STARTS_WITH_CI(osExt.c_str(), "x"))
881 : {
882 0 : eDataType = GDT_CFloat32;
883 : }
884 : else
885 : {
886 0 : CPLError(CE_Warning, CPLE_OpenFailed,
887 : "Unable to open band %d because extension %s is not "
888 : "handled. Skipping.",
889 : nRawBand + 1, osExt.c_str());
890 0 : nSkipped++;
891 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
892 0 : continue;
893 : }
894 :
895 8 : const int nBand = poDS->GetRasterCount() + 1;
896 :
897 8 : const int nPixelOffset = GDALGetDataTypeSizeBytes(eDataType);
898 0 : std::unique_ptr<GDALRasterBand> poBand;
899 :
900 8 : if (bTiled)
901 : {
902 3 : if (nTileXSize > INT_MAX / nTileYSize / nPixelOffset)
903 : {
904 1 : CPLError(CE_Failure, CPLE_AppDefined, "Too large tile");
905 1 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
906 1 : return nullptr;
907 : }
908 :
909 4 : poBand = std::make_unique<MFFTiledBand>(poDS.get(), nBand, fpRaw,
910 : nTileXSize, nTileYSize,
911 2 : eDataType, eByteOrder);
912 : }
913 : else
914 : {
915 10 : if (nPixelOffset != 0 &&
916 5 : poDS->GetRasterXSize() > INT_MAX / nPixelOffset)
917 : {
918 0 : CPLError(CE_Warning, CPLE_AppDefined,
919 : "Int overflow occurred... skipping");
920 0 : nSkipped++;
921 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
922 0 : continue;
923 : }
924 :
925 10 : poBand = RawRasterBand::Create(
926 5 : poDS.get(), nBand, fpRaw, 0, nPixelOffset,
927 5 : nPixelOffset * poDS->GetRasterXSize(), eDataType, eByteOrder,
928 5 : RawRasterBand::OwnFP::YES);
929 : }
930 :
931 7 : poDS->SetBand(nBand, std::move(poBand));
932 7 : }
933 :
934 : /* -------------------------------------------------------------------- */
935 : /* Check if we have bands. */
936 : /* -------------------------------------------------------------------- */
937 7 : if (poDS->GetRasterCount() == 0)
938 : {
939 0 : if (nSkipped > 0 && poOpenInfo->eAccess)
940 : {
941 0 : CPLError(CE_Failure, CPLE_OpenFailed,
942 : "Failed to open %d files that were apparently bands. "
943 : "Perhaps this dataset is readonly?",
944 : nSkipped);
945 0 : return nullptr;
946 : }
947 : else
948 : {
949 0 : CPLError(CE_Failure, CPLE_OpenFailed,
950 : "MFF header file read successfully, but no bands "
951 : "were successfully found and opened.");
952 0 : return nullptr;
953 : }
954 : }
955 :
956 : /* -------------------------------------------------------------------- */
957 : /* Set all information from the .hdr that isn't well know to be */
958 : /* metadata. */
959 : /* -------------------------------------------------------------------- */
960 81 : for (int i = 0; papszHdrLines[i] != nullptr; i++)
961 : {
962 74 : char *pszName = nullptr;
963 :
964 74 : const char *pszValue = CPLParseNameValue(papszHdrLines[i], &pszName);
965 74 : if (pszName == nullptr || pszValue == nullptr)
966 3 : continue;
967 :
968 71 : if (!EQUAL(pszName, "END") && !EQUAL(pszName, "FILE_TYPE") &&
969 68 : !EQUAL(pszName, "BYTE_ORDER") && !EQUAL(pszName, "no_columns") &&
970 65 : !EQUAL(pszName, "no_rows") && !EQUAL(pszName, "type") &&
971 63 : !EQUAL(pszName, "tile_size_rows") &&
972 61 : !EQUAL(pszName, "tile_size_columns") &&
973 59 : !EQUAL(pszName, "IMAGE_FILE_FORMAT") &&
974 52 : !EQUAL(pszName, "IMAGE_LINES") && !EQUAL(pszName, "LINE_SAMPLES"))
975 : {
976 38 : poDS->SetMetadataItem(pszName, pszValue);
977 : }
978 :
979 71 : CPLFree(pszName);
980 : }
981 :
982 : /* -------------------------------------------------------------------- */
983 : /* Any GCPs in header file? */
984 : /* -------------------------------------------------------------------- */
985 7 : poDS->ScanForGCPs();
986 7 : poDS->ScanForProjectionInfo();
987 7 : if (poDS->nGCPCount == 0)
988 1 : poDS->m_oGCPSRS.Clear();
989 :
990 : /* -------------------------------------------------------------------- */
991 : /* Initialize any PAM information. */
992 : /* -------------------------------------------------------------------- */
993 7 : poDS->SetDescription(poOpenInfo->pszFilename);
994 7 : poDS->TryLoadXML();
995 :
996 : /* -------------------------------------------------------------------- */
997 : /* Check for overviews. */
998 : /* -------------------------------------------------------------------- */
999 7 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1000 :
1001 7 : return poDS.release();
1002 : }
1003 :
1004 : /************************************************************************/
1005 : /* GDALRegister_MFF() */
1006 : /************************************************************************/
1007 :
1008 2038 : void GDALRegister_MFF()
1009 :
1010 : {
1011 2038 : if (GDALGetDriverByName("MFF") != nullptr)
1012 283 : return;
1013 :
1014 1755 : GDALDriver *poDriver = new GDALDriver();
1015 :
1016 1755 : poDriver->SetDescription("MFF");
1017 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1018 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF Raster");
1019 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff.html");
1020 1755 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
1021 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1022 :
1023 1755 : poDriver->pfnOpen = MFFDataset::Open;
1024 :
1025 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
1026 : }
|
