Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: FITS Driver
4 : * Purpose: Implement FITS raster read/write support
5 : * Author: Simon Perkins, s.perkins@lanl.gov
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2001, Simon Perkins
9 : * Copyright (c) 2008-2020, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2018, Chiara Marmo <chiara dot marmo at u-psud dot fr>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : // So that OFF_T is 64 bits
16 : #define _FILE_OFFSET_BITS 64
17 :
18 : #include "cpl_string.h"
19 : #include "gdal_frmts.h"
20 : #include "gdal_pam.h"
21 : #include "ogr_spatialref.h"
22 : #include "ogrsf_frmts.h"
23 : #include "fitsdrivercore.h"
24 :
25 : #include <string.h>
26 : #include <fitsio.h>
27 :
28 : #include <algorithm>
29 : #include <string>
30 : #include <cstring>
31 : #include <set>
32 : #include <vector>
33 :
34 : /************************************************************************/
35 : /* ==================================================================== */
36 : /* FITSDataset */
37 : /* ==================================================================== */
38 : /************************************************************************/
39 :
40 : class FITSRasterBand;
41 : class FITSLayer;
42 :
43 : class FITSDataset final : public GDALPamDataset
44 : {
45 :
46 : friend class FITSRasterBand;
47 : friend class FITSLayer;
48 :
49 : fitsfile *m_hFITS = nullptr;
50 :
51 : int m_hduNum = 0;
52 : GDALDataType m_gdalDataType = GDT_Unknown; // GDAL code for the image type
53 : int m_fitsDataType = 0; // FITS code for the image type
54 :
55 : bool m_isExistingFile = false;
56 : LONGLONG m_highestOffsetWritten = 0; // How much of image has been written
57 :
58 : bool m_bNoDataChanged = false;
59 : bool m_bNoDataSet = false;
60 : double m_dfNoDataValue = -9999.0;
61 :
62 : bool m_bMetadataChanged = false;
63 :
64 : CPLStringList m_aosSubdatasets{};
65 :
66 : OGRSpatialReference m_oSRS{};
67 :
68 : double m_adfGeoTransform[6];
69 : bool m_bGeoTransformValid = false;
70 :
71 : bool m_bFITSInfoChanged = false;
72 :
73 : std::vector<std::unique_ptr<FITSLayer>> m_apoLayers{};
74 :
75 : CPLErr Init(fitsfile *hFITS, bool isExistingFile, int hduNum);
76 :
77 : void LoadGeoreferencing();
78 : void LoadFITSInfo();
79 : void WriteFITSInfo();
80 : void LoadMetadata(GDALMajorObject *poTarget);
81 :
82 : public:
83 : FITSDataset(); // Others should not call this constructor explicitly
84 : ~FITSDataset();
85 :
86 : static GDALDataset *Open(GDALOpenInfo *);
87 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
88 : int nBandsIn, GDALDataType eType,
89 : char **papszParamList);
90 : static CPLErr Delete(const char *pszFilename);
91 :
92 : const OGRSpatialReference *GetSpatialRef() const override;
93 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
94 : virtual CPLErr GetGeoTransform(double *) override;
95 : virtual CPLErr SetGeoTransform(double *) override;
96 : char **GetMetadata(const char *papszDomain = nullptr) override;
97 :
98 75 : int GetLayerCount() override
99 : {
100 75 : return static_cast<int>(m_apoLayers.size());
101 : }
102 :
103 : OGRLayer *GetLayer(int) override;
104 :
105 : OGRLayer *ICreateLayer(const char *pszName,
106 : const OGRGeomFieldDefn *poGeomFieldDefn,
107 : CSLConstList papszOptions) override;
108 :
109 : int TestCapability(const char *pszCap) override;
110 :
111 : bool GetRawBinaryLayout(GDALDataset::RawBinaryLayout &) override;
112 : };
113 :
114 : /************************************************************************/
115 : /* ==================================================================== */
116 : /* FITSRasterBand */
117 : /* ==================================================================== */
118 : /************************************************************************/
119 :
120 : class FITSRasterBand final : public GDALPamRasterBand
121 : {
122 :
123 : friend class FITSDataset;
124 :
125 : bool m_bHaveOffsetScale = false;
126 : double m_dfOffset = 0.0;
127 : double m_dfScale = 1.0;
128 :
129 : protected:
130 : FITSDataset *m_poFDS = nullptr;
131 :
132 : bool m_bNoDataSet = false;
133 : double m_dfNoDataValue = -9999.0;
134 :
135 : public:
136 : FITSRasterBand(FITSDataset *, int);
137 : virtual ~FITSRasterBand();
138 :
139 : virtual CPLErr IReadBlock(int, int, void *) override;
140 : virtual CPLErr IWriteBlock(int, int, void *) override;
141 :
142 : virtual double GetNoDataValue(int *) override final;
143 : virtual CPLErr SetNoDataValue(double) override final;
144 : virtual CPLErr DeleteNoDataValue() override final;
145 :
146 : virtual double GetOffset(int *pbSuccess = nullptr) override final;
147 : virtual CPLErr SetOffset(double dfNewValue) override final;
148 : virtual double GetScale(int *pbSuccess = nullptr) override final;
149 : virtual CPLErr SetScale(double dfNewValue) override final;
150 : };
151 :
152 : /************************************************************************/
153 : /* ==================================================================== */
154 : /* FITSLayer */
155 : /* ==================================================================== */
156 : /************************************************************************/
157 : namespace gdal::FITS
158 : {
159 : struct ColDesc
160 : {
161 : std::string typechar{};
162 : int iCol = 0; // numbering starting at 1
163 : int iBit = 0; // numbering starting at 1
164 : int nRepeat = 0;
165 : int nItems = 1;
166 : double dfOffset = 0;
167 : double dfScale = 1;
168 : bool bHasNull = false;
169 : LONGLONG nNullValue = 0;
170 : int nTypeCode = 0; // unset
171 : };
172 : } // namespace gdal::FITS
173 :
174 : using namespace gdal::FITS;
175 :
176 : class FITSLayer final : public OGRLayer,
177 : public OGRGetNextFeatureThroughRaw<FITSLayer>
178 : {
179 : friend class FITSDataset;
180 : FITSDataset *m_poDS = nullptr;
181 : int m_hduNum = 0;
182 : OGRFeatureDefn *m_poFeatureDefn = nullptr;
183 : LONGLONG m_nCurRow = 1;
184 : LONGLONG m_nRows = 0;
185 :
186 : std::vector<ColDesc> m_aoColDescs;
187 :
188 : CPLStringList m_aosCreationOptions{};
189 :
190 : std::vector<int> m_anDeferredFieldsIndices{};
191 :
192 : OGRFeature *GetNextRawFeature();
193 : void SetActiveHDU();
194 : void RunDeferredFieldCreation(const OGRFeature *poFeature = nullptr);
195 : bool SetOrCreateFeature(const OGRFeature *poFeature, LONGLONG nRow);
196 :
197 : public:
198 : FITSLayer(FITSDataset *poDS, int hduNum, const char *pszExtName);
199 : ~FITSLayer();
200 :
201 63 : OGRFeatureDefn *GetLayerDefn() override
202 : {
203 63 : return m_poFeatureDefn;
204 : }
205 :
206 : void ResetReading() override;
207 : int TestCapability(const char *) override;
208 : OGRFeature *GetFeature(GIntBig) override;
209 : GIntBig GetFeatureCount(int bForce) override;
210 : OGRErr CreateField(const OGRFieldDefn *poField, int bApproxOK) override;
211 : OGRErr ICreateFeature(OGRFeature *poFeature) override;
212 : OGRErr ISetFeature(OGRFeature *poFeature) override;
213 : OGRErr DeleteFeature(GIntBig nFID) override;
214 :
215 23 : DEFINE_GET_NEXT_FEATURE_THROUGH_RAW(FITSLayer)
216 :
217 4 : void SetCreationOptions(CSLConstList papszOptions)
218 : {
219 4 : m_aosCreationOptions = papszOptions;
220 4 : }
221 : };
222 :
223 : /************************************************************************/
224 : /* FITSLayer() */
225 : /************************************************************************/
226 :
227 16 : FITSLayer::FITSLayer(FITSDataset *poDS, int hduNum, const char *pszExtName)
228 16 : : m_poDS(poDS), m_hduNum(hduNum)
229 : {
230 16 : if (pszExtName[0] != 0)
231 14 : m_poFeatureDefn = new OGRFeatureDefn(pszExtName);
232 : else
233 2 : m_poFeatureDefn =
234 2 : new OGRFeatureDefn(CPLSPrintf("Table HDU %d", hduNum));
235 16 : m_poFeatureDefn->Reference();
236 16 : m_poFeatureDefn->SetGeomType(wkbNone);
237 16 : SetDescription(m_poFeatureDefn->GetName());
238 :
239 16 : SetActiveHDU();
240 :
241 16 : m_poDS->LoadMetadata(this);
242 :
243 16 : int status = 0;
244 16 : fits_get_num_rowsll(m_poDS->m_hFITS, &m_nRows, &status);
245 16 : if (status)
246 : {
247 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_get_num_rowsll() failed");
248 : }
249 16 : int nCols = 0;
250 16 : status = 0;
251 16 : fits_get_num_cols(m_poDS->m_hFITS, &nCols, &status);
252 16 : if (status)
253 : {
254 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_get_num_cols() failed");
255 : }
256 :
257 32 : std::vector<std::string> aosNames(nCols);
258 32 : std::vector<char *> apszNames(nCols);
259 460 : for (int i = 0; i < nCols; i++)
260 : {
261 444 : aosNames[i].resize(80);
262 444 : apszNames[i] = &aosNames[i][0];
263 : }
264 :
265 16 : status = 0;
266 16 : fits_read_btblhdrll(m_poDS->m_hFITS, nCols, nullptr, nullptr,
267 : apszNames.data(), nullptr, nullptr, nullptr, nullptr,
268 : &status);
269 16 : if (status)
270 : {
271 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_read_btblhdrll() failed");
272 : }
273 :
274 460 : for (int i = 0; i < nCols; i++)
275 : {
276 444 : aosNames[i].resize(strlen(aosNames[i].c_str()));
277 :
278 : char typechar[80];
279 444 : LONGLONG nRepeat = 0;
280 444 : double dfScale = 0;
281 444 : double dfOffset = 0;
282 444 : status = 0;
283 444 : fits_get_bcolparmsll(m_poDS->m_hFITS, i + 1,
284 : nullptr, // column name
285 : nullptr, // unit
286 : typechar, &nRepeat, &dfScale, &dfOffset,
287 : nullptr, // nulval
288 : nullptr, // tdisp
289 : &status);
290 444 : if (status)
291 : {
292 0 : CPLError(CE_Failure, CPLE_AppDefined,
293 : "fits_get_bcolparmsll() failed");
294 : }
295 :
296 444 : ColDesc oCol;
297 :
298 444 : status = 0;
299 444 : fits_read_key(m_poDS->m_hFITS, TLONGLONG, CPLSPrintf("TNULL%d", i + 1),
300 : &oCol.nNullValue, nullptr, &status);
301 444 : oCol.bHasNull = status == 0;
302 :
303 444 : OGRFieldType eType = OFTString;
304 444 : OGRFieldSubType eSubType = OFSTNone;
305 444 : if (typechar[0] == 'L') // Logical
306 : {
307 12 : eType = OFTInteger;
308 12 : eSubType = OFSTBoolean;
309 : }
310 432 : else if (typechar[0] == 'X') // Bit array
311 : {
312 20 : if (nRepeat > 128)
313 : {
314 0 : CPLDebug("FITS", "Too large repetition count for column %s",
315 0 : aosNames[i].c_str());
316 0 : continue;
317 : }
318 360 : for (int j = 1; j <= nRepeat; j++)
319 : {
320 : OGRFieldDefn oFieldDefn(
321 340 : (aosNames[i] + CPLSPrintf("_bit%d", j)).c_str(),
322 680 : OFTInteger);
323 : // cppcheck-suppress danglingTemporaryLifetime
324 340 : m_poFeatureDefn->AddFieldDefn(&oFieldDefn);
325 :
326 680 : ColDesc oColBit;
327 340 : oColBit.typechar = typechar;
328 340 : oColBit.iCol = i + 1;
329 340 : oColBit.iBit = j;
330 340 : m_aoColDescs.emplace_back(oColBit);
331 : }
332 20 : continue;
333 : }
334 412 : else if (typechar[0] == 'B') // Unsigned byte
335 : {
336 30 : if (dfOffset == -128 && dfScale == 1)
337 : {
338 6 : eType = OFTInteger; // signed byte
339 6 : oCol.nTypeCode = TSBYTE;
340 : // fits_read_col() automatically offsets numeric values
341 6 : dfOffset = 0;
342 : }
343 24 : else if (dfOffset != 0 || dfScale != 1)
344 6 : eType = OFTReal;
345 : else
346 18 : eType = OFTInteger;
347 : }
348 382 : else if (typechar[0] == 'I') // 16-bit signed integer
349 : {
350 31 : if (dfOffset == 32768.0 && dfScale == 1)
351 : {
352 6 : eType = OFTInteger; // unsigned 16-bit integer
353 6 : oCol.nTypeCode = TUSHORT;
354 : // fits_read_col() automatically offsets numeric values
355 6 : dfOffset = 0;
356 : }
357 25 : else if (dfOffset != 0 || dfScale != 1)
358 6 : eType = OFTReal;
359 : else
360 : {
361 19 : eType = OFTInteger;
362 19 : eSubType = OFSTInt16;
363 : }
364 : }
365 351 : else if (typechar[0] == 'J') // 32-bit signed integer
366 : {
367 58 : if (dfOffset == 2147483648.0 && dfScale == 1)
368 : {
369 6 : eType = OFTInteger64; // unsigned 32-bit integer --> needs to
370 : // promote to 64 bits
371 6 : oCol.nTypeCode = TUINT;
372 : // fits_read_col() automatically offsets numeric values
373 6 : dfOffset = 0;
374 : }
375 52 : else if (dfOffset != 0 || dfScale != 1)
376 6 : eType = OFTReal;
377 : else
378 46 : eType = OFTInteger;
379 : }
380 293 : else if (typechar[0] == 'K') // 64-bit signed integer
381 : {
382 29 : if (dfOffset != 0 || dfScale != 1)
383 6 : eType = OFTReal;
384 : else
385 23 : eType = OFTInteger64;
386 : }
387 264 : else if (typechar[0] == 'A') // Character
388 : {
389 :
390 26 : status = 0;
391 26 : LONGLONG nWidth = 0;
392 26 : fits_get_coltypell(m_poDS->m_hFITS, i + 1,
393 : nullptr, // typecode
394 : nullptr, // repeat
395 : &nWidth, &status);
396 26 : if (status)
397 : {
398 0 : CPLError(CE_Failure, CPLE_AppDefined,
399 : "fits_get_coltypell() failed");
400 : }
401 26 : if (nRepeat >= 2 * nWidth && nWidth != 0)
402 : {
403 6 : oCol.nItems = static_cast<int>(nRepeat / nWidth);
404 6 : eType = OFTStringList;
405 6 : nRepeat = nWidth;
406 : }
407 : else
408 : {
409 20 : eType = OFTString;
410 : }
411 : }
412 238 : else if (typechar[0] == 'E') // IEEE754 32bit
413 : {
414 25 : eType = OFTReal;
415 25 : if (dfOffset == 0 && dfScale == 1)
416 19 : eSubType = OFSTFloat32;
417 : // fits_read_col() automatically scales numeric values
418 25 : dfOffset = 0;
419 25 : dfScale = 1;
420 : }
421 213 : else if (typechar[0] == 'D') // IEEE754 64bit
422 : {
423 51 : eType = OFTReal;
424 : // fits_read_col() automatically scales numeric values
425 51 : dfOffset = 0;
426 51 : dfScale = 1;
427 : }
428 162 : else if (typechar[0] == 'C') // IEEE754 32bit complex
429 : {
430 20 : eType = OFTString;
431 : // fits_read_col() automatically scales numeric values
432 20 : dfOffset = 0;
433 20 : dfScale = 1;
434 : }
435 142 : else if (typechar[0] == 'M') // IEEE754 64bit complex
436 : {
437 20 : eType = OFTString;
438 : // fits_read_col() automatically scales numeric values
439 20 : dfOffset = 0;
440 20 : dfScale = 1;
441 : }
442 122 : else if (typechar[0] == 'P' || typechar[0] == 'Q') // Array
443 : {
444 122 : if (typechar[1] == 'L')
445 : {
446 12 : nRepeat = 0;
447 12 : eType = OFTIntegerList;
448 12 : eSubType = OFSTBoolean;
449 : }
450 110 : else if (typechar[1] == 'B')
451 : {
452 6 : nRepeat = 0;
453 6 : eType = OFTIntegerList;
454 : }
455 104 : else if (typechar[1] == 'I')
456 : {
457 11 : nRepeat = 0;
458 11 : eType = OFTIntegerList;
459 11 : eSubType = OFSTInt16;
460 : }
461 93 : else if (typechar[1] == 'J')
462 : {
463 26 : nRepeat = 0;
464 26 : eType = OFTIntegerList;
465 : }
466 67 : else if (typechar[1] == 'K')
467 : {
468 11 : nRepeat = 0;
469 11 : eType = OFTInteger64List;
470 : }
471 56 : else if (typechar[1] == 'A')
472 : {
473 22 : eType = OFTString;
474 : }
475 34 : else if (typechar[1] == 'E')
476 : {
477 11 : nRepeat = 0;
478 11 : eType = OFTRealList;
479 11 : if (dfOffset == 0 && dfScale == 1)
480 11 : eSubType = OFSTFloat32;
481 : // fits_read_col() automatically scales numeric values
482 11 : dfOffset = 0;
483 11 : dfScale = 1;
484 : }
485 23 : else if (typechar[1] == 'D')
486 : {
487 11 : nRepeat = 0;
488 11 : eType = OFTRealList;
489 : // fits_read_col() automatically scales numeric values
490 11 : dfOffset = 0;
491 11 : dfScale = 1;
492 : }
493 12 : else if (typechar[1] == 'C')
494 : {
495 6 : nRepeat = 0;
496 6 : eType = OFTStringList;
497 : // fits_read_col() automatically scales numeric values
498 6 : dfOffset = 0;
499 6 : dfScale = 1;
500 : }
501 6 : else if (typechar[1] == 'M')
502 : {
503 6 : nRepeat = 0;
504 6 : eType = OFTStringList;
505 : // fits_read_col() automatically scales numeric values
506 6 : dfOffset = 0;
507 6 : dfScale = 1;
508 : }
509 : else
510 : {
511 0 : CPLDebug("FITS", "Unhandled type %s", typechar);
512 0 : continue;
513 : }
514 : }
515 : else
516 : {
517 0 : CPLDebug("FITS", "Unhandled type %s", typechar);
518 0 : continue;
519 : }
520 :
521 424 : if (nRepeat > 1 && typechar[0] != 'A')
522 : {
523 84 : if (eType == OFTInteger)
524 45 : eType = OFTIntegerList;
525 39 : else if (eType == OFTInteger64)
526 9 : eType = OFTInteger64List;
527 30 : else if (eType == OFTReal)
528 18 : eType = OFTRealList;
529 12 : else if (eType == OFTString)
530 12 : eType = OFTStringList;
531 : }
532 :
533 848 : OGRFieldDefn oFieldDefn(aosNames[i].c_str(), eType);
534 424 : oFieldDefn.SetSubType(eSubType);
535 424 : if (typechar[0] == 'A')
536 26 : oFieldDefn.SetWidth(static_cast<int>(nRepeat));
537 424 : m_poFeatureDefn->AddFieldDefn(&oFieldDefn);
538 :
539 424 : oCol.typechar = typechar;
540 424 : oCol.iCol = i + 1;
541 424 : oCol.nRepeat = static_cast<int>(nRepeat);
542 424 : oCol.dfOffset = dfOffset;
543 424 : oCol.dfScale = dfScale;
544 424 : m_aoColDescs.emplace_back(oCol);
545 : }
546 16 : }
547 :
548 : /************************************************************************/
549 : /* ~FITSLayer() */
550 : /************************************************************************/
551 :
552 32 : FITSLayer::~FITSLayer()
553 : {
554 16 : RunDeferredFieldCreation();
555 :
556 18 : for (int i = 0; i < m_aosCreationOptions.size(); i++)
557 : {
558 2 : if (STARTS_WITH_CI(m_aosCreationOptions[i], "REPEAT_"))
559 : {
560 1 : char *pszKey = nullptr;
561 1 : CPL_IGNORE_RET_VAL(
562 1 : CPLParseNameValue(m_aosCreationOptions[i], &pszKey));
563 1 : if (pszKey &&
564 1 : m_poFeatureDefn->GetFieldIndex(pszKey + strlen("REPEAT_")) < 0)
565 : {
566 0 : CPLError(CE_Warning, CPLE_AppDefined,
567 : "Creation option %s ignored as field does not exist",
568 : m_aosCreationOptions[i]);
569 : }
570 1 : CPLFree(pszKey);
571 : }
572 : }
573 :
574 16 : m_poFeatureDefn->Release();
575 32 : }
576 :
577 : /************************************************************************/
578 : /* SetActiveHDU() */
579 : /************************************************************************/
580 :
581 73 : void FITSLayer::SetActiveHDU()
582 : {
583 73 : int status = 0;
584 73 : fits_movabs_hdu(m_poDS->m_hFITS, m_hduNum, nullptr, &status);
585 73 : if (status != 0)
586 : {
587 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_movabs_hdu() failed: %d",
588 : status);
589 : }
590 73 : }
591 :
592 : /************************************************************************/
593 : /* GetFeatureCount() */
594 : /************************************************************************/
595 :
596 7 : GIntBig FITSLayer::GetFeatureCount(int bForce)
597 : {
598 7 : if (m_poAttrQuery == nullptr && m_poFilterGeom == nullptr)
599 7 : return m_nRows;
600 0 : return OGRLayer::GetFeatureCount(bForce);
601 : }
602 :
603 : /************************************************************************/
604 : /* ResetReading() */
605 : /************************************************************************/
606 :
607 4 : void FITSLayer::ResetReading()
608 : {
609 4 : m_nCurRow = 1;
610 4 : }
611 :
612 : /************************************************************************/
613 : /* ReadCol */
614 : /************************************************************************/
615 :
616 : template <typename T_FITS, typename T_GDAL, int TYPECODE> struct ReadCol
617 : {
618 545 : static void Read(fitsfile *hFITS, const ColDesc &colDesc, int iField,
619 : LONGLONG irow, OGRFeature *poFeature, int nRepeat)
620 : {
621 545 : int status = 0;
622 1090 : std::vector<T_FITS> x(nRepeat);
623 545 : fits_read_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat, nullptr,
624 : x.data(), nullptr, &status);
625 561 : if (nRepeat == 1 && colDesc.bHasNull &&
626 16 : x[0] == static_cast<T_FITS>(colDesc.nNullValue))
627 : {
628 5 : poFeature->SetFieldNull(iField);
629 : }
630 540 : else if (colDesc.dfScale != 1.0 || colDesc.dfOffset != 0.0)
631 : {
632 128 : std::vector<double> scaled;
633 64 : scaled.reserve(nRepeat);
634 128 : for (int i = 0; i < nRepeat; ++i)
635 : {
636 128 : scaled.push_back(static_cast<double>(x[i]) * colDesc.dfScale +
637 64 : colDesc.dfOffset);
638 : }
639 128 : poFeature->SetField(iField, nRepeat, scaled.data());
640 : }
641 476 : else if (nRepeat == 1)
642 : {
643 202 : poFeature->SetField(iField, static_cast<T_GDAL>(x[0]));
644 : }
645 : else
646 : {
647 548 : std::vector<T_GDAL> xGDAL;
648 274 : xGDAL.reserve(nRepeat);
649 890 : for (int i = 0; i < nRepeat; i++)
650 616 : xGDAL.push_back(x[i]);
651 274 : poFeature->SetField(iField, nRepeat, xGDAL.data());
652 : }
653 545 : }
654 : };
655 :
656 : /************************************************************************/
657 : /* GetNextRawFeature() */
658 : /************************************************************************/
659 :
660 23 : OGRFeature *FITSLayer::GetNextRawFeature()
661 : {
662 23 : auto poFeature = GetFeature(m_nCurRow);
663 23 : if (poFeature)
664 19 : m_nCurRow++;
665 23 : return poFeature;
666 : }
667 :
668 : /************************************************************************/
669 : /* GetFeature() */
670 : /************************************************************************/
671 :
672 29 : OGRFeature *FITSLayer::GetFeature(GIntBig nFID)
673 : {
674 29 : LONGLONG nRow = static_cast<LONGLONG>(nFID);
675 29 : if (nRow <= 0 || nRow > m_nRows)
676 6 : return nullptr;
677 :
678 23 : RunDeferredFieldCreation();
679 :
680 23 : OGRFeature *poFeature = new OGRFeature(m_poFeatureDefn);
681 :
682 23 : SetActiveHDU();
683 :
684 1761 : const auto ReadField = [this, &poFeature, nRow](const ColDesc &colDesc,
685 : int iField, char typechar,
686 3009 : int nRepeat)
687 : {
688 1761 : int status = 0;
689 1761 : if (typechar == 'L')
690 : {
691 128 : std::vector<char> x(nRepeat);
692 64 : fits_read_col(m_poDS->m_hFITS, TLOGICAL, colDesc.iCol, nRow, 1,
693 64 : nRepeat, nullptr, x.data(), nullptr, &status);
694 64 : if (nRepeat == 1)
695 : {
696 16 : poFeature->SetField(iField, x[0] == '1' ? 1 : 0);
697 : }
698 : else
699 : {
700 96 : std::vector<int> intValues;
701 48 : intValues.reserve(nRepeat);
702 134 : for (int i = 0; i < nRepeat; ++i)
703 : {
704 86 : intValues.push_back(x[i] == '1' ? 1 : 0);
705 : }
706 48 : poFeature->SetField(iField, nRepeat, intValues.data());
707 : }
708 : }
709 1697 : else if (typechar == 'X')
710 : {
711 782 : char x = 0;
712 782 : fits_read_col_bit(m_poDS->m_hFITS, colDesc.iCol, nRow, colDesc.iBit,
713 : 1, &x, &status);
714 782 : poFeature->SetField(iField, x);
715 : }
716 915 : else if (typechar == 'B')
717 : {
718 96 : if (colDesc.nTypeCode == TSBYTE)
719 : {
720 16 : ReadCol<signed char, int, TSBYTE>::Read(
721 16 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
722 : }
723 : else
724 : {
725 80 : ReadCol<GByte, int, TBYTE>::Read(
726 80 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
727 : }
728 : }
729 819 : else if (typechar == 'I')
730 : {
731 101 : if (colDesc.nTypeCode == TUSHORT)
732 : {
733 16 : ReadCol<GUInt16, int, TUSHORT>::Read(
734 16 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
735 : }
736 : else
737 : {
738 85 : ReadCol<GInt16, int, TSHORT>::Read(
739 85 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
740 : }
741 : }
742 718 : else if (typechar == 'J')
743 : {
744 171 : if (colDesc.nTypeCode == TUINT)
745 : {
746 16 : ReadCol<GUInt32, GIntBig, TUINT>::Read(
747 16 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
748 : }
749 : else
750 : {
751 155 : ReadCol<GInt32, int, TINT>::Read(
752 155 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
753 : }
754 : }
755 547 : else if (typechar == 'K')
756 : {
757 92 : ReadCol<GInt64, GIntBig, TLONGLONG>::Read(
758 92 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
759 : }
760 455 : else if (typechar == 'A') // Character
761 : {
762 109 : if (colDesc.nItems > 1)
763 : {
764 32 : CPLStringList aosList;
765 64 : for (int iItem = 1; iItem <= colDesc.nItems; iItem++)
766 : {
767 96 : std::string osStr;
768 48 : if (nRepeat)
769 : {
770 48 : osStr.resize(nRepeat);
771 48 : char *pszStr = &osStr[0];
772 48 : fits_read_col_str(m_poDS->m_hFITS, colDesc.iCol, nRow,
773 : iItem, 1, nullptr, &pszStr, nullptr,
774 : &status);
775 : }
776 48 : aosList.AddString(osStr.c_str());
777 : }
778 16 : poFeature->SetField(iField, aosList.List());
779 : }
780 : else
781 : {
782 186 : std::string osStr;
783 93 : if (nRepeat)
784 : {
785 93 : osStr.resize(nRepeat);
786 93 : char *pszStr = &osStr[0];
787 93 : fits_read_col_str(m_poDS->m_hFITS, colDesc.iCol, nRow, 1, 1,
788 : nullptr, &pszStr, nullptr, &status);
789 : }
790 93 : poFeature->SetField(iField, osStr.c_str());
791 : }
792 : }
793 346 : else if (typechar == 'E') // IEEE754 32bit
794 : {
795 85 : ReadCol<float, double, TFLOAT>::Read(
796 85 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature, nRepeat);
797 : }
798 261 : else if (typechar == 'D') // IEEE754 64bit
799 : {
800 258 : std::vector<double> x(nRepeat);
801 129 : fits_read_col(m_poDS->m_hFITS, TDOUBLE, colDesc.iCol, nRow, 1,
802 129 : nRepeat, nullptr, x.data(), nullptr, &status);
803 129 : if (nRepeat == 1)
804 83 : poFeature->SetField(iField, x[0]);
805 : else
806 46 : poFeature->SetField(iField, nRepeat, x.data());
807 : }
808 132 : else if (typechar == 'C') // IEEE754 32bit complex
809 : {
810 132 : std::vector<float> x(2 * nRepeat);
811 66 : fits_read_col(m_poDS->m_hFITS, TCOMPLEX, colDesc.iCol, nRow, 1,
812 66 : nRepeat, nullptr, x.data(), nullptr, &status);
813 132 : CPLStringList aosList;
814 159 : for (int i = 0; i < nRepeat; ++i)
815 : aosList.AddString(
816 93 : CPLSPrintf("%.17g + %.17gj", x[2 * i + 0], x[2 * i + 1]));
817 66 : if (nRepeat == 1)
818 34 : poFeature->SetField(iField, aosList[0]);
819 : else
820 32 : poFeature->SetField(iField, aosList.List());
821 : }
822 66 : else if (typechar == 'M') // IEEE754 64bit complex
823 : {
824 132 : std::vector<double> x(2 * nRepeat);
825 66 : fits_read_col(m_poDS->m_hFITS, TDBLCOMPLEX, colDesc.iCol, nRow, 1,
826 66 : nRepeat, nullptr, x.data(), nullptr, &status);
827 132 : CPLStringList aosList;
828 159 : for (int i = 0; i < nRepeat; ++i)
829 : {
830 : aosList.AddString(
831 93 : CPLSPrintf("%.17g + %.17gj", x[2 * i + 0], x[2 * i + 1]));
832 : }
833 66 : if (nRepeat == 1)
834 34 : poFeature->SetField(iField, aosList[0]);
835 : else
836 32 : poFeature->SetField(iField, aosList.List());
837 : }
838 : else
839 : {
840 0 : CPLDebug("FITS", "Unhandled typechar %c", typechar);
841 : }
842 1761 : if (status)
843 : {
844 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_read_col() failed");
845 : }
846 1761 : };
847 :
848 23 : const int nFieldCount = poFeature->GetFieldCount();
849 1784 : for (int iField = 0; iField < nFieldCount; iField++)
850 : {
851 1761 : const auto &colDesc = m_aoColDescs[iField];
852 1761 : if (colDesc.typechar[0] == 'P' || colDesc.typechar[0] == 'Q')
853 : {
854 295 : int status = 0;
855 295 : LONGLONG nRepeat = 0;
856 295 : fits_read_descriptll(m_poDS->m_hFITS, colDesc.iCol, nRow, &nRepeat,
857 : nullptr, &status);
858 295 : ReadField(colDesc, iField, colDesc.typechar[1],
859 : static_cast<int>(nRepeat));
860 : }
861 : else
862 : {
863 1466 : ReadField(colDesc, iField, colDesc.typechar[0], colDesc.nRepeat);
864 : }
865 : }
866 23 : poFeature->SetFID(nRow);
867 23 : return poFeature;
868 : }
869 :
870 : /************************************************************************/
871 : /* TestCapability() */
872 : /************************************************************************/
873 :
874 353 : int FITSLayer::TestCapability(const char *pszCap)
875 : {
876 353 : if (EQUAL(pszCap, OLCFastFeatureCount))
877 1 : return m_poAttrQuery == nullptr && m_poFilterGeom == nullptr;
878 :
879 352 : if (EQUAL(pszCap, OLCRandomRead))
880 1 : return true;
881 :
882 351 : if (EQUAL(pszCap, OLCCreateField) || EQUAL(pszCap, OLCSequentialWrite) ||
883 25 : EQUAL(pszCap, OLCRandomWrite) || EQUAL(pszCap, OLCDeleteFeature))
884 : {
885 332 : return m_poDS->GetAccess() == GA_Update;
886 : }
887 :
888 19 : return false;
889 : }
890 :
891 : /************************************************************************/
892 : /* RunDeferredFieldCreation() */
893 : /************************************************************************/
894 :
895 55 : void FITSLayer::RunDeferredFieldCreation(const OGRFeature *poFeature)
896 : {
897 55 : if (m_anDeferredFieldsIndices.empty())
898 50 : return;
899 :
900 5 : SetActiveHDU();
901 :
902 10 : CPLString osPendingBitFieldName{};
903 5 : int nPendingBitFieldSize = 0;
904 10 : std::set<CPLString> oSetBitFieldNames{};
905 :
906 161 : const auto FlushCreationPendingBitField = [this, &osPendingBitFieldName,
907 : &nPendingBitFieldSize,
908 505 : &oSetBitFieldNames]()
909 : {
910 161 : if (osPendingBitFieldName.empty())
911 153 : return;
912 :
913 : const int iCol =
914 8 : m_aoColDescs.empty() ? 1 : m_aoColDescs.back().iCol + 1;
915 144 : for (int iBit = 1; iBit <= nPendingBitFieldSize; iBit++)
916 : {
917 272 : ColDesc oCol;
918 136 : oCol.iCol = iCol;
919 136 : oCol.iBit = iBit;
920 136 : oCol.typechar = 'X';
921 136 : m_aoColDescs.emplace_back(oCol);
922 : }
923 :
924 8 : int status = 0;
925 8 : CPLString osTForm;
926 8 : osTForm.Printf("%dX", nPendingBitFieldSize);
927 8 : fits_insert_col(m_poDS->m_hFITS, iCol, &osPendingBitFieldName[0],
928 8 : &osTForm[0], &status);
929 8 : if (status != 0)
930 : {
931 0 : CPLError(CE_Failure, CPLE_AppDefined,
932 : "fits_insert_col() failed: %d", status);
933 : }
934 :
935 8 : oSetBitFieldNames.insert(osPendingBitFieldName);
936 8 : osPendingBitFieldName.clear();
937 8 : nPendingBitFieldSize = 0;
938 5 : };
939 :
940 : const bool bRepeatFromFirstFeature =
941 9 : poFeature != nullptr &&
942 4 : EQUAL(m_aosCreationOptions.FetchNameValueDef("COMPUTE_REPEAT",
943 : "AT_FIELD_CREATION"),
944 : "AT_FIRST_FEATURE_CREATION");
945 :
946 5 : char **papszMD = GetMetadata();
947 5 : bool bFirstMD = true;
948 :
949 10 : std::map<CPLString, std::map<CPLString, CPLString>> oMapColNameToMetadata;
950 :
951 : // Remap column related metadata (likely coming from source FITS) to
952 : // actual column numbers
953 10 : std::map<int, CPLString> oMapFITSMDColToName;
954 211 : for (auto papszIter = papszMD; papszIter && *papszIter; ++papszIter)
955 : {
956 206 : char *pszKey = nullptr;
957 206 : const char *pszValue = CPLParseNameValue(*papszIter, &pszKey);
958 206 : if (pszKey && pszValue)
959 : {
960 206 : bool bIgnore = false;
961 341 : for (const char *pszPrefix :
962 : {"TTYPE", "TFORM", "TUNIT", "TNULL", "TSCAL", "TZERO", "TDISP",
963 : "TDIM", "TBCOL", "TCTYP", "TCUNI", "TCRPX", "TCRVL", "TCDLT",
964 547 : "TRPOS"})
965 : {
966 537 : if (STARTS_WITH(pszKey, pszPrefix))
967 : {
968 196 : const char *pszCol = pszKey + strlen(pszPrefix);
969 196 : const int nCol = atoi(pszCol);
970 196 : if (!EQUAL(pszPrefix, "TTYPE"))
971 : {
972 109 : auto oIter = oMapFITSMDColToName.find(nCol);
973 218 : CPLString osColName;
974 109 : if (oIter == oMapFITSMDColToName.end())
975 : {
976 87 : const char *pszColName = CSLFetchNameValue(
977 : papszMD,
978 174 : (std::string("TTYPE") + pszCol).c_str());
979 87 : if (pszColName)
980 : {
981 87 : osColName = pszColName;
982 87 : osColName.Trim();
983 87 : oMapFITSMDColToName[nCol] = osColName;
984 : }
985 : }
986 : else
987 : {
988 22 : osColName = oIter->second;
989 : }
990 109 : if (!osColName.empty())
991 : {
992 218 : oMapColNameToMetadata[osColName][pszPrefix] =
993 327 : CPLString(pszValue).Trim();
994 : }
995 : }
996 196 : bIgnore = true;
997 196 : break;
998 : }
999 : }
1000 :
1001 206 : if (!bIgnore && strlen(pszKey) <= 8 && !EQUAL(pszKey, "TFIELDS") &&
1002 5 : !EQUAL(pszKey, "EXTNAME"))
1003 : {
1004 0 : if (bFirstMD)
1005 : {
1006 0 : int status = 0;
1007 0 : fits_write_key_longwarn(m_poDS->m_hFITS, &status);
1008 0 : bFirstMD = false;
1009 : }
1010 :
1011 0 : char *pszValueNonConst = const_cast<char *>(pszValue);
1012 0 : int status = 0;
1013 0 : fits_update_key_longstr(m_poDS->m_hFITS, pszKey,
1014 : pszValueNonConst, nullptr, &status);
1015 : }
1016 : }
1017 206 : CPLFree(pszKey);
1018 : }
1019 :
1020 298 : for (const int nFieldIdx : m_anDeferredFieldsIndices)
1021 : {
1022 293 : const auto poFieldDefn = m_poFeatureDefn->GetFieldDefn(nFieldIdx);
1023 293 : const auto &oFieldDefn = *poFieldDefn;
1024 293 : const char *pszFieldName = oFieldDefn.GetNameRef();
1025 293 : const OGRFieldType eType = oFieldDefn.GetType();
1026 293 : const OGRFieldSubType eSubType = oFieldDefn.GetSubType();
1027 :
1028 293 : if (eType == OFTInteger)
1029 : {
1030 160 : const char *pszBit = strstr(pszFieldName, "_bit");
1031 160 : long iBit = 0;
1032 160 : char *pszEndPtr = nullptr;
1033 136 : if (pszBit &&
1034 136 : (iBit = strtol(pszBit + strlen("_bit"), &pszEndPtr, 10)) > 0 &&
1035 296 : pszEndPtr && *pszEndPtr == '\0')
1036 : {
1037 136 : std::string osName;
1038 136 : osName.assign(pszFieldName, pszBit - pszFieldName);
1039 136 : if (oSetBitFieldNames.find(osName) == oSetBitFieldNames.end())
1040 : {
1041 268 : if (!osPendingBitFieldName.empty() &&
1042 132 : osPendingBitFieldName != osName)
1043 : {
1044 4 : FlushCreationPendingBitField();
1045 : }
1046 :
1047 136 : if (osPendingBitFieldName.empty())
1048 : {
1049 8 : osPendingBitFieldName = std::move(osName);
1050 8 : nPendingBitFieldSize = 1;
1051 8 : continue;
1052 : }
1053 128 : else if (iBit == nPendingBitFieldSize + 1)
1054 : {
1055 128 : nPendingBitFieldSize++;
1056 128 : continue;
1057 : }
1058 : }
1059 : }
1060 : }
1061 :
1062 157 : FlushCreationPendingBitField();
1063 :
1064 314 : CPLString osTForm;
1065 314 : ColDesc oCol;
1066 157 : oCol.iCol = m_aoColDescs.empty() ? 1 : m_aoColDescs.back().iCol + 1;
1067 157 : oCol.nRepeat = 1;
1068 :
1069 314 : CPLString osRepeat;
1070 157 : const char *pszRepeat = m_aosCreationOptions.FetchNameValue(
1071 314 : (CPLString("REPEAT_") + pszFieldName).c_str());
1072 :
1073 : const auto &osTFormFromMD =
1074 157 : oMapColNameToMetadata[pszFieldName]["TFORM"];
1075 :
1076 : // For fields of type list, determine if we can know if it has a fixed
1077 : // number of elements
1078 157 : if (eType == OFTIntegerList || eType == OFTInteger64List ||
1079 : eType == OFTRealList)
1080 : {
1081 : // First take into account the REPEAT_{FIELD_NAME} creatin option
1082 64 : if (pszRepeat)
1083 : {
1084 1 : osRepeat = pszRepeat;
1085 1 : oCol.nRepeat = atoi(pszRepeat);
1086 : }
1087 : // Then if COMPUTE_REPEAT=AT_FIRST_FEATURE_CREATION was specified
1088 : // and we have a feature, then look at the number of items in the
1089 : // field
1090 78 : else if (bRepeatFromFirstFeature &&
1091 15 : poFeature->IsFieldSetAndNotNull(nFieldIdx))
1092 : {
1093 15 : int nCount = 0;
1094 15 : if (eType == OFTIntegerList)
1095 : {
1096 10 : CPL_IGNORE_RET_VAL(
1097 10 : poFeature->GetFieldAsIntegerList(nFieldIdx, &nCount));
1098 : }
1099 5 : else if (eType == OFTInteger64List)
1100 : {
1101 2 : CPL_IGNORE_RET_VAL(
1102 2 : poFeature->GetFieldAsInteger64List(nFieldIdx, &nCount));
1103 : }
1104 3 : else if (eType == OFTRealList)
1105 : {
1106 3 : CPL_IGNORE_RET_VAL(
1107 3 : poFeature->GetFieldAsDoubleList(nFieldIdx, &nCount));
1108 : }
1109 : else
1110 : {
1111 0 : CPLAssert(false);
1112 : }
1113 15 : osRepeat.Printf("%d", nCount);
1114 15 : oCol.nRepeat = nCount;
1115 : }
1116 : else
1117 : {
1118 64 : if (!osTFormFromMD.empty() && osTFormFromMD[0] >= '1' &&
1119 16 : osTFormFromMD[0] <= '9')
1120 : {
1121 8 : oCol.nRepeat = atoi(osTFormFromMD.c_str());
1122 8 : osRepeat.Printf("%d", oCol.nRepeat);
1123 : }
1124 : else
1125 : {
1126 40 : oCol.nRepeat = 0;
1127 40 : oCol.typechar = 'P';
1128 : }
1129 64 : }
1130 : }
1131 93 : else if (pszRepeat)
1132 : {
1133 0 : CPLError(CE_Warning, CPLE_AppDefined,
1134 : "%s ignored on a non-List data type",
1135 0 : (CPLString("REPEAT_") + pszFieldName).c_str());
1136 : }
1137 :
1138 157 : switch (eType)
1139 : {
1140 64 : case OFTIntegerList:
1141 : case OFTInteger:
1142 64 : if (eSubType == OFSTInt16)
1143 : {
1144 12 : oCol.typechar += 'I';
1145 12 : oCol.nTypeCode = TSHORT;
1146 : }
1147 : else
1148 : {
1149 52 : oCol.typechar += 'J';
1150 52 : oCol.nTypeCode = TINT;
1151 : }
1152 64 : break;
1153 :
1154 16 : case OFTInteger64List:
1155 : case OFTInteger64:
1156 16 : oCol.typechar += 'K';
1157 16 : oCol.nTypeCode = TLONGLONG;
1158 16 : break;
1159 :
1160 49 : case OFTRealList:
1161 : case OFTReal:
1162 49 : if (eSubType == OFSTFloat32)
1163 : {
1164 12 : oCol.typechar += 'E';
1165 12 : oCol.nTypeCode = TFLOAT;
1166 : }
1167 : else
1168 : {
1169 37 : oCol.typechar += 'D';
1170 37 : oCol.nTypeCode = TDOUBLE;
1171 : }
1172 49 : break;
1173 :
1174 28 : case OFTString:
1175 28 : if (osTFormFromMD == "C")
1176 : {
1177 2 : oCol.typechar = "C";
1178 2 : oCol.nTypeCode = TCOMPLEX;
1179 : }
1180 26 : else if (osTFormFromMD == "M")
1181 : {
1182 2 : oCol.typechar = "M";
1183 2 : oCol.nTypeCode = TDBLCOMPLEX;
1184 : }
1185 : else
1186 : {
1187 24 : if (oFieldDefn.GetWidth() == 0)
1188 : {
1189 16 : oCol.typechar = "PA";
1190 : }
1191 : else
1192 : {
1193 8 : oCol.typechar = "A";
1194 8 : oCol.nRepeat = oFieldDefn.GetWidth();
1195 8 : osTForm = CPLSPrintf("%dA", oCol.nRepeat);
1196 : }
1197 24 : oCol.nTypeCode = TSTRING;
1198 : }
1199 28 : break;
1200 :
1201 0 : default:
1202 0 : CPLError(CE_Failure, CPLE_NotSupported,
1203 : "Unsupported field type: should not happen");
1204 0 : break;
1205 : }
1206 :
1207 314 : CPLString osTType(pszFieldName);
1208 157 : if (osTForm.empty())
1209 : {
1210 109 : if ((eType == OFTIntegerList || eType == OFTInteger64List ||
1211 298 : eType == OFTRealList) &&
1212 64 : !osRepeat.empty())
1213 : {
1214 24 : osTForm = osRepeat;
1215 24 : osTForm += oCol.typechar;
1216 : }
1217 : else
1218 : {
1219 125 : osTForm = oCol.typechar;
1220 : }
1221 : }
1222 157 : CPL_IGNORE_RET_VAL(osRepeat);
1223 157 : int status = 0;
1224 157 : fits_insert_col(m_poDS->m_hFITS, oCol.iCol, &osTType[0], &osTForm[0],
1225 : &status);
1226 157 : if (status != 0)
1227 : {
1228 0 : CPLError(CE_Failure, CPLE_AppDefined,
1229 : "fits_insert_col() failed: %d", status);
1230 : }
1231 :
1232 : // Set unit from metadata
1233 471 : auto osUnit = oMapColNameToMetadata[pszFieldName]["TUNIT"];
1234 157 : if (!osUnit.empty())
1235 : {
1236 0 : CPLString osKey;
1237 0 : osKey.Printf("TUNIT%d", oCol.iCol);
1238 0 : fits_update_key_longstr(m_poDS->m_hFITS, osKey.data(),
1239 0 : osUnit.data(), nullptr, &status);
1240 : }
1241 :
1242 157 : m_aoColDescs.emplace_back(oCol);
1243 : }
1244 :
1245 5 : m_anDeferredFieldsIndices.clear();
1246 5 : CPLAssert(static_cast<int>(m_aoColDescs.size()) ==
1247 : m_poFeatureDefn->GetFieldCount());
1248 : }
1249 :
1250 : /************************************************************************/
1251 : /* CreateField() */
1252 : /************************************************************************/
1253 :
1254 313 : OGRErr FITSLayer::CreateField(const OGRFieldDefn *poField, int /* bApproxOK */)
1255 : {
1256 313 : if (!TestCapability(OLCCreateField))
1257 0 : return OGRERR_FAILURE;
1258 313 : if (m_poFeatureDefn->GetFieldIndex(poField->GetNameRef()) >= 0)
1259 : {
1260 0 : CPLError(CE_Failure, CPLE_NotSupported,
1261 : "A field with name %s already exists", poField->GetNameRef());
1262 0 : return OGRERR_FAILURE;
1263 : }
1264 313 : if (poField->GetType() == OFTStringList)
1265 : {
1266 20 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported field type");
1267 20 : return OGRERR_FAILURE;
1268 : }
1269 :
1270 293 : m_anDeferredFieldsIndices.emplace_back(m_poFeatureDefn->GetFieldCount());
1271 293 : m_poFeatureDefn->AddFieldDefn(poField);
1272 293 : return OGRERR_NONE;
1273 : }
1274 :
1275 : /************************************************************************/
1276 : /* WriteCol */
1277 : /************************************************************************/
1278 :
1279 0 : template <class T> static T Round(double dfValue)
1280 : {
1281 0 : return static_cast<T>(floor(dfValue + 0.5));
1282 : }
1283 :
1284 0 : template <> double Round<double>(double dfValue)
1285 : {
1286 0 : return dfValue;
1287 : }
1288 :
1289 0 : template <> float Round<float>(double dfValue)
1290 : {
1291 0 : return static_cast<float>(dfValue);
1292 : }
1293 :
1294 : template <typename T_FITS, typename T_GDAL, int TYPECODE,
1295 : T_GDAL (OGRFeature::*GetFieldMethod)(int) const,
1296 : const T_GDAL *(OGRFeature::*GetFieldListMethod)(int, int *) const>
1297 : struct WriteCol
1298 : {
1299 449 : static int Write(fitsfile *hFITS, const ColDesc &colDesc, int iField,
1300 : LONGLONG irow, const OGRFeature *poFeature)
1301 : {
1302 449 : int status = 0;
1303 449 : int nRepeat = colDesc.nRepeat;
1304 449 : const auto poFieldDefn = poFeature->GetFieldDefnRef(iField);
1305 449 : const auto eOGRType = poFieldDefn->GetType();
1306 449 : int nCount = 0;
1307 : // cppcheck-suppress constStatement
1308 449 : const T_GDAL *panList =
1309 309 : (eOGRType == OFTIntegerList || eOGRType == OFTInteger64List ||
1310 : eOGRType == OFTRealList)
1311 758 : ? (poFeature->*GetFieldListMethod)(iField, &nCount)
1312 : : nullptr;
1313 449 : if (panList)
1314 : {
1315 224 : nRepeat = nRepeat == 0 ? nCount : std::min(nRepeat, nCount);
1316 224 : if (nCount > nRepeat)
1317 : {
1318 8 : CPLError(CE_Warning, CPLE_AppDefined,
1319 : "Field %s of feature " CPL_FRMT_GIB " had %d "
1320 : "elements, but had to be truncated to %d",
1321 : poFieldDefn->GetNameRef(), static_cast<GIntBig>(irow),
1322 : nCount, nRepeat);
1323 : }
1324 : }
1325 : else
1326 : {
1327 225 : nRepeat = 1;
1328 : }
1329 :
1330 449 : if (nRepeat == 0)
1331 32 : return 0;
1332 :
1333 417 : if (colDesc.bHasNull && nRepeat == 1 && poFeature->IsFieldNull(iField))
1334 : {
1335 0 : T_FITS x = static_cast<T_FITS>(colDesc.nNullValue);
1336 0 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat, &x,
1337 : &status);
1338 : }
1339 417 : else if (nRepeat == 1)
1340 : {
1341 225 : const auto val =
1342 225 : panList ? panList[0] : (poFeature->*GetFieldMethod)(iField);
1343 225 : if (colDesc.dfScale != 1.0 || colDesc.dfOffset != 0.0)
1344 : {
1345 0 : T_FITS x =
1346 0 : Round<T_FITS>((val - colDesc.dfOffset) / colDesc.dfScale);
1347 0 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat,
1348 0 : &x, &status);
1349 : }
1350 : else
1351 : {
1352 225 : T_FITS x = static_cast<T_FITS>(val);
1353 225 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat,
1354 : &x, &status);
1355 : }
1356 : }
1357 : else
1358 : {
1359 192 : CPLAssert(panList);
1360 :
1361 384 : std::vector<T_FITS> x;
1362 192 : x.reserve(nRepeat);
1363 192 : if (colDesc.dfScale != 1.0 || colDesc.dfOffset != 0.0)
1364 : {
1365 0 : for (int i = 0; i < nRepeat; i++)
1366 : {
1367 0 : x.push_back(Round<T_FITS>((panList[i] - colDesc.dfOffset) /
1368 0 : colDesc.dfScale));
1369 : }
1370 0 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat,
1371 0 : x.data(), &status);
1372 : }
1373 : else
1374 : {
1375 656 : for (int i = 0; i < nRepeat; i++)
1376 : {
1377 464 : x.push_back(static_cast<T_FITS>(panList[i]));
1378 : }
1379 192 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat,
1380 : x.data(), &status);
1381 : }
1382 : }
1383 417 : return status;
1384 : }
1385 : };
1386 :
1387 : /************************************************************************/
1388 : /* WriteComplex */
1389 : /************************************************************************/
1390 :
1391 : template <typename T, int TYPECODE> struct WriteComplex
1392 : {
1393 12 : static int Write(fitsfile *hFITS, const ColDesc &colDesc, int iField,
1394 : LONGLONG irow, const OGRFeature *poFeature)
1395 : {
1396 12 : int status = 0;
1397 12 : const auto poFieldDefn = poFeature->GetFieldDefnRef(iField);
1398 12 : if (poFieldDefn->GetType() == OFTStringList)
1399 : {
1400 0 : auto papszStrings = poFeature->GetFieldAsStringList(iField);
1401 0 : const int nCount = CSLCount(papszStrings);
1402 0 : int nRepeat = colDesc.nRepeat;
1403 0 : nRepeat = nRepeat == 0 ? nCount : std::min(nRepeat, nCount);
1404 0 : if (nRepeat > nCount)
1405 : {
1406 0 : CPLError(CE_Warning, CPLE_AppDefined,
1407 : "Field %s of feature " CPL_FRMT_GIB " had %d "
1408 : "elements, but had to be truncated to %d",
1409 : poFieldDefn->GetNameRef(), static_cast<GIntBig>(irow),
1410 : nRepeat, nCount);
1411 : }
1412 0 : std::vector<T> x(2 * nRepeat);
1413 0 : for (int i = 0; i < nRepeat; i++)
1414 : {
1415 0 : double re = 0;
1416 0 : double im = 0;
1417 0 : CPLsscanf(papszStrings[i], "%lf + %lfj", &re, &im);
1418 0 : x[2 * i] = static_cast<T>(re);
1419 0 : x[2 * i + 1] = static_cast<T>(im);
1420 : }
1421 0 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, nRepeat,
1422 : x.data(), &status);
1423 : }
1424 : else
1425 : {
1426 24 : std::vector<T> x(2);
1427 12 : double re = 0;
1428 12 : double im = 0;
1429 12 : CPLsscanf(poFeature->GetFieldAsString(iField), "%lf + %lfj", &re,
1430 : &im);
1431 12 : x[0] = static_cast<T>(re);
1432 12 : x[1] = static_cast<T>(im);
1433 12 : fits_write_col(hFITS, TYPECODE, colDesc.iCol, irow, 1, 1, x.data(),
1434 : &status);
1435 : }
1436 12 : return status;
1437 : }
1438 : };
1439 :
1440 : /************************************************************************/
1441 : /* SetOrCreateFeature() */
1442 : /************************************************************************/
1443 :
1444 14 : bool FITSLayer::SetOrCreateFeature(const OGRFeature *poFeature, LONGLONG nRow)
1445 : {
1446 14 : SetActiveHDU();
1447 :
1448 3631 : const auto WriteField = [this, &poFeature, nRow](int iField)
1449 : {
1450 1023 : const auto poFieldDefn = poFeature->GetFieldDefnRef(iField);
1451 1023 : const auto &colDesc = m_aoColDescs[iField];
1452 : const char typechar =
1453 1836 : (colDesc.typechar[0] == 'P' || colDesc.typechar[0] == 'Q')
1454 1836 : ? colDesc.typechar[1]
1455 813 : : colDesc.typechar[0];
1456 1023 : int nRepeat = colDesc.nRepeat;
1457 1023 : int status = 0;
1458 1023 : if (typechar == 'L')
1459 : {
1460 0 : const auto ToLogical = [](int x) -> char { return x ? '1' : '0'; };
1461 :
1462 0 : if (poFieldDefn->GetType() == OFTIntegerList)
1463 : {
1464 0 : int nCount = 0;
1465 : const int *panVals =
1466 0 : poFeature->GetFieldAsIntegerList(iField, &nCount);
1467 0 : nRepeat = nRepeat == 0 ? nCount : std::min(nRepeat, nCount);
1468 0 : if (nRepeat > nCount)
1469 : {
1470 0 : CPLError(CE_Warning, CPLE_AppDefined,
1471 : "Field %s of feature " CPL_FRMT_GIB " had %d "
1472 : "elements, but had to be truncated to %d",
1473 : poFieldDefn->GetNameRef(),
1474 : static_cast<GIntBig>(nRow), nRepeat, nCount);
1475 : }
1476 0 : std::vector<char> x(nRepeat);
1477 0 : for (int i = 0; i < nRepeat; i++)
1478 : {
1479 0 : x[i] = ToLogical(panVals[i]);
1480 : }
1481 0 : fits_write_col(m_poDS->m_hFITS, TLOGICAL, colDesc.iCol, nRow, 1,
1482 0 : nRepeat, x.data(), &status);
1483 : }
1484 : else
1485 : {
1486 0 : char x = ToLogical(poFeature->GetFieldAsInteger(iField));
1487 0 : fits_write_col(m_poDS->m_hFITS, TLOGICAL, colDesc.iCol, nRow, 1,
1488 : nRepeat, &x, &status);
1489 : }
1490 : }
1491 1023 : else if (typechar == 'X') // bit array
1492 : {
1493 476 : char flag = poFeature->GetFieldAsInteger(iField) ? 0x80 : 0;
1494 476 : fits_write_col_bit(m_poDS->m_hFITS, colDesc.iCol, nRow,
1495 476 : colDesc.iBit, 1, &flag, &status);
1496 : }
1497 547 : else if (typechar == 'B')
1498 : {
1499 0 : if (colDesc.nTypeCode == TSBYTE)
1500 : {
1501 0 : status = WriteCol<
1502 : signed char, int, TSBYTE, &OGRFeature::GetFieldAsInteger,
1503 0 : &OGRFeature::GetFieldAsIntegerList>::Write(m_poDS->m_hFITS,
1504 : colDesc, iField,
1505 : nRow, poFeature);
1506 : }
1507 : else
1508 : {
1509 0 : status = WriteCol<
1510 : GByte, int, TBYTE, &OGRFeature::GetFieldAsInteger,
1511 0 : &OGRFeature::GetFieldAsIntegerList>::Write(m_poDS->m_hFITS,
1512 : colDesc, iField,
1513 : nRow, poFeature);
1514 : }
1515 : }
1516 547 : else if (typechar == 'I')
1517 : {
1518 42 : if (colDesc.nTypeCode == TUSHORT)
1519 : {
1520 0 : status = WriteCol<
1521 : GUInt16, int, TUSHORT, &OGRFeature::GetFieldAsInteger,
1522 0 : &OGRFeature::GetFieldAsIntegerList>::Write(m_poDS->m_hFITS,
1523 : colDesc, iField,
1524 : nRow, poFeature);
1525 : }
1526 : else
1527 : {
1528 42 : status = WriteCol<
1529 : GInt16, int, TSHORT, &OGRFeature::GetFieldAsInteger,
1530 42 : &OGRFeature::GetFieldAsIntegerList>::Write(m_poDS->m_hFITS,
1531 : colDesc, iField,
1532 : nRow, poFeature);
1533 : }
1534 : }
1535 505 : else if (typechar == 'J')
1536 : {
1537 182 : if (colDesc.nTypeCode == TUINT)
1538 : {
1539 0 : status = WriteCol<
1540 : GUInt32, GIntBig, TUINT, &OGRFeature::GetFieldAsInteger64,
1541 0 : &OGRFeature::GetFieldAsInteger64List>::Write(m_poDS
1542 : ->m_hFITS,
1543 : colDesc,
1544 : iField, nRow,
1545 : poFeature);
1546 : }
1547 : else
1548 : {
1549 182 : status = WriteCol<
1550 : GInt32, int, TINT, &OGRFeature::GetFieldAsInteger,
1551 182 : &OGRFeature::GetFieldAsIntegerList>::Write(m_poDS->m_hFITS,
1552 : colDesc, iField,
1553 : nRow, poFeature);
1554 : }
1555 : }
1556 323 : else if (typechar == 'K')
1557 : {
1558 56 : status = WriteCol<
1559 : GInt64, GIntBig, TLONGLONG, &OGRFeature::GetFieldAsInteger64,
1560 56 : &OGRFeature::GetFieldAsInteger64List>::Write(m_poDS->m_hFITS,
1561 : colDesc, iField,
1562 : nRow, poFeature);
1563 : }
1564 267 : else if (typechar == 'A') // Character
1565 : {
1566 86 : if (poFieldDefn->GetType() == OFTStringList)
1567 : {
1568 0 : auto papszStrings = poFeature->GetFieldAsStringList(iField);
1569 0 : const int nStringCount = CSLCount(papszStrings);
1570 0 : const int nItems = std::min(colDesc.nItems, nStringCount);
1571 0 : if (nItems > nStringCount)
1572 : {
1573 0 : CPLError(CE_Warning, CPLE_AppDefined,
1574 : "Field %s of feature " CPL_FRMT_GIB " had %d "
1575 : "elements, but had to be truncated to %d",
1576 : poFieldDefn->GetNameRef(),
1577 : static_cast<GIntBig>(nRow), nItems, nStringCount);
1578 : }
1579 0 : fits_write_col_str(m_poDS->m_hFITS, colDesc.iCol, nRow, 1,
1580 : nItems, papszStrings, &status);
1581 : }
1582 : else
1583 : {
1584 : char *pszStr =
1585 86 : const_cast<char *>(poFeature->GetFieldAsString(iField));
1586 86 : fits_write_col_str(m_poDS->m_hFITS, colDesc.iCol, nRow, 1, 1,
1587 : &pszStr, &status);
1588 : }
1589 : }
1590 181 : else if (typechar == 'E') // IEEE754 32bit
1591 : {
1592 42 : status = WriteCol<
1593 : float, double, TFLOAT, &OGRFeature::GetFieldAsDouble,
1594 42 : &OGRFeature::GetFieldAsDoubleList>::Write(m_poDS->m_hFITS,
1595 : colDesc, iField, nRow,
1596 : poFeature);
1597 : }
1598 139 : else if (typechar == 'D') // IEEE754 64bit
1599 : {
1600 127 : status = WriteCol<
1601 : double, double, TDOUBLE, &OGRFeature::GetFieldAsDouble,
1602 127 : &OGRFeature::GetFieldAsDoubleList>::Write(m_poDS->m_hFITS,
1603 : colDesc, iField, nRow,
1604 : poFeature);
1605 : }
1606 12 : else if (typechar == 'C') // IEEE754 32bit complex
1607 : {
1608 6 : status = WriteComplex<float, TCOMPLEX>::Write(
1609 6 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature);
1610 : }
1611 6 : else if (typechar == 'M') // IEEE754 64bit complex
1612 : {
1613 6 : status = WriteComplex<double, TDBLCOMPLEX>::Write(
1614 6 : m_poDS->m_hFITS, colDesc, iField, nRow, poFeature);
1615 : }
1616 : else
1617 : {
1618 0 : CPLDebug("FITS", "Unhandled typechar %c", typechar);
1619 : }
1620 1023 : if (status)
1621 : {
1622 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_write_col() failed");
1623 : }
1624 1023 : return status == 0;
1625 14 : };
1626 :
1627 14 : bool bOK = true;
1628 14 : const int nFieldCount = poFeature->GetFieldCount();
1629 1037 : for (int iField = 0; iField < nFieldCount; iField++)
1630 : {
1631 1023 : if (!WriteField(iField))
1632 0 : bOK = false;
1633 : }
1634 14 : return bOK;
1635 : }
1636 :
1637 : /************************************************************************/
1638 : /* ICreateFeature() */
1639 : /************************************************************************/
1640 :
1641 13 : OGRErr FITSLayer::ICreateFeature(OGRFeature *poFeature)
1642 : {
1643 13 : if (!TestCapability(OLCSequentialWrite))
1644 0 : return OGRERR_FAILURE;
1645 :
1646 13 : RunDeferredFieldCreation(poFeature);
1647 :
1648 13 : m_nRows++;
1649 13 : SetActiveHDU();
1650 13 : const bool bOK = SetOrCreateFeature(poFeature, m_nRows);
1651 13 : poFeature->SetFID(m_nRows);
1652 :
1653 13 : return bOK ? OGRERR_NONE : OGRERR_FAILURE;
1654 : }
1655 :
1656 : /************************************************************************/
1657 : /* ISetFeature() */
1658 : /************************************************************************/
1659 :
1660 3 : OGRErr FITSLayer::ISetFeature(OGRFeature *poFeature)
1661 : {
1662 3 : if (!TestCapability(OLCRandomWrite))
1663 0 : return OGRERR_FAILURE;
1664 :
1665 3 : RunDeferredFieldCreation();
1666 :
1667 3 : const GIntBig nRow = poFeature->GetFID();
1668 3 : if (nRow <= 0 || nRow > m_nRows)
1669 2 : return OGRERR_NON_EXISTING_FEATURE;
1670 :
1671 1 : SetActiveHDU();
1672 1 : const bool bOK = SetOrCreateFeature(poFeature, nRow);
1673 1 : return bOK ? OGRERR_NONE : OGRERR_FAILURE;
1674 : }
1675 :
1676 : /************************************************************************/
1677 : /* DeleteFeature() */
1678 : /************************************************************************/
1679 :
1680 3 : OGRErr FITSLayer::DeleteFeature(GIntBig nFID)
1681 : {
1682 3 : if (!TestCapability(OLCDeleteFeature))
1683 0 : return OGRERR_FAILURE;
1684 :
1685 3 : if (nFID <= 0 || nFID > m_nRows)
1686 2 : return OGRERR_NON_EXISTING_FEATURE;
1687 :
1688 1 : SetActiveHDU();
1689 :
1690 1 : int status = 0;
1691 1 : fits_delete_rows(m_poDS->m_hFITS, nFID, 1, &status);
1692 1 : m_nRows--;
1693 1 : return status == 0 ? OGRERR_NONE : OGRERR_FAILURE;
1694 : }
1695 :
1696 : /************************************************************************/
1697 : /* FITSRasterBand() */
1698 : /************************************************************************/
1699 :
1700 133 : FITSRasterBand::FITSRasterBand(FITSDataset *poDSIn, int nBandIn)
1701 133 : : m_poFDS(poDSIn)
1702 : {
1703 133 : poDS = poDSIn;
1704 133 : nBand = nBandIn;
1705 133 : eDataType = poDSIn->m_gdalDataType;
1706 133 : nBlockXSize = poDSIn->nRasterXSize;
1707 133 : nBlockYSize = 1;
1708 133 : }
1709 :
1710 : /************************************************************************/
1711 : /* ~FITSRasterBand() */
1712 : /************************************************************************/
1713 :
1714 266 : FITSRasterBand::~FITSRasterBand()
1715 : {
1716 133 : FlushCache(true);
1717 266 : }
1718 :
1719 : /************************************************************************/
1720 : /* IReadBlock() */
1721 : /************************************************************************/
1722 :
1723 160 : CPLErr FITSRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
1724 : void *pImage)
1725 : {
1726 : // A FITS block is one row (we assume BSQ formatted data)
1727 160 : FITSDataset *dataset = m_poFDS;
1728 160 : fitsfile *hFITS = dataset->m_hFITS;
1729 160 : int status = 0;
1730 :
1731 : // Since a FITS block is a whole row, nBlockXOff must be zero
1732 : // and the row number equals the y block offset. Also, nBlockYOff
1733 : // cannot be greater than the number of rows
1734 160 : CPLAssert(nBlockXOff == 0);
1735 160 : CPLAssert(nBlockYOff < nRasterYSize);
1736 :
1737 : // Calculate offsets and read in the data. Note that FITS array offsets
1738 : // start at 1 at the bottom left...
1739 160 : LONGLONG offset =
1740 160 : static_cast<LONGLONG>(nBand - 1) * nRasterXSize * nRasterYSize +
1741 160 : (static_cast<LONGLONG>(nRasterYSize - 1 - nBlockYOff) * nRasterXSize +
1742 : 1);
1743 160 : long nElements = nRasterXSize;
1744 :
1745 : // If we haven't written this block to the file yet, then attempting
1746 : // to read causes an error, so in this case, just return zeros.
1747 160 : if (!dataset->m_isExistingFile && offset > dataset->m_highestOffsetWritten)
1748 : {
1749 0 : memset(pImage, 0,
1750 0 : nBlockXSize * nBlockYSize * GDALGetDataTypeSize(eDataType) / 8);
1751 0 : return CE_None;
1752 : }
1753 :
1754 : // Otherwise read in the image data
1755 160 : fits_read_img(hFITS, dataset->m_fitsDataType, offset, nElements, nullptr,
1756 : pImage, nullptr, &status);
1757 :
1758 : // Capture special case of non-zero status due to data range
1759 : // overflow Standard GDAL policy is to silently truncate, which is
1760 : // what CFITSIO does, in addition to returning NUM_OVERFLOW (412) as
1761 : // the status.
1762 160 : if (status == NUM_OVERFLOW)
1763 0 : status = 0;
1764 :
1765 160 : if (status)
1766 : {
1767 0 : CPLError(CE_Failure, CPLE_AppDefined,
1768 : "Couldn't read image data from FITS file (%d).", status);
1769 0 : return CE_Failure;
1770 : }
1771 :
1772 160 : return CE_None;
1773 : }
1774 :
1775 : /************************************************************************/
1776 : /* IWriteBlock() */
1777 : /* */
1778 : /************************************************************************/
1779 :
1780 450 : CPLErr FITSRasterBand::IWriteBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
1781 : void *pImage)
1782 : {
1783 450 : FITSDataset *dataset = m_poFDS;
1784 450 : fitsfile *hFITS = dataset->m_hFITS;
1785 450 : int status = 0;
1786 :
1787 : // Since a FITS block is a whole row, nBlockXOff must be zero
1788 : // and the row number equals the y block offset. Also, nBlockYOff
1789 : // cannot be greater than the number of rows
1790 :
1791 : // Calculate offsets and read in the data. Note that FITS array offsets
1792 : // start at 1 at the bottom left...
1793 450 : LONGLONG offset =
1794 450 : static_cast<LONGLONG>(nBand - 1) * nRasterXSize * nRasterYSize +
1795 450 : (static_cast<LONGLONG>(nRasterYSize - 1 - nBlockYOff) * nRasterXSize +
1796 : 1);
1797 450 : long nElements = nRasterXSize;
1798 450 : fits_write_img(hFITS, dataset->m_fitsDataType, offset, nElements, pImage,
1799 : &status);
1800 :
1801 : // Capture special case of non-zero status due to data range
1802 : // overflow Standard GDAL policy is to silently truncate, which is
1803 : // what CFITSIO does, in addition to returning NUM_OVERFLOW (412) as
1804 : // the status.
1805 450 : if (status == NUM_OVERFLOW)
1806 0 : status = 0;
1807 :
1808 : // Check for other errors
1809 450 : if (status)
1810 : {
1811 0 : CPLError(CE_Failure, CPLE_AppDefined,
1812 : "Error writing image data to FITS file (%d).", status);
1813 0 : return CE_Failure;
1814 : }
1815 :
1816 : // When we write a block, update the offset counter that we've written
1817 450 : if (offset > dataset->m_highestOffsetWritten)
1818 33 : dataset->m_highestOffsetWritten = offset;
1819 :
1820 450 : return CE_None;
1821 : }
1822 :
1823 : /************************************************************************/
1824 : /* ==================================================================== */
1825 : /* FITSDataset */
1826 : /* ==================================================================== */
1827 : /************************************************************************/
1828 :
1829 : // Some useful utility functions
1830 :
1831 : // Simple static function to determine if FITS header keyword should
1832 : // be saved in meta data.
1833 : static const char *const ignorableFITSHeaders[] = {
1834 : "SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "NAXIS3", "END",
1835 : "XTENSION", "PCOUNT", "GCOUNT", "EXTEND", "CONTINUE", "COMMENT", "",
1836 : "LONGSTRN", "BZERO", "BSCALE", "BLANK", "CHECKSUM", "DATASUM",
1837 : };
1838 :
1839 2269 : static bool isIgnorableFITSHeader(const char *name)
1840 : {
1841 39230 : for (const char *keyword : ignorableFITSHeaders)
1842 : {
1843 37667 : if (strcmp(name, keyword) == 0)
1844 706 : return true;
1845 : }
1846 1563 : return false;
1847 : }
1848 :
1849 : /************************************************************************/
1850 : /* FITSDataset() */
1851 : /************************************************************************/
1852 :
1853 103 : FITSDataset::FITSDataset()
1854 : {
1855 103 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1856 103 : m_adfGeoTransform[0] = 0;
1857 103 : m_adfGeoTransform[1] = 1;
1858 103 : m_adfGeoTransform[2] = 0;
1859 103 : m_adfGeoTransform[3] = 0;
1860 103 : m_adfGeoTransform[4] = 0;
1861 103 : m_adfGeoTransform[5] = 1;
1862 103 : }
1863 :
1864 : /************************************************************************/
1865 : /* ~FITSDataset() */
1866 : /************************************************************************/
1867 :
1868 206 : FITSDataset::~FITSDataset()
1869 : {
1870 :
1871 103 : int status = 0;
1872 103 : if (m_hFITS)
1873 : {
1874 103 : m_apoLayers.clear();
1875 :
1876 103 : if (m_hduNum > 0 && eAccess == GA_Update)
1877 : {
1878 : // Only do this if we've successfully opened the file and update
1879 : // capability. Write any meta data to the file that's compatible
1880 : // with FITS.
1881 41 : fits_movabs_hdu(m_hFITS, m_hduNum, nullptr, &status);
1882 41 : fits_write_key_longwarn(m_hFITS, &status);
1883 41 : if (status)
1884 : {
1885 0 : CPLError(CE_Warning, CPLE_AppDefined,
1886 : "Couldn't move to HDU %d in FITS file %s (%d).\n",
1887 : m_hduNum, GetDescription(), status);
1888 : }
1889 41 : char **metaData = FITSDataset::GetMetadata();
1890 41 : int count = CSLCount(metaData);
1891 65 : for (int i = 0; i < count; ++i)
1892 : {
1893 24 : const char *field = CSLGetField(metaData, i);
1894 24 : if (strlen(field) == 0)
1895 0 : continue;
1896 : else
1897 : {
1898 24 : char *key = nullptr;
1899 24 : const char *value = CPLParseNameValue(field, &key);
1900 : // FITS keys must be less than 8 chars
1901 26 : if (key != nullptr && strlen(key) <= 8 &&
1902 2 : !isIgnorableFITSHeader(key))
1903 : {
1904 : // Although FITS provides support for different value
1905 : // types, the GDAL Metadata mechanism works only with
1906 : // string values. Prior to about 2003-05-02, this driver
1907 : // would attempt to guess the value type from the
1908 : // metadata value string amd then would use the
1909 : // appropriate type-specific FITS keyword update
1910 : // routine. This was found to be troublesome (e.g. a
1911 : // numeric version string with leading zeros would be
1912 : // interpreted as a number and might get those leading
1913 : // zeros stripped), and so now the driver writes every
1914 : // value as a string. In practice this is not a problem
1915 : // since most FITS reading routines will convert from
1916 : // strings to numbers automatically, but if you want
1917 : // finer control, use the underlying FITS handle. Note:
1918 : // to avoid a compiler warning we copy the const value
1919 : // string to a non const one.
1920 2 : char *valueCpy = CPLStrdup(value);
1921 2 : fits_update_key_longstr(m_hFITS, key, valueCpy, nullptr,
1922 : &status);
1923 2 : CPLFree(valueCpy);
1924 :
1925 : // Check for errors.
1926 2 : if (status)
1927 : {
1928 : // Throw a warning with CFITSIO error status, then
1929 : // ignore status
1930 0 : CPLError(
1931 : CE_Warning, CPLE_AppDefined,
1932 : "Couldn't update key %s in FITS file %s (%d).",
1933 : key, GetDescription(), status);
1934 0 : status = 0;
1935 0 : return;
1936 : }
1937 : }
1938 : // Must free up key
1939 24 : CPLFree(key);
1940 : }
1941 : }
1942 :
1943 : // Writing nodata value
1944 41 : if (m_gdalDataType != GDT_Float32 && m_gdalDataType != GDT_Float64)
1945 : {
1946 33 : fits_update_key(m_hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue,
1947 : nullptr, &status);
1948 33 : if (status)
1949 : {
1950 : // Throw a warning with CFITSIO error status, then ignore
1951 : // status
1952 0 : CPLError(CE_Warning, CPLE_AppDefined,
1953 : "Couldn't update key BLANK in FITS file %s (%d).",
1954 : GetDescription(), status);
1955 0 : status = 0;
1956 0 : return;
1957 : }
1958 : }
1959 :
1960 : // Writing Scale and offset if defined
1961 : int pbSuccess;
1962 41 : GDALRasterBand *poSrcBand = GDALPamDataset::GetRasterBand(1);
1963 41 : double dfScale = poSrcBand->GetScale(&pbSuccess);
1964 41 : double dfOffset = poSrcBand->GetOffset(&pbSuccess);
1965 41 : if (m_bMetadataChanged)
1966 : {
1967 1 : fits_update_key(m_hFITS, TDOUBLE, "BSCALE", &dfScale, nullptr,
1968 : &status);
1969 1 : if (status)
1970 : {
1971 : // Throw a warning with CFITSIO error status, then ignore
1972 : // status
1973 0 : CPLError(CE_Warning, CPLE_AppDefined,
1974 : "Couldn't update key BSCALE in FITS file %s (%d).",
1975 : GetDescription(), status);
1976 0 : status = 0;
1977 0 : return;
1978 : }
1979 1 : fits_update_key(m_hFITS, TDOUBLE, "BZERO", &dfOffset, nullptr,
1980 : &status);
1981 1 : if (status)
1982 : {
1983 : // Throw a warning with CFITSIO error status, then ignore
1984 : // status
1985 0 : CPLError(CE_Warning, CPLE_AppDefined,
1986 : "Couldn't update key BZERO in FITS file %s (%d).",
1987 : GetDescription(), status);
1988 0 : status = 0;
1989 0 : return;
1990 : }
1991 : }
1992 :
1993 : // Copy georeferencing info to PAM if the profile is not FITS
1994 41 : GDALPamDataset::SetSpatialRef(GDALPamDataset::GetSpatialRef());
1995 :
1996 : // Write geographic info
1997 41 : if (m_bFITSInfoChanged)
1998 : {
1999 40 : WriteFITSInfo();
2000 : }
2001 :
2002 : // Make sure we flush the raster cache before we close the file!
2003 41 : FlushCache(true);
2004 : }
2005 :
2006 : // Close the FITS handle
2007 103 : fits_close_file(m_hFITS, &status);
2008 103 : if (status != 0)
2009 : {
2010 0 : CPLError(CE_Failure, CPLE_AppDefined,
2011 : "fits_close_file() failed with %d", status);
2012 : }
2013 : }
2014 206 : }
2015 :
2016 : /************************************************************************/
2017 : /* GetRawBinaryLayout() */
2018 : /************************************************************************/
2019 :
2020 2 : bool FITSDataset::GetRawBinaryLayout(GDALDataset::RawBinaryLayout &sLayout)
2021 : {
2022 2 : if (m_hduNum == 0)
2023 0 : return false;
2024 2 : int status = 0;
2025 2 : if (fits_is_compressed_image(m_hFITS, &status))
2026 0 : return false;
2027 2 : GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
2028 2 : if (eDT == GDT_UInt16 || eDT == GDT_UInt32)
2029 0 : return false; // are supported as native signed with offset
2030 :
2031 2 : sLayout.osRawFilename = GetDescription();
2032 : static_assert(sizeof(OFF_T) == 8, "OFF_T should be 64 bits !");
2033 2 : OFF_T headerstart = 0;
2034 2 : OFF_T datastart = 0;
2035 2 : OFF_T dataend = 0;
2036 2 : fits_get_hduoff(m_hFITS, &headerstart, &datastart, &dataend, &status);
2037 2 : if (nBands > 1)
2038 0 : sLayout.eInterleaving = RawBinaryLayout::Interleaving::BSQ;
2039 2 : sLayout.eDataType = eDT;
2040 2 : sLayout.bLittleEndianOrder = false;
2041 2 : sLayout.nImageOffset = static_cast<GIntBig>(datastart);
2042 2 : sLayout.nPixelOffset = GDALGetDataTypeSizeBytes(eDT);
2043 2 : sLayout.nLineOffset = sLayout.nPixelOffset * nRasterXSize;
2044 2 : sLayout.nBandOffset = sLayout.nLineOffset * nRasterYSize;
2045 2 : return true;
2046 : }
2047 :
2048 : /************************************************************************/
2049 : /* Init() */
2050 : /************************************************************************/
2051 :
2052 80 : CPLErr FITSDataset::Init(fitsfile *hFITS, bool isExistingFile, int hduNum)
2053 : {
2054 :
2055 80 : m_hFITS = hFITS;
2056 80 : m_isExistingFile = isExistingFile;
2057 :
2058 80 : int status = 0;
2059 : double offset;
2060 :
2061 80 : int hduType = 0;
2062 80 : fits_movabs_hdu(hFITS, hduNum, &hduType, &status);
2063 80 : if (status)
2064 : {
2065 1 : CPLError(CE_Failure, CPLE_AppDefined,
2066 : "Couldn't move to HDU %d in FITS file %s (%d).", hduNum,
2067 1 : GetDescription(), status);
2068 1 : return CE_Failure;
2069 : }
2070 :
2071 79 : if (hduType != IMAGE_HDU)
2072 : {
2073 0 : CPLError(CE_Failure, CPLE_AppDefined, "HDU %d is not an image.",
2074 : hduNum);
2075 0 : return CE_Failure;
2076 : }
2077 :
2078 : // Get the image info for this dataset (note that all bands in a FITS
2079 : // dataset have the same type)
2080 79 : int bitpix = 0;
2081 79 : int naxis = 0;
2082 79 : const int maxdim = 3;
2083 79 : long naxes[maxdim] = {0, 0, 0};
2084 79 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes, &status);
2085 79 : if (status)
2086 : {
2087 0 : CPLError(CE_Failure, CPLE_AppDefined,
2088 : "Couldn't determine image parameters of FITS file %s (%d)",
2089 0 : GetDescription(), status);
2090 0 : return CE_Failure;
2091 : }
2092 :
2093 79 : m_hduNum = hduNum;
2094 :
2095 79 : fits_read_key(hFITS, TDOUBLE, "BZERO", &offset, nullptr, &status);
2096 79 : if (status)
2097 : {
2098 : // BZERO is not mandatory offset defaulted to 0 if BZERO is missing
2099 63 : status = 0;
2100 63 : offset = 0.;
2101 : }
2102 :
2103 79 : fits_read_key(hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue, nullptr, &status);
2104 79 : m_bNoDataSet = !status;
2105 79 : status = 0;
2106 :
2107 : // Determine data type and nodata value if BLANK keyword is absent
2108 79 : if (bitpix == BYTE_IMG)
2109 : {
2110 33 : m_gdalDataType = GDT_Byte;
2111 33 : m_fitsDataType = TBYTE;
2112 : }
2113 46 : else if (bitpix == SHORT_IMG)
2114 : {
2115 18 : if (offset == 32768.)
2116 : {
2117 7 : m_gdalDataType = GDT_UInt16;
2118 7 : m_fitsDataType = TUSHORT;
2119 : }
2120 : else
2121 : {
2122 11 : m_gdalDataType = GDT_Int16;
2123 11 : m_fitsDataType = TSHORT;
2124 : }
2125 : }
2126 28 : else if (bitpix == LONG_IMG)
2127 : {
2128 14 : if (offset == 2147483648.)
2129 : {
2130 7 : m_gdalDataType = GDT_UInt32;
2131 7 : m_fitsDataType = TUINT;
2132 : }
2133 : else
2134 : {
2135 7 : m_gdalDataType = GDT_Int32;
2136 7 : m_fitsDataType = TINT;
2137 : }
2138 : }
2139 14 : else if (bitpix == FLOAT_IMG)
2140 : {
2141 7 : m_gdalDataType = GDT_Float32;
2142 7 : m_fitsDataType = TFLOAT;
2143 : }
2144 7 : else if (bitpix == DOUBLE_IMG)
2145 : {
2146 7 : m_gdalDataType = GDT_Float64;
2147 7 : m_fitsDataType = TDOUBLE;
2148 : }
2149 : else
2150 : {
2151 0 : CPLError(CE_Failure, CPLE_AppDefined,
2152 0 : "FITS file %s has unknown data type: %d.", GetDescription(),
2153 : bitpix);
2154 0 : return CE_Failure;
2155 : }
2156 :
2157 : // Determine image dimensions - we assume BSQ ordering
2158 79 : if (naxis == 2)
2159 : {
2160 55 : nRasterXSize = static_cast<int>(naxes[0]);
2161 55 : nRasterYSize = static_cast<int>(naxes[1]);
2162 55 : nBands = 1;
2163 : }
2164 24 : else if (naxis == 3)
2165 : {
2166 24 : nRasterXSize = static_cast<int>(naxes[0]);
2167 24 : nRasterYSize = static_cast<int>(naxes[1]);
2168 24 : nBands = static_cast<int>(naxes[2]);
2169 : }
2170 : else
2171 : {
2172 0 : CPLError(CE_Failure, CPLE_AppDefined,
2173 : "FITS file %s does not have 2 or 3 dimensions.",
2174 0 : GetDescription());
2175 0 : return CE_Failure;
2176 : }
2177 :
2178 : // Create the bands
2179 212 : for (int i = 0; i < nBands; ++i)
2180 133 : SetBand(i + 1, new FITSRasterBand(this, i + 1));
2181 :
2182 79 : return CE_None;
2183 : }
2184 :
2185 : /************************************************************************/
2186 : /* LoadMetadata() */
2187 : /************************************************************************/
2188 :
2189 62 : void FITSDataset::LoadMetadata(GDALMajorObject *poTarget)
2190 : {
2191 : // Read header information from file and use it to set metadata
2192 : // This process understands the CONTINUE standard for long strings.
2193 : // We don't bother to capture header names that duplicate information
2194 : // already captured elsewhere (e.g. image dimensions and type)
2195 : int keyNum;
2196 : char key[100];
2197 : char value[100];
2198 62 : CPLStringList aosMD;
2199 :
2200 62 : int nKeys = 0;
2201 62 : int nMoreKeys = 0;
2202 62 : int status = 0;
2203 62 : fits_get_hdrspace(m_hFITS, &nKeys, &nMoreKeys, &status);
2204 2329 : for (keyNum = 1; keyNum <= nKeys; keyNum++)
2205 : {
2206 2267 : fits_read_keyn(m_hFITS, keyNum, key, value, nullptr, &status);
2207 2267 : if (status)
2208 : {
2209 0 : CPLError(CE_Failure, CPLE_AppDefined,
2210 : "Error while reading key %d from FITS file %s (%d)",
2211 0 : keyNum, GetDescription(), status);
2212 0 : return;
2213 : }
2214 2267 : if (strcmp(key, "END") == 0)
2215 : {
2216 : // We should not get here in principle since the END
2217 : // keyword shouldn't be counted in nKeys, but who knows.
2218 0 : break;
2219 : }
2220 2267 : else if (isIgnorableFITSHeader(key))
2221 : {
2222 : // Ignore it
2223 : }
2224 : else
2225 : { // Going to store something, but check for long strings etc
2226 : // Strip off leading and trailing quote if present
2227 1561 : char *newValue = value;
2228 1561 : if (value[0] == '\'' && value[strlen(value) - 1] == '\'')
2229 : {
2230 996 : newValue = value + 1;
2231 996 : value[strlen(value) - 1] = '\0';
2232 : }
2233 : // Check for long string
2234 1561 : if (strrchr(newValue, '&') == newValue + strlen(newValue) - 1)
2235 : {
2236 : // Value string ends in "&", so use long string conventions
2237 0 : char *longString = nullptr;
2238 0 : fits_read_key_longstr(m_hFITS, key, &longString, nullptr,
2239 : &status);
2240 : // Note that read_key_longstr already strips quotes
2241 0 : if (status)
2242 : {
2243 0 : CPLError(CE_Failure, CPLE_AppDefined,
2244 : "Error while reading long string for key %s from "
2245 : "FITS file %s (%d)",
2246 0 : key, GetDescription(), status);
2247 0 : return;
2248 : }
2249 0 : poTarget->SetMetadataItem(key, longString);
2250 0 : free(longString);
2251 : }
2252 : else
2253 : { // Normal keyword
2254 1561 : poTarget->SetMetadataItem(key, newValue);
2255 : }
2256 : }
2257 : }
2258 : }
2259 :
2260 : /************************************************************************/
2261 : /* GetMetadata() */
2262 : /************************************************************************/
2263 :
2264 70 : char **FITSDataset::GetMetadata(const char *pszDomain)
2265 :
2266 : {
2267 70 : if (pszDomain != nullptr && EQUAL(pszDomain, "SUBDATASETS"))
2268 : {
2269 8 : return m_aosSubdatasets.List();
2270 : }
2271 :
2272 62 : return GDALPamDataset::GetMetadata(pszDomain);
2273 : }
2274 :
2275 : /************************************************************************/
2276 : /* GetLayer() */
2277 : /************************************************************************/
2278 :
2279 15 : OGRLayer *FITSDataset::GetLayer(int idx)
2280 : {
2281 15 : if (idx < 0 || idx >= GetLayerCount())
2282 2 : return nullptr;
2283 13 : return m_apoLayers[idx].get();
2284 : }
2285 :
2286 : /************************************************************************/
2287 : /* ICreateLayer() */
2288 : /************************************************************************/
2289 :
2290 4 : OGRLayer *FITSDataset::ICreateLayer(const char *pszName,
2291 : const OGRGeomFieldDefn *poGeomFieldDefn,
2292 : CSLConstList papszOptions)
2293 : {
2294 4 : if (!TestCapability(ODsCCreateLayer))
2295 0 : return nullptr;
2296 :
2297 4 : const auto eGType = poGeomFieldDefn ? poGeomFieldDefn->GetType() : wkbNone;
2298 4 : if (eGType != wkbNone)
2299 : {
2300 0 : CPLError(CE_Failure, CPLE_NotSupported, "Spatial tables not supported");
2301 0 : return nullptr;
2302 : }
2303 :
2304 4 : int status = 0;
2305 4 : int numHDUs = 0;
2306 4 : fits_get_num_hdus(m_hFITS, &numHDUs, &status);
2307 4 : if (status != 0)
2308 : {
2309 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_get_num_hdus() failed: %d",
2310 : status);
2311 0 : return nullptr;
2312 : }
2313 :
2314 4 : fits_create_tbl(m_hFITS, BINARY_TBL,
2315 : 0, // number of initial rows
2316 : 0, // nfields,
2317 : nullptr, // ttype,
2318 : nullptr, // tform
2319 : nullptr, // tunits
2320 : pszName, &status);
2321 4 : if (status != 0)
2322 : {
2323 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create layer");
2324 0 : return nullptr;
2325 : }
2326 :
2327 : // If calling fits_get_num_hdus() here on a freshly new created file,
2328 : // it reports only one HDU, missing the initial dummy HDU
2329 4 : if (numHDUs == 0)
2330 : {
2331 4 : numHDUs = 2;
2332 : }
2333 : else
2334 : {
2335 0 : numHDUs++;
2336 : }
2337 :
2338 4 : auto poLayer = new FITSLayer(this, numHDUs, pszName);
2339 4 : poLayer->SetCreationOptions(papszOptions);
2340 4 : m_apoLayers.emplace_back(std::unique_ptr<FITSLayer>(poLayer));
2341 4 : return m_apoLayers.back().get();
2342 : }
2343 :
2344 : /************************************************************************/
2345 : /* TestCapability() */
2346 : /************************************************************************/
2347 :
2348 20 : int FITSDataset::TestCapability(const char *pszCap)
2349 : {
2350 20 : if (EQUAL(pszCap, ODsCCreateLayer))
2351 8 : return eAccess == GA_Update;
2352 12 : return false;
2353 : }
2354 :
2355 : /************************************************************************/
2356 : /* Open() */
2357 : /************************************************************************/
2358 :
2359 63 : GDALDataset *FITSDataset::Open(GDALOpenInfo *poOpenInfo)
2360 : {
2361 :
2362 63 : if (!FITSDriverIdentify(poOpenInfo))
2363 0 : return nullptr;
2364 :
2365 126 : CPLString osFilename(poOpenInfo->pszFilename);
2366 63 : int iSelectedHDU = 0;
2367 63 : if (STARTS_WITH(poOpenInfo->pszFilename, "FITS:"))
2368 : {
2369 : const CPLStringList aosTokens(
2370 9 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
2371 9 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
2372 9 : if (aosTokens.size() != 3)
2373 : {
2374 2 : return nullptr;
2375 : }
2376 7 : osFilename = aosTokens[1];
2377 7 : iSelectedHDU = atoi(aosTokens[2]);
2378 7 : if (iSelectedHDU <= 0)
2379 : {
2380 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid HDU number");
2381 1 : return nullptr;
2382 : }
2383 : }
2384 :
2385 : // Get access mode and attempt to open the file
2386 60 : int status = 0;
2387 60 : fitsfile *hFITS = nullptr;
2388 60 : if (poOpenInfo->eAccess == GA_ReadOnly)
2389 58 : fits_open_file(&hFITS, osFilename.c_str(), READONLY, &status);
2390 : else
2391 2 : fits_open_file(&hFITS, osFilename.c_str(), READWRITE, &status);
2392 60 : if (status)
2393 : {
2394 1 : CPLError(CE_Failure, CPLE_AppDefined,
2395 : "Error while opening FITS file %s (%d).\n", osFilename.c_str(),
2396 : status);
2397 1 : fits_close_file(hFITS, &status);
2398 1 : return nullptr;
2399 : }
2400 : // Create a FITSDataset object
2401 118 : auto dataset = std::make_unique<FITSDataset>();
2402 59 : dataset->m_isExistingFile = true;
2403 59 : dataset->m_hFITS = hFITS;
2404 59 : dataset->eAccess = poOpenInfo->eAccess;
2405 59 : dataset->SetPhysicalFilename(osFilename);
2406 :
2407 : /* -------------------------------------------------------------------- */
2408 : /* Iterate over HDUs */
2409 : /* -------------------------------------------------------------------- */
2410 59 : bool firstHDUIsDummy = false;
2411 59 : int firstValidHDU = 0;
2412 118 : CPLStringList aosSubdatasets;
2413 59 : bool hasVector = false;
2414 59 : if (iSelectedHDU == 0)
2415 : {
2416 54 : int numHDUs = 0;
2417 54 : fits_get_num_hdus(hFITS, &numHDUs, &status);
2418 54 : if (numHDUs <= 0)
2419 : {
2420 0 : return nullptr;
2421 : }
2422 :
2423 131 : for (int iHDU = 1; iHDU <= numHDUs; iHDU++)
2424 : {
2425 77 : int hduType = 0;
2426 77 : fits_movabs_hdu(hFITS, iHDU, &hduType, &status);
2427 77 : if (status)
2428 : {
2429 31 : continue;
2430 : }
2431 :
2432 77 : char szExtname[81] = {0};
2433 77 : fits_read_key(hFITS, TSTRING, "EXTNAME", szExtname, nullptr,
2434 : &status);
2435 77 : status = 0;
2436 77 : int nExtVer = 0;
2437 77 : fits_read_key(hFITS, TINT, "EXTVER", &nExtVer, nullptr, &status);
2438 77 : status = 0;
2439 77 : CPLString osExtname(szExtname);
2440 77 : if (nExtVer > 0)
2441 0 : osExtname += CPLSPrintf(" %d", nExtVer);
2442 :
2443 77 : if (hduType == BINARY_TBL)
2444 : {
2445 14 : hasVector = true;
2446 14 : if ((poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0)
2447 : {
2448 24 : dataset->m_apoLayers.push_back(std::unique_ptr<FITSLayer>(
2449 12 : new FITSLayer(dataset.get(), iHDU, osExtname.c_str())));
2450 : }
2451 : }
2452 :
2453 77 : if (hduType != IMAGE_HDU)
2454 : {
2455 14 : continue;
2456 : }
2457 :
2458 63 : int bitpix = 0;
2459 63 : int naxis = 0;
2460 63 : const int maxdim = 3;
2461 63 : long naxes[maxdim] = {0, 0, 0};
2462 63 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes, &status);
2463 63 : if (status)
2464 : {
2465 0 : continue;
2466 : }
2467 :
2468 63 : if (naxis != 2 && naxis != 3)
2469 : {
2470 17 : if (naxis == 0 && iHDU == 1)
2471 : {
2472 17 : firstHDUIsDummy = true;
2473 : }
2474 17 : continue;
2475 : }
2476 :
2477 46 : if ((poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0)
2478 : {
2479 41 : const int nIdx = aosSubdatasets.size() / 2 + 1;
2480 : aosSubdatasets.AddNameValue(
2481 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
2482 : CPLSPrintf("FITS:\"%s\":%d", poOpenInfo->pszFilename,
2483 41 : iHDU));
2484 : CPLString osDesc(CPLSPrintf(
2485 : "HDU %d (%dx%d, %d band%s)", iHDU,
2486 41 : static_cast<int>(naxes[0]), static_cast<int>(naxes[1]),
2487 51 : naxis == 3 ? static_cast<int>(naxes[2]) : 1,
2488 82 : (naxis == 3 && naxes[2] > 1) ? "s" : ""));
2489 41 : if (!osExtname.empty())
2490 : {
2491 5 : osDesc += ", ";
2492 5 : osDesc += osExtname;
2493 : }
2494 : aosSubdatasets.AddNameValue(
2495 41 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx), osDesc);
2496 : }
2497 :
2498 46 : if (firstValidHDU == 0)
2499 : {
2500 41 : firstValidHDU = iHDU;
2501 : }
2502 : }
2503 54 : if (aosSubdatasets.size() == 2)
2504 : {
2505 35 : aosSubdatasets.Clear();
2506 : }
2507 : }
2508 : else
2509 : {
2510 5 : if (iSelectedHDU != 1)
2511 : {
2512 4 : int hduType = 0;
2513 4 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2514 4 : if (status == 0)
2515 : {
2516 4 : int bitpix = 0;
2517 4 : int naxis = 0;
2518 4 : const int maxdim = 3;
2519 4 : long naxes[maxdim] = {0, 0, 0};
2520 4 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes,
2521 : &status);
2522 4 : if (status == 0 && naxis == 0)
2523 : {
2524 2 : firstHDUIsDummy = true;
2525 : }
2526 : }
2527 4 : status = 0;
2528 : }
2529 5 : firstValidHDU = iSelectedHDU;
2530 : }
2531 :
2532 59 : const bool hasRaster = firstValidHDU > 0;
2533 59 : const bool hasRasterAndIsAllowed =
2534 59 : hasRaster && (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0;
2535 :
2536 59 : if (!hasRasterAndIsAllowed &&
2537 16 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0 &&
2538 3 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) == 0)
2539 : {
2540 2 : if (hasVector)
2541 : {
2542 2 : std::string osPath;
2543 1 : osPath.resize(1024);
2544 1 : if (CPLGetExecPath(&osPath[0], static_cast<int>(osPath.size())))
2545 : {
2546 1 : osPath = CPLGetBasename(osPath.c_str());
2547 : }
2548 1 : if (osPath == "gdalinfo")
2549 : {
2550 0 : CPLError(CE_Failure, CPLE_AppDefined,
2551 : "This FITS dataset does not contain any image, but "
2552 : "contains binary table(s) that could be opened "
2553 : "in vector mode with ogrinfo.");
2554 : }
2555 : else
2556 : {
2557 1 : CPLError(CE_Failure, CPLE_AppDefined,
2558 : "This FITS dataset does not contain any image, but "
2559 : "contains binary table(s) that could be opened "
2560 : "in vector mode.");
2561 : }
2562 : }
2563 : else
2564 : {
2565 1 : CPLError(CE_Failure, CPLE_AppDefined,
2566 : "Cannot find HDU of image type with 2 or 3 axes.");
2567 : }
2568 2 : return nullptr;
2569 : }
2570 :
2571 57 : if (dataset->m_apoLayers.empty() &&
2572 60 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0 &&
2573 3 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0)
2574 : {
2575 3 : if (hasRaster)
2576 : {
2577 4 : std::string osPath;
2578 2 : osPath.resize(1024);
2579 2 : if (CPLGetExecPath(&osPath[0], static_cast<int>(osPath.size())))
2580 : {
2581 2 : osPath = CPLGetBasename(osPath.c_str());
2582 : }
2583 2 : if (osPath == "ogrinfo")
2584 : {
2585 0 : CPLError(CE_Failure, CPLE_AppDefined,
2586 : "This FITS dataset does not contain any binary "
2587 : "table, but contains image(s) that could be opened "
2588 : "in raster mode with gdalinfo.");
2589 : }
2590 : else
2591 : {
2592 2 : CPLError(CE_Failure, CPLE_AppDefined,
2593 : "This FITS dataset does not contain any binary "
2594 : "table, but contains image(s) that could be opened "
2595 : "in raster mode.");
2596 : }
2597 : }
2598 : else
2599 : {
2600 1 : CPLError(CE_Failure, CPLE_AppDefined,
2601 : "Cannot find binary table(s).");
2602 : }
2603 3 : return nullptr;
2604 : }
2605 :
2606 54 : dataset->m_aosSubdatasets = aosSubdatasets;
2607 :
2608 : // Set up the description and initialize the dataset
2609 54 : dataset->SetDescription(poOpenInfo->pszFilename);
2610 54 : if (hasRasterAndIsAllowed)
2611 : {
2612 43 : if (aosSubdatasets.size() > 2)
2613 : {
2614 3 : firstValidHDU = 0;
2615 3 : int hduType = 0;
2616 3 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2617 : }
2618 : else
2619 : {
2620 80 : if (firstValidHDU != 0 &&
2621 40 : dataset->Init(hFITS, true, firstValidHDU) != CE_None)
2622 : {
2623 1 : return nullptr;
2624 : }
2625 : }
2626 : }
2627 :
2628 : // If the first HDU is a dummy one, load its metadata first, and then
2629 : // add/override it by the one of the image HDU
2630 53 : if (firstHDUIsDummy && firstValidHDU > 1)
2631 : {
2632 4 : int hduType = 0;
2633 4 : status = 0;
2634 4 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2635 4 : if (status == 0)
2636 : {
2637 4 : dataset->LoadMetadata(dataset.get());
2638 : }
2639 4 : status = 0;
2640 4 : fits_movabs_hdu(hFITS, firstValidHDU, &hduType, &status);
2641 4 : if (status)
2642 : {
2643 0 : return nullptr;
2644 : }
2645 : }
2646 53 : if (hasRasterAndIsAllowed)
2647 : {
2648 42 : dataset->LoadMetadata(dataset.get());
2649 42 : dataset->LoadFITSInfo();
2650 : }
2651 :
2652 : /* -------------------------------------------------------------------- */
2653 : /* Initialize any information. */
2654 : /* -------------------------------------------------------------------- */
2655 53 : dataset->SetDescription(poOpenInfo->pszFilename);
2656 53 : dataset->TryLoadXML();
2657 :
2658 : /* -------------------------------------------------------------------- */
2659 : /* Check for external overviews. */
2660 : /* -------------------------------------------------------------------- */
2661 106 : dataset->oOvManager.Initialize(dataset.get(), poOpenInfo->pszFilename,
2662 53 : poOpenInfo->GetSiblingFiles());
2663 :
2664 53 : return dataset.release();
2665 : }
2666 :
2667 : /************************************************************************/
2668 : /* Create() */
2669 : /* */
2670 : /* Create a new FITS file. */
2671 : /************************************************************************/
2672 :
2673 67 : GDALDataset *FITSDataset::Create(const char *pszFilename, int nXSize,
2674 : int nYSize, int nBandsIn, GDALDataType eType,
2675 : CPL_UNUSED char **papszParamList)
2676 : {
2677 67 : int status = 0;
2678 :
2679 67 : if (nXSize == 0 && nYSize == 0 && nBandsIn == 0 && eType == GDT_Unknown)
2680 : {
2681 : // Create the file - to force creation, we prepend the name with '!'
2682 8 : CPLString extFilename("!");
2683 4 : extFilename += pszFilename;
2684 4 : fitsfile *hFITS = nullptr;
2685 4 : fits_create_file(&hFITS, extFilename, &status);
2686 4 : if (status)
2687 : {
2688 0 : CPLError(CE_Failure, CPLE_AppDefined,
2689 : "Couldn't create FITS file %s (%d).\n", pszFilename,
2690 : status);
2691 0 : return nullptr;
2692 : }
2693 :
2694 : // Likely vector creation mode
2695 4 : FITSDataset *dataset = new FITSDataset();
2696 4 : dataset->m_hFITS = hFITS;
2697 4 : dataset->eAccess = GA_Update;
2698 4 : dataset->SetDescription(pszFilename);
2699 4 : return dataset;
2700 : }
2701 :
2702 : // No creation options are defined. The BSCALE/BZERO options were
2703 : // removed on 2002-07-02 by Simon Perkins because they introduced
2704 : // excessive complications and didn't really fit into the GDAL
2705 : // paradigm.
2706 : // 2018 - BZERO BSCALE keywords are now set using SetScale() and
2707 : // SetOffset() functions
2708 :
2709 63 : if (nXSize < 1 || nYSize < 1 || nBandsIn < 1)
2710 : {
2711 2 : CPLError(CE_Failure, CPLE_AppDefined,
2712 : "Attempt to create %dx%dx%d raster FITS file, but width, "
2713 : "height and bands"
2714 : " must be positive.",
2715 : nXSize, nYSize, nBandsIn);
2716 :
2717 2 : return nullptr;
2718 : }
2719 :
2720 : // Determine FITS type of image
2721 : int bitpix;
2722 61 : if (eType == GDT_Byte)
2723 : {
2724 17 : bitpix = BYTE_IMG;
2725 : }
2726 44 : else if (eType == GDT_UInt16)
2727 : {
2728 4 : bitpix = USHORT_IMG;
2729 : }
2730 40 : else if (eType == GDT_Int16)
2731 : {
2732 6 : bitpix = SHORT_IMG;
2733 : }
2734 34 : else if (eType == GDT_UInt32)
2735 : {
2736 4 : bitpix = ULONG_IMG;
2737 : }
2738 30 : else if (eType == GDT_Int32)
2739 : {
2740 4 : bitpix = LONG_IMG;
2741 : }
2742 26 : else if (eType == GDT_Float32)
2743 4 : bitpix = FLOAT_IMG;
2744 22 : else if (eType == GDT_Float64)
2745 4 : bitpix = DOUBLE_IMG;
2746 : else
2747 : {
2748 18 : CPLError(CE_Failure, CPLE_AppDefined,
2749 : "GDALDataType (%d) unsupported for FITS", eType);
2750 18 : return nullptr;
2751 : }
2752 :
2753 : // Create the file - to force creation, we prepend the name with '!'
2754 86 : CPLString extFilename("!");
2755 43 : extFilename += pszFilename;
2756 43 : fitsfile *hFITS = nullptr;
2757 43 : fits_create_file(&hFITS, extFilename, &status);
2758 43 : if (status)
2759 : {
2760 3 : CPLError(CE_Failure, CPLE_AppDefined,
2761 : "Couldn't create FITS file %s (%d).\n", pszFilename, status);
2762 3 : return nullptr;
2763 : }
2764 :
2765 : // Now create an image of appropriate size and type
2766 40 : long naxes[3] = {nXSize, nYSize, nBandsIn};
2767 40 : int naxis = (nBandsIn == 1) ? 2 : 3;
2768 40 : fits_create_img(hFITS, bitpix, naxis, naxes, &status);
2769 :
2770 : // Check the status
2771 40 : if (status)
2772 : {
2773 0 : CPLError(CE_Failure, CPLE_AppDefined,
2774 : "Couldn't create image within FITS file %s (%d).", pszFilename,
2775 : status);
2776 0 : fits_close_file(hFITS, &status);
2777 0 : return nullptr;
2778 : }
2779 :
2780 40 : FITSDataset *dataset = new FITSDataset();
2781 40 : dataset->nRasterXSize = nXSize;
2782 40 : dataset->nRasterYSize = nYSize;
2783 40 : dataset->eAccess = GA_Update;
2784 40 : dataset->SetDescription(pszFilename);
2785 :
2786 : // Init recalculates a lot of stuff we already know, but...
2787 40 : if (dataset->Init(hFITS, false, 1) != CE_None)
2788 : {
2789 0 : delete dataset;
2790 0 : return nullptr;
2791 : }
2792 : else
2793 : {
2794 40 : return dataset;
2795 : }
2796 : }
2797 :
2798 : /************************************************************************/
2799 : /* Delete() */
2800 : /************************************************************************/
2801 :
2802 13 : CPLErr FITSDataset::Delete(const char *pszFilename)
2803 : {
2804 13 : return VSIUnlink(pszFilename) == 0 ? CE_None : CE_Failure;
2805 : }
2806 :
2807 : /************************************************************************/
2808 : /* WriteFITSInfo() */
2809 : /************************************************************************/
2810 :
2811 40 : void FITSDataset::WriteFITSInfo()
2812 :
2813 : {
2814 40 : int status = 0;
2815 :
2816 40 : const double PI = std::atan(1.0) * 4;
2817 40 : const double DEG2RAD = PI / 180.;
2818 :
2819 40 : double falseEast = 0;
2820 40 : double falseNorth = 0;
2821 :
2822 : double cfactor, mres, mapres, UpperLeftCornerX, UpperLeftCornerY;
2823 : double crpix1, crpix2;
2824 :
2825 : /* -------------------------------------------------------------------- */
2826 : /* Write out projection definition. */
2827 : /* -------------------------------------------------------------------- */
2828 40 : const bool bHasProjection = !m_oSRS.IsEmpty();
2829 40 : if (bHasProjection)
2830 : {
2831 :
2832 : // Set according to coordinate system (thanks to Trent Hare - USGS)
2833 :
2834 40 : std::string object, ctype1, ctype2;
2835 :
2836 40 : const char *target = m_oSRS.GetAttrValue("DATUM", 0);
2837 40 : if (target)
2838 : {
2839 40 : if (strstr(target, "Moon"))
2840 : {
2841 0 : object.assign("Moon");
2842 0 : ctype1.assign("SE");
2843 0 : ctype2.assign("SE");
2844 : }
2845 40 : else if (strstr(target, "Mercury"))
2846 : {
2847 0 : object.assign("Mercury");
2848 0 : ctype1.assign("ME");
2849 0 : ctype2.assign("ME");
2850 : }
2851 40 : else if (strstr(target, "Venus"))
2852 : {
2853 0 : object.assign("Venus");
2854 0 : ctype1.assign("VE");
2855 0 : ctype2.assign("VE");
2856 : }
2857 40 : else if (strstr(target, "Mars"))
2858 : {
2859 0 : object.assign("Mars");
2860 0 : ctype1.assign("MA");
2861 0 : ctype2.assign("MA");
2862 : }
2863 40 : else if (strstr(target, "Jupiter"))
2864 : {
2865 0 : object.assign("Jupiter");
2866 0 : ctype1.assign("JU");
2867 0 : ctype2.assign("JU");
2868 : }
2869 40 : else if (strstr(target, "Saturn"))
2870 : {
2871 0 : object.assign("Saturn");
2872 0 : ctype1.assign("SA");
2873 0 : ctype2.assign("SA");
2874 : }
2875 40 : else if (strstr(target, "Uranus"))
2876 : {
2877 0 : object.assign("Uranus");
2878 0 : ctype1.assign("UR");
2879 0 : ctype2.assign("UR");
2880 : }
2881 40 : else if (strstr(target, "Neptune"))
2882 : {
2883 0 : object.assign("Neptune");
2884 0 : ctype1.assign("NE");
2885 0 : ctype2.assign("NE");
2886 : }
2887 : else
2888 : {
2889 40 : object.assign("Earth");
2890 40 : ctype1.assign("EA");
2891 40 : ctype2.assign("EA");
2892 : }
2893 :
2894 40 : fits_update_key(
2895 : m_hFITS, TSTRING, "OBJECT",
2896 40 : const_cast<void *>(static_cast<const void *>(object.c_str())),
2897 : nullptr, &status);
2898 : }
2899 :
2900 40 : double aradius = m_oSRS.GetSemiMajor();
2901 40 : double bradius = aradius;
2902 40 : double cradius = m_oSRS.GetSemiMinor();
2903 :
2904 40 : cfactor = aradius * DEG2RAD;
2905 :
2906 40 : fits_update_key(m_hFITS, TDOUBLE, "A_RADIUS", &aradius, nullptr,
2907 : &status);
2908 40 : if (status)
2909 : {
2910 : // Throw a warning with CFITSIO error status, then ignore status
2911 0 : CPLError(CE_Warning, CPLE_AppDefined,
2912 : "Couldn't update key A_RADIUS in FITS file %s (%d).",
2913 0 : GetDescription(), status);
2914 0 : status = 0;
2915 0 : return;
2916 : }
2917 40 : fits_update_key(m_hFITS, TDOUBLE, "B_RADIUS", &bradius, nullptr,
2918 : &status);
2919 40 : if (status)
2920 : {
2921 : // Throw a warning with CFITSIO error status, then ignore status
2922 0 : CPLError(CE_Warning, CPLE_AppDefined,
2923 : "Couldn't update key B_RADIUS in FITS file %s (%d).",
2924 0 : GetDescription(), status);
2925 0 : status = 0;
2926 0 : return;
2927 : }
2928 40 : fits_update_key(m_hFITS, TDOUBLE, "C_RADIUS", &cradius, nullptr,
2929 : &status);
2930 40 : if (status)
2931 : {
2932 : // Throw a warning with CFITSIO error status, then ignore status
2933 0 : CPLError(CE_Warning, CPLE_AppDefined,
2934 : "Couldn't update key C_RADIUS in FITS file %s (%d).",
2935 0 : GetDescription(), status);
2936 0 : status = 0;
2937 0 : return;
2938 : }
2939 :
2940 40 : const char *unit = m_oSRS.GetAttrValue("UNIT", 0);
2941 :
2942 40 : ctype1.append("LN-");
2943 40 : ctype2.append("LT-");
2944 :
2945 : // strcat(ctype1a, "PX-");
2946 : // strcat(ctype2a, "PY-");
2947 :
2948 40 : std::string fitsproj;
2949 40 : const char *projection = m_oSRS.GetAttrValue("PROJECTION", 0);
2950 40 : double centlon = 0, centlat = 0;
2951 :
2952 40 : if (projection)
2953 : {
2954 12 : if (strstr(projection, "Sinusoidal"))
2955 : {
2956 0 : fitsproj.assign("SFL");
2957 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2958 : }
2959 12 : else if (strstr(projection, "Equirectangular"))
2960 : {
2961 0 : fitsproj.assign("CAR");
2962 0 : centlat = m_oSRS.GetProjParm("standard_parallel_1", 0, nullptr);
2963 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2964 : }
2965 12 : else if (strstr(projection, "Orthographic"))
2966 : {
2967 0 : fitsproj.assign("SIN");
2968 0 : centlat = m_oSRS.GetProjParm("standard_parallel_1", 0, nullptr);
2969 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2970 : }
2971 12 : else if (strstr(projection, "Mercator_1SP") ||
2972 12 : strstr(projection, "Mercator"))
2973 : {
2974 12 : fitsproj.assign("MER");
2975 12 : centlat = m_oSRS.GetProjParm("standard_parallel_1", 0, nullptr);
2976 12 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2977 : }
2978 0 : else if (strstr(projection, "Polar_Stereographic") ||
2979 0 : strstr(projection, "Stereographic_South_Pole") ||
2980 0 : strstr(projection, "Stereographic_North_Pole"))
2981 : {
2982 0 : fitsproj.assign("STG");
2983 0 : centlat = m_oSRS.GetProjParm("latitude_of_origin", 0, nullptr);
2984 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2985 : }
2986 :
2987 : /*
2988 : #Transverse Mercator is supported in FITS via
2989 : specific MER parameters. # need some more testing... #if
2990 : EQUAL(mapProjection,"Transverse_Mercator"): # mapProjection =
2991 : "MER" # centLat = hSRS.GetProjParm('standard_parallel_1') #
2992 : centLon = hSRS.GetProjParm('central_meridian') # TMscale =
2993 : hSRS.GetProjParm('scale_factor') # #Need to research when TM
2994 : actually applies false values # #but planetary is almost
2995 : always 0.0 # falseEast = hSRS.GetProjParm('false_easting') #
2996 : falseNorth = hSRS.GetProjParm('false_northing')
2997 : */
2998 :
2999 12 : ctype1.append(fitsproj);
3000 12 : ctype2.append(fitsproj);
3001 :
3002 12 : fits_update_key(
3003 : m_hFITS, TSTRING, "CTYPE1",
3004 12 : const_cast<void *>(static_cast<const void *>(ctype1.c_str())),
3005 : nullptr, &status);
3006 12 : if (status)
3007 : {
3008 : // Throw a warning with CFITSIO error status, then ignore status
3009 0 : CPLError(CE_Warning, CPLE_AppDefined,
3010 : "Couldn't update key CTYPE1 in FITS file %s (%d).",
3011 0 : GetDescription(), status);
3012 0 : status = 0;
3013 0 : return;
3014 : }
3015 :
3016 12 : fits_update_key(
3017 : m_hFITS, TSTRING, "CTYPE2",
3018 12 : const_cast<void *>(static_cast<const void *>(ctype2.c_str())),
3019 : nullptr, &status);
3020 12 : if (status)
3021 : {
3022 : // Throw a warning with CFITSIO error status, then ignore status
3023 0 : CPLError(CE_Warning, CPLE_AppDefined,
3024 : "Couldn't update key CTYPE2 in FITS file %s (%d).",
3025 0 : GetDescription(), status);
3026 0 : status = 0;
3027 0 : return;
3028 : }
3029 : }
3030 :
3031 40 : UpperLeftCornerX = m_adfGeoTransform[0] - falseEast;
3032 40 : UpperLeftCornerY = m_adfGeoTransform[3] - falseNorth;
3033 :
3034 40 : if (centlon > 180.)
3035 : {
3036 0 : centlon = centlon - 180.;
3037 : }
3038 40 : if (strstr(unit, "metre"))
3039 : {
3040 : // convert degrees/pixel to m/pixel
3041 12 : mapres = 1. / m_adfGeoTransform[1]; // mapres is pixel/meters
3042 12 : mres = m_adfGeoTransform[1] / cfactor; // mres is deg/pixel
3043 12 : crpix1 = -(UpperLeftCornerX * mapres) + centlon / mres + 0.5;
3044 : // assuming that center latitude is also the origin of the
3045 : // coordinate system: this is not always true. More generic
3046 : // implementation coming soon
3047 12 : crpix2 = (UpperLeftCornerY * mapres) + 0.5; // - (centlat / mres);
3048 : }
3049 28 : else if (strstr(unit, "degree"))
3050 : {
3051 : // convert m/pixel to pixel/degree
3052 28 : mapres =
3053 28 : 1. / m_adfGeoTransform[1] / cfactor; // mapres is pixel/deg
3054 28 : mres = m_adfGeoTransform[1]; // mres is meters/pixel
3055 28 : crpix1 = -(UpperLeftCornerX * mres) + centlon / mapres + 0.5;
3056 : // assuming that center latitude is also the origin of the
3057 : // coordinate system: this is not always true. More generic
3058 : // implementation coming soon
3059 28 : crpix2 = (UpperLeftCornerY * mres) + 0.5; // - (centlat / mapres);
3060 : }
3061 :
3062 : /// Write WCS CRPIXia CRVALia CTYPEia here
3063 :
3064 40 : fits_update_key(m_hFITS, TDOUBLE, "CRVAL1", ¢lon, nullptr, &status);
3065 40 : if (status)
3066 : {
3067 : // Throw a warning with CFITSIO error status, then ignore status
3068 0 : CPLError(CE_Warning, CPLE_AppDefined,
3069 : "Couldn't update key CRVAL1 in FITS file %s (%d).",
3070 0 : GetDescription(), status);
3071 0 : status = 0;
3072 0 : return;
3073 : }
3074 40 : fits_update_key(m_hFITS, TDOUBLE, "CRVAL2", ¢lat, nullptr, &status);
3075 40 : if (status)
3076 : {
3077 : // Throw a warning with CFITSIO error status, then ignore status
3078 0 : CPLError(CE_Warning, CPLE_AppDefined,
3079 : "Couldn't update key CRVAL2 in FITS file %s (%d).",
3080 0 : GetDescription(), status);
3081 0 : status = 0;
3082 0 : return;
3083 : }
3084 40 : fits_update_key(m_hFITS, TDOUBLE, "CRPIX1", &crpix1, nullptr, &status);
3085 40 : if (status)
3086 : {
3087 : // Throw a warning with CFITSIO error status, then ignore status
3088 0 : CPLError(CE_Warning, CPLE_AppDefined,
3089 : "Couldn't update key CRPIX1 in FITS file %s (%d).",
3090 0 : GetDescription(), status);
3091 0 : status = 0;
3092 0 : return;
3093 : }
3094 40 : fits_update_key(m_hFITS, TDOUBLE, "CRPIX2", &crpix2, nullptr, &status);
3095 40 : if (status)
3096 : {
3097 : // Throw a warning with CFITSIO error status, then ignore status
3098 0 : CPLError(CE_Warning, CPLE_AppDefined,
3099 : "Couldn't update key CRPIX2 in FITS file %s (%d).",
3100 0 : GetDescription(), status);
3101 0 : status = 0;
3102 0 : return;
3103 : }
3104 :
3105 : /* --------------------------------------------------------------------
3106 : */
3107 : /* Write geotransform if valid. */
3108 : /* --------------------------------------------------------------------
3109 : */
3110 40 : if (m_bGeoTransformValid)
3111 : {
3112 :
3113 : /* --------------------------------------------------------------------
3114 : */
3115 : /* Write the transform. */
3116 : /* --------------------------------------------------------------------
3117 : */
3118 :
3119 : /// Write WCS CDELTia and PCi_ja here
3120 :
3121 : double cd[4];
3122 40 : cd[0] = m_adfGeoTransform[1] / cfactor;
3123 40 : cd[1] = m_adfGeoTransform[2] / cfactor;
3124 40 : cd[2] = m_adfGeoTransform[4] / cfactor;
3125 40 : cd[3] = m_adfGeoTransform[5] / cfactor;
3126 :
3127 : double pc[4];
3128 40 : pc[0] = 1.;
3129 40 : pc[1] = cd[1] / cd[0];
3130 40 : pc[2] = cd[2] / cd[3];
3131 40 : pc[3] = -1.;
3132 :
3133 40 : fits_update_key(m_hFITS, TDOUBLE, "CDELT1", &cd[0], nullptr,
3134 : &status);
3135 40 : if (status)
3136 : {
3137 : // Throw a warning with CFITSIO error status, then ignore status
3138 0 : CPLError(CE_Warning, CPLE_AppDefined,
3139 : "Couldn't update key CDELT1 in FITS file %s (%d).",
3140 0 : GetDescription(), status);
3141 0 : status = 0;
3142 0 : return;
3143 : }
3144 :
3145 40 : fits_update_key(m_hFITS, TDOUBLE, "CDELT2", &cd[3], nullptr,
3146 : &status);
3147 40 : if (status)
3148 : {
3149 : // Throw a warning with CFITSIO error status, then ignore status
3150 0 : CPLError(CE_Warning, CPLE_AppDefined,
3151 : "Couldn't update key CDELT2 in FITS file %s (%d).",
3152 0 : GetDescription(), status);
3153 0 : status = 0;
3154 0 : return;
3155 : }
3156 :
3157 40 : fits_update_key(m_hFITS, TDOUBLE, "PC1_1", &pc[0], nullptr,
3158 : &status);
3159 40 : if (status)
3160 : {
3161 : // Throw a warning with CFITSIO error status, then ignore status
3162 0 : CPLError(CE_Warning, CPLE_AppDefined,
3163 : "Couldn't update key PC1_1 in FITS file %s (%d).",
3164 0 : GetDescription(), status);
3165 0 : status = 0;
3166 0 : return;
3167 : }
3168 :
3169 40 : fits_update_key(m_hFITS, TDOUBLE, "PC1_2", &pc[1], nullptr,
3170 : &status);
3171 40 : if (status)
3172 : {
3173 : // Throw a warning with CFITSIO error status, then ignore status
3174 0 : CPLError(CE_Warning, CPLE_AppDefined,
3175 : "Couldn't update key PC1_2 in FITS file %s (%d).",
3176 0 : GetDescription(), status);
3177 0 : status = 0;
3178 0 : return;
3179 : }
3180 :
3181 40 : fits_update_key(m_hFITS, TDOUBLE, "PC2_1", &pc[2], nullptr,
3182 : &status);
3183 40 : if (status)
3184 : {
3185 : // Throw a warning with CFITSIO error status, then ignore status
3186 0 : CPLError(CE_Warning, CPLE_AppDefined,
3187 : "Couldn't update key PC2_1 in FITS file %s (%d).",
3188 0 : GetDescription(), status);
3189 0 : status = 0;
3190 0 : return;
3191 : }
3192 :
3193 40 : fits_update_key(m_hFITS, TDOUBLE, "PC2_2", &pc[3], nullptr,
3194 : &status);
3195 40 : if (status)
3196 : {
3197 : // Throw a warning with CFITSIO error status, then ignore status
3198 0 : CPLError(CE_Warning, CPLE_AppDefined,
3199 : "Couldn't update key PC2_2 in FITS file %s (%d).",
3200 0 : GetDescription(), status);
3201 0 : status = 0;
3202 0 : return;
3203 : }
3204 : }
3205 : }
3206 : }
3207 :
3208 : /************************************************************************/
3209 : /* GetSpatialRef() */
3210 : /************************************************************************/
3211 :
3212 4 : const OGRSpatialReference *FITSDataset::GetSpatialRef() const
3213 :
3214 : {
3215 4 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
3216 : }
3217 :
3218 : /************************************************************************/
3219 : /* SetSpatialRef() */
3220 : /************************************************************************/
3221 :
3222 40 : CPLErr FITSDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
3223 :
3224 : {
3225 40 : if (poSRS == nullptr || poSRS->IsEmpty())
3226 : {
3227 0 : m_oSRS.Clear();
3228 : }
3229 : else
3230 : {
3231 40 : m_oSRS = *poSRS;
3232 40 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3233 : }
3234 :
3235 40 : m_bFITSInfoChanged = true;
3236 :
3237 40 : return CE_None;
3238 : }
3239 :
3240 : /************************************************************************/
3241 : /* GetGeoTransform() */
3242 : /************************************************************************/
3243 :
3244 19 : CPLErr FITSDataset::GetGeoTransform(double *padfTransform)
3245 :
3246 : {
3247 19 : memcpy(padfTransform, m_adfGeoTransform, sizeof(double) * 6);
3248 :
3249 19 : if (!m_bGeoTransformValid)
3250 0 : return CE_Failure;
3251 :
3252 19 : return CE_None;
3253 : }
3254 :
3255 : /************************************************************************/
3256 : /* SetGeoTransform() */
3257 : /************************************************************************/
3258 :
3259 40 : CPLErr FITSDataset::SetGeoTransform(double *padfTransform)
3260 :
3261 : {
3262 40 : memcpy(m_adfGeoTransform, padfTransform, sizeof(double) * 6);
3263 40 : m_bGeoTransformValid = true;
3264 :
3265 40 : return CE_None;
3266 : }
3267 :
3268 : /************************************************************************/
3269 : /* GetOffset() */
3270 : /************************************************************************/
3271 :
3272 44 : double FITSRasterBand::GetOffset(int *pbSuccess)
3273 :
3274 : {
3275 44 : if (pbSuccess)
3276 43 : *pbSuccess = m_bHaveOffsetScale;
3277 44 : return m_dfOffset;
3278 : }
3279 :
3280 : /************************************************************************/
3281 : /* SetOffset() */
3282 : /************************************************************************/
3283 :
3284 1 : CPLErr FITSRasterBand::SetOffset(double dfNewValue)
3285 :
3286 : {
3287 1 : if (!m_bHaveOffsetScale || dfNewValue != m_dfOffset)
3288 1 : m_poFDS->m_bMetadataChanged = true;
3289 :
3290 1 : m_bHaveOffsetScale = true;
3291 1 : m_dfOffset = dfNewValue;
3292 1 : return CE_None;
3293 : }
3294 :
3295 : /************************************************************************/
3296 : /* GetScale() */
3297 : /************************************************************************/
3298 :
3299 44 : double FITSRasterBand::GetScale(int *pbSuccess)
3300 :
3301 : {
3302 44 : if (pbSuccess)
3303 43 : *pbSuccess = m_bHaveOffsetScale;
3304 44 : return m_dfScale;
3305 : }
3306 :
3307 : /************************************************************************/
3308 : /* SetScale() */
3309 : /************************************************************************/
3310 :
3311 1 : CPLErr FITSRasterBand::SetScale(double dfNewValue)
3312 :
3313 : {
3314 1 : if (!m_bHaveOffsetScale || dfNewValue != m_dfScale)
3315 1 : m_poFDS->m_bMetadataChanged = true;
3316 :
3317 1 : m_bHaveOffsetScale = true;
3318 1 : m_dfScale = dfNewValue;
3319 1 : return CE_None;
3320 : }
3321 :
3322 : /************************************************************************/
3323 : /* GetNoDataValue() */
3324 : /************************************************************************/
3325 :
3326 2 : double FITSRasterBand::GetNoDataValue(int *pbSuccess)
3327 :
3328 : {
3329 2 : if (m_bNoDataSet)
3330 : {
3331 0 : if (pbSuccess)
3332 0 : *pbSuccess = TRUE;
3333 :
3334 0 : return m_dfNoDataValue;
3335 : }
3336 :
3337 2 : if (m_poFDS->m_bNoDataSet)
3338 : {
3339 2 : if (pbSuccess)
3340 2 : *pbSuccess = TRUE;
3341 :
3342 2 : return m_poFDS->m_dfNoDataValue;
3343 : }
3344 :
3345 0 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
3346 : }
3347 :
3348 : /************************************************************************/
3349 : /* SetNoDataValue() */
3350 : /************************************************************************/
3351 :
3352 1 : CPLErr FITSRasterBand::SetNoDataValue(double dfNoData)
3353 :
3354 : {
3355 1 : if (m_poFDS->m_bNoDataSet && m_poFDS->m_dfNoDataValue == dfNoData)
3356 : {
3357 0 : m_bNoDataSet = true;
3358 0 : m_dfNoDataValue = dfNoData;
3359 0 : return CE_None;
3360 : }
3361 :
3362 1 : m_poFDS->m_bNoDataSet = true;
3363 1 : m_poFDS->m_dfNoDataValue = dfNoData;
3364 :
3365 1 : m_poFDS->m_bNoDataChanged = true;
3366 :
3367 1 : m_bNoDataSet = true;
3368 1 : m_dfNoDataValue = dfNoData;
3369 1 : return CE_None;
3370 : }
3371 :
3372 : /************************************************************************/
3373 : /* DeleteNoDataValue() */
3374 : /************************************************************************/
3375 :
3376 0 : CPLErr FITSRasterBand::DeleteNoDataValue()
3377 :
3378 : {
3379 0 : if (!m_poFDS->m_bNoDataSet)
3380 0 : return CE_None;
3381 :
3382 0 : m_poFDS->m_bNoDataSet = false;
3383 0 : m_poFDS->m_dfNoDataValue = -9999.0;
3384 :
3385 0 : m_poFDS->m_bNoDataChanged = true;
3386 :
3387 0 : m_bNoDataSet = false;
3388 0 : m_dfNoDataValue = -9999.0;
3389 0 : return CE_None;
3390 : }
3391 :
3392 : /************************************************************************/
3393 : /* LoadGeoreferencing() */
3394 : /************************************************************************/
3395 :
3396 42 : void FITSDataset::LoadGeoreferencing()
3397 : {
3398 42 : int status = 0;
3399 : double crpix1, crpix2, crval1, crval2, cdelt1, cdelt2, pc[4], cd[4];
3400 42 : double aRadius, cRadius, invFlattening = 0.0;
3401 42 : double falseEast = 0.0, falseNorth = 0.0, scale = 1.0;
3402 : char target[81], ctype[81];
3403 42 : std::string GeogName, DatumName, projName;
3404 :
3405 42 : const double PI = std::atan(1.0) * 4;
3406 42 : const double DEG2RAD = PI / 180.;
3407 :
3408 : /* -------------------------------------------------------------------- */
3409 : /* Get the transform from the FITS file. */
3410 : /* -------------------------------------------------------------------- */
3411 :
3412 42 : fits_read_key(m_hFITS, TSTRING, "OBJECT", target, nullptr, &status);
3413 42 : if (status)
3414 : {
3415 9 : strncpy(target, "Undefined", 10);
3416 9 : CPLDebug("FITS", "OBJECT keyword is missing");
3417 9 : status = 0;
3418 : }
3419 :
3420 42 : GeogName.assign("GCS_");
3421 42 : GeogName.append(target);
3422 42 : DatumName.assign("D_");
3423 42 : DatumName.append(target);
3424 :
3425 42 : fits_read_key(m_hFITS, TDOUBLE, "A_RADIUS", &aRadius, nullptr, &status);
3426 42 : if (status)
3427 : {
3428 9 : CPLDebug("FITS", "No Radii keyword available, metadata will not "
3429 : "contain DATUM information.");
3430 9 : return;
3431 : }
3432 : else
3433 : {
3434 33 : fits_read_key(m_hFITS, TDOUBLE, "C_RADIUS", &cRadius, nullptr, &status);
3435 33 : if (status)
3436 : {
3437 0 : CPLError(CE_Warning, CPLE_AppDefined,
3438 : "No polar radius keyword available, setting C_RADIUS = "
3439 : "A_RADIUS");
3440 0 : cRadius = aRadius;
3441 0 : status = 0;
3442 : }
3443 33 : if (aRadius != cRadius)
3444 : {
3445 33 : invFlattening = aRadius / (aRadius - cRadius);
3446 : }
3447 : }
3448 :
3449 : /* Waiting for linear keywords standardization only deg ctype are used */
3450 : /* Check if WCS are there */
3451 33 : fits_read_key(m_hFITS, TSTRING, "CTYPE1", ctype, nullptr, &status);
3452 33 : if (!status)
3453 : {
3454 : /* Check if angular WCS are there */
3455 16 : if (strstr(ctype, "LN"))
3456 : {
3457 : /* Reading reference points */
3458 16 : fits_read_key(m_hFITS, TDOUBLE, "CRPIX1", &crpix1, nullptr,
3459 : &status);
3460 16 : fits_read_key(m_hFITS, TDOUBLE, "CRPIX2", &crpix2, nullptr,
3461 : &status);
3462 16 : fits_read_key(m_hFITS, TDOUBLE, "CRVAL1", &crval1, nullptr,
3463 : &status);
3464 16 : fits_read_key(m_hFITS, TDOUBLE, "CRVAL2", &crval2, nullptr,
3465 : &status);
3466 16 : if (status)
3467 : {
3468 0 : CPLError(CE_Failure, CPLE_AppDefined,
3469 : "No CRPIX / CRVAL keyword available, the raster "
3470 : "cannot be georeferenced.");
3471 0 : status = 0;
3472 : }
3473 : else
3474 : {
3475 : /* Check for CDELT and PC matrix representation */
3476 16 : fits_read_key(m_hFITS, TDOUBLE, "CDELT1", &cdelt1, nullptr,
3477 : &status);
3478 16 : if (!status)
3479 : {
3480 16 : fits_read_key(m_hFITS, TDOUBLE, "CDELT2", &cdelt2, nullptr,
3481 : &status);
3482 16 : fits_read_key(m_hFITS, TDOUBLE, "PC1_1", &pc[0], nullptr,
3483 : &status);
3484 16 : fits_read_key(m_hFITS, TDOUBLE, "PC1_2", &pc[1], nullptr,
3485 : &status);
3486 16 : fits_read_key(m_hFITS, TDOUBLE, "PC2_1", &pc[2], nullptr,
3487 : &status);
3488 16 : fits_read_key(m_hFITS, TDOUBLE, "PC2_2", &pc[3], nullptr,
3489 : &status);
3490 16 : cd[0] = cdelt1 * pc[0];
3491 16 : cd[1] = cdelt1 * pc[1];
3492 16 : cd[2] = cdelt2 * pc[2];
3493 16 : cd[3] = cdelt2 * pc[3];
3494 16 : status = 0;
3495 : }
3496 : else
3497 : {
3498 : /* Look for CD matrix representation */
3499 0 : fits_read_key(m_hFITS, TDOUBLE, "CD1_1", &cd[0], nullptr,
3500 : &status);
3501 0 : fits_read_key(m_hFITS, TDOUBLE, "CD1_2", &cd[1], nullptr,
3502 : &status);
3503 0 : fits_read_key(m_hFITS, TDOUBLE, "CD2_1", &cd[2], nullptr,
3504 : &status);
3505 0 : fits_read_key(m_hFITS, TDOUBLE, "CD2_2", &cd[3], nullptr,
3506 : &status);
3507 : }
3508 :
3509 16 : double radfac = DEG2RAD * aRadius;
3510 :
3511 16 : m_adfGeoTransform[1] = cd[0] * radfac;
3512 16 : m_adfGeoTransform[2] = cd[1] * radfac;
3513 16 : m_adfGeoTransform[4] = cd[2] * radfac;
3514 16 : m_adfGeoTransform[5] = -cd[3] * radfac;
3515 16 : if (crval1 > 180.)
3516 : {
3517 0 : crval1 = crval1 - 180.;
3518 : }
3519 :
3520 : /* NOTA BENE: FITS standard define pixel integers at the center
3521 : of the pixel, 0.5 must be subtract to have UpperLeft corner
3522 : */
3523 16 : m_adfGeoTransform[0] =
3524 16 : crval1 * radfac - m_adfGeoTransform[1] * (crpix1 - 0.5);
3525 : // assuming that center latitude is also the origin of the
3526 : // coordinate system: this is not always true. More generic
3527 : // implementation coming soon
3528 16 : m_adfGeoTransform[3] = -m_adfGeoTransform[5] * (crpix2 - 0.5);
3529 : //+ crval2 * radfac;
3530 16 : m_bGeoTransformValid = true;
3531 : }
3532 :
3533 16 : char *pstr = strrchr(ctype, '-');
3534 16 : if (pstr)
3535 : {
3536 16 : pstr += 1;
3537 :
3538 : /* Defining projection type
3539 : Following http://www.gdal.org/ogr__srs__api_8h.html (GDAL)
3540 : and
3541 : http://www.aanda.org/component/article?access=bibcode&bibcode=&bibcode=2002A%2526A...395.1077CFUL
3542 : (FITS)
3543 : */
3544 :
3545 : /* Sinusoidal / SFL projection */
3546 16 : if (strcmp(pstr, "SFL") == 0)
3547 : {
3548 0 : projName.assign("Sinusoidal_");
3549 0 : m_oSRS.SetSinusoidal(crval1, falseEast, falseNorth);
3550 :
3551 : /* Mercator, Oblique (Hotine) Mercator, Transverse Mercator
3552 : */
3553 : /* Mercator / MER projection */
3554 : }
3555 16 : else if (strcmp(pstr, "MER") == 0)
3556 : {
3557 16 : projName.assign("Mercator_");
3558 16 : m_oSRS.SetMercator(crval2, crval1, scale, falseEast,
3559 : falseNorth);
3560 :
3561 : /* Equirectangular / CAR projection */
3562 : }
3563 0 : else if (strcmp(pstr, "CAR") == 0)
3564 : {
3565 0 : projName.assign("Equirectangular_");
3566 : /*
3567 : The standard_parallel_1 defines where the local radius is
3568 : calculated not the center of Y Cartesian system (which is
3569 : latitude_of_origin) But FITS WCS only supports projections
3570 : on the sphere we assume here that the local radius is the
3571 : one computed at the projection center
3572 : */
3573 0 : m_oSRS.SetEquirectangular2(crval2, crval1, crval2,
3574 : falseEast, falseNorth);
3575 : /* Lambert Azimuthal Equal Area / ZEA projection */
3576 : }
3577 0 : else if (strcmp(pstr, "ZEA") == 0)
3578 : {
3579 0 : projName.assign("Lambert_Azimuthal_Equal_Area_");
3580 0 : m_oSRS.SetLAEA(crval2, crval1, falseEast, falseNorth);
3581 :
3582 : /* Lambert Conformal Conic 1SP / COO projection */
3583 : }
3584 0 : else if (strcmp(pstr, "COO") == 0)
3585 : {
3586 0 : projName.assign("Lambert_Conformal_Conic_1SP_");
3587 0 : m_oSRS.SetLCC1SP(crval2, crval1, scale, falseEast,
3588 : falseNorth);
3589 :
3590 : /* Orthographic / SIN projection */
3591 : }
3592 0 : else if (strcmp(pstr, "SIN") == 0)
3593 : {
3594 0 : projName.assign("Orthographic_");
3595 0 : m_oSRS.SetOrthographic(crval2, crval1, falseEast,
3596 : falseNorth);
3597 :
3598 : /* Point Perspective / AZP projection */
3599 : }
3600 0 : else if (strcmp(pstr, "AZP") == 0)
3601 : {
3602 0 : projName.assign("perspective_point_height_");
3603 0 : m_oSRS.SetProjection(SRS_PP_PERSPECTIVE_POINT_HEIGHT);
3604 : /* # appears to need height... maybe center lon/lat */
3605 :
3606 : /* Polar Stereographic / STG projection */
3607 : }
3608 0 : else if (strcmp(pstr, "STG") == 0)
3609 : {
3610 0 : projName.assign("Polar_Stereographic_");
3611 0 : m_oSRS.SetStereographic(crval2, crval1, scale, falseEast,
3612 : falseNorth);
3613 : }
3614 : else
3615 : {
3616 0 : CPLError(CE_Failure, CPLE_AppDefined,
3617 : "Unknown projection.");
3618 : }
3619 :
3620 16 : projName.append(target);
3621 16 : m_oSRS.SetProjParm(SRS_PP_FALSE_EASTING, 0.0);
3622 16 : m_oSRS.SetProjParm(SRS_PP_FALSE_NORTHING, 0.0);
3623 :
3624 16 : m_oSRS.SetNode("PROJCS", projName.c_str());
3625 :
3626 16 : m_oSRS.SetGeogCS(GeogName.c_str(), DatumName.c_str(), target,
3627 : aRadius, invFlattening, "Reference_Meridian",
3628 : 0.0, "degree", 0.0174532925199433);
3629 : }
3630 : else
3631 : {
3632 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown projection.");
3633 : }
3634 : }
3635 : }
3636 : else
3637 : {
3638 17 : CPLError(CE_Warning, CPLE_AppDefined,
3639 : "No CTYPE keywords: no geospatial information available.");
3640 : }
3641 : }
3642 :
3643 : /************************************************************************/
3644 : /* LoadFITSInfo() */
3645 : /************************************************************************/
3646 :
3647 42 : void FITSDataset::LoadFITSInfo()
3648 :
3649 : {
3650 42 : int status = 0;
3651 : int bitpix;
3652 : double dfScale, dfOffset;
3653 :
3654 42 : LoadGeoreferencing();
3655 :
3656 42 : CPLAssert(!m_bMetadataChanged);
3657 42 : CPLAssert(!m_bNoDataChanged);
3658 :
3659 42 : m_bMetadataChanged = false;
3660 42 : m_bNoDataChanged = false;
3661 :
3662 42 : bitpix = this->m_fitsDataType;
3663 42 : FITSRasterBand *poBand = cpl::down_cast<FITSRasterBand *>(GetRasterBand(1));
3664 :
3665 42 : if (bitpix != TUSHORT && bitpix != TUINT)
3666 : {
3667 36 : fits_read_key(m_hFITS, TDOUBLE, "BSCALE", &dfScale, nullptr, &status);
3668 36 : if (status)
3669 : {
3670 34 : status = 0;
3671 34 : dfScale = 1.;
3672 : }
3673 36 : fits_read_key(m_hFITS, TDOUBLE, "BZERO", &dfOffset, nullptr, &status);
3674 36 : if (status)
3675 : {
3676 34 : status = 0;
3677 34 : dfOffset = 0.;
3678 : }
3679 36 : if (dfScale != 1. || dfOffset != 0.)
3680 : {
3681 2 : poBand->m_bHaveOffsetScale = true;
3682 2 : poBand->m_dfScale = dfScale;
3683 2 : poBand->m_dfOffset = dfOffset;
3684 : }
3685 : }
3686 :
3687 42 : fits_read_key(m_hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue, nullptr,
3688 : &status);
3689 42 : m_bNoDataSet = !status;
3690 42 : }
3691 :
3692 : /************************************************************************/
3693 : /* GDALRegister_FITS() */
3694 : /************************************************************************/
3695 :
3696 11 : void GDALRegister_FITS()
3697 :
3698 : {
3699 11 : if (GDALGetDriverByName(DRIVER_NAME) != nullptr)
3700 0 : return;
3701 :
3702 11 : GDALDriver *poDriver = new GDALDriver();
3703 11 : FITSDriverSetCommonMetadata(poDriver);
3704 :
3705 11 : poDriver->pfnOpen = FITSDataset::Open;
3706 11 : poDriver->pfnCreate = FITSDataset::Create;
3707 11 : poDriver->pfnDelete = FITSDataset::Delete;
3708 :
3709 11 : GetGDALDriverManager()->RegisterDriver(poDriver);
3710 : }
|
