Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: FITS Driver
4 : * Purpose: Implement FITS raster read/write support
5 : * Author: Simon Perkins, s.perkins@lanl.gov
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2001, Simon Perkins
9 : * Copyright (c) 2008-2020, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2018, Chiara Marmo <chiara dot marmo at u-psud dot fr>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : // So that OFF_T is 64 bits
16 : #define _FILE_OFFSET_BITS 64
17 :
18 : #include "cpl_string.h"
19 : #include "gdal_frmts.h"
20 : #include "gdal_pam.h"
21 : #include "ogr_spatialref.h"
22 : #include "ogrsf_frmts.h"
23 : #include "fitsdrivercore.h"
24 :
25 : #include <string.h>
26 : #include <fitsio.h>
27 :
28 : #include <algorithm>
29 : #include <string>
30 : #include <cstring>
31 : #include <set>
32 : #include <vector>
33 :
34 : /************************************************************************/
35 : /* ==================================================================== */
36 : /* FITSDataset */
37 : /* ==================================================================== */
38 : /************************************************************************/
39 :
40 : class FITSRasterBand;
41 : class FITSLayer;
42 :
43 : class FITSDataset final : public GDALPamDataset
44 : {
45 :
46 : friend class FITSRasterBand;
47 : friend class FITSLayer;
48 :
49 : fitsfile *m_hFITS = nullptr;
50 :
51 : int m_hduNum = 0;
52 : GDALDataType m_gdalDataType = GDT_Unknown; // GDAL code for the image type
53 : int m_fitsDataType = 0; // FITS code for the image type
54 :
55 : bool m_isExistingFile = false;
56 : LONGLONG m_highestOffsetWritten = 0; // How much of image has been written
57 :
58 : bool m_bNoDataChanged = false;
59 : bool m_bNoDataSet = false;
60 : double m_dfNoDataValue = -9999.0;
61 :
62 : bool m_bMetadataChanged = false;
63 :
64 : CPLStringList m_aosSubdatasets{};
65 :
66 : OGRSpatialReference m_oSRS{};
67 :
68 : double m_adfGeoTransform[6];
69 : bool m_bGeoTransformValid = false;
70 :
71 : bool m_bFITSInfoChanged = false;
72 :
73 : std::vector<std::unique_ptr<FITSLayer>> m_apoLayers{};
74 :
75 : CPLErr Init(fitsfile *hFITS, bool isExistingFile, int hduNum);
76 :
77 : void LoadGeoreferencing();
78 : void LoadFITSInfo();
79 : void WriteFITSInfo();
80 : void LoadMetadata(GDALMajorObject *poTarget);
81 :
82 : 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(double *) override;
97 : virtual CPLErr SetGeoTransform(double *) 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 : m_adfGeoTransform[0] = 0;
1864 103 : m_adfGeoTransform[1] = 1;
1865 103 : m_adfGeoTransform[2] = 0;
1866 103 : m_adfGeoTransform[3] = 0;
1867 103 : m_adfGeoTransform[4] = 0;
1868 103 : m_adfGeoTransform[5] = 1;
1869 103 : }
1870 :
1871 : /************************************************************************/
1872 : /* ~FITSDataset() */
1873 : /************************************************************************/
1874 :
1875 206 : FITSDataset::~FITSDataset()
1876 : {
1877 :
1878 103 : int status = 0;
1879 103 : if (m_hFITS)
1880 : {
1881 103 : m_apoLayers.clear();
1882 :
1883 103 : if (m_hduNum > 0 && eAccess == GA_Update)
1884 : {
1885 : // Only do this if we've successfully opened the file and update
1886 : // capability. Write any meta data to the file that's compatible
1887 : // with FITS.
1888 41 : fits_movabs_hdu(m_hFITS, m_hduNum, nullptr, &status);
1889 41 : fits_write_key_longwarn(m_hFITS, &status);
1890 41 : if (status)
1891 : {
1892 0 : CPLError(CE_Warning, CPLE_AppDefined,
1893 : "Couldn't move to HDU %d in FITS file %s (%d).\n",
1894 : m_hduNum, GetDescription(), status);
1895 : }
1896 41 : char **metaData = FITSDataset::GetMetadata();
1897 41 : int count = CSLCount(metaData);
1898 65 : for (int i = 0; i < count; ++i)
1899 : {
1900 24 : const char *field = CSLGetField(metaData, i);
1901 24 : if (strlen(field) == 0)
1902 0 : continue;
1903 : else
1904 : {
1905 24 : char *key = nullptr;
1906 24 : const char *value = CPLParseNameValue(field, &key);
1907 : // FITS keys must be less than 8 chars
1908 26 : if (key != nullptr && strlen(key) <= 8 &&
1909 2 : !isIgnorableFITSHeader(key))
1910 : {
1911 : // Although FITS provides support for different value
1912 : // types, the GDAL Metadata mechanism works only with
1913 : // string values. Prior to about 2003-05-02, this driver
1914 : // would attempt to guess the value type from the
1915 : // metadata value string amd then would use the
1916 : // appropriate type-specific FITS keyword update
1917 : // routine. This was found to be troublesome (e.g. a
1918 : // numeric version string with leading zeros would be
1919 : // interpreted as a number and might get those leading
1920 : // zeros stripped), and so now the driver writes every
1921 : // value as a string. In practice this is not a problem
1922 : // since most FITS reading routines will convert from
1923 : // strings to numbers automatically, but if you want
1924 : // finer control, use the underlying FITS handle. Note:
1925 : // to avoid a compiler warning we copy the const value
1926 : // string to a non const one.
1927 2 : char *valueCpy = CPLStrdup(value);
1928 2 : fits_update_key_longstr(m_hFITS, key, valueCpy, nullptr,
1929 : &status);
1930 2 : CPLFree(valueCpy);
1931 :
1932 : // Check for errors.
1933 2 : if (status)
1934 : {
1935 : // Throw a warning with CFITSIO error status, then
1936 : // ignore status
1937 0 : CPLError(
1938 : CE_Warning, CPLE_AppDefined,
1939 : "Couldn't update key %s in FITS file %s (%d).",
1940 : key, GetDescription(), status);
1941 0 : status = 0;
1942 0 : return;
1943 : }
1944 : }
1945 : // Must free up key
1946 24 : CPLFree(key);
1947 : }
1948 : }
1949 :
1950 : // Writing nodata value
1951 41 : if (m_gdalDataType != GDT_Float32 && m_gdalDataType != GDT_Float64)
1952 : {
1953 33 : fits_update_key(m_hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue,
1954 : nullptr, &status);
1955 33 : if (status)
1956 : {
1957 : // Throw a warning with CFITSIO error status, then ignore
1958 : // status
1959 0 : CPLError(CE_Warning, CPLE_AppDefined,
1960 : "Couldn't update key BLANK in FITS file %s (%d).",
1961 : GetDescription(), status);
1962 0 : status = 0;
1963 0 : return;
1964 : }
1965 : }
1966 :
1967 : // Writing Scale and offset if defined
1968 : int pbSuccess;
1969 41 : GDALRasterBand *poSrcBand = GDALPamDataset::GetRasterBand(1);
1970 41 : double dfScale = poSrcBand->GetScale(&pbSuccess);
1971 41 : double dfOffset = poSrcBand->GetOffset(&pbSuccess);
1972 41 : if (m_bMetadataChanged)
1973 : {
1974 1 : fits_update_key(m_hFITS, TDOUBLE, "BSCALE", &dfScale, nullptr,
1975 : &status);
1976 1 : if (status)
1977 : {
1978 : // Throw a warning with CFITSIO error status, then ignore
1979 : // status
1980 0 : CPLError(CE_Warning, CPLE_AppDefined,
1981 : "Couldn't update key BSCALE in FITS file %s (%d).",
1982 : GetDescription(), status);
1983 0 : status = 0;
1984 0 : return;
1985 : }
1986 1 : fits_update_key(m_hFITS, TDOUBLE, "BZERO", &dfOffset, nullptr,
1987 : &status);
1988 1 : if (status)
1989 : {
1990 : // Throw a warning with CFITSIO error status, then ignore
1991 : // status
1992 0 : CPLError(CE_Warning, CPLE_AppDefined,
1993 : "Couldn't update key BZERO in FITS file %s (%d).",
1994 : GetDescription(), status);
1995 0 : status = 0;
1996 0 : return;
1997 : }
1998 : }
1999 :
2000 : // Copy georeferencing info to PAM if the profile is not FITS
2001 41 : GDALPamDataset::SetSpatialRef(GDALPamDataset::GetSpatialRef());
2002 :
2003 : // Write geographic info
2004 41 : if (m_bFITSInfoChanged)
2005 : {
2006 40 : WriteFITSInfo();
2007 : }
2008 :
2009 : // Make sure we flush the raster cache before we close the file!
2010 41 : FlushCache(true);
2011 : }
2012 :
2013 : // Close the FITS handle
2014 103 : fits_close_file(m_hFITS, &status);
2015 103 : if (status != 0)
2016 : {
2017 0 : CPLError(CE_Failure, CPLE_AppDefined,
2018 : "fits_close_file() failed with %d", status);
2019 : }
2020 : }
2021 206 : }
2022 :
2023 : /************************************************************************/
2024 : /* GetRawBinaryLayout() */
2025 : /************************************************************************/
2026 :
2027 2 : bool FITSDataset::GetRawBinaryLayout(GDALDataset::RawBinaryLayout &sLayout)
2028 : {
2029 2 : if (m_hduNum == 0)
2030 0 : return false;
2031 2 : int status = 0;
2032 2 : if (fits_is_compressed_image(m_hFITS, &status))
2033 0 : return false;
2034 2 : GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
2035 2 : if (eDT == GDT_UInt16 || eDT == GDT_UInt32)
2036 0 : return false; // are supported as native signed with offset
2037 :
2038 2 : sLayout.osRawFilename = GetDescription();
2039 : static_assert(sizeof(OFF_T) == 8, "OFF_T should be 64 bits !");
2040 2 : OFF_T headerstart = 0;
2041 2 : OFF_T datastart = 0;
2042 2 : OFF_T dataend = 0;
2043 2 : fits_get_hduoff(m_hFITS, &headerstart, &datastart, &dataend, &status);
2044 2 : if (nBands > 1)
2045 0 : sLayout.eInterleaving = RawBinaryLayout::Interleaving::BSQ;
2046 2 : sLayout.eDataType = eDT;
2047 2 : sLayout.bLittleEndianOrder = false;
2048 2 : sLayout.nImageOffset = static_cast<GIntBig>(datastart);
2049 2 : sLayout.nPixelOffset = GDALGetDataTypeSizeBytes(eDT);
2050 2 : sLayout.nLineOffset = sLayout.nPixelOffset * nRasterXSize;
2051 2 : sLayout.nBandOffset = sLayout.nLineOffset * nRasterYSize;
2052 2 : return true;
2053 : }
2054 :
2055 : /************************************************************************/
2056 : /* Init() */
2057 : /************************************************************************/
2058 :
2059 80 : CPLErr FITSDataset::Init(fitsfile *hFITS, bool isExistingFile, int hduNum)
2060 : {
2061 :
2062 80 : m_hFITS = hFITS;
2063 80 : m_isExistingFile = isExistingFile;
2064 :
2065 80 : int status = 0;
2066 : double offset;
2067 :
2068 80 : int hduType = 0;
2069 80 : fits_movabs_hdu(hFITS, hduNum, &hduType, &status);
2070 80 : if (status)
2071 : {
2072 1 : CPLError(CE_Failure, CPLE_AppDefined,
2073 : "Couldn't move to HDU %d in FITS file %s (%d).", hduNum,
2074 1 : GetDescription(), status);
2075 1 : return CE_Failure;
2076 : }
2077 :
2078 79 : if (hduType != IMAGE_HDU)
2079 : {
2080 0 : CPLError(CE_Failure, CPLE_AppDefined, "HDU %d is not an image.",
2081 : hduNum);
2082 0 : return CE_Failure;
2083 : }
2084 :
2085 : // Get the image info for this dataset (note that all bands in a FITS
2086 : // dataset have the same type)
2087 79 : int bitpix = 0;
2088 79 : int naxis = 0;
2089 79 : const int maxdim = 3;
2090 79 : long naxes[maxdim] = {0, 0, 0};
2091 79 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes, &status);
2092 79 : if (status)
2093 : {
2094 0 : CPLError(CE_Failure, CPLE_AppDefined,
2095 : "Couldn't determine image parameters of FITS file %s (%d)",
2096 0 : GetDescription(), status);
2097 0 : return CE_Failure;
2098 : }
2099 :
2100 79 : m_hduNum = hduNum;
2101 :
2102 79 : fits_read_key(hFITS, TDOUBLE, "BZERO", &offset, nullptr, &status);
2103 79 : if (status)
2104 : {
2105 : // BZERO is not mandatory offset defaulted to 0 if BZERO is missing
2106 63 : status = 0;
2107 63 : offset = 0.;
2108 : }
2109 :
2110 79 : fits_read_key(hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue, nullptr, &status);
2111 79 : m_bNoDataSet = !status;
2112 79 : status = 0;
2113 :
2114 : // Determine data type and nodata value if BLANK keyword is absent
2115 79 : if (bitpix == BYTE_IMG)
2116 : {
2117 33 : m_gdalDataType = GDT_Byte;
2118 33 : m_fitsDataType = TBYTE;
2119 : }
2120 46 : else if (bitpix == SHORT_IMG)
2121 : {
2122 18 : if (offset == 32768.)
2123 : {
2124 7 : m_gdalDataType = GDT_UInt16;
2125 7 : m_fitsDataType = TUSHORT;
2126 : }
2127 : else
2128 : {
2129 11 : m_gdalDataType = GDT_Int16;
2130 11 : m_fitsDataType = TSHORT;
2131 : }
2132 : }
2133 28 : else if (bitpix == LONG_IMG)
2134 : {
2135 14 : if (offset == 2147483648.)
2136 : {
2137 7 : m_gdalDataType = GDT_UInt32;
2138 7 : m_fitsDataType = TUINT;
2139 : }
2140 : else
2141 : {
2142 7 : m_gdalDataType = GDT_Int32;
2143 7 : m_fitsDataType = TINT;
2144 : }
2145 : }
2146 14 : else if (bitpix == FLOAT_IMG)
2147 : {
2148 7 : m_gdalDataType = GDT_Float32;
2149 7 : m_fitsDataType = TFLOAT;
2150 : }
2151 7 : else if (bitpix == DOUBLE_IMG)
2152 : {
2153 7 : m_gdalDataType = GDT_Float64;
2154 7 : m_fitsDataType = TDOUBLE;
2155 : }
2156 : else
2157 : {
2158 0 : CPLError(CE_Failure, CPLE_AppDefined,
2159 0 : "FITS file %s has unknown data type: %d.", GetDescription(),
2160 : bitpix);
2161 0 : return CE_Failure;
2162 : }
2163 :
2164 : // Determine image dimensions - we assume BSQ ordering
2165 79 : if (naxis == 2)
2166 : {
2167 55 : nRasterXSize = static_cast<int>(naxes[0]);
2168 55 : nRasterYSize = static_cast<int>(naxes[1]);
2169 55 : nBands = 1;
2170 : }
2171 24 : else if (naxis == 3)
2172 : {
2173 24 : nRasterXSize = static_cast<int>(naxes[0]);
2174 24 : nRasterYSize = static_cast<int>(naxes[1]);
2175 24 : nBands = static_cast<int>(naxes[2]);
2176 : }
2177 : else
2178 : {
2179 0 : CPLError(CE_Failure, CPLE_AppDefined,
2180 : "FITS file %s does not have 2 or 3 dimensions.",
2181 0 : GetDescription());
2182 0 : return CE_Failure;
2183 : }
2184 :
2185 : // Create the bands
2186 212 : for (int i = 0; i < nBands; ++i)
2187 133 : SetBand(i + 1, new FITSRasterBand(this, i + 1));
2188 :
2189 79 : return CE_None;
2190 : }
2191 :
2192 : /************************************************************************/
2193 : /* LoadMetadata() */
2194 : /************************************************************************/
2195 :
2196 62 : void FITSDataset::LoadMetadata(GDALMajorObject *poTarget)
2197 : {
2198 : // Read header information from file and use it to set metadata
2199 : // This process understands the CONTINUE standard for long strings.
2200 : // We don't bother to capture header names that duplicate information
2201 : // already captured elsewhere (e.g. image dimensions and type)
2202 : int keyNum;
2203 : char key[100];
2204 : char value[100];
2205 62 : CPLStringList aosMD;
2206 :
2207 62 : int nKeys = 0;
2208 62 : int nMoreKeys = 0;
2209 62 : int status = 0;
2210 62 : fits_get_hdrspace(m_hFITS, &nKeys, &nMoreKeys, &status);
2211 2329 : for (keyNum = 1; keyNum <= nKeys; keyNum++)
2212 : {
2213 2267 : fits_read_keyn(m_hFITS, keyNum, key, value, nullptr, &status);
2214 2267 : if (status)
2215 : {
2216 0 : CPLError(CE_Failure, CPLE_AppDefined,
2217 : "Error while reading key %d from FITS file %s (%d)",
2218 0 : keyNum, GetDescription(), status);
2219 0 : return;
2220 : }
2221 2267 : if (strcmp(key, "END") == 0)
2222 : {
2223 : // We should not get here in principle since the END
2224 : // keyword shouldn't be counted in nKeys, but who knows.
2225 0 : break;
2226 : }
2227 2267 : else if (isIgnorableFITSHeader(key))
2228 : {
2229 : // Ignore it
2230 : }
2231 : else
2232 : { // Going to store something, but check for long strings etc
2233 : // Strip off leading and trailing quote if present
2234 1561 : char *newValue = value;
2235 1561 : if (value[0] == '\'' && value[strlen(value) - 1] == '\'')
2236 : {
2237 996 : newValue = value + 1;
2238 996 : value[strlen(value) - 1] = '\0';
2239 : }
2240 : // Check for long string
2241 1561 : if (strrchr(newValue, '&') == newValue + strlen(newValue) - 1)
2242 : {
2243 : // Value string ends in "&", so use long string conventions
2244 0 : char *longString = nullptr;
2245 0 : fits_read_key_longstr(m_hFITS, key, &longString, nullptr,
2246 : &status);
2247 : // Note that read_key_longstr already strips quotes
2248 0 : if (status)
2249 : {
2250 0 : CPLError(CE_Failure, CPLE_AppDefined,
2251 : "Error while reading long string for key %s from "
2252 : "FITS file %s (%d)",
2253 0 : key, GetDescription(), status);
2254 0 : return;
2255 : }
2256 0 : poTarget->SetMetadataItem(key, longString);
2257 0 : free(longString);
2258 : }
2259 : else
2260 : { // Normal keyword
2261 1561 : poTarget->SetMetadataItem(key, newValue);
2262 : }
2263 : }
2264 : }
2265 : }
2266 :
2267 : /************************************************************************/
2268 : /* GetMetadata() */
2269 : /************************************************************************/
2270 :
2271 70 : char **FITSDataset::GetMetadata(const char *pszDomain)
2272 :
2273 : {
2274 70 : if (pszDomain != nullptr && EQUAL(pszDomain, "SUBDATASETS"))
2275 : {
2276 8 : return m_aosSubdatasets.List();
2277 : }
2278 :
2279 62 : return GDALPamDataset::GetMetadata(pszDomain);
2280 : }
2281 :
2282 : /************************************************************************/
2283 : /* GetLayer() */
2284 : /************************************************************************/
2285 :
2286 15 : OGRLayer *FITSDataset::GetLayer(int idx)
2287 : {
2288 15 : if (idx < 0 || idx >= GetLayerCount())
2289 2 : return nullptr;
2290 13 : return m_apoLayers[idx].get();
2291 : }
2292 :
2293 : /************************************************************************/
2294 : /* ICreateLayer() */
2295 : /************************************************************************/
2296 :
2297 4 : OGRLayer *FITSDataset::ICreateLayer(const char *pszName,
2298 : const OGRGeomFieldDefn *poGeomFieldDefn,
2299 : CSLConstList papszOptions)
2300 : {
2301 4 : if (!TestCapability(ODsCCreateLayer))
2302 0 : return nullptr;
2303 :
2304 4 : const auto eGType = poGeomFieldDefn ? poGeomFieldDefn->GetType() : wkbNone;
2305 4 : if (eGType != wkbNone)
2306 : {
2307 0 : CPLError(CE_Failure, CPLE_NotSupported, "Spatial tables not supported");
2308 0 : return nullptr;
2309 : }
2310 :
2311 4 : int status = 0;
2312 4 : int numHDUs = 0;
2313 4 : fits_get_num_hdus(m_hFITS, &numHDUs, &status);
2314 4 : if (status != 0)
2315 : {
2316 0 : CPLError(CE_Failure, CPLE_AppDefined, "fits_get_num_hdus() failed: %d",
2317 : status);
2318 0 : return nullptr;
2319 : }
2320 :
2321 4 : fits_create_tbl(m_hFITS, BINARY_TBL,
2322 : 0, // number of initial rows
2323 : 0, // nfields,
2324 : nullptr, // ttype,
2325 : nullptr, // tform
2326 : nullptr, // tunits
2327 : pszName, &status);
2328 4 : if (status != 0)
2329 : {
2330 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create layer");
2331 0 : return nullptr;
2332 : }
2333 :
2334 : // If calling fits_get_num_hdus() here on a freshly new created file,
2335 : // it reports only one HDU, missing the initial dummy HDU
2336 4 : if (numHDUs == 0)
2337 : {
2338 4 : numHDUs = 2;
2339 : }
2340 : else
2341 : {
2342 0 : numHDUs++;
2343 : }
2344 :
2345 4 : auto poLayer = new FITSLayer(this, numHDUs, pszName);
2346 4 : poLayer->SetCreationOptions(papszOptions);
2347 4 : m_apoLayers.emplace_back(std::unique_ptr<FITSLayer>(poLayer));
2348 4 : return m_apoLayers.back().get();
2349 : }
2350 :
2351 : /************************************************************************/
2352 : /* TestCapability() */
2353 : /************************************************************************/
2354 :
2355 20 : int FITSDataset::TestCapability(const char *pszCap)
2356 : {
2357 20 : if (EQUAL(pszCap, ODsCCreateLayer))
2358 8 : return eAccess == GA_Update;
2359 12 : return false;
2360 : }
2361 :
2362 : /************************************************************************/
2363 : /* Open() */
2364 : /************************************************************************/
2365 :
2366 63 : GDALDataset *FITSDataset::Open(GDALOpenInfo *poOpenInfo)
2367 : {
2368 :
2369 63 : if (!FITSDriverIdentify(poOpenInfo))
2370 0 : return nullptr;
2371 :
2372 126 : CPLString osFilename(poOpenInfo->pszFilename);
2373 63 : int iSelectedHDU = 0;
2374 63 : if (STARTS_WITH(poOpenInfo->pszFilename, "FITS:"))
2375 : {
2376 : const CPLStringList aosTokens(
2377 9 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
2378 9 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
2379 9 : if (aosTokens.size() != 3)
2380 : {
2381 2 : return nullptr;
2382 : }
2383 7 : osFilename = aosTokens[1];
2384 7 : iSelectedHDU = atoi(aosTokens[2]);
2385 7 : if (iSelectedHDU <= 0)
2386 : {
2387 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid HDU number");
2388 1 : return nullptr;
2389 : }
2390 : }
2391 :
2392 : // Get access mode and attempt to open the file
2393 60 : int status = 0;
2394 60 : fitsfile *hFITS = nullptr;
2395 60 : if (poOpenInfo->eAccess == GA_ReadOnly)
2396 58 : fits_open_file(&hFITS, osFilename.c_str(), READONLY, &status);
2397 : else
2398 2 : fits_open_file(&hFITS, osFilename.c_str(), READWRITE, &status);
2399 60 : if (status)
2400 : {
2401 1 : CPLError(CE_Failure, CPLE_AppDefined,
2402 : "Error while opening FITS file %s (%d).\n", osFilename.c_str(),
2403 : status);
2404 1 : fits_close_file(hFITS, &status);
2405 1 : return nullptr;
2406 : }
2407 : // Create a FITSDataset object
2408 118 : auto dataset = std::make_unique<FITSDataset>();
2409 59 : dataset->m_isExistingFile = true;
2410 59 : dataset->m_hFITS = hFITS;
2411 59 : dataset->eAccess = poOpenInfo->eAccess;
2412 59 : dataset->SetPhysicalFilename(osFilename);
2413 :
2414 : /* -------------------------------------------------------------------- */
2415 : /* Iterate over HDUs */
2416 : /* -------------------------------------------------------------------- */
2417 59 : bool firstHDUIsDummy = false;
2418 59 : int firstValidHDU = 0;
2419 118 : CPLStringList aosSubdatasets;
2420 59 : bool hasVector = false;
2421 59 : if (iSelectedHDU == 0)
2422 : {
2423 54 : int numHDUs = 0;
2424 54 : fits_get_num_hdus(hFITS, &numHDUs, &status);
2425 54 : if (numHDUs <= 0)
2426 : {
2427 0 : return nullptr;
2428 : }
2429 :
2430 131 : for (int iHDU = 1; iHDU <= numHDUs; iHDU++)
2431 : {
2432 77 : int hduType = 0;
2433 77 : fits_movabs_hdu(hFITS, iHDU, &hduType, &status);
2434 77 : if (status)
2435 : {
2436 31 : continue;
2437 : }
2438 :
2439 77 : char szExtname[81] = {0};
2440 77 : fits_read_key(hFITS, TSTRING, "EXTNAME", szExtname, nullptr,
2441 : &status);
2442 77 : status = 0;
2443 77 : int nExtVer = 0;
2444 77 : fits_read_key(hFITS, TINT, "EXTVER", &nExtVer, nullptr, &status);
2445 77 : status = 0;
2446 77 : CPLString osExtname(szExtname);
2447 77 : if (nExtVer > 0)
2448 0 : osExtname += CPLSPrintf(" %d", nExtVer);
2449 :
2450 77 : if (hduType == BINARY_TBL)
2451 : {
2452 14 : hasVector = true;
2453 14 : if ((poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0)
2454 : {
2455 24 : dataset->m_apoLayers.push_back(std::unique_ptr<FITSLayer>(
2456 12 : new FITSLayer(dataset.get(), iHDU, osExtname.c_str())));
2457 : }
2458 : }
2459 :
2460 77 : if (hduType != IMAGE_HDU)
2461 : {
2462 14 : continue;
2463 : }
2464 :
2465 63 : int bitpix = 0;
2466 63 : int naxis = 0;
2467 63 : const int maxdim = 3;
2468 63 : long naxes[maxdim] = {0, 0, 0};
2469 63 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes, &status);
2470 63 : if (status)
2471 : {
2472 0 : continue;
2473 : }
2474 :
2475 63 : if (naxis != 2 && naxis != 3)
2476 : {
2477 17 : if (naxis == 0 && iHDU == 1)
2478 : {
2479 17 : firstHDUIsDummy = true;
2480 : }
2481 17 : continue;
2482 : }
2483 :
2484 46 : if ((poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0)
2485 : {
2486 41 : const int nIdx = aosSubdatasets.size() / 2 + 1;
2487 : aosSubdatasets.AddNameValue(
2488 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
2489 : CPLSPrintf("FITS:\"%s\":%d", poOpenInfo->pszFilename,
2490 41 : iHDU));
2491 : CPLString osDesc(CPLSPrintf(
2492 : "HDU %d (%dx%d, %d band%s)", iHDU,
2493 41 : static_cast<int>(naxes[0]), static_cast<int>(naxes[1]),
2494 51 : naxis == 3 ? static_cast<int>(naxes[2]) : 1,
2495 82 : (naxis == 3 && naxes[2] > 1) ? "s" : ""));
2496 41 : if (!osExtname.empty())
2497 : {
2498 5 : osDesc += ", ";
2499 5 : osDesc += osExtname;
2500 : }
2501 : aosSubdatasets.AddNameValue(
2502 41 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx), osDesc);
2503 : }
2504 :
2505 46 : if (firstValidHDU == 0)
2506 : {
2507 41 : firstValidHDU = iHDU;
2508 : }
2509 : }
2510 54 : if (aosSubdatasets.size() == 2)
2511 : {
2512 35 : aosSubdatasets.Clear();
2513 : }
2514 : }
2515 : else
2516 : {
2517 5 : if (iSelectedHDU != 1)
2518 : {
2519 4 : int hduType = 0;
2520 4 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2521 4 : if (status == 0)
2522 : {
2523 4 : int bitpix = 0;
2524 4 : int naxis = 0;
2525 4 : const int maxdim = 3;
2526 4 : long naxes[maxdim] = {0, 0, 0};
2527 4 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes,
2528 : &status);
2529 4 : if (status == 0 && naxis == 0)
2530 : {
2531 2 : firstHDUIsDummy = true;
2532 : }
2533 : }
2534 4 : status = 0;
2535 : }
2536 5 : firstValidHDU = iSelectedHDU;
2537 : }
2538 :
2539 59 : const bool hasRaster = firstValidHDU > 0;
2540 59 : const bool hasRasterAndIsAllowed =
2541 59 : hasRaster && (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0;
2542 :
2543 59 : if (!hasRasterAndIsAllowed &&
2544 16 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0 &&
2545 3 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) == 0)
2546 : {
2547 2 : if (hasVector)
2548 : {
2549 2 : std::string osPath;
2550 1 : osPath.resize(1024);
2551 1 : if (CPLGetExecPath(&osPath[0], static_cast<int>(osPath.size())))
2552 : {
2553 1 : osPath = CPLGetBasenameSafe(osPath.c_str());
2554 : }
2555 1 : if (osPath == "gdalinfo")
2556 : {
2557 0 : CPLError(CE_Failure, CPLE_AppDefined,
2558 : "This FITS dataset does not contain any image, but "
2559 : "contains binary table(s) that could be opened "
2560 : "in vector mode with ogrinfo.");
2561 : }
2562 : else
2563 : {
2564 1 : CPLError(CE_Failure, CPLE_AppDefined,
2565 : "This FITS dataset does not contain any image, but "
2566 : "contains binary table(s) that could be opened "
2567 : "in vector mode.");
2568 : }
2569 : }
2570 : else
2571 : {
2572 1 : CPLError(CE_Failure, CPLE_AppDefined,
2573 : "Cannot find HDU of image type with 2 or 3 axes.");
2574 : }
2575 2 : return nullptr;
2576 : }
2577 :
2578 57 : if (dataset->m_apoLayers.empty() &&
2579 60 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0 &&
2580 3 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0)
2581 : {
2582 3 : if (hasRaster)
2583 : {
2584 4 : std::string osPath;
2585 2 : osPath.resize(1024);
2586 2 : if (CPLGetExecPath(&osPath[0], static_cast<int>(osPath.size())))
2587 : {
2588 2 : osPath = CPLGetBasenameSafe(osPath.c_str());
2589 : }
2590 2 : if (osPath == "ogrinfo")
2591 : {
2592 0 : CPLError(CE_Failure, CPLE_AppDefined,
2593 : "This FITS dataset does not contain any binary "
2594 : "table, but contains image(s) that could be opened "
2595 : "in raster mode with gdalinfo.");
2596 : }
2597 : else
2598 : {
2599 2 : CPLError(CE_Failure, CPLE_AppDefined,
2600 : "This FITS dataset does not contain any binary "
2601 : "table, but contains image(s) that could be opened "
2602 : "in raster mode.");
2603 : }
2604 : }
2605 : else
2606 : {
2607 1 : CPLError(CE_Failure, CPLE_AppDefined,
2608 : "Cannot find binary table(s).");
2609 : }
2610 3 : return nullptr;
2611 : }
2612 :
2613 54 : dataset->m_aosSubdatasets = aosSubdatasets;
2614 :
2615 : // Set up the description and initialize the dataset
2616 54 : dataset->SetDescription(poOpenInfo->pszFilename);
2617 54 : if (hasRasterAndIsAllowed)
2618 : {
2619 43 : if (aosSubdatasets.size() > 2)
2620 : {
2621 3 : firstValidHDU = 0;
2622 3 : int hduType = 0;
2623 3 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2624 : }
2625 : else
2626 : {
2627 80 : if (firstValidHDU != 0 &&
2628 40 : dataset->Init(hFITS, true, firstValidHDU) != CE_None)
2629 : {
2630 1 : return nullptr;
2631 : }
2632 : }
2633 : }
2634 :
2635 : // If the first HDU is a dummy one, load its metadata first, and then
2636 : // add/override it by the one of the image HDU
2637 53 : if (firstHDUIsDummy && firstValidHDU > 1)
2638 : {
2639 4 : int hduType = 0;
2640 4 : status = 0;
2641 4 : fits_movabs_hdu(hFITS, 1, &hduType, &status);
2642 4 : if (status == 0)
2643 : {
2644 4 : dataset->LoadMetadata(dataset.get());
2645 : }
2646 4 : status = 0;
2647 4 : fits_movabs_hdu(hFITS, firstValidHDU, &hduType, &status);
2648 4 : if (status)
2649 : {
2650 0 : return nullptr;
2651 : }
2652 : }
2653 53 : if (hasRasterAndIsAllowed)
2654 : {
2655 42 : dataset->LoadMetadata(dataset.get());
2656 42 : dataset->LoadFITSInfo();
2657 : }
2658 :
2659 : /* -------------------------------------------------------------------- */
2660 : /* Initialize any information. */
2661 : /* -------------------------------------------------------------------- */
2662 53 : dataset->SetDescription(poOpenInfo->pszFilename);
2663 53 : dataset->TryLoadXML();
2664 :
2665 : /* -------------------------------------------------------------------- */
2666 : /* Check for external overviews. */
2667 : /* -------------------------------------------------------------------- */
2668 106 : dataset->oOvManager.Initialize(dataset.get(), poOpenInfo->pszFilename,
2669 53 : poOpenInfo->GetSiblingFiles());
2670 :
2671 53 : return dataset.release();
2672 : }
2673 :
2674 : /************************************************************************/
2675 : /* Create() */
2676 : /* */
2677 : /* Create a new FITS file. */
2678 : /************************************************************************/
2679 :
2680 67 : GDALDataset *FITSDataset::Create(const char *pszFilename, int nXSize,
2681 : int nYSize, int nBandsIn, GDALDataType eType,
2682 : CPL_UNUSED char **papszParamList)
2683 : {
2684 67 : int status = 0;
2685 :
2686 67 : if (nXSize == 0 && nYSize == 0 && nBandsIn == 0 && eType == GDT_Unknown)
2687 : {
2688 : // Create the file - to force creation, we prepend the name with '!'
2689 8 : CPLString extFilename("!");
2690 4 : extFilename += pszFilename;
2691 4 : fitsfile *hFITS = nullptr;
2692 4 : fits_create_file(&hFITS, extFilename, &status);
2693 4 : if (status)
2694 : {
2695 0 : CPLError(CE_Failure, CPLE_AppDefined,
2696 : "Couldn't create FITS file %s (%d).\n", pszFilename,
2697 : status);
2698 0 : return nullptr;
2699 : }
2700 :
2701 : // Likely vector creation mode
2702 4 : FITSDataset *dataset = new FITSDataset();
2703 4 : dataset->m_hFITS = hFITS;
2704 4 : dataset->eAccess = GA_Update;
2705 4 : dataset->SetDescription(pszFilename);
2706 4 : return dataset;
2707 : }
2708 :
2709 : // No creation options are defined. The BSCALE/BZERO options were
2710 : // removed on 2002-07-02 by Simon Perkins because they introduced
2711 : // excessive complications and didn't really fit into the GDAL
2712 : // paradigm.
2713 : // 2018 - BZERO BSCALE keywords are now set using SetScale() and
2714 : // SetOffset() functions
2715 :
2716 63 : if (nXSize < 1 || nYSize < 1 || nBandsIn < 1)
2717 : {
2718 2 : CPLError(CE_Failure, CPLE_AppDefined,
2719 : "Attempt to create %dx%dx%d raster FITS file, but width, "
2720 : "height and bands"
2721 : " must be positive.",
2722 : nXSize, nYSize, nBandsIn);
2723 :
2724 2 : return nullptr;
2725 : }
2726 :
2727 : // Determine FITS type of image
2728 : int bitpix;
2729 61 : if (eType == GDT_Byte)
2730 : {
2731 17 : bitpix = BYTE_IMG;
2732 : }
2733 44 : else if (eType == GDT_UInt16)
2734 : {
2735 4 : bitpix = USHORT_IMG;
2736 : }
2737 40 : else if (eType == GDT_Int16)
2738 : {
2739 6 : bitpix = SHORT_IMG;
2740 : }
2741 34 : else if (eType == GDT_UInt32)
2742 : {
2743 4 : bitpix = ULONG_IMG;
2744 : }
2745 30 : else if (eType == GDT_Int32)
2746 : {
2747 4 : bitpix = LONG_IMG;
2748 : }
2749 26 : else if (eType == GDT_Float32)
2750 4 : bitpix = FLOAT_IMG;
2751 22 : else if (eType == GDT_Float64)
2752 4 : bitpix = DOUBLE_IMG;
2753 : else
2754 : {
2755 18 : CPLError(CE_Failure, CPLE_AppDefined,
2756 : "GDALDataType (%d) unsupported for FITS", eType);
2757 18 : return nullptr;
2758 : }
2759 :
2760 : // Create the file - to force creation, we prepend the name with '!'
2761 86 : CPLString extFilename("!");
2762 43 : extFilename += pszFilename;
2763 43 : fitsfile *hFITS = nullptr;
2764 43 : fits_create_file(&hFITS, extFilename, &status);
2765 43 : if (status)
2766 : {
2767 3 : CPLError(CE_Failure, CPLE_AppDefined,
2768 : "Couldn't create FITS file %s (%d).\n", pszFilename, status);
2769 3 : return nullptr;
2770 : }
2771 :
2772 : // Now create an image of appropriate size and type
2773 40 : long naxes[3] = {nXSize, nYSize, nBandsIn};
2774 40 : int naxis = (nBandsIn == 1) ? 2 : 3;
2775 40 : fits_create_img(hFITS, bitpix, naxis, naxes, &status);
2776 :
2777 : // Check the status
2778 40 : if (status)
2779 : {
2780 0 : CPLError(CE_Failure, CPLE_AppDefined,
2781 : "Couldn't create image within FITS file %s (%d).", pszFilename,
2782 : status);
2783 0 : fits_close_file(hFITS, &status);
2784 0 : return nullptr;
2785 : }
2786 :
2787 40 : FITSDataset *dataset = new FITSDataset();
2788 40 : dataset->nRasterXSize = nXSize;
2789 40 : dataset->nRasterYSize = nYSize;
2790 40 : dataset->eAccess = GA_Update;
2791 40 : dataset->SetDescription(pszFilename);
2792 :
2793 : // Init recalculates a lot of stuff we already know, but...
2794 40 : if (dataset->Init(hFITS, false, 1) != CE_None)
2795 : {
2796 0 : delete dataset;
2797 0 : return nullptr;
2798 : }
2799 : else
2800 : {
2801 40 : return dataset;
2802 : }
2803 : }
2804 :
2805 : /************************************************************************/
2806 : /* Delete() */
2807 : /************************************************************************/
2808 :
2809 13 : CPLErr FITSDataset::Delete(const char *pszFilename)
2810 : {
2811 13 : return VSIUnlink(pszFilename) == 0 ? CE_None : CE_Failure;
2812 : }
2813 :
2814 : /************************************************************************/
2815 : /* WriteFITSInfo() */
2816 : /************************************************************************/
2817 :
2818 40 : void FITSDataset::WriteFITSInfo()
2819 :
2820 : {
2821 40 : int status = 0;
2822 :
2823 40 : const double PI = std::atan(1.0) * 4;
2824 40 : const double DEG2RAD = PI / 180.;
2825 :
2826 40 : double falseEast = 0;
2827 40 : double falseNorth = 0;
2828 :
2829 : double cfactor, mres, mapres, UpperLeftCornerX, UpperLeftCornerY;
2830 : double crpix1, crpix2;
2831 :
2832 : /* -------------------------------------------------------------------- */
2833 : /* Write out projection definition. */
2834 : /* -------------------------------------------------------------------- */
2835 40 : const bool bHasProjection = !m_oSRS.IsEmpty();
2836 40 : if (bHasProjection)
2837 : {
2838 :
2839 : // Set according to coordinate system (thanks to Trent Hare - USGS)
2840 :
2841 40 : std::string object, ctype1, ctype2;
2842 :
2843 40 : const char *target = m_oSRS.GetAttrValue("DATUM", 0);
2844 40 : if (target)
2845 : {
2846 40 : if (strstr(target, "Moon"))
2847 : {
2848 0 : object.assign("Moon");
2849 0 : ctype1.assign("SE");
2850 0 : ctype2.assign("SE");
2851 : }
2852 40 : else if (strstr(target, "Mercury"))
2853 : {
2854 0 : object.assign("Mercury");
2855 0 : ctype1.assign("ME");
2856 0 : ctype2.assign("ME");
2857 : }
2858 40 : else if (strstr(target, "Venus"))
2859 : {
2860 0 : object.assign("Venus");
2861 0 : ctype1.assign("VE");
2862 0 : ctype2.assign("VE");
2863 : }
2864 40 : else if (strstr(target, "Mars"))
2865 : {
2866 0 : object.assign("Mars");
2867 0 : ctype1.assign("MA");
2868 0 : ctype2.assign("MA");
2869 : }
2870 40 : else if (strstr(target, "Jupiter"))
2871 : {
2872 0 : object.assign("Jupiter");
2873 0 : ctype1.assign("JU");
2874 0 : ctype2.assign("JU");
2875 : }
2876 40 : else if (strstr(target, "Saturn"))
2877 : {
2878 0 : object.assign("Saturn");
2879 0 : ctype1.assign("SA");
2880 0 : ctype2.assign("SA");
2881 : }
2882 40 : else if (strstr(target, "Uranus"))
2883 : {
2884 0 : object.assign("Uranus");
2885 0 : ctype1.assign("UR");
2886 0 : ctype2.assign("UR");
2887 : }
2888 40 : else if (strstr(target, "Neptune"))
2889 : {
2890 0 : object.assign("Neptune");
2891 0 : ctype1.assign("NE");
2892 0 : ctype2.assign("NE");
2893 : }
2894 : else
2895 : {
2896 40 : object.assign("Earth");
2897 40 : ctype1.assign("EA");
2898 40 : ctype2.assign("EA");
2899 : }
2900 :
2901 40 : fits_update_key(
2902 : m_hFITS, TSTRING, "OBJECT",
2903 40 : const_cast<void *>(static_cast<const void *>(object.c_str())),
2904 : nullptr, &status);
2905 : }
2906 :
2907 40 : double aradius = m_oSRS.GetSemiMajor();
2908 40 : double bradius = aradius;
2909 40 : double cradius = m_oSRS.GetSemiMinor();
2910 :
2911 40 : cfactor = aradius * DEG2RAD;
2912 :
2913 40 : fits_update_key(m_hFITS, TDOUBLE, "A_RADIUS", &aradius, nullptr,
2914 : &status);
2915 40 : if (status)
2916 : {
2917 : // Throw a warning with CFITSIO error status, then ignore status
2918 0 : CPLError(CE_Warning, CPLE_AppDefined,
2919 : "Couldn't update key A_RADIUS in FITS file %s (%d).",
2920 0 : GetDescription(), status);
2921 0 : status = 0;
2922 0 : return;
2923 : }
2924 40 : fits_update_key(m_hFITS, TDOUBLE, "B_RADIUS", &bradius, nullptr,
2925 : &status);
2926 40 : if (status)
2927 : {
2928 : // Throw a warning with CFITSIO error status, then ignore status
2929 0 : CPLError(CE_Warning, CPLE_AppDefined,
2930 : "Couldn't update key B_RADIUS in FITS file %s (%d).",
2931 0 : GetDescription(), status);
2932 0 : status = 0;
2933 0 : return;
2934 : }
2935 40 : fits_update_key(m_hFITS, TDOUBLE, "C_RADIUS", &cradius, nullptr,
2936 : &status);
2937 40 : if (status)
2938 : {
2939 : // Throw a warning with CFITSIO error status, then ignore status
2940 0 : CPLError(CE_Warning, CPLE_AppDefined,
2941 : "Couldn't update key C_RADIUS in FITS file %s (%d).",
2942 0 : GetDescription(), status);
2943 0 : status = 0;
2944 0 : return;
2945 : }
2946 :
2947 40 : const char *unit = m_oSRS.GetAttrValue("UNIT", 0);
2948 :
2949 40 : ctype1.append("LN-");
2950 40 : ctype2.append("LT-");
2951 :
2952 : // strcat(ctype1a, "PX-");
2953 : // strcat(ctype2a, "PY-");
2954 :
2955 40 : std::string fitsproj;
2956 40 : const char *projection = m_oSRS.GetAttrValue("PROJECTION", 0);
2957 40 : double centlon = 0, centlat = 0;
2958 :
2959 40 : if (projection)
2960 : {
2961 12 : if (strstr(projection, "Sinusoidal"))
2962 : {
2963 0 : fitsproj.assign("SFL");
2964 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2965 : }
2966 12 : else if (strstr(projection, "Equirectangular"))
2967 : {
2968 0 : fitsproj.assign("CAR");
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, "Orthographic"))
2973 : {
2974 0 : fitsproj.assign("SIN");
2975 0 : centlat = m_oSRS.GetProjParm("standard_parallel_1", 0, nullptr);
2976 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2977 : }
2978 12 : else if (strstr(projection, "Mercator_1SP") ||
2979 12 : strstr(projection, "Mercator"))
2980 : {
2981 12 : fitsproj.assign("MER");
2982 12 : centlat = m_oSRS.GetProjParm("standard_parallel_1", 0, nullptr);
2983 12 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2984 : }
2985 0 : else if (strstr(projection, "Polar_Stereographic") ||
2986 0 : strstr(projection, "Stereographic_South_Pole") ||
2987 0 : strstr(projection, "Stereographic_North_Pole"))
2988 : {
2989 0 : fitsproj.assign("STG");
2990 0 : centlat = m_oSRS.GetProjParm("latitude_of_origin", 0, nullptr);
2991 0 : centlon = m_oSRS.GetProjParm("central_meridian", 0, nullptr);
2992 : }
2993 :
2994 : /*
2995 : #Transverse Mercator is supported in FITS via
2996 : specific MER parameters. # need some more testing... #if
2997 : EQUAL(mapProjection,"Transverse_Mercator"): # mapProjection =
2998 : "MER" # centLat = hSRS.GetProjParm('standard_parallel_1') #
2999 : centLon = hSRS.GetProjParm('central_meridian') # TMscale =
3000 : hSRS.GetProjParm('scale_factor') # #Need to research when TM
3001 : actually applies false values # #but planetary is almost
3002 : always 0.0 # falseEast = hSRS.GetProjParm('false_easting') #
3003 : falseNorth = hSRS.GetProjParm('false_northing')
3004 : */
3005 :
3006 12 : ctype1.append(fitsproj);
3007 12 : ctype2.append(fitsproj);
3008 :
3009 12 : fits_update_key(
3010 : m_hFITS, TSTRING, "CTYPE1",
3011 12 : const_cast<void *>(static_cast<const void *>(ctype1.c_str())),
3012 : nullptr, &status);
3013 12 : if (status)
3014 : {
3015 : // Throw a warning with CFITSIO error status, then ignore status
3016 0 : CPLError(CE_Warning, CPLE_AppDefined,
3017 : "Couldn't update key CTYPE1 in FITS file %s (%d).",
3018 0 : GetDescription(), status);
3019 0 : status = 0;
3020 0 : return;
3021 : }
3022 :
3023 12 : fits_update_key(
3024 : m_hFITS, TSTRING, "CTYPE2",
3025 12 : const_cast<void *>(static_cast<const void *>(ctype2.c_str())),
3026 : nullptr, &status);
3027 12 : if (status)
3028 : {
3029 : // Throw a warning with CFITSIO error status, then ignore status
3030 0 : CPLError(CE_Warning, CPLE_AppDefined,
3031 : "Couldn't update key CTYPE2 in FITS file %s (%d).",
3032 0 : GetDescription(), status);
3033 0 : status = 0;
3034 0 : return;
3035 : }
3036 : }
3037 :
3038 40 : UpperLeftCornerX = m_adfGeoTransform[0] - falseEast;
3039 40 : UpperLeftCornerY = m_adfGeoTransform[3] - falseNorth;
3040 :
3041 40 : if (centlon > 180.)
3042 : {
3043 0 : centlon = centlon - 180.;
3044 : }
3045 40 : if (strstr(unit, "metre"))
3046 : {
3047 : // convert degrees/pixel to m/pixel
3048 12 : mapres = 1. / m_adfGeoTransform[1]; // mapres is pixel/meters
3049 12 : mres = m_adfGeoTransform[1] / cfactor; // mres is deg/pixel
3050 12 : crpix1 = -(UpperLeftCornerX * mapres) + centlon / mres + 0.5;
3051 : // assuming that center latitude is also the origin of the
3052 : // coordinate system: this is not always true. More generic
3053 : // implementation coming soon
3054 12 : crpix2 = (UpperLeftCornerY * mapres) + 0.5; // - (centlat / mres);
3055 : }
3056 28 : else if (strstr(unit, "degree"))
3057 : {
3058 : // convert m/pixel to pixel/degree
3059 28 : mapres =
3060 28 : 1. / m_adfGeoTransform[1] / cfactor; // mapres is pixel/deg
3061 28 : mres = m_adfGeoTransform[1]; // mres is meters/pixel
3062 28 : crpix1 = -(UpperLeftCornerX * mres) + centlon / mapres + 0.5;
3063 : // assuming that center latitude is also the origin of the
3064 : // coordinate system: this is not always true. More generic
3065 : // implementation coming soon
3066 28 : crpix2 = (UpperLeftCornerY * mres) + 0.5; // - (centlat / mapres);
3067 : }
3068 :
3069 : /// Write WCS CRPIXia CRVALia CTYPEia here
3070 :
3071 40 : fits_update_key(m_hFITS, TDOUBLE, "CRVAL1", ¢lon, nullptr, &status);
3072 40 : if (status)
3073 : {
3074 : // Throw a warning with CFITSIO error status, then ignore status
3075 0 : CPLError(CE_Warning, CPLE_AppDefined,
3076 : "Couldn't update key CRVAL1 in FITS file %s (%d).",
3077 0 : GetDescription(), status);
3078 0 : status = 0;
3079 0 : return;
3080 : }
3081 40 : fits_update_key(m_hFITS, TDOUBLE, "CRVAL2", ¢lat, nullptr, &status);
3082 40 : if (status)
3083 : {
3084 : // Throw a warning with CFITSIO error status, then ignore status
3085 0 : CPLError(CE_Warning, CPLE_AppDefined,
3086 : "Couldn't update key CRVAL2 in FITS file %s (%d).",
3087 0 : GetDescription(), status);
3088 0 : status = 0;
3089 0 : return;
3090 : }
3091 40 : fits_update_key(m_hFITS, TDOUBLE, "CRPIX1", &crpix1, nullptr, &status);
3092 40 : if (status)
3093 : {
3094 : // Throw a warning with CFITSIO error status, then ignore status
3095 0 : CPLError(CE_Warning, CPLE_AppDefined,
3096 : "Couldn't update key CRPIX1 in FITS file %s (%d).",
3097 0 : GetDescription(), status);
3098 0 : status = 0;
3099 0 : return;
3100 : }
3101 40 : fits_update_key(m_hFITS, TDOUBLE, "CRPIX2", &crpix2, nullptr, &status);
3102 40 : if (status)
3103 : {
3104 : // Throw a warning with CFITSIO error status, then ignore status
3105 0 : CPLError(CE_Warning, CPLE_AppDefined,
3106 : "Couldn't update key CRPIX2 in FITS file %s (%d).",
3107 0 : GetDescription(), status);
3108 0 : status = 0;
3109 0 : return;
3110 : }
3111 :
3112 : /* --------------------------------------------------------------------
3113 : */
3114 : /* Write geotransform if valid. */
3115 : /* --------------------------------------------------------------------
3116 : */
3117 40 : if (m_bGeoTransformValid)
3118 : {
3119 :
3120 : /* --------------------------------------------------------------------
3121 : */
3122 : /* Write the transform. */
3123 : /* --------------------------------------------------------------------
3124 : */
3125 :
3126 : /// Write WCS CDELTia and PCi_ja here
3127 :
3128 : double cd[4];
3129 40 : cd[0] = m_adfGeoTransform[1] / cfactor;
3130 40 : cd[1] = m_adfGeoTransform[2] / cfactor;
3131 40 : cd[2] = m_adfGeoTransform[4] / cfactor;
3132 40 : cd[3] = m_adfGeoTransform[5] / cfactor;
3133 :
3134 : double pc[4];
3135 40 : pc[0] = 1.;
3136 40 : pc[1] = cd[1] / cd[0];
3137 40 : pc[2] = cd[2] / cd[3];
3138 40 : pc[3] = -1.;
3139 :
3140 40 : fits_update_key(m_hFITS, TDOUBLE, "CDELT1", &cd[0], nullptr,
3141 : &status);
3142 40 : if (status)
3143 : {
3144 : // Throw a warning with CFITSIO error status, then ignore status
3145 0 : CPLError(CE_Warning, CPLE_AppDefined,
3146 : "Couldn't update key CDELT1 in FITS file %s (%d).",
3147 0 : GetDescription(), status);
3148 0 : status = 0;
3149 0 : return;
3150 : }
3151 :
3152 40 : fits_update_key(m_hFITS, TDOUBLE, "CDELT2", &cd[3], nullptr,
3153 : &status);
3154 40 : if (status)
3155 : {
3156 : // Throw a warning with CFITSIO error status, then ignore status
3157 0 : CPLError(CE_Warning, CPLE_AppDefined,
3158 : "Couldn't update key CDELT2 in FITS file %s (%d).",
3159 0 : GetDescription(), status);
3160 0 : status = 0;
3161 0 : return;
3162 : }
3163 :
3164 40 : fits_update_key(m_hFITS, TDOUBLE, "PC1_1", &pc[0], nullptr,
3165 : &status);
3166 40 : if (status)
3167 : {
3168 : // Throw a warning with CFITSIO error status, then ignore status
3169 0 : CPLError(CE_Warning, CPLE_AppDefined,
3170 : "Couldn't update key PC1_1 in FITS file %s (%d).",
3171 0 : GetDescription(), status);
3172 0 : status = 0;
3173 0 : return;
3174 : }
3175 :
3176 40 : fits_update_key(m_hFITS, TDOUBLE, "PC1_2", &pc[1], nullptr,
3177 : &status);
3178 40 : if (status)
3179 : {
3180 : // Throw a warning with CFITSIO error status, then ignore status
3181 0 : CPLError(CE_Warning, CPLE_AppDefined,
3182 : "Couldn't update key PC1_2 in FITS file %s (%d).",
3183 0 : GetDescription(), status);
3184 0 : status = 0;
3185 0 : return;
3186 : }
3187 :
3188 40 : fits_update_key(m_hFITS, TDOUBLE, "PC2_1", &pc[2], nullptr,
3189 : &status);
3190 40 : if (status)
3191 : {
3192 : // Throw a warning with CFITSIO error status, then ignore status
3193 0 : CPLError(CE_Warning, CPLE_AppDefined,
3194 : "Couldn't update key PC2_1 in FITS file %s (%d).",
3195 0 : GetDescription(), status);
3196 0 : status = 0;
3197 0 : return;
3198 : }
3199 :
3200 40 : fits_update_key(m_hFITS, TDOUBLE, "PC2_2", &pc[3], nullptr,
3201 : &status);
3202 40 : if (status)
3203 : {
3204 : // Throw a warning with CFITSIO error status, then ignore status
3205 0 : CPLError(CE_Warning, CPLE_AppDefined,
3206 : "Couldn't update key PC2_2 in FITS file %s (%d).",
3207 0 : GetDescription(), status);
3208 0 : status = 0;
3209 0 : return;
3210 : }
3211 : }
3212 : }
3213 : }
3214 :
3215 : /************************************************************************/
3216 : /* GetSpatialRef() */
3217 : /************************************************************************/
3218 :
3219 4 : const OGRSpatialReference *FITSDataset::GetSpatialRef() const
3220 :
3221 : {
3222 4 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
3223 : }
3224 :
3225 : /************************************************************************/
3226 : /* SetSpatialRef() */
3227 : /************************************************************************/
3228 :
3229 40 : CPLErr FITSDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
3230 :
3231 : {
3232 40 : if (poSRS == nullptr || poSRS->IsEmpty())
3233 : {
3234 0 : m_oSRS.Clear();
3235 : }
3236 : else
3237 : {
3238 40 : m_oSRS = *poSRS;
3239 40 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3240 : }
3241 :
3242 40 : m_bFITSInfoChanged = true;
3243 :
3244 40 : return CE_None;
3245 : }
3246 :
3247 : /************************************************************************/
3248 : /* GetGeoTransform() */
3249 : /************************************************************************/
3250 :
3251 19 : CPLErr FITSDataset::GetGeoTransform(double *padfTransform)
3252 :
3253 : {
3254 19 : memcpy(padfTransform, m_adfGeoTransform, sizeof(double) * 6);
3255 :
3256 19 : if (!m_bGeoTransformValid)
3257 0 : return CE_Failure;
3258 :
3259 19 : return CE_None;
3260 : }
3261 :
3262 : /************************************************************************/
3263 : /* SetGeoTransform() */
3264 : /************************************************************************/
3265 :
3266 40 : CPLErr FITSDataset::SetGeoTransform(double *padfTransform)
3267 :
3268 : {
3269 40 : memcpy(m_adfGeoTransform, padfTransform, sizeof(double) * 6);
3270 40 : m_bGeoTransformValid = true;
3271 :
3272 40 : return CE_None;
3273 : }
3274 :
3275 : /************************************************************************/
3276 : /* GetOffset() */
3277 : /************************************************************************/
3278 :
3279 44 : double FITSRasterBand::GetOffset(int *pbSuccess)
3280 :
3281 : {
3282 44 : if (pbSuccess)
3283 43 : *pbSuccess = m_bHaveOffsetScale;
3284 44 : return m_dfOffset;
3285 : }
3286 :
3287 : /************************************************************************/
3288 : /* SetOffset() */
3289 : /************************************************************************/
3290 :
3291 1 : CPLErr FITSRasterBand::SetOffset(double dfNewValue)
3292 :
3293 : {
3294 1 : if (!m_bHaveOffsetScale || dfNewValue != m_dfOffset)
3295 1 : m_poFDS->m_bMetadataChanged = true;
3296 :
3297 1 : m_bHaveOffsetScale = true;
3298 1 : m_dfOffset = dfNewValue;
3299 1 : return CE_None;
3300 : }
3301 :
3302 : /************************************************************************/
3303 : /* GetScale() */
3304 : /************************************************************************/
3305 :
3306 44 : double FITSRasterBand::GetScale(int *pbSuccess)
3307 :
3308 : {
3309 44 : if (pbSuccess)
3310 43 : *pbSuccess = m_bHaveOffsetScale;
3311 44 : return m_dfScale;
3312 : }
3313 :
3314 : /************************************************************************/
3315 : /* SetScale() */
3316 : /************************************************************************/
3317 :
3318 1 : CPLErr FITSRasterBand::SetScale(double dfNewValue)
3319 :
3320 : {
3321 1 : if (!m_bHaveOffsetScale || dfNewValue != m_dfScale)
3322 1 : m_poFDS->m_bMetadataChanged = true;
3323 :
3324 1 : m_bHaveOffsetScale = true;
3325 1 : m_dfScale = dfNewValue;
3326 1 : return CE_None;
3327 : }
3328 :
3329 : /************************************************************************/
3330 : /* GetNoDataValue() */
3331 : /************************************************************************/
3332 :
3333 2 : double FITSRasterBand::GetNoDataValue(int *pbSuccess)
3334 :
3335 : {
3336 2 : if (m_bNoDataSet)
3337 : {
3338 0 : if (pbSuccess)
3339 0 : *pbSuccess = TRUE;
3340 :
3341 0 : return m_dfNoDataValue;
3342 : }
3343 :
3344 2 : if (m_poFDS->m_bNoDataSet)
3345 : {
3346 2 : if (pbSuccess)
3347 2 : *pbSuccess = TRUE;
3348 :
3349 2 : return m_poFDS->m_dfNoDataValue;
3350 : }
3351 :
3352 0 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
3353 : }
3354 :
3355 : /************************************************************************/
3356 : /* SetNoDataValue() */
3357 : /************************************************************************/
3358 :
3359 1 : CPLErr FITSRasterBand::SetNoDataValue(double dfNoData)
3360 :
3361 : {
3362 1 : if (m_poFDS->m_bNoDataSet && m_poFDS->m_dfNoDataValue == dfNoData)
3363 : {
3364 0 : m_bNoDataSet = true;
3365 0 : m_dfNoDataValue = dfNoData;
3366 0 : return CE_None;
3367 : }
3368 :
3369 1 : m_poFDS->m_bNoDataSet = true;
3370 1 : m_poFDS->m_dfNoDataValue = dfNoData;
3371 :
3372 1 : m_poFDS->m_bNoDataChanged = true;
3373 :
3374 1 : m_bNoDataSet = true;
3375 1 : m_dfNoDataValue = dfNoData;
3376 1 : return CE_None;
3377 : }
3378 :
3379 : /************************************************************************/
3380 : /* DeleteNoDataValue() */
3381 : /************************************************************************/
3382 :
3383 0 : CPLErr FITSRasterBand::DeleteNoDataValue()
3384 :
3385 : {
3386 0 : if (!m_poFDS->m_bNoDataSet)
3387 0 : return CE_None;
3388 :
3389 0 : m_poFDS->m_bNoDataSet = false;
3390 0 : m_poFDS->m_dfNoDataValue = -9999.0;
3391 :
3392 0 : m_poFDS->m_bNoDataChanged = true;
3393 :
3394 0 : m_bNoDataSet = false;
3395 0 : m_dfNoDataValue = -9999.0;
3396 0 : return CE_None;
3397 : }
3398 :
3399 : /************************************************************************/
3400 : /* LoadGeoreferencing() */
3401 : /************************************************************************/
3402 :
3403 42 : void FITSDataset::LoadGeoreferencing()
3404 : {
3405 42 : int status = 0;
3406 : double crpix1, crpix2, crval1, crval2, cdelt1, cdelt2, pc[4], cd[4];
3407 42 : double aRadius, cRadius, invFlattening = 0.0;
3408 42 : double falseEast = 0.0, falseNorth = 0.0, scale = 1.0;
3409 : char target[81], ctype[81];
3410 42 : std::string GeogName, DatumName, projName;
3411 :
3412 42 : const double PI = std::atan(1.0) * 4;
3413 42 : const double DEG2RAD = PI / 180.;
3414 :
3415 : /* -------------------------------------------------------------------- */
3416 : /* Get the transform from the FITS file. */
3417 : /* -------------------------------------------------------------------- */
3418 :
3419 42 : fits_read_key(m_hFITS, TSTRING, "OBJECT", target, nullptr, &status);
3420 42 : if (status)
3421 : {
3422 9 : strncpy(target, "Undefined", 10);
3423 9 : CPLDebug("FITS", "OBJECT keyword is missing");
3424 9 : status = 0;
3425 : }
3426 :
3427 42 : GeogName.assign("GCS_");
3428 42 : GeogName.append(target);
3429 42 : DatumName.assign("D_");
3430 42 : DatumName.append(target);
3431 :
3432 42 : fits_read_key(m_hFITS, TDOUBLE, "A_RADIUS", &aRadius, nullptr, &status);
3433 42 : if (status)
3434 : {
3435 9 : CPLDebug("FITS", "No Radii keyword available, metadata will not "
3436 : "contain DATUM information.");
3437 9 : return;
3438 : }
3439 : else
3440 : {
3441 33 : fits_read_key(m_hFITS, TDOUBLE, "C_RADIUS", &cRadius, nullptr, &status);
3442 33 : if (status)
3443 : {
3444 0 : CPLError(CE_Warning, CPLE_AppDefined,
3445 : "No polar radius keyword available, setting C_RADIUS = "
3446 : "A_RADIUS");
3447 0 : cRadius = aRadius;
3448 0 : status = 0;
3449 : }
3450 33 : if (aRadius != cRadius)
3451 : {
3452 33 : invFlattening = aRadius / (aRadius - cRadius);
3453 : }
3454 : }
3455 :
3456 : /* Waiting for linear keywords standardization only deg ctype are used */
3457 : /* Check if WCS are there */
3458 33 : fits_read_key(m_hFITS, TSTRING, "CTYPE1", ctype, nullptr, &status);
3459 33 : if (!status)
3460 : {
3461 : /* Check if angular WCS are there */
3462 16 : if (strstr(ctype, "LN"))
3463 : {
3464 : /* Reading reference points */
3465 16 : fits_read_key(m_hFITS, TDOUBLE, "CRPIX1", &crpix1, nullptr,
3466 : &status);
3467 16 : fits_read_key(m_hFITS, TDOUBLE, "CRPIX2", &crpix2, nullptr,
3468 : &status);
3469 16 : fits_read_key(m_hFITS, TDOUBLE, "CRVAL1", &crval1, nullptr,
3470 : &status);
3471 16 : fits_read_key(m_hFITS, TDOUBLE, "CRVAL2", &crval2, nullptr,
3472 : &status);
3473 16 : if (status)
3474 : {
3475 0 : CPLError(CE_Failure, CPLE_AppDefined,
3476 : "No CRPIX / CRVAL keyword available, the raster "
3477 : "cannot be georeferenced.");
3478 0 : status = 0;
3479 : }
3480 : else
3481 : {
3482 : /* Check for CDELT and PC matrix representation */
3483 16 : fits_read_key(m_hFITS, TDOUBLE, "CDELT1", &cdelt1, nullptr,
3484 : &status);
3485 16 : if (!status)
3486 : {
3487 16 : fits_read_key(m_hFITS, TDOUBLE, "CDELT2", &cdelt2, nullptr,
3488 : &status);
3489 16 : fits_read_key(m_hFITS, TDOUBLE, "PC1_1", &pc[0], nullptr,
3490 : &status);
3491 16 : fits_read_key(m_hFITS, TDOUBLE, "PC1_2", &pc[1], nullptr,
3492 : &status);
3493 16 : fits_read_key(m_hFITS, TDOUBLE, "PC2_1", &pc[2], nullptr,
3494 : &status);
3495 16 : fits_read_key(m_hFITS, TDOUBLE, "PC2_2", &pc[3], nullptr,
3496 : &status);
3497 16 : cd[0] = cdelt1 * pc[0];
3498 16 : cd[1] = cdelt1 * pc[1];
3499 16 : cd[2] = cdelt2 * pc[2];
3500 16 : cd[3] = cdelt2 * pc[3];
3501 16 : status = 0;
3502 : }
3503 : else
3504 : {
3505 : /* Look for CD matrix representation */
3506 0 : fits_read_key(m_hFITS, TDOUBLE, "CD1_1", &cd[0], nullptr,
3507 : &status);
3508 0 : fits_read_key(m_hFITS, TDOUBLE, "CD1_2", &cd[1], nullptr,
3509 : &status);
3510 0 : fits_read_key(m_hFITS, TDOUBLE, "CD2_1", &cd[2], nullptr,
3511 : &status);
3512 0 : fits_read_key(m_hFITS, TDOUBLE, "CD2_2", &cd[3], nullptr,
3513 : &status);
3514 : }
3515 :
3516 16 : double radfac = DEG2RAD * aRadius;
3517 :
3518 16 : m_adfGeoTransform[1] = cd[0] * radfac;
3519 16 : m_adfGeoTransform[2] = cd[1] * radfac;
3520 16 : m_adfGeoTransform[4] = cd[2] * radfac;
3521 16 : m_adfGeoTransform[5] = -cd[3] * radfac;
3522 16 : if (crval1 > 180.)
3523 : {
3524 0 : crval1 = crval1 - 180.;
3525 : }
3526 :
3527 : /* NOTA BENE: FITS standard define pixel integers at the center
3528 : of the pixel, 0.5 must be subtract to have UpperLeft corner
3529 : */
3530 16 : m_adfGeoTransform[0] =
3531 16 : crval1 * radfac - m_adfGeoTransform[1] * (crpix1 - 0.5);
3532 : // assuming that center latitude is also the origin of the
3533 : // coordinate system: this is not always true. More generic
3534 : // implementation coming soon
3535 16 : m_adfGeoTransform[3] = -m_adfGeoTransform[5] * (crpix2 - 0.5);
3536 : //+ crval2 * radfac;
3537 16 : m_bGeoTransformValid = true;
3538 : }
3539 :
3540 16 : char *pstr = strrchr(ctype, '-');
3541 16 : if (pstr)
3542 : {
3543 16 : pstr += 1;
3544 :
3545 : /* Defining projection type
3546 : Following http://www.gdal.org/ogr__srs__api_8h.html (GDAL)
3547 : and
3548 : http://www.aanda.org/component/article?access=bibcode&bibcode=&bibcode=2002A%2526A...395.1077CFUL
3549 : (FITS)
3550 : */
3551 :
3552 : /* Sinusoidal / SFL projection */
3553 16 : if (strcmp(pstr, "SFL") == 0)
3554 : {
3555 0 : projName.assign("Sinusoidal_");
3556 0 : m_oSRS.SetSinusoidal(crval1, falseEast, falseNorth);
3557 :
3558 : /* Mercator, Oblique (Hotine) Mercator, Transverse Mercator
3559 : */
3560 : /* Mercator / MER projection */
3561 : }
3562 16 : else if (strcmp(pstr, "MER") == 0)
3563 : {
3564 16 : projName.assign("Mercator_");
3565 16 : m_oSRS.SetMercator(crval2, crval1, scale, falseEast,
3566 : falseNorth);
3567 :
3568 : /* Equirectangular / CAR projection */
3569 : }
3570 0 : else if (strcmp(pstr, "CAR") == 0)
3571 : {
3572 0 : projName.assign("Equirectangular_");
3573 : /*
3574 : The standard_parallel_1 defines where the local radius is
3575 : calculated not the center of Y Cartesian system (which is
3576 : latitude_of_origin) But FITS WCS only supports projections
3577 : on the sphere we assume here that the local radius is the
3578 : one computed at the projection center
3579 : */
3580 0 : m_oSRS.SetEquirectangular2(crval2, crval1, crval2,
3581 : falseEast, falseNorth);
3582 : /* Lambert Azimuthal Equal Area / ZEA projection */
3583 : }
3584 0 : else if (strcmp(pstr, "ZEA") == 0)
3585 : {
3586 0 : projName.assign("Lambert_Azimuthal_Equal_Area_");
3587 0 : m_oSRS.SetLAEA(crval2, crval1, falseEast, falseNorth);
3588 :
3589 : /* Lambert Conformal Conic 1SP / COO projection */
3590 : }
3591 0 : else if (strcmp(pstr, "COO") == 0)
3592 : {
3593 0 : projName.assign("Lambert_Conformal_Conic_1SP_");
3594 0 : m_oSRS.SetLCC1SP(crval2, crval1, scale, falseEast,
3595 : falseNorth);
3596 :
3597 : /* Orthographic / SIN projection */
3598 : }
3599 0 : else if (strcmp(pstr, "SIN") == 0)
3600 : {
3601 0 : projName.assign("Orthographic_");
3602 0 : m_oSRS.SetOrthographic(crval2, crval1, falseEast,
3603 : falseNorth);
3604 :
3605 : /* Point Perspective / AZP projection */
3606 : }
3607 0 : else if (strcmp(pstr, "AZP") == 0)
3608 : {
3609 0 : projName.assign("perspective_point_height_");
3610 0 : m_oSRS.SetProjection(SRS_PP_PERSPECTIVE_POINT_HEIGHT);
3611 : /* # appears to need height... maybe center lon/lat */
3612 :
3613 : /* Polar Stereographic / STG projection */
3614 : }
3615 0 : else if (strcmp(pstr, "STG") == 0)
3616 : {
3617 0 : projName.assign("Polar_Stereographic_");
3618 0 : m_oSRS.SetStereographic(crval2, crval1, scale, falseEast,
3619 : falseNorth);
3620 : }
3621 : else
3622 : {
3623 0 : CPLError(CE_Failure, CPLE_AppDefined,
3624 : "Unknown projection.");
3625 : }
3626 :
3627 16 : projName.append(target);
3628 16 : m_oSRS.SetProjParm(SRS_PP_FALSE_EASTING, 0.0);
3629 16 : m_oSRS.SetProjParm(SRS_PP_FALSE_NORTHING, 0.0);
3630 :
3631 16 : m_oSRS.SetNode("PROJCS", projName.c_str());
3632 :
3633 16 : m_oSRS.SetGeogCS(GeogName.c_str(), DatumName.c_str(), target,
3634 : aRadius, invFlattening, "Reference_Meridian",
3635 : 0.0, "degree", 0.0174532925199433);
3636 : }
3637 : else
3638 : {
3639 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown projection.");
3640 : }
3641 : }
3642 : }
3643 : else
3644 : {
3645 17 : CPLError(CE_Warning, CPLE_AppDefined,
3646 : "No CTYPE keywords: no geospatial information available.");
3647 : }
3648 : }
3649 :
3650 : /************************************************************************/
3651 : /* LoadFITSInfo() */
3652 : /************************************************************************/
3653 :
3654 42 : void FITSDataset::LoadFITSInfo()
3655 :
3656 : {
3657 42 : int status = 0;
3658 : int bitpix;
3659 : double dfScale, dfOffset;
3660 :
3661 42 : LoadGeoreferencing();
3662 :
3663 42 : CPLAssert(!m_bMetadataChanged);
3664 42 : CPLAssert(!m_bNoDataChanged);
3665 :
3666 42 : m_bMetadataChanged = false;
3667 42 : m_bNoDataChanged = false;
3668 :
3669 42 : bitpix = this->m_fitsDataType;
3670 42 : FITSRasterBand *poBand = cpl::down_cast<FITSRasterBand *>(GetRasterBand(1));
3671 :
3672 42 : if (bitpix != TUSHORT && bitpix != TUINT)
3673 : {
3674 36 : fits_read_key(m_hFITS, TDOUBLE, "BSCALE", &dfScale, nullptr, &status);
3675 36 : if (status)
3676 : {
3677 34 : status = 0;
3678 34 : dfScale = 1.;
3679 : }
3680 36 : fits_read_key(m_hFITS, TDOUBLE, "BZERO", &dfOffset, nullptr, &status);
3681 36 : if (status)
3682 : {
3683 34 : status = 0;
3684 34 : dfOffset = 0.;
3685 : }
3686 36 : if (dfScale != 1. || dfOffset != 0.)
3687 : {
3688 2 : poBand->m_bHaveOffsetScale = true;
3689 2 : poBand->m_dfScale = dfScale;
3690 2 : poBand->m_dfOffset = dfOffset;
3691 : }
3692 : }
3693 :
3694 42 : fits_read_key(m_hFITS, TDOUBLE, "BLANK", &m_dfNoDataValue, nullptr,
3695 : &status);
3696 42 : m_bNoDataSet = !status;
3697 42 : }
3698 :
3699 : /************************************************************************/
3700 : /* GDALRegister_FITS() */
3701 : /************************************************************************/
3702 :
3703 12 : void GDALRegister_FITS()
3704 :
3705 : {
3706 12 : if (GDALGetDriverByName(DRIVER_NAME) != nullptr)
3707 0 : return;
3708 :
3709 12 : GDALDriver *poDriver = new GDALDriver();
3710 12 : FITSDriverSetCommonMetadata(poDriver);
3711 :
3712 12 : poDriver->pfnOpen = FITSDataset::Open;
3713 12 : poDriver->pfnCreate = FITSDataset::Create;
3714 12 : poDriver->pfnDelete = FITSDataset::Delete;
3715 :
3716 12 : GetGDALDriverManager()->RegisterDriver(poDriver);
3717 : }
|
