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 : static int GetMFFProjectionType(const OGRSpatialReference *poSRS);
33 :
34 : /************************************************************************/
35 : /* ==================================================================== */
36 : /* MFFDataset */
37 : /* ==================================================================== */
38 : /************************************************************************/
39 :
40 : class MFFDataset final : public RawDataset
41 : {
42 : int nGCPCount;
43 : GDAL_GCP *pasGCPList;
44 :
45 : OGRSpatialReference m_oSRS{};
46 : OGRSpatialReference m_oGCPSRS{};
47 : double adfGeoTransform[6];
48 : char **m_papszFileList;
49 :
50 : void ScanForGCPs();
51 : void ScanForProjectionInfo();
52 :
53 : CPL_DISALLOW_COPY_ASSIGN(MFFDataset)
54 :
55 : CPLErr Close() override;
56 :
57 : public:
58 : MFFDataset();
59 : ~MFFDataset() override;
60 :
61 : char **papszHdrLines;
62 :
63 : VSILFILE **pafpBandFiles;
64 :
65 : char **GetFileList() override;
66 :
67 : int GetGCPCount() override;
68 :
69 0 : const OGRSpatialReference *GetGCPSpatialRef() const override
70 : {
71 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
72 : }
73 :
74 : const GDAL_GCP *GetGCPs() override;
75 :
76 9 : const OGRSpatialReference *GetSpatialRef() const override
77 : {
78 9 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
79 : }
80 :
81 : CPLErr GetGeoTransform(double *) override;
82 :
83 : static GDALDataset *Open(GDALOpenInfo *);
84 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
85 : int nBandsIn, GDALDataType eType,
86 : char **papszParamList);
87 : static GDALDataset *CreateCopy(const char *pszFilename,
88 : GDALDataset *poSrcDS, int bStrict,
89 : char **papszOptions,
90 : GDALProgressFunc pfnProgress,
91 : void *pProgressData);
92 : };
93 :
94 : /************************************************************************/
95 : /* ==================================================================== */
96 : /* MFFTiledBand */
97 : /* ==================================================================== */
98 : /************************************************************************/
99 :
100 : class MFFTiledBand final : public GDALRasterBand
101 : {
102 : friend class MFFDataset;
103 :
104 : VSILFILE *fpRaw;
105 : RawRasterBand::ByteOrder eByteOrder;
106 :
107 : CPL_DISALLOW_COPY_ASSIGN(MFFTiledBand)
108 :
109 : public:
110 : MFFTiledBand(MFFDataset *, int, VSILFILE *, int, int, GDALDataType,
111 : RawRasterBand::ByteOrder);
112 : ~MFFTiledBand() override;
113 :
114 : CPLErr IReadBlock(int, int, void *) override;
115 : };
116 :
117 : /************************************************************************/
118 : /* MFFTiledBand() */
119 : /************************************************************************/
120 :
121 2 : MFFTiledBand::MFFTiledBand(MFFDataset *poDSIn, int nBandIn, VSILFILE *fp,
122 : int nTileXSize, int nTileYSize,
123 : GDALDataType eDataTypeIn,
124 2 : RawRasterBand::ByteOrder eByteOrderIn)
125 2 : : fpRaw(fp), eByteOrder(eByteOrderIn)
126 : {
127 2 : poDS = poDSIn;
128 2 : nBand = nBandIn;
129 :
130 2 : eDataType = eDataTypeIn;
131 :
132 2 : nBlockXSize = nTileXSize;
133 2 : nBlockYSize = nTileYSize;
134 2 : }
135 :
136 : /************************************************************************/
137 : /* ~MFFTiledBand() */
138 : /************************************************************************/
139 :
140 4 : MFFTiledBand::~MFFTiledBand()
141 :
142 : {
143 2 : if (VSIFCloseL(fpRaw) != 0)
144 : {
145 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
146 : }
147 4 : }
148 :
149 : /************************************************************************/
150 : /* IReadBlock() */
151 : /************************************************************************/
152 :
153 1 : CPLErr MFFTiledBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
154 :
155 : {
156 1 : const int nTilesPerRow = (nRasterXSize + nBlockXSize - 1) / nBlockXSize;
157 1 : const int nWordSize = GDALGetDataTypeSize(eDataType) / 8;
158 1 : const int nBlockSize = nWordSize * nBlockXSize * nBlockYSize;
159 :
160 1 : const vsi_l_offset nOffset =
161 1 : nBlockSize *
162 1 : (nBlockXOff + static_cast<vsi_l_offset>(nBlockYOff) * nTilesPerRow);
163 :
164 2 : if (VSIFSeekL(fpRaw, nOffset, SEEK_SET) == -1 ||
165 1 : VSIFReadL(pImage, 1, nBlockSize, fpRaw) < 1)
166 : {
167 0 : CPLError(CE_Failure, CPLE_FileIO,
168 : "Read of tile %d/%d failed with fseek or fread error.",
169 : nBlockXOff, nBlockYOff);
170 0 : return CE_Failure;
171 : }
172 :
173 1 : if (eByteOrder != RawRasterBand::NATIVE_BYTE_ORDER && nWordSize > 1)
174 : {
175 0 : if (GDALDataTypeIsComplex(eDataType))
176 : {
177 0 : GDALSwapWords(pImage, nWordSize / 2, nBlockXSize * nBlockYSize,
178 : nWordSize);
179 0 : GDALSwapWords(reinterpret_cast<GByte *>(pImage) + nWordSize / 2,
180 0 : nWordSize / 2, nBlockXSize * nBlockYSize, nWordSize);
181 : }
182 : else
183 0 : GDALSwapWords(pImage, nWordSize, nBlockXSize * nBlockYSize,
184 : nWordSize);
185 : }
186 :
187 1 : return CE_None;
188 : }
189 :
190 : /************************************************************************/
191 : /* MFF Spheroids */
192 : /************************************************************************/
193 :
194 : class MFFSpheroidList : public SpheroidList
195 : {
196 : public:
197 : MFFSpheroidList();
198 :
199 15 : ~MFFSpheroidList()
200 15 : {
201 15 : }
202 : };
203 :
204 15 : MFFSpheroidList ::MFFSpheroidList()
205 : {
206 15 : num_spheroids = 18;
207 :
208 15 : epsilonR = 0.1;
209 15 : epsilonI = 0.000001;
210 :
211 15 : spheroids[0].SetValuesByRadii("SPHERE", 6371007.0, 6371007.0);
212 15 : spheroids[1].SetValuesByRadii("EVEREST", 6377304.0, 6356103.0);
213 15 : spheroids[2].SetValuesByRadii("BESSEL", 6377397.0, 6356082.0);
214 15 : spheroids[3].SetValuesByRadii("AIRY", 6377563.0, 6356300.0);
215 15 : spheroids[4].SetValuesByRadii("CLARKE_1858", 6378294.0, 6356621.0);
216 15 : spheroids[5].SetValuesByRadii("CLARKE_1866", 6378206.4, 6356583.8);
217 15 : spheroids[6].SetValuesByRadii("CLARKE_1880", 6378249.0, 6356517.0);
218 15 : spheroids[7].SetValuesByRadii("HAYFORD", 6378388.0, 6356915.0);
219 15 : spheroids[8].SetValuesByRadii("KRASOVSKI", 6378245.0, 6356863.0);
220 15 : spheroids[9].SetValuesByRadii("HOUGH", 6378270.0, 6356794.0);
221 15 : spheroids[10].SetValuesByRadii("FISHER_60", 6378166.0, 6356784.0);
222 15 : spheroids[11].SetValuesByRadii("KAULA", 6378165.0, 6356345.0);
223 15 : spheroids[12].SetValuesByRadii("IUGG_67", 6378160.0, 6356775.0);
224 15 : spheroids[13].SetValuesByRadii("FISHER_68", 6378150.0, 6356330.0);
225 15 : spheroids[14].SetValuesByRadii("WGS_72", 6378135.0, 6356751.0);
226 15 : spheroids[15].SetValuesByRadii("IUGG_75", 6378140.0, 6356755.0);
227 15 : spheroids[16].SetValuesByRadii("WGS_84", 6378137.0, 6356752.0);
228 15 : spheroids[17].SetValuesByRadii("HUGHES", 6378273.0, 6356889.4);
229 15 : }
230 :
231 : /************************************************************************/
232 : /* ==================================================================== */
233 : /* MFFDataset */
234 : /* ==================================================================== */
235 : /************************************************************************/
236 :
237 : /************************************************************************/
238 : /* MFFDataset() */
239 : /************************************************************************/
240 :
241 29 : MFFDataset::MFFDataset()
242 : : nGCPCount(0), pasGCPList(nullptr), m_papszFileList(nullptr),
243 29 : papszHdrLines(nullptr), pafpBandFiles(nullptr)
244 : {
245 29 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
246 29 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
247 29 : adfGeoTransform[0] = 0.0;
248 29 : adfGeoTransform[1] = 1.0;
249 29 : adfGeoTransform[2] = 0.0;
250 29 : adfGeoTransform[3] = 0.0;
251 29 : adfGeoTransform[4] = 0.0;
252 29 : adfGeoTransform[5] = 1.0;
253 29 : }
254 :
255 : /************************************************************************/
256 : /* ~MFFDataset() */
257 : /************************************************************************/
258 :
259 58 : MFFDataset::~MFFDataset()
260 :
261 : {
262 29 : MFFDataset::Close();
263 58 : }
264 :
265 : /************************************************************************/
266 : /* Close() */
267 : /************************************************************************/
268 :
269 58 : CPLErr MFFDataset::Close()
270 : {
271 58 : CPLErr eErr = CE_None;
272 58 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
273 : {
274 29 : if (MFFDataset::FlushCache(true) != CE_None)
275 0 : eErr = CE_Failure;
276 :
277 29 : CSLDestroy(papszHdrLines);
278 29 : if (pafpBandFiles)
279 : {
280 0 : for (int i = 0; i < GetRasterCount(); i++)
281 : {
282 0 : if (pafpBandFiles[i])
283 : {
284 0 : if (VSIFCloseL(pafpBandFiles[i]) != 0)
285 : {
286 0 : eErr = CE_Failure;
287 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
288 : }
289 : }
290 : }
291 0 : CPLFree(pafpBandFiles);
292 : }
293 :
294 29 : if (nGCPCount > 0)
295 : {
296 6 : GDALDeinitGCPs(nGCPCount, pasGCPList);
297 : }
298 29 : CPLFree(pasGCPList);
299 29 : CSLDestroy(m_papszFileList);
300 :
301 29 : if (GDALPamDataset::Close() != CE_None)
302 0 : eErr = CE_Failure;
303 : }
304 58 : return eErr;
305 : }
306 :
307 : /************************************************************************/
308 : /* GetFileList() */
309 : /************************************************************************/
310 :
311 3 : char **MFFDataset::GetFileList()
312 : {
313 3 : char **papszFileList = RawDataset::GetFileList();
314 3 : papszFileList = CSLInsertStrings(papszFileList, -1, m_papszFileList);
315 3 : return papszFileList;
316 : }
317 :
318 : /************************************************************************/
319 : /* GetGCPCount() */
320 : /************************************************************************/
321 :
322 0 : int MFFDataset::GetGCPCount()
323 :
324 : {
325 0 : return nGCPCount;
326 : }
327 :
328 : /************************************************************************/
329 : /* GetGeoTransform() */
330 : /************************************************************************/
331 :
332 9 : CPLErr MFFDataset::GetGeoTransform(double *padfTransform)
333 :
334 : {
335 9 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
336 9 : return CE_None;
337 : }
338 :
339 : /************************************************************************/
340 : /* GetGCP() */
341 : /************************************************************************/
342 :
343 0 : const GDAL_GCP *MFFDataset::GetGCPs()
344 :
345 : {
346 0 : return pasGCPList;
347 : }
348 :
349 : /************************************************************************/
350 : /* ScanForGCPs() */
351 : /************************************************************************/
352 :
353 29 : void MFFDataset::ScanForGCPs()
354 :
355 : {
356 29 : int NUM_GCPS = 0;
357 :
358 29 : if (CSLFetchNameValue(papszHdrLines, "NUM_GCPS") != nullptr)
359 4 : NUM_GCPS = atoi(CSLFetchNameValue(papszHdrLines, "NUM_GCPS"));
360 29 : if (NUM_GCPS < 0)
361 0 : return;
362 :
363 29 : nGCPCount = 0;
364 29 : pasGCPList =
365 29 : static_cast<GDAL_GCP *>(VSICalloc(sizeof(GDAL_GCP), 5 + NUM_GCPS));
366 29 : if (pasGCPList == nullptr)
367 0 : return;
368 :
369 174 : for (int nCorner = 0; nCorner < 5; nCorner++)
370 : {
371 145 : const char *pszBase = nullptr;
372 145 : double dfRasterX = 0.0;
373 145 : double dfRasterY = 0.0;
374 :
375 145 : if (nCorner == 0)
376 : {
377 29 : dfRasterX = 0.5;
378 29 : dfRasterY = 0.5;
379 29 : pszBase = "TOP_LEFT_CORNER";
380 : }
381 116 : else if (nCorner == 1)
382 : {
383 29 : dfRasterX = GetRasterXSize() - 0.5;
384 29 : dfRasterY = 0.5;
385 29 : pszBase = "TOP_RIGHT_CORNER";
386 : }
387 87 : else if (nCorner == 2)
388 : {
389 29 : dfRasterX = GetRasterXSize() - 0.5;
390 29 : dfRasterY = GetRasterYSize() - 0.5;
391 29 : pszBase = "BOTTOM_RIGHT_CORNER";
392 : }
393 58 : else if (nCorner == 3)
394 : {
395 29 : dfRasterX = 0.5;
396 29 : dfRasterY = GetRasterYSize() - 0.5;
397 29 : pszBase = "BOTTOM_LEFT_CORNER";
398 : }
399 : else /* if( nCorner == 4 ) */
400 : {
401 29 : dfRasterX = GetRasterXSize() / 2.0;
402 29 : dfRasterY = GetRasterYSize() / 2.0;
403 29 : pszBase = "CENTRE";
404 : }
405 :
406 145 : char szLatName[40] = {'\0'};
407 145 : char szLongName[40] = {'\0'};
408 145 : snprintf(szLatName, sizeof(szLatName), "%s_LATITUDE", pszBase);
409 145 : snprintf(szLongName, sizeof(szLongName), "%s_LONGITUDE", pszBase);
410 :
411 155 : if (CSLFetchNameValue(papszHdrLines, szLatName) != nullptr &&
412 10 : CSLFetchNameValue(papszHdrLines, szLongName) != nullptr)
413 : {
414 10 : GDALInitGCPs(1, pasGCPList + nGCPCount);
415 :
416 10 : CPLFree(pasGCPList[nGCPCount].pszId);
417 :
418 10 : pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
419 :
420 20 : pasGCPList[nGCPCount].dfGCPX =
421 10 : CPLAtof(CSLFetchNameValue(papszHdrLines, szLongName));
422 20 : pasGCPList[nGCPCount].dfGCPY =
423 10 : CPLAtof(CSLFetchNameValue(papszHdrLines, szLatName));
424 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
425 :
426 10 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
427 10 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
428 :
429 10 : nGCPCount++;
430 : }
431 : }
432 :
433 : /* -------------------------------------------------------------------- */
434 : /* Collect standalone GCPs. They look like: */
435 : /* */
436 : /* GCPn = row, col, lat, long */
437 : /* GCP1 = 1, 1, 45.0, -75.0 */
438 : /* -------------------------------------------------------------------- */
439 33 : for (int i = 0; i < NUM_GCPS; i++)
440 : {
441 4 : char szName[25] = {'\0'};
442 4 : snprintf(szName, sizeof(szName), "GCP%d", i + 1);
443 4 : if (CSLFetchNameValue(papszHdrLines, szName) == nullptr)
444 0 : continue;
445 :
446 4 : char **papszTokens = CSLTokenizeStringComplex(
447 4 : CSLFetchNameValue(papszHdrLines, szName), ",", FALSE, FALSE);
448 4 : if (CSLCount(papszTokens) == 4)
449 : {
450 4 : GDALInitGCPs(1, pasGCPList + nGCPCount);
451 :
452 4 : CPLFree(pasGCPList[nGCPCount].pszId);
453 4 : pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
454 :
455 4 : pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[3]);
456 4 : pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[2]);
457 4 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
458 4 : pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]) + 0.5;
459 4 : pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[0]) + 0.5;
460 :
461 4 : nGCPCount++;
462 : }
463 :
464 4 : CSLDestroy(papszTokens);
465 : }
466 : }
467 :
468 : /************************************************************************/
469 : /* ScanForProjectionInfo */
470 : /************************************************************************/
471 :
472 29 : void MFFDataset::ScanForProjectionInfo()
473 : {
474 : const char *pszProjName =
475 29 : CSLFetchNameValue(papszHdrLines, "PROJECTION_NAME");
476 : const char *pszOriginLong =
477 29 : CSLFetchNameValue(papszHdrLines, "PROJECTION_ORIGIN_LONGITUDE");
478 : const char *pszSpheroidName =
479 29 : CSLFetchNameValue(papszHdrLines, "SPHEROID_NAME");
480 :
481 29 : if (pszProjName == nullptr)
482 : {
483 23 : m_oSRS.Clear();
484 23 : m_oGCPSRS.Clear();
485 23 : return;
486 : }
487 6 : else if ((!EQUAL(pszProjName, "utm")) && (!EQUAL(pszProjName, "ll")))
488 : {
489 0 : CPLError(CE_Warning, CPLE_AppDefined,
490 : "Only utm and lat/long projections are currently supported.");
491 0 : m_oSRS.Clear();
492 0 : m_oGCPSRS.Clear();
493 0 : return;
494 : }
495 6 : MFFSpheroidList *mffEllipsoids = new MFFSpheroidList;
496 :
497 12 : OGRSpatialReference oProj;
498 6 : oProj.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
499 6 : if (EQUAL(pszProjName, "utm"))
500 : {
501 : int nZone;
502 :
503 6 : if (pszOriginLong == nullptr)
504 : {
505 : // If origin not specified, assume 0.0.
506 4 : CPLError(
507 : CE_Warning, CPLE_AppDefined,
508 : "No projection origin longitude specified. Assuming 0.0.");
509 4 : nZone = 31;
510 : }
511 : else
512 2 : nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
513 :
514 6 : if (nGCPCount >= 5 && pasGCPList[4].dfGCPY < 0)
515 0 : oProj.SetUTM(nZone, 0);
516 : else
517 6 : oProj.SetUTM(nZone, 1);
518 :
519 6 : if (pszOriginLong != nullptr)
520 2 : oProj.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
521 : }
522 :
523 12 : OGRSpatialReference oLL;
524 6 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
525 6 : if (pszOriginLong != nullptr)
526 2 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
527 :
528 6 : if (pszSpheroidName == nullptr)
529 : {
530 4 : CPLError(CE_Warning, CPLE_AppDefined,
531 : "Unspecified ellipsoid. Using wgs-84 parameters.\n");
532 :
533 4 : oProj.SetWellKnownGeogCS("WGS84");
534 4 : oLL.SetWellKnownGeogCS("WGS84");
535 : }
536 : else
537 : {
538 2 : if (mffEllipsoids->SpheroidInList(pszSpheroidName))
539 : {
540 2 : oProj.SetGeogCS(
541 : "unknown", "unknown", pszSpheroidName,
542 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
543 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName));
544 2 : oLL.SetGeogCS(
545 : "unknown", "unknown", pszSpheroidName,
546 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
547 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName));
548 : }
549 0 : else if (EQUAL(pszSpheroidName, "USER_DEFINED"))
550 : {
551 : const char *pszSpheroidEqRadius =
552 0 : CSLFetchNameValue(papszHdrLines, "SPHEROID_EQUATORIAL_RADIUS");
553 : const char *pszSpheroidPolarRadius =
554 0 : CSLFetchNameValue(papszHdrLines, "SPHEROID_POLAR_RADIUS");
555 0 : if ((pszSpheroidEqRadius != nullptr) &&
556 : (pszSpheroidPolarRadius != nullptr))
557 : {
558 0 : const double eq_radius = CPLAtof(pszSpheroidEqRadius);
559 0 : const double polar_radius = CPLAtof(pszSpheroidPolarRadius);
560 0 : oProj.SetGeogCS("unknown", "unknown", "unknown", eq_radius,
561 0 : eq_radius / (eq_radius - polar_radius));
562 0 : oLL.SetGeogCS("unknown", "unknown", "unknown", eq_radius,
563 0 : eq_radius / (eq_radius - polar_radius));
564 : }
565 : else
566 : {
567 0 : CPLError(CE_Warning, CPLE_AppDefined,
568 : "Radii not specified for user-defined ellipsoid. "
569 : "Using wgs-84 parameters.");
570 0 : oProj.SetWellKnownGeogCS("WGS84");
571 0 : oLL.SetWellKnownGeogCS("WGS84");
572 : }
573 : }
574 : else
575 : {
576 0 : CPLError(CE_Warning, CPLE_AppDefined,
577 : "Unrecognized ellipsoid. Using wgs-84 parameters.");
578 0 : oProj.SetWellKnownGeogCS("WGS84");
579 0 : oLL.SetWellKnownGeogCS("WGS84");
580 : }
581 : }
582 :
583 : /* If a geotransform is sufficient to represent the GCP's (i.e. each */
584 : /* estimated gcp is within 0.25*pixel size of the actual value- this */
585 : /* is the test applied by GDALGCPsToGeoTransform), store the */
586 : /* geotransform. */
587 6 : bool transform_ok = false;
588 :
589 6 : if (EQUAL(pszProjName, "LL"))
590 : {
591 0 : transform_ok = CPL_TO_BOOL(
592 0 : GDALGCPsToGeoTransform(nGCPCount, pasGCPList, adfGeoTransform, 0));
593 : }
594 : else
595 : {
596 : OGRCoordinateTransformation *poTransform =
597 6 : OGRCreateCoordinateTransformation(&oLL, &oProj);
598 6 : bool bSuccess = true;
599 6 : if (poTransform == nullptr)
600 : {
601 0 : CPLErrorReset();
602 0 : bSuccess = FALSE;
603 : }
604 :
605 : double *dfPrjX =
606 6 : static_cast<double *>(CPLMalloc(nGCPCount * sizeof(double)));
607 : double *dfPrjY =
608 6 : static_cast<double *>(CPLMalloc(nGCPCount * sizeof(double)));
609 :
610 20 : for (int gcp_index = 0; gcp_index < nGCPCount; gcp_index++)
611 : {
612 14 : dfPrjX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
613 14 : dfPrjY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
614 :
615 28 : if (bSuccess && !poTransform->Transform(1, &(dfPrjX[gcp_index]),
616 14 : &(dfPrjY[gcp_index])))
617 0 : bSuccess = FALSE;
618 : }
619 :
620 6 : if (bSuccess)
621 : {
622 20 : for (int gcp_index = 0; gcp_index < nGCPCount; gcp_index++)
623 : {
624 14 : pasGCPList[gcp_index].dfGCPX = dfPrjX[gcp_index];
625 14 : pasGCPList[gcp_index].dfGCPY = dfPrjY[gcp_index];
626 : }
627 6 : transform_ok = CPL_TO_BOOL(GDALGCPsToGeoTransform(
628 6 : nGCPCount, pasGCPList, adfGeoTransform, 0));
629 : }
630 :
631 6 : if (poTransform)
632 6 : delete poTransform;
633 :
634 6 : CPLFree(dfPrjX);
635 6 : CPLFree(dfPrjY);
636 : }
637 :
638 6 : m_oSRS = oProj;
639 6 : m_oGCPSRS = std::move(oProj);
640 :
641 6 : if (!transform_ok)
642 : {
643 : /* transform is sufficient in some cases (slant range, standalone gcps)
644 : */
645 4 : adfGeoTransform[0] = 0.0;
646 4 : adfGeoTransform[1] = 1.0;
647 4 : adfGeoTransform[2] = 0.0;
648 4 : adfGeoTransform[3] = 0.0;
649 4 : adfGeoTransform[4] = 0.0;
650 4 : adfGeoTransform[5] = 1.0;
651 4 : m_oSRS.Clear();
652 : }
653 :
654 6 : delete mffEllipsoids;
655 : }
656 :
657 : /************************************************************************/
658 : /* Open() */
659 : /************************************************************************/
660 :
661 31066 : GDALDataset *MFFDataset::Open(GDALOpenInfo *poOpenInfo)
662 :
663 : {
664 : /* -------------------------------------------------------------------- */
665 : /* We assume the user is pointing to the header file. */
666 : /* -------------------------------------------------------------------- */
667 31066 : if (poOpenInfo->nHeaderBytes < 17 || poOpenInfo->fpL == nullptr)
668 27687 : return nullptr;
669 :
670 3379 : if (!poOpenInfo->IsExtensionEqualToCI("hdr"))
671 3340 : return nullptr;
672 :
673 : /* -------------------------------------------------------------------- */
674 : /* Load the .hdr file, and compress white space out around the */
675 : /* equal sign. */
676 : /* -------------------------------------------------------------------- */
677 39 : char **papszHdrLines = CSLLoad(poOpenInfo->pszFilename);
678 39 : if (papszHdrLines == nullptr)
679 0 : return nullptr;
680 :
681 : // Remove spaces. e.g.
682 : // SPHEROID_NAME = CLARKE_1866 -> SPHEROID_NAME=CLARKE_1866
683 712 : for (int i = 0; papszHdrLines[i] != nullptr; i++)
684 : {
685 673 : int iDst = 0;
686 673 : char *pszLine = papszHdrLines[i];
687 :
688 18285 : for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
689 : {
690 17612 : if (pszLine[iSrc] != ' ')
691 : {
692 15438 : pszLine[iDst++] = pszLine[iSrc];
693 : }
694 : }
695 673 : pszLine[iDst] = '\0';
696 : }
697 :
698 : /* -------------------------------------------------------------------- */
699 : /* Verify it is an MFF file. */
700 : /* -------------------------------------------------------------------- */
701 68 : if (CSLFetchNameValue(papszHdrLines, "IMAGE_FILE_FORMAT") != nullptr &&
702 29 : !EQUAL(CSLFetchNameValue(papszHdrLines, "IMAGE_FILE_FORMAT"), "MFF"))
703 : {
704 0 : CSLDestroy(papszHdrLines);
705 0 : return nullptr;
706 : }
707 :
708 39 : if ((CSLFetchNameValue(papszHdrLines, "IMAGE_LINES") == nullptr ||
709 49 : CSLFetchNameValue(papszHdrLines, "LINE_SAMPLES") == nullptr) &&
710 10 : (CSLFetchNameValue(papszHdrLines, "no_rows") == nullptr ||
711 0 : CSLFetchNameValue(papszHdrLines, "no_columns") == nullptr))
712 : {
713 10 : CSLDestroy(papszHdrLines);
714 10 : return nullptr;
715 : }
716 :
717 : /* -------------------------------------------------------------------- */
718 : /* Create a corresponding GDALDataset. */
719 : /* -------------------------------------------------------------------- */
720 58 : auto poDS = std::make_unique<MFFDataset>();
721 :
722 29 : poDS->papszHdrLines = papszHdrLines;
723 :
724 29 : poDS->eAccess = poOpenInfo->eAccess;
725 :
726 : /* -------------------------------------------------------------------- */
727 : /* Set some dataset wide information. */
728 : /* -------------------------------------------------------------------- */
729 31 : if (CSLFetchNameValue(papszHdrLines, "no_rows") != nullptr &&
730 2 : CSLFetchNameValue(papszHdrLines, "no_columns") != nullptr)
731 : {
732 0 : poDS->nRasterXSize =
733 0 : atoi(CSLFetchNameValue(papszHdrLines, "no_columns"));
734 0 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszHdrLines, "no_rows"));
735 : }
736 : else
737 : {
738 58 : poDS->nRasterXSize =
739 29 : atoi(CSLFetchNameValue(papszHdrLines, "LINE_SAMPLES"));
740 29 : poDS->nRasterYSize =
741 29 : atoi(CSLFetchNameValue(papszHdrLines, "IMAGE_LINES"));
742 : }
743 :
744 29 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
745 : {
746 0 : return nullptr;
747 : }
748 :
749 29 : RawRasterBand::ByteOrder eByteOrder = RawRasterBand::NATIVE_BYTE_ORDER;
750 :
751 29 : const char *pszByteOrder = CSLFetchNameValue(papszHdrLines, "BYTE_ORDER");
752 29 : if (pszByteOrder)
753 : {
754 25 : eByteOrder = EQUAL(pszByteOrder, "LSB")
755 25 : ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
756 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN;
757 : }
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Get some information specific to APP tiled files. */
761 : /* -------------------------------------------------------------------- */
762 29 : int nTileXSize = 0;
763 29 : int nTileYSize = 0;
764 29 : const char *pszRefinedType = CSLFetchNameValue(papszHdrLines, "type");
765 29 : const bool bTiled = CSLFetchNameValue(papszHdrLines, "no_rows") != nullptr;
766 :
767 29 : if (bTiled)
768 : {
769 2 : if (CSLFetchNameValue(papszHdrLines, "tile_size_rows"))
770 2 : nTileYSize =
771 2 : atoi(CSLFetchNameValue(papszHdrLines, "tile_size_rows"));
772 2 : if (CSLFetchNameValue(papszHdrLines, "tile_size_columns"))
773 2 : nTileXSize =
774 2 : atoi(CSLFetchNameValue(papszHdrLines, "tile_size_columns"));
775 :
776 2 : if (nTileXSize <= 0 || nTileYSize <= 0 ||
777 6 : poDS->nRasterXSize - 1 > INT_MAX - nTileXSize ||
778 2 : poDS->nRasterYSize - 1 > INT_MAX - nTileYSize)
779 : {
780 0 : return nullptr;
781 : }
782 : }
783 :
784 : /* -------------------------------------------------------------------- */
785 : /* Read the directory to find matching band files. */
786 : /* -------------------------------------------------------------------- */
787 : char *const pszTargetPath =
788 29 : CPLStrdup(CPLGetPathSafe(poOpenInfo->pszFilename).c_str());
789 : char *const pszTargetBase =
790 29 : CPLStrdup(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str());
791 : char **papszDirFiles =
792 29 : VSIReadDir(CPLGetPathSafe(poOpenInfo->pszFilename).c_str());
793 29 : if (papszDirFiles == nullptr)
794 : {
795 0 : CPLFree(pszTargetPath);
796 0 : CPLFree(pszTargetBase);
797 0 : return nullptr;
798 : }
799 :
800 29 : int nSkipped = 0;
801 29 : for (int nRawBand = 0; true; nRawBand++)
802 : {
803 : /* Find the next raw band file. */
804 :
805 86 : int i = 0; // Used after for.
806 352 : for (; papszDirFiles[i] != nullptr; i++)
807 : {
808 323 : if (!EQUAL(CPLGetBasenameSafe(papszDirFiles[i]).c_str(),
809 : pszTargetBase))
810 84 : continue;
811 :
812 : const std::string osExtension =
813 239 : CPLGetExtensionSafe(papszDirFiles[i]);
814 239 : if (osExtension.size() >= 2 &&
815 239 : isdigit(static_cast<unsigned char>(osExtension[1])) &&
816 535 : atoi(osExtension.c_str() + 1) == nRawBand &&
817 57 : strchr("bBcCiIjJrRxXzZ", osExtension[0]) != nullptr)
818 57 : break;
819 : }
820 :
821 86 : if (papszDirFiles[i] == nullptr)
822 29 : break;
823 :
824 : /* open the file for required level of access */
825 : const std::string osRawFilename =
826 57 : CPLFormFilenameSafe(pszTargetPath, papszDirFiles[i], nullptr);
827 :
828 57 : VSILFILE *fpRaw = nullptr;
829 57 : if (poOpenInfo->eAccess == GA_Update)
830 50 : fpRaw = VSIFOpenL(osRawFilename.c_str(), "rb+");
831 : else
832 7 : fpRaw = VSIFOpenL(osRawFilename.c_str(), "rb");
833 :
834 57 : if (fpRaw == nullptr)
835 : {
836 0 : CPLError(CE_Warning, CPLE_OpenFailed,
837 : "Unable to open %s ... skipping.", osRawFilename.c_str());
838 0 : nSkipped++;
839 0 : continue;
840 : }
841 114 : poDS->m_papszFileList =
842 57 : CSLAddString(poDS->m_papszFileList, osRawFilename.c_str());
843 :
844 57 : GDALDataType eDataType = GDT_Unknown;
845 57 : const std::string osExt = CPLGetExtensionSafe(papszDirFiles[i]);
846 57 : if (pszRefinedType != nullptr)
847 : {
848 0 : if (EQUAL(pszRefinedType, "C*4"))
849 0 : eDataType = GDT_CFloat32;
850 0 : else if (EQUAL(pszRefinedType, "C*8"))
851 0 : eDataType = GDT_CFloat64;
852 0 : else if (EQUAL(pszRefinedType, "R*4"))
853 0 : eDataType = GDT_Float32;
854 0 : else if (EQUAL(pszRefinedType, "R*8"))
855 0 : eDataType = GDT_Float64;
856 0 : else if (EQUAL(pszRefinedType, "I*1"))
857 0 : eDataType = GDT_Byte;
858 0 : else if (EQUAL(pszRefinedType, "I*2"))
859 0 : eDataType = GDT_Int16;
860 0 : else if (EQUAL(pszRefinedType, "I*4"))
861 0 : eDataType = GDT_Int32;
862 0 : else if (EQUAL(pszRefinedType, "U*2"))
863 0 : eDataType = GDT_UInt16;
864 0 : else if (EQUAL(pszRefinedType, "U*4"))
865 0 : eDataType = GDT_UInt32;
866 0 : else if (EQUAL(pszRefinedType, "J*1"))
867 : {
868 0 : CPLError(
869 : CE_Warning, CPLE_OpenFailed,
870 : "Unable to open band %d because type J*1 is not handled. "
871 : "Skipping.",
872 : nRawBand + 1);
873 0 : nSkipped++;
874 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
875 0 : continue; // Does not support 1 byte complex.
876 : }
877 0 : else if (EQUAL(pszRefinedType, "J*2"))
878 0 : eDataType = GDT_CInt16;
879 0 : else if (EQUAL(pszRefinedType, "K*4"))
880 0 : eDataType = GDT_CInt32;
881 : else
882 : {
883 0 : CPLError(
884 : CE_Warning, CPLE_OpenFailed,
885 : "Unable to open band %d because type %s is not handled. "
886 : "Skipping.\n",
887 : nRawBand + 1, pszRefinedType);
888 0 : nSkipped++;
889 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
890 0 : continue;
891 : }
892 : }
893 57 : else if (STARTS_WITH_CI(osExt.c_str(), "b"))
894 : {
895 36 : eDataType = GDT_Byte;
896 : }
897 21 : else if (STARTS_WITH_CI(osExt.c_str(), "i"))
898 : {
899 5 : eDataType = GDT_UInt16;
900 : }
901 16 : else if (STARTS_WITH_CI(osExt.c_str(), "j"))
902 : {
903 5 : eDataType = GDT_CInt16;
904 : }
905 11 : else if (STARTS_WITH_CI(osExt.c_str(), "r"))
906 : {
907 6 : eDataType = GDT_Float32;
908 : }
909 5 : else if (STARTS_WITH_CI(osExt.c_str(), "x"))
910 : {
911 5 : eDataType = GDT_CFloat32;
912 : }
913 : else
914 : {
915 0 : CPLError(CE_Warning, CPLE_OpenFailed,
916 : "Unable to open band %d because extension %s is not "
917 : "handled. Skipping.",
918 : nRawBand + 1, osExt.c_str());
919 0 : nSkipped++;
920 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
921 0 : continue;
922 : }
923 :
924 57 : const int nBand = poDS->GetRasterCount() + 1;
925 :
926 57 : const int nPixelOffset = GDALGetDataTypeSizeBytes(eDataType);
927 0 : std::unique_ptr<GDALRasterBand> poBand;
928 :
929 57 : if (bTiled)
930 : {
931 4 : poBand = std::make_unique<MFFTiledBand>(poDS.get(), nBand, fpRaw,
932 : nTileXSize, nTileYSize,
933 2 : eDataType, eByteOrder);
934 : }
935 : else
936 : {
937 110 : if (nPixelOffset != 0 &&
938 55 : poDS->GetRasterXSize() > INT_MAX / nPixelOffset)
939 : {
940 0 : CPLError(CE_Warning, CPLE_AppDefined,
941 : "Int overflow occurred... skipping");
942 0 : nSkipped++;
943 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
944 0 : continue;
945 : }
946 :
947 110 : poBand = RawRasterBand::Create(
948 55 : poDS.get(), nBand, fpRaw, 0, nPixelOffset,
949 55 : nPixelOffset * poDS->GetRasterXSize(), eDataType, eByteOrder,
950 55 : RawRasterBand::OwnFP::YES);
951 : }
952 :
953 57 : poDS->SetBand(nBand, std::move(poBand));
954 57 : }
955 :
956 29 : CPLFree(pszTargetPath);
957 29 : CPLFree(pszTargetBase);
958 29 : CSLDestroy(papszDirFiles);
959 :
960 : /* -------------------------------------------------------------------- */
961 : /* Check if we have bands. */
962 : /* -------------------------------------------------------------------- */
963 29 : if (poDS->GetRasterCount() == 0)
964 : {
965 0 : if (nSkipped > 0 && poOpenInfo->eAccess)
966 : {
967 0 : CPLError(CE_Failure, CPLE_OpenFailed,
968 : "Failed to open %d files that were apparently bands. "
969 : "Perhaps this dataset is readonly?",
970 : nSkipped);
971 0 : return nullptr;
972 : }
973 : else
974 : {
975 0 : CPLError(CE_Failure, CPLE_OpenFailed,
976 : "MFF header file read successfully, but no bands "
977 : "were successfully found and opened.");
978 0 : return nullptr;
979 : }
980 : }
981 :
982 : /* -------------------------------------------------------------------- */
983 : /* Set all information from the .hdr that isn't well know to be */
984 : /* metadata. */
985 : /* -------------------------------------------------------------------- */
986 226 : for (int i = 0; papszHdrLines[i] != nullptr; i++)
987 : {
988 197 : char *pszName = nullptr;
989 :
990 197 : const char *pszValue = CPLParseNameValue(papszHdrLines[i], &pszName);
991 197 : if (pszName == nullptr || pszValue == nullptr)
992 16 : continue;
993 :
994 181 : if (!EQUAL(pszName, "END") && !EQUAL(pszName, "FILE_TYPE") &&
995 156 : !EQUAL(pszName, "BYTE_ORDER") && !EQUAL(pszName, "no_columns") &&
996 131 : !EQUAL(pszName, "no_rows") && !EQUAL(pszName, "type") &&
997 129 : !EQUAL(pszName, "tile_size_rows") &&
998 127 : !EQUAL(pszName, "tile_size_columns") &&
999 125 : !EQUAL(pszName, "IMAGE_FILE_FORMAT") &&
1000 96 : !EQUAL(pszName, "IMAGE_LINES") && !EQUAL(pszName, "LINE_SAMPLES"))
1001 : {
1002 38 : poDS->SetMetadataItem(pszName, pszValue);
1003 : }
1004 :
1005 181 : CPLFree(pszName);
1006 : }
1007 :
1008 : /* -------------------------------------------------------------------- */
1009 : /* Any GCPs in header file? */
1010 : /* -------------------------------------------------------------------- */
1011 29 : poDS->ScanForGCPs();
1012 29 : poDS->ScanForProjectionInfo();
1013 29 : if (poDS->nGCPCount == 0)
1014 23 : poDS->m_oGCPSRS.Clear();
1015 :
1016 : /* -------------------------------------------------------------------- */
1017 : /* Initialize any PAM information. */
1018 : /* -------------------------------------------------------------------- */
1019 29 : poDS->SetDescription(poOpenInfo->pszFilename);
1020 29 : poDS->TryLoadXML();
1021 :
1022 : /* -------------------------------------------------------------------- */
1023 : /* Check for overviews. */
1024 : /* -------------------------------------------------------------------- */
1025 29 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1026 :
1027 29 : return poDS.release();
1028 : }
1029 :
1030 9 : int GetMFFProjectionType(const OGRSpatialReference *poSRS)
1031 : {
1032 9 : if (poSRS == nullptr)
1033 : {
1034 0 : return MFFPRJ_NONE;
1035 : }
1036 9 : if (poSRS->IsProjected() && poSRS->GetAttrValue("PROJECTION") &&
1037 0 : EQUAL(poSRS->GetAttrValue("PROJECTION"), SRS_PT_TRANSVERSE_MERCATOR))
1038 : {
1039 0 : return MFFPRJ_UTM;
1040 : }
1041 9 : else if (poSRS->IsGeographic())
1042 : {
1043 9 : return MFFPRJ_LL;
1044 : }
1045 : else
1046 : {
1047 0 : return MFFPRJ_UNRECOGNIZED;
1048 : }
1049 : }
1050 :
1051 : /************************************************************************/
1052 : /* Create() */
1053 : /************************************************************************/
1054 :
1055 50 : GDALDataset *MFFDataset::Create(const char *pszFilenameIn, int nXSize,
1056 : int nYSize, int nBandsIn, GDALDataType eType,
1057 : char **papszParamList)
1058 :
1059 : {
1060 : /* -------------------------------------------------------------------- */
1061 : /* Verify input options. */
1062 : /* -------------------------------------------------------------------- */
1063 50 : if (nBandsIn <= 0)
1064 : {
1065 1 : CPLError(CE_Failure, CPLE_NotSupported,
1066 : "MFF driver does not support %d bands.", nBandsIn);
1067 1 : return nullptr;
1068 : }
1069 :
1070 49 : if (eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 &&
1071 27 : eType != GDT_CInt16 && eType != GDT_CFloat32)
1072 : {
1073 24 : CPLError(CE_Failure, CPLE_AppDefined,
1074 : "Attempt to create MFF file with currently unsupported\n"
1075 : "data type (%s).\n",
1076 : GDALGetDataTypeName(eType));
1077 :
1078 24 : return nullptr;
1079 : }
1080 :
1081 : /* -------------------------------------------------------------------- */
1082 : /* Establish the base filename (path+filename, less extension). */
1083 : /* -------------------------------------------------------------------- */
1084 : char *pszBaseFilename =
1085 25 : static_cast<char *>(CPLMalloc(strlen(pszFilenameIn) + 5));
1086 25 : strcpy(pszBaseFilename, pszFilenameIn);
1087 :
1088 100 : for (int i = static_cast<int>(strlen(pszBaseFilename)) - 1; i > 0; i--)
1089 : {
1090 100 : if (pszBaseFilename[i] == '.')
1091 : {
1092 0 : pszBaseFilename[i] = '\0';
1093 0 : break;
1094 : }
1095 :
1096 100 : if (pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\')
1097 : break;
1098 : }
1099 :
1100 : /* -------------------------------------------------------------------- */
1101 : /* Create the header file. */
1102 : /* -------------------------------------------------------------------- */
1103 : std::string osFilename =
1104 50 : CPLFormFilenameSafe(nullptr, pszBaseFilename, "hdr");
1105 :
1106 25 : VSILFILE *fp = VSIFOpenL(osFilename.c_str(), "wt");
1107 25 : if (fp == nullptr)
1108 : {
1109 3 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.\n",
1110 : osFilename.c_str());
1111 3 : CPLFree(pszBaseFilename);
1112 3 : return nullptr;
1113 : }
1114 :
1115 22 : bool bOK = VSIFPrintfL(fp, "IMAGE_FILE_FORMAT = MFF\n") >= 0;
1116 22 : bOK &= VSIFPrintfL(fp, "FILE_TYPE = IMAGE\n") >= 0;
1117 22 : bOK &= VSIFPrintfL(fp, "IMAGE_LINES = %d\n", nYSize) >= 0;
1118 22 : bOK &= VSIFPrintfL(fp, "LINE_SAMPLES = %d\n", nXSize) >= 0;
1119 : #ifdef CPL_MSB
1120 : bOK &= VSIFPrintfL(fp, "BYTE_ORDER = MSB\n") >= 0;
1121 : #else
1122 22 : bOK &= VSIFPrintfL(fp, "BYTE_ORDER = LSB\n") >= 0;
1123 : #endif
1124 :
1125 22 : if (CSLFetchNameValue(papszParamList, "NO_END") == nullptr)
1126 13 : bOK &= VSIFPrintfL(fp, "END\n") >= 0;
1127 :
1128 22 : if (VSIFCloseL(fp) != 0)
1129 0 : bOK = false;
1130 :
1131 : /* -------------------------------------------------------------------- */
1132 : /* Create the data files, but don't bother writing any data to them.*/
1133 : /* -------------------------------------------------------------------- */
1134 72 : for (int iBand = 0; bOK && iBand < nBandsIn; iBand++)
1135 : {
1136 50 : char szExtension[4] = {'\0'};
1137 :
1138 50 : if (eType == GDT_Byte)
1139 30 : CPLsnprintf(szExtension, sizeof(szExtension), "b%02d", iBand);
1140 20 : else if (eType == GDT_UInt16)
1141 5 : CPLsnprintf(szExtension, sizeof(szExtension), "i%02d", iBand);
1142 15 : else if (eType == GDT_Float32)
1143 5 : CPLsnprintf(szExtension, sizeof(szExtension), "r%02d", iBand);
1144 10 : else if (eType == GDT_CInt16)
1145 5 : CPLsnprintf(szExtension, sizeof(szExtension), "j%02d", iBand);
1146 5 : else if (eType == GDT_CFloat32)
1147 5 : CPLsnprintf(szExtension, sizeof(szExtension), "x%02d", iBand);
1148 :
1149 50 : osFilename = CPLFormFilenameSafe(nullptr, pszBaseFilename, szExtension);
1150 50 : fp = VSIFOpenL(osFilename.c_str(), "wb");
1151 50 : if (fp == nullptr)
1152 : {
1153 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.\n",
1154 : osFilename.c_str());
1155 0 : CPLFree(pszBaseFilename);
1156 0 : return nullptr;
1157 : }
1158 :
1159 50 : bOK &= VSIFWriteL("", 1, 1, fp) == 1;
1160 50 : if (VSIFCloseL(fp) != 0)
1161 0 : bOK = false;
1162 : }
1163 :
1164 22 : if (!bOK)
1165 : {
1166 0 : CPLFree(pszBaseFilename);
1167 0 : return nullptr;
1168 : }
1169 :
1170 : /* -------------------------------------------------------------------- */
1171 : /* Open the dataset normally. */
1172 : /* -------------------------------------------------------------------- */
1173 22 : strcat(pszBaseFilename, ".hdr");
1174 : GDALDataset *poDS =
1175 22 : GDALDataset::Open(pszBaseFilename, GDAL_OF_RASTER | GDAL_OF_UPDATE);
1176 22 : CPLFree(pszBaseFilename);
1177 :
1178 22 : return poDS;
1179 : }
1180 :
1181 : /************************************************************************/
1182 : /* CreateCopy() */
1183 : /************************************************************************/
1184 :
1185 19 : GDALDataset *MFFDataset::CreateCopy(const char *pszFilename,
1186 : GDALDataset *poSrcDS, int /* bStrict */,
1187 : char **papszOptions,
1188 : GDALProgressFunc pfnProgress,
1189 : void *pProgressData)
1190 : {
1191 19 : const int nBands = poSrcDS->GetRasterCount();
1192 19 : if (nBands == 0)
1193 : {
1194 1 : CPLError(CE_Failure, CPLE_NotSupported,
1195 : "MFF driver does not support source dataset with zero band.");
1196 1 : return nullptr;
1197 : }
1198 :
1199 18 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1200 18 : if (!pfnProgress(0.0, nullptr, pProgressData))
1201 0 : return nullptr;
1202 :
1203 : // Check that other bands match type- sets type
1204 : // to unknown if they differ.
1205 28 : for (int iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++)
1206 : {
1207 10 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
1208 10 : eType = GDALDataTypeUnion(eType, poBand->GetRasterDataType());
1209 : }
1210 :
1211 18 : char **newpapszOptions = CSLDuplicate(papszOptions);
1212 18 : newpapszOptions = CSLSetNameValue(newpapszOptions, "NO_END", "TRUE");
1213 :
1214 18 : MFFDataset *poDS = reinterpret_cast<MFFDataset *>(Create(
1215 : pszFilename, poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
1216 : poSrcDS->GetRasterCount(), eType, newpapszOptions));
1217 :
1218 18 : CSLDestroy(newpapszOptions);
1219 :
1220 18 : if (poDS == nullptr)
1221 9 : return nullptr;
1222 :
1223 : /* -------------------------------------------------------------------- */
1224 : /* Copy the image data. */
1225 : /* -------------------------------------------------------------------- */
1226 9 : const int nXSize = poDS->GetRasterXSize();
1227 9 : const int nYSize = poDS->GetRasterYSize();
1228 :
1229 9 : int nBlockXSize = 0;
1230 9 : int nBlockYSize = 0;
1231 9 : poDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
1232 :
1233 9 : const int nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize) *
1234 9 : ((nYSize + nBlockYSize - 1) / nBlockYSize) *
1235 9 : poSrcDS->GetRasterCount();
1236 :
1237 9 : int nBlocksDone = 0;
1238 28 : for (int iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++)
1239 : {
1240 19 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand + 1);
1241 19 : GDALRasterBand *poDstBand = poDS->GetRasterBand(iBand + 1);
1242 :
1243 57 : void *pData = CPLMalloc(static_cast<size_t>(nBlockXSize) * nBlockYSize *
1244 19 : GDALGetDataTypeSizeBytes(eType));
1245 :
1246 209 : for (int iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize)
1247 : {
1248 380 : for (int iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize)
1249 : {
1250 190 : if (!pfnProgress((nBlocksDone++) /
1251 190 : static_cast<float>(nBlockTotal),
1252 : nullptr, pProgressData))
1253 : {
1254 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1255 0 : delete poDS;
1256 0 : CPLFree(pData);
1257 :
1258 : GDALDriver *poMFFDriver =
1259 0 : static_cast<GDALDriver *>(GDALGetDriverByName("MFF"));
1260 0 : poMFFDriver->Delete(pszFilename);
1261 0 : return nullptr;
1262 : }
1263 :
1264 190 : const int nTBXSize = std::min(nBlockXSize, nXSize - iXOffset);
1265 190 : const int nTBYSize = std::min(nBlockYSize, nYSize - iYOffset);
1266 :
1267 190 : CPLErr eErr = poSrcBand->RasterIO(
1268 : GF_Read, iXOffset, iYOffset, nTBXSize, nTBYSize, pData,
1269 : nTBXSize, nTBYSize, eType, 0, 0, nullptr);
1270 :
1271 190 : if (eErr != CE_None)
1272 : {
1273 0 : delete poDS;
1274 0 : CPLFree(pData);
1275 0 : return nullptr;
1276 : }
1277 :
1278 190 : eErr = poDstBand->RasterIO(GF_Write, iXOffset, iYOffset,
1279 : nTBXSize, nTBYSize, pData, nTBXSize,
1280 : nTBYSize, eType, 0, 0, nullptr);
1281 :
1282 190 : if (eErr != CE_None)
1283 : {
1284 0 : delete poDS;
1285 0 : CPLFree(pData);
1286 0 : return nullptr;
1287 : }
1288 : }
1289 : }
1290 :
1291 19 : CPLFree(pData);
1292 : }
1293 :
1294 : /* -------------------------------------------------------------------- */
1295 : /* Copy georeferencing information, if enough is available. */
1296 : /* -------------------------------------------------------------------- */
1297 :
1298 : /* -------------------------------------------------------------------- */
1299 : /* Establish the base filename (path+filename, less extension). */
1300 : /* -------------------------------------------------------------------- */
1301 : char *pszBaseFilename =
1302 9 : static_cast<char *>(CPLMalloc(strlen(pszFilename) + 5));
1303 9 : strcpy(pszBaseFilename, pszFilename);
1304 :
1305 36 : for (int i = static_cast<int>(strlen(pszBaseFilename)) - 1; i > 0; i--)
1306 : {
1307 36 : if (pszBaseFilename[i] == '.')
1308 : {
1309 0 : pszBaseFilename[i] = '\0';
1310 0 : break;
1311 : }
1312 :
1313 36 : if (pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\')
1314 : break;
1315 : }
1316 :
1317 : const std::string osFilenameGEO =
1318 18 : CPLFormFilenameSafe(nullptr, pszBaseFilename, "hdr");
1319 :
1320 9 : VSILFILE *fp = VSIFOpenL(osFilenameGEO.c_str(), "at");
1321 9 : if (fp == nullptr)
1322 : {
1323 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1324 : "Couldn't open %s for appending.\n", osFilenameGEO.c_str());
1325 0 : CPLFree(pszBaseFilename);
1326 0 : return nullptr;
1327 : }
1328 :
1329 : /* MFF requires corner and center gcps */
1330 9 : bool georef_created = false;
1331 :
1332 : double *padfTiepoints =
1333 9 : static_cast<double *>(CPLMalloc(2 * sizeof(double) * 5));
1334 :
1335 9 : const int src_prj = GetMFFProjectionType(poSrcDS->GetSpatialRef());
1336 :
1337 9 : if ((src_prj != MFFPRJ_NONE) && (src_prj != MFFPRJ_UNRECOGNIZED))
1338 : {
1339 : double *tempGeoTransform =
1340 9 : static_cast<double *>(CPLMalloc(6 * sizeof(double)));
1341 :
1342 18 : if ((poSrcDS->GetGeoTransform(tempGeoTransform) == CE_None) &&
1343 9 : (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0 ||
1344 0 : tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0 ||
1345 0 : tempGeoTransform[4] != 0.0 ||
1346 0 : std::abs(tempGeoTransform[5]) != 1.0))
1347 : {
1348 9 : padfTiepoints[0] = tempGeoTransform[0] + tempGeoTransform[1] * 0.5 +
1349 9 : tempGeoTransform[2] * 0.5;
1350 :
1351 9 : padfTiepoints[1] = tempGeoTransform[3] + tempGeoTransform[4] * 0.5 +
1352 9 : tempGeoTransform[5] * 0.5;
1353 :
1354 9 : padfTiepoints[2] =
1355 18 : tempGeoTransform[0] + tempGeoTransform[2] * 0.5 +
1356 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize() - 0.5);
1357 :
1358 9 : padfTiepoints[3] =
1359 18 : tempGeoTransform[3] + tempGeoTransform[5] * 0.5 +
1360 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize() - 0.5);
1361 :
1362 9 : padfTiepoints[4] =
1363 18 : tempGeoTransform[0] + tempGeoTransform[1] * 0.5 +
1364 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize() - 0.5);
1365 :
1366 9 : padfTiepoints[5] =
1367 18 : tempGeoTransform[3] + tempGeoTransform[4] * 0.5 +
1368 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize() - 0.5);
1369 :
1370 9 : padfTiepoints[6] =
1371 18 : tempGeoTransform[0] +
1372 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize() - 0.5) +
1373 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize() - 0.5);
1374 :
1375 9 : padfTiepoints[7] =
1376 18 : tempGeoTransform[3] +
1377 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize() - 0.5) +
1378 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize() - 0.5);
1379 :
1380 9 : padfTiepoints[8] =
1381 18 : tempGeoTransform[0] +
1382 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize()) / 2.0 +
1383 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize()) / 2.0;
1384 :
1385 9 : padfTiepoints[9] =
1386 18 : tempGeoTransform[3] +
1387 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize()) / 2.0 +
1388 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize()) / 2.0;
1389 :
1390 18 : OGRSpatialReference oUTMorLL;
1391 9 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
1392 9 : if (poSrcSRS)
1393 9 : oUTMorLL = *poSrcSRS;
1394 9 : auto poLLSRS = oUTMorLL.CloneGeogCS();
1395 9 : if (poLLSRS && oUTMorLL.IsProjected())
1396 : {
1397 0 : poLLSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1398 : OGRCoordinateTransformation *poTransform =
1399 0 : OGRCreateCoordinateTransformation(&oUTMorLL, poLLSRS);
1400 :
1401 : // projected coordinate system- need to translate gcps */
1402 0 : bool bSuccess = poTransform != nullptr;
1403 :
1404 0 : for (int index = 0; index < 5; index++)
1405 : {
1406 : // TODO: If bSuccess is false, set it to false?
1407 0 : if (!bSuccess || !poTransform->Transform(
1408 : 1, &(padfTiepoints[index * 2]),
1409 : &(padfTiepoints[index * 2 + 1])))
1410 0 : bSuccess = false;
1411 : }
1412 0 : if (bSuccess)
1413 0 : georef_created = true;
1414 : }
1415 : else
1416 : {
1417 9 : georef_created = true;
1418 : }
1419 9 : delete poLLSRS;
1420 : }
1421 9 : CPLFree(tempGeoTransform);
1422 : }
1423 :
1424 9 : bool bOK = true;
1425 9 : if (georef_created)
1426 : {
1427 : /* --------------------------------------------------------------------
1428 : */
1429 : /* top left */
1430 : /* --------------------------------------------------------------------
1431 : */
1432 18 : bOK &= VSIFPrintfL(fp, "TOP_LEFT_CORNER_LATITUDE = %.10f\n",
1433 9 : padfTiepoints[1]) >= 0;
1434 9 : bOK &= VSIFPrintfL(fp, "TOP_LEFT_CORNER_LONGITUDE = %.10f\n",
1435 9 : padfTiepoints[0]) >= 0;
1436 : /* --------------------------------------------------------------------
1437 : */
1438 : /* top_right */
1439 : /* --------------------------------------------------------------------
1440 : */
1441 18 : bOK &= VSIFPrintfL(fp, "TOP_RIGHT_CORNER_LATITUDE = %.10f\n",
1442 9 : padfTiepoints[3]) >= 0;
1443 18 : bOK &= VSIFPrintfL(fp, "TOP_RIGHT_CORNER_LONGITUDE = %.10f\n",
1444 9 : padfTiepoints[2]) >= 0;
1445 : /* --------------------------------------------------------------------
1446 : */
1447 : /* bottom_left */
1448 : /* --------------------------------------------------------------------
1449 : */
1450 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_LEFT_CORNER_LATITUDE = %.10f\n",
1451 9 : padfTiepoints[5]) >= 0;
1452 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_LEFT_CORNER_LONGITUDE = %.10f\n",
1453 9 : padfTiepoints[4]) >= 0;
1454 : /* --------------------------------------------------------------------
1455 : */
1456 : /* bottom_right */
1457 : /* --------------------------------------------------------------------
1458 : */
1459 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_RIGHT_CORNER_LATITUDE = %.10f\n",
1460 9 : padfTiepoints[7]) >= 0;
1461 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_RIGHT_CORNER_LONGITUDE = %.10f\n",
1462 9 : padfTiepoints[6]) >= 0;
1463 : /* --------------------------------------------------------------------
1464 : */
1465 : /* Center */
1466 : /* --------------------------------------------------------------------
1467 : */
1468 9 : bOK &=
1469 9 : VSIFPrintfL(fp, "CENTRE_LATITUDE = %.10f\n", padfTiepoints[9]) >= 0;
1470 18 : bOK &= VSIFPrintfL(fp, "CENTRE_LONGITUDE = %.10f\n",
1471 9 : padfTiepoints[8]) >= 0;
1472 : /* -------------------------------------------------------------------
1473 : */
1474 : /* Ellipsoid/projection */
1475 : /* --------------------------------------------------------------------*/
1476 :
1477 9 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
1478 9 : char *spheroid_name = nullptr;
1479 :
1480 9 : if (poSrcSRS != nullptr)
1481 : {
1482 9 : if (poSrcSRS->IsProjected() &&
1483 9 : poSrcSRS->GetAttrValue("PROJECTION") != nullptr &&
1484 0 : EQUAL(poSrcSRS->GetAttrValue("PROJECTION"),
1485 : SRS_PT_TRANSVERSE_MERCATOR))
1486 : {
1487 0 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = UTM\n") >= 0;
1488 0 : OGRErr ogrerrorOl = OGRERR_NONE;
1489 0 : bOK &=
1490 0 : VSIFPrintfL(fp, "PROJECTION_ORIGIN_LONGITUDE = %f\n",
1491 : poSrcSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,
1492 0 : 0.0, &ogrerrorOl)) >= 0;
1493 : }
1494 9 : else if (poSrcSRS->IsGeographic())
1495 : {
1496 9 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = LL\n") >= 0;
1497 : }
1498 : else
1499 : {
1500 0 : CPLError(CE_Warning, CPLE_AppDefined,
1501 : "Unrecognized projection- no georeferencing "
1502 : "information transferred.");
1503 0 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = LL\n") >= 0;
1504 : }
1505 9 : OGRErr ogrerrorEq = OGRERR_NONE;
1506 9 : const double eq_radius = poSrcSRS->GetSemiMajor(&ogrerrorEq);
1507 9 : OGRErr ogrerrorInvf = OGRERR_NONE;
1508 : const double inv_flattening =
1509 9 : poSrcSRS->GetInvFlattening(&ogrerrorInvf);
1510 9 : if (ogrerrorEq == OGRERR_NONE && ogrerrorInvf == OGRERR_NONE)
1511 : {
1512 9 : MFFSpheroidList *mffEllipsoids = new MFFSpheroidList;
1513 : spheroid_name =
1514 9 : mffEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(
1515 : eq_radius, inv_flattening);
1516 9 : if (spheroid_name != nullptr)
1517 : {
1518 0 : bOK &= VSIFPrintfL(fp, "SPHEROID_NAME = %s\n",
1519 0 : spheroid_name) >= 0;
1520 : }
1521 : else
1522 : {
1523 18 : bOK &= VSIFPrintfL(fp,
1524 : "SPHEROID_NAME = USER_DEFINED\n"
1525 : "SPHEROID_EQUATORIAL_RADIUS = %.10f\n"
1526 : "SPHEROID_POLAR_RADIUS = %.10f\n",
1527 : eq_radius,
1528 : eq_radius *
1529 9 : (1 - 1.0 / inv_flattening)) >= 0;
1530 : }
1531 9 : delete mffEllipsoids;
1532 9 : CPLFree(spheroid_name);
1533 : }
1534 : }
1535 : }
1536 :
1537 9 : CPLFree(padfTiepoints);
1538 9 : bOK &= VSIFPrintfL(fp, "END\n") >= 0;
1539 9 : if (VSIFCloseL(fp) != 0)
1540 0 : bOK = false;
1541 :
1542 9 : if (!bOK)
1543 : {
1544 0 : delete poDS;
1545 0 : CPLFree(pszBaseFilename);
1546 0 : return nullptr;
1547 : }
1548 :
1549 : /* End of georeferencing stuff */
1550 :
1551 : /* Make sure image data gets flushed */
1552 28 : for (int iBand = 0; iBand < poDS->GetRasterCount(); iBand++)
1553 : {
1554 : RawRasterBand *poDstBand =
1555 19 : reinterpret_cast<RawRasterBand *>(poDS->GetRasterBand(iBand + 1));
1556 19 : poDstBand->FlushCache(false);
1557 : }
1558 :
1559 9 : if (!pfnProgress(1.0, nullptr, pProgressData))
1560 : {
1561 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1562 0 : delete poDS;
1563 :
1564 : GDALDriver *poMFFDriver =
1565 0 : static_cast<GDALDriver *>(GDALGetDriverByName("MFF"));
1566 0 : poMFFDriver->Delete(pszFilename);
1567 0 : CPLFree(pszBaseFilename);
1568 0 : return nullptr;
1569 : }
1570 :
1571 9 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1572 9 : CPLFree(pszBaseFilename);
1573 :
1574 9 : return poDS;
1575 : }
1576 :
1577 : /************************************************************************/
1578 : /* GDALRegister_MFF() */
1579 : /************************************************************************/
1580 :
1581 1682 : void GDALRegister_MFF()
1582 :
1583 : {
1584 1682 : if (GDALGetDriverByName("MFF") != nullptr)
1585 301 : return;
1586 :
1587 1381 : GDALDriver *poDriver = new GDALDriver();
1588 :
1589 1381 : poDriver->SetDescription("MFF");
1590 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1591 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF Raster");
1592 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff.html");
1593 1381 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
1594 1381 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1595 1381 : "Byte UInt16 Float32 CInt16 CFloat32");
1596 :
1597 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1598 :
1599 1381 : poDriver->pfnOpen = MFFDataset::Open;
1600 1381 : poDriver->pfnCreate = MFFDataset::Create;
1601 1381 : poDriver->pfnCreateCopy = MFFDataset::CreateCopy;
1602 :
1603 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
1604 : }
|
