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