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