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