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