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 30723 : GDALDataset *MFFDataset::Open(GDALOpenInfo *poOpenInfo)
662 :
663 : {
664 : /* -------------------------------------------------------------------- */
665 : /* We assume the user is pointing to the header file. */
666 : /* -------------------------------------------------------------------- */
667 30723 : if (poOpenInfo->nHeaderBytes < 17 || poOpenInfo->fpL == nullptr)
668 27422 : return nullptr;
669 :
670 3301 : if (!EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hdr"))
671 3262 : 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 29 : char *const pszTargetPath = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
788 : char *const pszTargetBase =
789 29 : CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
790 29 : char **papszDirFiles = VSIReadDir(CPLGetPath(poOpenInfo->pszFilename));
791 29 : if (papszDirFiles == nullptr)
792 : {
793 0 : CPLFree(pszTargetPath);
794 0 : CPLFree(pszTargetBase);
795 0 : return nullptr;
796 : }
797 :
798 29 : int nSkipped = 0;
799 29 : for (int nRawBand = 0; true; nRawBand++)
800 : {
801 86 : const char *pszExtension = nullptr;
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(CPLGetBasename(papszDirFiles[i]), pszTargetBase))
809 84 : continue;
810 :
811 239 : pszExtension = CPLGetExtension(papszDirFiles[i]);
812 239 : if (strlen(pszExtension) >= 2 &&
813 239 : isdigit(static_cast<unsigned char>(pszExtension[1])) &&
814 166 : atoi(pszExtension + 1) == nRawBand &&
815 57 : strchr("bBcCiIjJrRxXzZ", pszExtension[0]) != nullptr)
816 57 : break;
817 : }
818 :
819 86 : if (papszDirFiles[i] == nullptr)
820 29 : break;
821 :
822 : /* open the file for required level of access */
823 : const char *pszRawFilename =
824 57 : CPLFormFilename(pszTargetPath, papszDirFiles[i], nullptr);
825 :
826 57 : VSILFILE *fpRaw = nullptr;
827 57 : if (poOpenInfo->eAccess == GA_Update)
828 50 : fpRaw = VSIFOpenL(pszRawFilename, "rb+");
829 : else
830 7 : fpRaw = VSIFOpenL(pszRawFilename, "rb");
831 :
832 57 : if (fpRaw == nullptr)
833 : {
834 0 : CPLError(CE_Warning, CPLE_OpenFailed,
835 : "Unable to open %s ... skipping.", pszRawFilename);
836 0 : nSkipped++;
837 0 : continue;
838 : }
839 114 : poDS->m_papszFileList =
840 57 : CSLAddString(poDS->m_papszFileList, pszRawFilename);
841 :
842 57 : GDALDataType eDataType = GDT_Unknown;
843 57 : pszExtension = CPLGetExtension(papszDirFiles[i]);
844 57 : if (pszRefinedType != nullptr)
845 : {
846 0 : if (EQUAL(pszRefinedType, "C*4"))
847 0 : eDataType = GDT_CFloat32;
848 0 : else if (EQUAL(pszRefinedType, "C*8"))
849 0 : eDataType = GDT_CFloat64;
850 0 : else if (EQUAL(pszRefinedType, "R*4"))
851 0 : eDataType = GDT_Float32;
852 0 : else if (EQUAL(pszRefinedType, "R*8"))
853 0 : eDataType = GDT_Float64;
854 0 : else if (EQUAL(pszRefinedType, "I*1"))
855 0 : eDataType = GDT_Byte;
856 0 : else if (EQUAL(pszRefinedType, "I*2"))
857 0 : eDataType = GDT_Int16;
858 0 : else if (EQUAL(pszRefinedType, "I*4"))
859 0 : eDataType = GDT_Int32;
860 0 : else if (EQUAL(pszRefinedType, "U*2"))
861 0 : eDataType = GDT_UInt16;
862 0 : else if (EQUAL(pszRefinedType, "U*4"))
863 0 : eDataType = GDT_UInt32;
864 0 : else if (EQUAL(pszRefinedType, "J*1"))
865 : {
866 0 : CPLError(
867 : CE_Warning, CPLE_OpenFailed,
868 : "Unable to open band %d because type J*1 is not handled. "
869 : "Skipping.",
870 : nRawBand + 1);
871 0 : nSkipped++;
872 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
873 0 : continue; // Does not support 1 byte complex.
874 : }
875 0 : else if (EQUAL(pszRefinedType, "J*2"))
876 0 : eDataType = GDT_CInt16;
877 0 : else if (EQUAL(pszRefinedType, "K*4"))
878 0 : eDataType = GDT_CInt32;
879 : else
880 : {
881 0 : CPLError(
882 : CE_Warning, CPLE_OpenFailed,
883 : "Unable to open band %d because type %s is not handled. "
884 : "Skipping.\n",
885 : nRawBand + 1, pszRefinedType);
886 0 : nSkipped++;
887 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
888 0 : continue;
889 : }
890 : }
891 57 : else if (STARTS_WITH_CI(pszExtension, "b"))
892 : {
893 36 : eDataType = GDT_Byte;
894 : }
895 21 : else if (STARTS_WITH_CI(pszExtension, "i"))
896 : {
897 5 : eDataType = GDT_UInt16;
898 : }
899 16 : else if (STARTS_WITH_CI(pszExtension, "j"))
900 : {
901 5 : eDataType = GDT_CInt16;
902 : }
903 11 : else if (STARTS_WITH_CI(pszExtension, "r"))
904 : {
905 6 : eDataType = GDT_Float32;
906 : }
907 5 : else if (STARTS_WITH_CI(pszExtension, "x"))
908 : {
909 5 : eDataType = GDT_CFloat32;
910 : }
911 : else
912 : {
913 0 : CPLError(CE_Warning, CPLE_OpenFailed,
914 : "Unable to open band %d because extension %s is not "
915 : "handled. Skipping.",
916 : nRawBand + 1, pszExtension);
917 0 : nSkipped++;
918 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
919 0 : continue;
920 : }
921 :
922 57 : const int nBand = poDS->GetRasterCount() + 1;
923 :
924 57 : const int nPixelOffset = GDALGetDataTypeSizeBytes(eDataType);
925 0 : std::unique_ptr<GDALRasterBand> poBand;
926 :
927 57 : if (bTiled)
928 : {
929 4 : poBand = std::make_unique<MFFTiledBand>(poDS.get(), nBand, fpRaw,
930 : nTileXSize, nTileYSize,
931 2 : eDataType, eByteOrder);
932 : }
933 : else
934 : {
935 110 : if (nPixelOffset != 0 &&
936 55 : poDS->GetRasterXSize() > INT_MAX / nPixelOffset)
937 : {
938 0 : CPLError(CE_Warning, CPLE_AppDefined,
939 : "Int overflow occurred... skipping");
940 0 : nSkipped++;
941 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpRaw));
942 0 : continue;
943 : }
944 :
945 110 : poBand = RawRasterBand::Create(
946 55 : poDS.get(), nBand, fpRaw, 0, nPixelOffset,
947 55 : nPixelOffset * poDS->GetRasterXSize(), eDataType, eByteOrder,
948 55 : RawRasterBand::OwnFP::YES);
949 : }
950 :
951 57 : poDS->SetBand(nBand, std::move(poBand));
952 57 : }
953 :
954 29 : CPLFree(pszTargetPath);
955 29 : CPLFree(pszTargetBase);
956 29 : CSLDestroy(papszDirFiles);
957 :
958 : /* -------------------------------------------------------------------- */
959 : /* Check if we have bands. */
960 : /* -------------------------------------------------------------------- */
961 29 : if (poDS->GetRasterCount() == 0)
962 : {
963 0 : if (nSkipped > 0 && poOpenInfo->eAccess)
964 : {
965 0 : CPLError(CE_Failure, CPLE_OpenFailed,
966 : "Failed to open %d files that were apparently bands. "
967 : "Perhaps this dataset is readonly?",
968 : nSkipped);
969 0 : return nullptr;
970 : }
971 : else
972 : {
973 0 : CPLError(CE_Failure, CPLE_OpenFailed,
974 : "MFF header file read successfully, but no bands "
975 : "were successfully found and opened.");
976 0 : return nullptr;
977 : }
978 : }
979 :
980 : /* -------------------------------------------------------------------- */
981 : /* Set all information from the .hdr that isn't well know to be */
982 : /* metadata. */
983 : /* -------------------------------------------------------------------- */
984 226 : for (int i = 0; papszHdrLines[i] != nullptr; i++)
985 : {
986 197 : char *pszName = nullptr;
987 :
988 197 : const char *pszValue = CPLParseNameValue(papszHdrLines[i], &pszName);
989 197 : if (pszName == nullptr || pszValue == nullptr)
990 16 : continue;
991 :
992 181 : if (!EQUAL(pszName, "END") && !EQUAL(pszName, "FILE_TYPE") &&
993 156 : !EQUAL(pszName, "BYTE_ORDER") && !EQUAL(pszName, "no_columns") &&
994 131 : !EQUAL(pszName, "no_rows") && !EQUAL(pszName, "type") &&
995 129 : !EQUAL(pszName, "tile_size_rows") &&
996 127 : !EQUAL(pszName, "tile_size_columns") &&
997 125 : !EQUAL(pszName, "IMAGE_FILE_FORMAT") &&
998 96 : !EQUAL(pszName, "IMAGE_LINES") && !EQUAL(pszName, "LINE_SAMPLES"))
999 : {
1000 38 : poDS->SetMetadataItem(pszName, pszValue);
1001 : }
1002 :
1003 181 : CPLFree(pszName);
1004 : }
1005 :
1006 : /* -------------------------------------------------------------------- */
1007 : /* Any GCPs in header file? */
1008 : /* -------------------------------------------------------------------- */
1009 29 : poDS->ScanForGCPs();
1010 29 : poDS->ScanForProjectionInfo();
1011 29 : if (poDS->nGCPCount == 0)
1012 23 : poDS->m_oGCPSRS.Clear();
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Initialize any PAM information. */
1016 : /* -------------------------------------------------------------------- */
1017 29 : poDS->SetDescription(poOpenInfo->pszFilename);
1018 29 : poDS->TryLoadXML();
1019 :
1020 : /* -------------------------------------------------------------------- */
1021 : /* Check for overviews. */
1022 : /* -------------------------------------------------------------------- */
1023 29 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1024 :
1025 29 : return poDS.release();
1026 : }
1027 :
1028 9 : int GetMFFProjectionType(const OGRSpatialReference *poSRS)
1029 : {
1030 9 : if (poSRS == nullptr)
1031 : {
1032 0 : return MFFPRJ_NONE;
1033 : }
1034 9 : if (poSRS->IsProjected() && poSRS->GetAttrValue("PROJECTION") &&
1035 0 : EQUAL(poSRS->GetAttrValue("PROJECTION"), SRS_PT_TRANSVERSE_MERCATOR))
1036 : {
1037 0 : return MFFPRJ_UTM;
1038 : }
1039 9 : else if (poSRS->IsGeographic())
1040 : {
1041 9 : return MFFPRJ_LL;
1042 : }
1043 : else
1044 : {
1045 0 : return MFFPRJ_UNRECOGNIZED;
1046 : }
1047 : }
1048 :
1049 : /************************************************************************/
1050 : /* Create() */
1051 : /************************************************************************/
1052 :
1053 50 : GDALDataset *MFFDataset::Create(const char *pszFilenameIn, int nXSize,
1054 : int nYSize, int nBandsIn, GDALDataType eType,
1055 : char **papszParamList)
1056 :
1057 : {
1058 : /* -------------------------------------------------------------------- */
1059 : /* Verify input options. */
1060 : /* -------------------------------------------------------------------- */
1061 50 : if (nBandsIn <= 0)
1062 : {
1063 1 : CPLError(CE_Failure, CPLE_NotSupported,
1064 : "MFF driver does not support %d bands.", nBandsIn);
1065 1 : return nullptr;
1066 : }
1067 :
1068 49 : if (eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 &&
1069 27 : eType != GDT_CInt16 && eType != GDT_CFloat32)
1070 : {
1071 24 : CPLError(CE_Failure, CPLE_AppDefined,
1072 : "Attempt to create MFF file with currently unsupported\n"
1073 : "data type (%s).\n",
1074 : GDALGetDataTypeName(eType));
1075 :
1076 24 : return nullptr;
1077 : }
1078 :
1079 : /* -------------------------------------------------------------------- */
1080 : /* Establish the base filename (path+filename, less extension). */
1081 : /* -------------------------------------------------------------------- */
1082 : char *pszBaseFilename =
1083 25 : static_cast<char *>(CPLMalloc(strlen(pszFilenameIn) + 5));
1084 25 : strcpy(pszBaseFilename, pszFilenameIn);
1085 :
1086 100 : for (int i = static_cast<int>(strlen(pszBaseFilename)) - 1; i > 0; i--)
1087 : {
1088 100 : if (pszBaseFilename[i] == '.')
1089 : {
1090 0 : pszBaseFilename[i] = '\0';
1091 0 : break;
1092 : }
1093 :
1094 100 : if (pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\')
1095 : break;
1096 : }
1097 :
1098 : /* -------------------------------------------------------------------- */
1099 : /* Create the header file. */
1100 : /* -------------------------------------------------------------------- */
1101 25 : const char *pszFilename = CPLFormFilename(nullptr, pszBaseFilename, "hdr");
1102 :
1103 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "wt");
1104 25 : if (fp == nullptr)
1105 : {
1106 3 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.\n",
1107 : pszFilename);
1108 3 : CPLFree(pszBaseFilename);
1109 3 : return nullptr;
1110 : }
1111 :
1112 22 : bool bOK = VSIFPrintfL(fp, "IMAGE_FILE_FORMAT = MFF\n") >= 0;
1113 22 : bOK &= VSIFPrintfL(fp, "FILE_TYPE = IMAGE\n") >= 0;
1114 22 : bOK &= VSIFPrintfL(fp, "IMAGE_LINES = %d\n", nYSize) >= 0;
1115 22 : bOK &= VSIFPrintfL(fp, "LINE_SAMPLES = %d\n", nXSize) >= 0;
1116 : #ifdef CPL_MSB
1117 : bOK &= VSIFPrintfL(fp, "BYTE_ORDER = MSB\n") >= 0;
1118 : #else
1119 22 : bOK &= VSIFPrintfL(fp, "BYTE_ORDER = LSB\n") >= 0;
1120 : #endif
1121 :
1122 22 : if (CSLFetchNameValue(papszParamList, "NO_END") == nullptr)
1123 13 : bOK &= VSIFPrintfL(fp, "END\n") >= 0;
1124 :
1125 22 : if (VSIFCloseL(fp) != 0)
1126 0 : bOK = false;
1127 :
1128 : /* -------------------------------------------------------------------- */
1129 : /* Create the data files, but don't bother writing any data to them.*/
1130 : /* -------------------------------------------------------------------- */
1131 72 : for (int iBand = 0; bOK && iBand < nBandsIn; iBand++)
1132 : {
1133 50 : char szExtension[4] = {'\0'};
1134 :
1135 50 : if (eType == GDT_Byte)
1136 30 : CPLsnprintf(szExtension, sizeof(szExtension), "b%02d", iBand);
1137 20 : else if (eType == GDT_UInt16)
1138 5 : CPLsnprintf(szExtension, sizeof(szExtension), "i%02d", iBand);
1139 15 : else if (eType == GDT_Float32)
1140 5 : CPLsnprintf(szExtension, sizeof(szExtension), "r%02d", iBand);
1141 10 : else if (eType == GDT_CInt16)
1142 5 : CPLsnprintf(szExtension, sizeof(szExtension), "j%02d", iBand);
1143 5 : else if (eType == GDT_CFloat32)
1144 5 : CPLsnprintf(szExtension, sizeof(szExtension), "x%02d", iBand);
1145 :
1146 50 : pszFilename = CPLFormFilename(nullptr, pszBaseFilename, szExtension);
1147 50 : fp = VSIFOpenL(pszFilename, "wb");
1148 50 : if (fp == nullptr)
1149 : {
1150 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Couldn't create %s.\n",
1151 : pszFilename);
1152 0 : CPLFree(pszBaseFilename);
1153 0 : return nullptr;
1154 : }
1155 :
1156 50 : bOK &= VSIFWriteL("", 1, 1, fp) == 1;
1157 50 : if (VSIFCloseL(fp) != 0)
1158 0 : bOK = false;
1159 : }
1160 :
1161 22 : if (!bOK)
1162 : {
1163 0 : CPLFree(pszBaseFilename);
1164 0 : return nullptr;
1165 : }
1166 :
1167 : /* -------------------------------------------------------------------- */
1168 : /* Open the dataset normally. */
1169 : /* -------------------------------------------------------------------- */
1170 22 : strcat(pszBaseFilename, ".hdr");
1171 : GDALDataset *poDS =
1172 22 : static_cast<GDALDataset *>(GDALOpen(pszBaseFilename, GA_Update));
1173 22 : CPLFree(pszBaseFilename);
1174 :
1175 22 : return poDS;
1176 : }
1177 :
1178 : /************************************************************************/
1179 : /* CreateCopy() */
1180 : /************************************************************************/
1181 :
1182 19 : GDALDataset *MFFDataset::CreateCopy(const char *pszFilename,
1183 : GDALDataset *poSrcDS, int /* bStrict */,
1184 : char **papszOptions,
1185 : GDALProgressFunc pfnProgress,
1186 : void *pProgressData)
1187 : {
1188 19 : const int nBands = poSrcDS->GetRasterCount();
1189 19 : if (nBands == 0)
1190 : {
1191 1 : CPLError(CE_Failure, CPLE_NotSupported,
1192 : "MFF driver does not support source dataset with zero band.");
1193 1 : return nullptr;
1194 : }
1195 :
1196 18 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1197 18 : if (!pfnProgress(0.0, nullptr, pProgressData))
1198 0 : return nullptr;
1199 :
1200 : // Check that other bands match type- sets type
1201 : // to unknown if they differ.
1202 28 : for (int iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++)
1203 : {
1204 10 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
1205 10 : eType = GDALDataTypeUnion(eType, poBand->GetRasterDataType());
1206 : }
1207 :
1208 18 : char **newpapszOptions = CSLDuplicate(papszOptions);
1209 18 : newpapszOptions = CSLSetNameValue(newpapszOptions, "NO_END", "TRUE");
1210 :
1211 18 : MFFDataset *poDS = reinterpret_cast<MFFDataset *>(Create(
1212 : pszFilename, poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
1213 : poSrcDS->GetRasterCount(), eType, newpapszOptions));
1214 :
1215 18 : CSLDestroy(newpapszOptions);
1216 :
1217 18 : if (poDS == nullptr)
1218 9 : return nullptr;
1219 :
1220 : /* -------------------------------------------------------------------- */
1221 : /* Copy the image data. */
1222 : /* -------------------------------------------------------------------- */
1223 9 : const int nXSize = poDS->GetRasterXSize();
1224 9 : const int nYSize = poDS->GetRasterYSize();
1225 :
1226 9 : int nBlockXSize = 0;
1227 9 : int nBlockYSize = 0;
1228 9 : poDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
1229 :
1230 9 : const int nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize) *
1231 9 : ((nYSize + nBlockYSize - 1) / nBlockYSize) *
1232 9 : poSrcDS->GetRasterCount();
1233 :
1234 9 : int nBlocksDone = 0;
1235 28 : for (int iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++)
1236 : {
1237 19 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand + 1);
1238 19 : GDALRasterBand *poDstBand = poDS->GetRasterBand(iBand + 1);
1239 :
1240 57 : void *pData = CPLMalloc(static_cast<size_t>(nBlockXSize) * nBlockYSize *
1241 19 : GDALGetDataTypeSizeBytes(eType));
1242 :
1243 209 : for (int iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize)
1244 : {
1245 380 : for (int iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize)
1246 : {
1247 190 : if (!pfnProgress((nBlocksDone++) /
1248 190 : static_cast<float>(nBlockTotal),
1249 : nullptr, pProgressData))
1250 : {
1251 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1252 0 : delete poDS;
1253 0 : CPLFree(pData);
1254 :
1255 : GDALDriver *poMFFDriver =
1256 0 : static_cast<GDALDriver *>(GDALGetDriverByName("MFF"));
1257 0 : poMFFDriver->Delete(pszFilename);
1258 0 : return nullptr;
1259 : }
1260 :
1261 190 : const int nTBXSize = std::min(nBlockXSize, nXSize - iXOffset);
1262 190 : const int nTBYSize = std::min(nBlockYSize, nYSize - iYOffset);
1263 :
1264 190 : CPLErr eErr = poSrcBand->RasterIO(
1265 : GF_Read, iXOffset, iYOffset, nTBXSize, nTBYSize, pData,
1266 : nTBXSize, nTBYSize, eType, 0, 0, nullptr);
1267 :
1268 190 : if (eErr != CE_None)
1269 : {
1270 0 : delete poDS;
1271 0 : CPLFree(pData);
1272 0 : return nullptr;
1273 : }
1274 :
1275 190 : eErr = poDstBand->RasterIO(GF_Write, iXOffset, iYOffset,
1276 : nTBXSize, nTBYSize, pData, nTBXSize,
1277 : nTBYSize, eType, 0, 0, nullptr);
1278 :
1279 190 : if (eErr != CE_None)
1280 : {
1281 0 : delete poDS;
1282 0 : CPLFree(pData);
1283 0 : return nullptr;
1284 : }
1285 : }
1286 : }
1287 :
1288 19 : CPLFree(pData);
1289 : }
1290 :
1291 : /* -------------------------------------------------------------------- */
1292 : /* Copy georeferencing information, if enough is available. */
1293 : /* -------------------------------------------------------------------- */
1294 :
1295 : /* -------------------------------------------------------------------- */
1296 : /* Establish the base filename (path+filename, less extension). */
1297 : /* -------------------------------------------------------------------- */
1298 : char *pszBaseFilename =
1299 9 : static_cast<char *>(CPLMalloc(strlen(pszFilename) + 5));
1300 9 : strcpy(pszBaseFilename, pszFilename);
1301 :
1302 36 : for (int i = static_cast<int>(strlen(pszBaseFilename)) - 1; i > 0; i--)
1303 : {
1304 36 : if (pszBaseFilename[i] == '.')
1305 : {
1306 0 : pszBaseFilename[i] = '\0';
1307 0 : break;
1308 : }
1309 :
1310 36 : if (pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\')
1311 : break;
1312 : }
1313 :
1314 : const char *pszFilenameGEO =
1315 9 : CPLFormFilename(nullptr, pszBaseFilename, "hdr");
1316 :
1317 9 : VSILFILE *fp = VSIFOpenL(pszFilenameGEO, "at");
1318 9 : if (fp == nullptr)
1319 : {
1320 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1321 : "Couldn't open %s for appending.\n", pszFilenameGEO);
1322 0 : CPLFree(pszBaseFilename);
1323 0 : return nullptr;
1324 : }
1325 :
1326 : /* MFF requires corner and center gcps */
1327 9 : bool georef_created = false;
1328 :
1329 : double *padfTiepoints =
1330 9 : static_cast<double *>(CPLMalloc(2 * sizeof(double) * 5));
1331 :
1332 9 : const int src_prj = GetMFFProjectionType(poSrcDS->GetSpatialRef());
1333 :
1334 9 : if ((src_prj != MFFPRJ_NONE) && (src_prj != MFFPRJ_UNRECOGNIZED))
1335 : {
1336 : double *tempGeoTransform =
1337 9 : static_cast<double *>(CPLMalloc(6 * sizeof(double)));
1338 :
1339 18 : if ((poSrcDS->GetGeoTransform(tempGeoTransform) == CE_None) &&
1340 9 : (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0 ||
1341 0 : tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0 ||
1342 0 : tempGeoTransform[4] != 0.0 ||
1343 0 : std::abs(tempGeoTransform[5]) != 1.0))
1344 : {
1345 9 : padfTiepoints[0] = tempGeoTransform[0] + tempGeoTransform[1] * 0.5 +
1346 9 : tempGeoTransform[2] * 0.5;
1347 :
1348 9 : padfTiepoints[1] = tempGeoTransform[3] + tempGeoTransform[4] * 0.5 +
1349 9 : tempGeoTransform[5] * 0.5;
1350 :
1351 9 : padfTiepoints[2] =
1352 18 : tempGeoTransform[0] + tempGeoTransform[2] * 0.5 +
1353 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize() - 0.5);
1354 :
1355 9 : padfTiepoints[3] =
1356 18 : tempGeoTransform[3] + tempGeoTransform[5] * 0.5 +
1357 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize() - 0.5);
1358 :
1359 9 : padfTiepoints[4] =
1360 18 : tempGeoTransform[0] + tempGeoTransform[1] * 0.5 +
1361 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize() - 0.5);
1362 :
1363 9 : padfTiepoints[5] =
1364 18 : tempGeoTransform[3] + tempGeoTransform[4] * 0.5 +
1365 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize() - 0.5);
1366 :
1367 9 : padfTiepoints[6] =
1368 18 : tempGeoTransform[0] +
1369 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize() - 0.5) +
1370 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize() - 0.5);
1371 :
1372 9 : padfTiepoints[7] =
1373 18 : tempGeoTransform[3] +
1374 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize() - 0.5) +
1375 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize() - 0.5);
1376 :
1377 9 : padfTiepoints[8] =
1378 18 : tempGeoTransform[0] +
1379 9 : tempGeoTransform[1] * (poSrcDS->GetRasterXSize()) / 2.0 +
1380 9 : tempGeoTransform[2] * (poSrcDS->GetRasterYSize()) / 2.0;
1381 :
1382 9 : padfTiepoints[9] =
1383 18 : tempGeoTransform[3] +
1384 9 : tempGeoTransform[4] * (poSrcDS->GetRasterXSize()) / 2.0 +
1385 9 : tempGeoTransform[5] * (poSrcDS->GetRasterYSize()) / 2.0;
1386 :
1387 18 : OGRSpatialReference oUTMorLL;
1388 9 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
1389 9 : if (poSrcSRS)
1390 9 : oUTMorLL = *poSrcSRS;
1391 9 : auto poLLSRS = oUTMorLL.CloneGeogCS();
1392 9 : if (poLLSRS && oUTMorLL.IsProjected())
1393 : {
1394 0 : poLLSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1395 : OGRCoordinateTransformation *poTransform =
1396 0 : OGRCreateCoordinateTransformation(&oUTMorLL, poLLSRS);
1397 :
1398 : // projected coordinate system- need to translate gcps */
1399 0 : bool bSuccess = poTransform != nullptr;
1400 :
1401 0 : for (int index = 0; index < 5; index++)
1402 : {
1403 : // TODO: If bSuccess is false, set it to false?
1404 0 : if (!bSuccess || !poTransform->Transform(
1405 : 1, &(padfTiepoints[index * 2]),
1406 : &(padfTiepoints[index * 2 + 1])))
1407 0 : bSuccess = false;
1408 : }
1409 0 : if (bSuccess)
1410 0 : georef_created = true;
1411 : }
1412 : else
1413 : {
1414 9 : georef_created = true;
1415 : }
1416 9 : delete poLLSRS;
1417 : }
1418 9 : CPLFree(tempGeoTransform);
1419 : }
1420 :
1421 9 : bool bOK = true;
1422 9 : if (georef_created)
1423 : {
1424 : /* --------------------------------------------------------------------
1425 : */
1426 : /* top left */
1427 : /* --------------------------------------------------------------------
1428 : */
1429 18 : bOK &= VSIFPrintfL(fp, "TOP_LEFT_CORNER_LATITUDE = %.10f\n",
1430 9 : padfTiepoints[1]) >= 0;
1431 9 : bOK &= VSIFPrintfL(fp, "TOP_LEFT_CORNER_LONGITUDE = %.10f\n",
1432 9 : padfTiepoints[0]) >= 0;
1433 : /* --------------------------------------------------------------------
1434 : */
1435 : /* top_right */
1436 : /* --------------------------------------------------------------------
1437 : */
1438 18 : bOK &= VSIFPrintfL(fp, "TOP_RIGHT_CORNER_LATITUDE = %.10f\n",
1439 9 : padfTiepoints[3]) >= 0;
1440 18 : bOK &= VSIFPrintfL(fp, "TOP_RIGHT_CORNER_LONGITUDE = %.10f\n",
1441 9 : padfTiepoints[2]) >= 0;
1442 : /* --------------------------------------------------------------------
1443 : */
1444 : /* bottom_left */
1445 : /* --------------------------------------------------------------------
1446 : */
1447 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_LEFT_CORNER_LATITUDE = %.10f\n",
1448 9 : padfTiepoints[5]) >= 0;
1449 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_LEFT_CORNER_LONGITUDE = %.10f\n",
1450 9 : padfTiepoints[4]) >= 0;
1451 : /* --------------------------------------------------------------------
1452 : */
1453 : /* bottom_right */
1454 : /* --------------------------------------------------------------------
1455 : */
1456 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_RIGHT_CORNER_LATITUDE = %.10f\n",
1457 9 : padfTiepoints[7]) >= 0;
1458 18 : bOK &= VSIFPrintfL(fp, "BOTTOM_RIGHT_CORNER_LONGITUDE = %.10f\n",
1459 9 : padfTiepoints[6]) >= 0;
1460 : /* --------------------------------------------------------------------
1461 : */
1462 : /* Center */
1463 : /* --------------------------------------------------------------------
1464 : */
1465 9 : bOK &=
1466 9 : VSIFPrintfL(fp, "CENTRE_LATITUDE = %.10f\n", padfTiepoints[9]) >= 0;
1467 18 : bOK &= VSIFPrintfL(fp, "CENTRE_LONGITUDE = %.10f\n",
1468 9 : padfTiepoints[8]) >= 0;
1469 : /* -------------------------------------------------------------------
1470 : */
1471 : /* Ellipsoid/projection */
1472 : /* --------------------------------------------------------------------*/
1473 :
1474 9 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
1475 9 : char *spheroid_name = nullptr;
1476 :
1477 9 : if (poSrcSRS != nullptr)
1478 : {
1479 9 : if (poSrcSRS->IsProjected() &&
1480 9 : poSrcSRS->GetAttrValue("PROJECTION") != nullptr &&
1481 0 : EQUAL(poSrcSRS->GetAttrValue("PROJECTION"),
1482 : SRS_PT_TRANSVERSE_MERCATOR))
1483 : {
1484 0 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = UTM\n") >= 0;
1485 0 : OGRErr ogrerrorOl = OGRERR_NONE;
1486 0 : bOK &=
1487 0 : VSIFPrintfL(fp, "PROJECTION_ORIGIN_LONGITUDE = %f\n",
1488 : poSrcSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,
1489 0 : 0.0, &ogrerrorOl)) >= 0;
1490 : }
1491 9 : else if (poSrcSRS->IsGeographic())
1492 : {
1493 9 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = LL\n") >= 0;
1494 : }
1495 : else
1496 : {
1497 0 : CPLError(CE_Warning, CPLE_AppDefined,
1498 : "Unrecognized projection- no georeferencing "
1499 : "information transferred.");
1500 0 : bOK &= VSIFPrintfL(fp, "PROJECTION_NAME = LL\n") >= 0;
1501 : }
1502 9 : OGRErr ogrerrorEq = OGRERR_NONE;
1503 9 : const double eq_radius = poSrcSRS->GetSemiMajor(&ogrerrorEq);
1504 9 : OGRErr ogrerrorInvf = OGRERR_NONE;
1505 : const double inv_flattening =
1506 9 : poSrcSRS->GetInvFlattening(&ogrerrorInvf);
1507 9 : if (ogrerrorEq == OGRERR_NONE && ogrerrorInvf == OGRERR_NONE)
1508 : {
1509 9 : MFFSpheroidList *mffEllipsoids = new MFFSpheroidList;
1510 : spheroid_name =
1511 9 : mffEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(
1512 : eq_radius, inv_flattening);
1513 9 : if (spheroid_name != nullptr)
1514 : {
1515 0 : bOK &= VSIFPrintfL(fp, "SPHEROID_NAME = %s\n",
1516 0 : spheroid_name) >= 0;
1517 : }
1518 : else
1519 : {
1520 18 : bOK &= VSIFPrintfL(fp,
1521 : "SPHEROID_NAME = USER_DEFINED\n"
1522 : "SPHEROID_EQUATORIAL_RADIUS = %.10f\n"
1523 : "SPHEROID_POLAR_RADIUS = %.10f\n",
1524 : eq_radius,
1525 : eq_radius *
1526 9 : (1 - 1.0 / inv_flattening)) >= 0;
1527 : }
1528 9 : delete mffEllipsoids;
1529 9 : CPLFree(spheroid_name);
1530 : }
1531 : }
1532 : }
1533 :
1534 9 : CPLFree(padfTiepoints);
1535 9 : bOK &= VSIFPrintfL(fp, "END\n") >= 0;
1536 9 : if (VSIFCloseL(fp) != 0)
1537 0 : bOK = false;
1538 :
1539 9 : if (!bOK)
1540 : {
1541 0 : delete poDS;
1542 0 : CPLFree(pszBaseFilename);
1543 0 : return nullptr;
1544 : }
1545 :
1546 : /* End of georeferencing stuff */
1547 :
1548 : /* Make sure image data gets flushed */
1549 28 : for (int iBand = 0; iBand < poDS->GetRasterCount(); iBand++)
1550 : {
1551 : RawRasterBand *poDstBand =
1552 19 : reinterpret_cast<RawRasterBand *>(poDS->GetRasterBand(iBand + 1));
1553 19 : poDstBand->FlushCache(false);
1554 : }
1555 :
1556 9 : if (!pfnProgress(1.0, nullptr, pProgressData))
1557 : {
1558 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1559 0 : delete poDS;
1560 :
1561 : GDALDriver *poMFFDriver =
1562 0 : static_cast<GDALDriver *>(GDALGetDriverByName("MFF"));
1563 0 : poMFFDriver->Delete(pszFilename);
1564 0 : CPLFree(pszBaseFilename);
1565 0 : return nullptr;
1566 : }
1567 :
1568 9 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1569 9 : CPLFree(pszBaseFilename);
1570 :
1571 9 : return poDS;
1572 : }
1573 :
1574 : /************************************************************************/
1575 : /* GDALRegister_MFF() */
1576 : /************************************************************************/
1577 :
1578 1595 : void GDALRegister_MFF()
1579 :
1580 : {
1581 1595 : if (GDALGetDriverByName("MFF") != nullptr)
1582 302 : return;
1583 :
1584 1293 : GDALDriver *poDriver = new GDALDriver();
1585 :
1586 1293 : poDriver->SetDescription("MFF");
1587 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1588 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF Raster");
1589 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff.html");
1590 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
1591 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1592 1293 : "Byte UInt16 Float32 CInt16 CFloat32");
1593 :
1594 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1595 :
1596 1293 : poDriver->pfnOpen = MFFDataset::Open;
1597 1293 : poDriver->pfnCreate = MFFDataset::Create;
1598 1293 : poDriver->pfnCreateCopy = MFFDataset::CreateCopy;
1599 :
1600 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1601 : }
|
