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