1 : /******************************************************************************
2 : * $Id: fitsdataset.cpp 15409 2008-09-22 19:08:15Z rouault $
3 : *
4 : * Project: FITS Driver
5 : * Purpose: Implement FITS raster read/write support
6 : * Author: Simon Perkins, s.perkins@lanl.gov
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, Simon Perkins
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 :
31 :
32 : #include "gdal_pam.h"
33 : #include "cpl_string.h"
34 : #include <string.h>
35 :
36 : CPL_CVSID("$Id: fitsdataset.cpp 15409 2008-09-22 19:08:15Z rouault $");
37 :
38 : CPL_C_START
39 : #include <fitsio.h>
40 : void GDALRegister_FITS(void);
41 : CPL_C_END
42 :
43 : /************************************************************************/
44 : /* ==================================================================== */
45 : /* FITSDataset */
46 : /* ==================================================================== */
47 : /************************************************************************/
48 :
49 : class FITSRasterBand;
50 :
51 : class FITSDataset : public GDALPamDataset {
52 :
53 : friend class FITSRasterBand;
54 :
55 : fitsfile* hFITS;
56 :
57 : GDALDataType gdalDataType; // GDAL code for the image type
58 : int fitsDataType; // FITS code for the image type
59 :
60 : bool isExistingFile;
61 : long highestOffsetWritten; // How much of image has been written
62 :
63 : FITSDataset(); // Others shouldn't call this constructor explicitly
64 :
65 : CPLErr Init(fitsfile* hFITS_, bool isExistingFile_);
66 :
67 : public:
68 : ~FITSDataset();
69 :
70 : static GDALDataset* Open(GDALOpenInfo* );
71 : static GDALDataset* Create(const char* pszFilename,
72 : int nXSize, int nYSize, int nBands,
73 : GDALDataType eType,
74 : char** papszParmList);
75 :
76 : };
77 :
78 : /************************************************************************/
79 : /* ==================================================================== */
80 : /* FITSRasterBand */
81 : /* ==================================================================== */
82 : /************************************************************************/
83 :
84 : class FITSRasterBand : public GDALPamRasterBand {
85 :
86 : friend class FITSDataset;
87 :
88 : public:
89 :
90 : FITSRasterBand(FITSDataset*, int);
91 : ~FITSRasterBand();
92 :
93 : virtual CPLErr IReadBlock( int, int, void * );
94 : virtual CPLErr IWriteBlock( int, int, void * );
95 : };
96 :
97 :
98 : /************************************************************************/
99 : /* FITSRasterBand() */
100 : /************************************************************************/
101 :
102 58 : FITSRasterBand::FITSRasterBand(FITSDataset *poDS, int nBand) {
103 :
104 58 : this->poDS = poDS;
105 58 : this->nBand = nBand;
106 58 : eDataType = poDS->gdalDataType;
107 58 : nBlockXSize = poDS->nRasterXSize;;
108 58 : nBlockYSize = 1;
109 58 : }
110 :
111 : /************************************************************************/
112 : /* ~FITSRasterBand() */
113 : /************************************************************************/
114 :
115 116 : FITSRasterBand::~FITSRasterBand() {
116 58 : FlushCache();
117 116 : }
118 :
119 : /************************************************************************/
120 : /* IReadBlock() */
121 : /************************************************************************/
122 :
123 100 : CPLErr FITSRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
124 : void* pImage ) {
125 :
126 : // A FITS block is one row (we assume BSQ formatted data)
127 100 : FITSDataset* dataset = (FITSDataset*) poDS;
128 100 : fitsfile* hFITS = dataset->hFITS;
129 100 : int status = 0;
130 :
131 : // Since a FITS block is a whole row, nBlockXOff must be zero
132 : // and the row number equals the y block offset. Also, nBlockYOff
133 : // cannot be greater than the number of rows
134 : CPLAssert(nBlockXOff == 0);
135 : CPLAssert(nBlockYOff < nRasterYSize);
136 :
137 : // Calculate offsets and read in the data. Note that FITS array offsets
138 : // start at 1...
139 : long offset = (nBand - 1) * nRasterXSize * nRasterYSize +
140 100 : nBlockYOff * nRasterXSize + 1;
141 100 : long nElements = nRasterXSize;
142 :
143 : // If we haven't written this block to the file yet, then attempting
144 : // to read causes an error, so in this case, just return zeros.
145 100 : if (!dataset->isExistingFile && offset > dataset->highestOffsetWritten) {
146 : memset(pImage, 0, nBlockXSize * nBlockYSize
147 0 : * GDALGetDataTypeSize(eDataType) / 8);
148 0 : return CE_None;
149 : }
150 :
151 : // Otherwise read in the image data
152 : fits_read_img(hFITS, dataset->fitsDataType, offset, nElements,
153 100 : 0, pImage, 0, &status);
154 100 : if (status) {
155 : CPLError(CE_Failure, CPLE_AppDefined,
156 0 : "Couldn't read image data from FITS file (%d).", status);
157 0 : return CE_Failure;
158 : }
159 :
160 100 : return CE_None;
161 : }
162 :
163 : /************************************************************************/
164 : /* IWriteBlock() */
165 : /* */
166 : /************************************************************************/
167 :
168 310 : CPLErr FITSRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
169 : void* pImage) {
170 :
171 310 : FITSDataset* dataset = (FITSDataset*) poDS;
172 310 : fitsfile* hFITS = dataset->hFITS;
173 310 : int status = 0;
174 :
175 : // Since a FITS block is a whole row, nBlockXOff must be zero
176 : // and the row number equals the y block offset. Also, nBlockYOff
177 : // cannot be greater than the number of rows
178 :
179 : // Calculate offsets and read in the data. Note that FITS array offsets
180 : // start at 1...
181 : long offset = (nBand - 1) * nRasterXSize * nRasterYSize +
182 310 : nBlockYOff * nRasterXSize + 1;
183 310 : long nElements = nRasterXSize;
184 : fits_write_img(hFITS, dataset->fitsDataType, offset, nElements,
185 310 : pImage, &status);
186 :
187 : // Capture special case of non-zero status due to data range
188 : // overflow Standard GDAL policy is to silently truncate, which is
189 : // what CFITSIO does, in addition to returning NUM_OVERFLOW (412) as
190 : // the status.
191 310 : if (status == NUM_OVERFLOW)
192 0 : status = 0;
193 :
194 : // Check for other errors
195 310 : if (status) {
196 : CPLError(CE_Failure, CPLE_AppDefined,
197 0 : "Error writing image data to FITS file (%d).", status);
198 0 : return CE_Failure;
199 : }
200 :
201 : // When we write a block, update the offset counter that we've written
202 310 : if (offset > dataset->highestOffsetWritten)
203 310 : dataset->highestOffsetWritten = offset;
204 :
205 310 : return CE_None;
206 : }
207 :
208 : /************************************************************************/
209 : /* ==================================================================== */
210 : /* FITSDataset */
211 : /* ==================================================================== */
212 : /************************************************************************/
213 :
214 : // Some useful utility functions
215 :
216 : // Simple static function to determine if FITS header keyword should
217 : // be saved in meta data.
218 : static const int ignorableHeaderCount = 15;
219 : static const char* ignorableFITSHeaders[ignorableHeaderCount] = {
220 : "SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "NAXIS3", "END",
221 : "XTENSION", "PCOUNT", "GCOUNT", "EXTEND", "CONTINUE",
222 : "COMMENT", "", "LONGSTRN"
223 : };
224 399 : static bool isIgnorableFITSHeader(const char* name) {
225 3148 : for (int i = 0; i < ignorableHeaderCount; ++i) {
226 3140 : if (strcmp(name, ignorableFITSHeaders[i]) == 0)
227 391 : return true;
228 : }
229 8 : return false;
230 : }
231 :
232 :
233 : /************************************************************************/
234 : /* FITSDataset() */
235 : /************************************************************************/
236 :
237 39 : FITSDataset::FITSDataset() {
238 39 : hFITS = 0;
239 39 : }
240 :
241 : /************************************************************************/
242 : /* ~FITSDataset() */
243 : /************************************************************************/
244 :
245 78 : FITSDataset::~FITSDataset() {
246 :
247 : int status;
248 39 : if (hFITS) {
249 39 : if(eAccess == GA_Update) { // Only do this if we've successfully opened the file and update capability
250 : // Write any meta data to the file that's compatible with FITS
251 26 : status = 0;
252 26 : fits_movabs_hdu(hFITS, 1, 0, &status);
253 26 : fits_write_key_longwarn(hFITS, &status);
254 26 : if (status) {
255 : CPLError(CE_Warning, CPLE_AppDefined,
256 : "Couldn't move to first HDU in FITS file %s (%d).\n",
257 0 : GetDescription(), status);
258 : }
259 26 : char** metaData = GetMetadata();
260 26 : int count = CSLCount(metaData);
261 42 : for (int i = 0; i < count; ++i) {
262 16 : const char* field = CSLGetField(metaData, i);
263 16 : if (strlen(field) == 0)
264 0 : continue;
265 : else {
266 : char* key;
267 16 : const char* value = CPLParseNameValue(field, &key);
268 : // FITS keys must be less than 8 chars
269 16 : if (strlen(key) <= 8 && !isIgnorableFITSHeader(key)) {
270 : // Although FITS provides support for different value
271 : // types, the GDAL Metadata mechanism works only with
272 : // string values. Prior to about 2003-05-02, this driver
273 : // would attempt to guess the value type from the metadata
274 : // value string amd then would use the appropriate
275 : // type-specific FITS keyword update routine. This was
276 : // found to be troublesome (e.g. a numeric version string
277 : // with leading zeros would be interpreted as a number
278 : // and might get those leading zeros stripped), and so now
279 : // the driver writes every value as a string. In practice
280 : // this is not a problem since most FITS reading routines
281 : // will convert from strings to numbers automatically, but
282 : // if you want finer control, use the underlying FITS
283 : // handle. Note: to avoid a compiler warning we copy the
284 : // const value string to a non const one...
285 2 : char* valueCpy = strdup(value);
286 2 : fits_update_key_longstr(hFITS, key, valueCpy, 0, &status);
287 2 : free(valueCpy);
288 :
289 : // Check for errors
290 2 : if (status) {
291 : CPLError(CE_Warning, CPLE_AppDefined,
292 : "Couldn't update key %s in FITS file %s (%d).",
293 0 : key, GetDescription(), status);
294 : return;
295 : }
296 : }
297 : // Must free up key
298 16 : CPLFree(key);
299 : }
300 : }
301 :
302 : // Make sure we flush the raster cache before we close the file!
303 26 : FlushCache();
304 : }
305 :
306 : // Close the FITS handle - ignore the error status
307 39 : status = 0;
308 39 : fits_close_file(hFITS, &status);
309 :
310 : }
311 39 : }
312 :
313 :
314 : /************************************************************************/
315 : /* Init() */
316 : /************************************************************************/
317 :
318 39 : CPLErr FITSDataset::Init(fitsfile* hFITS_, bool isExistingFile_) {
319 :
320 39 : hFITS = hFITS_;
321 39 : isExistingFile = isExistingFile_;
322 39 : highestOffsetWritten = 0;
323 39 : int status = 0;
324 :
325 : // Move to the primary HDU
326 39 : fits_movabs_hdu(hFITS, 1, 0, &status);
327 39 : if (status) {
328 : CPLError(CE_Failure, CPLE_AppDefined,
329 : "Couldn't move to first HDU in FITS file %s (%d).\n",
330 0 : GetDescription(), status);
331 0 : return CE_Failure;
332 : }
333 :
334 : // The cftsio driver automatically rescales data on read and write
335 : // if BSCALE and BZERO are defined as header keywords. This behavior
336 : // causes overflows with GDAL and is slightly mysterious, so we
337 : // disable this rescaling here.
338 39 : fits_set_bscale(hFITS, 1.0, 0.0, &status);
339 :
340 : // Get the image info for this dataset (note that all bands in a FITS dataset
341 : // have the same type)
342 : int bitpix;
343 : int naxis;
344 39 : const int maxdim = 3;
345 : long naxes[maxdim];
346 39 : fits_get_img_param(hFITS, maxdim, &bitpix, &naxis, naxes, &status);
347 39 : if (status) {
348 : CPLError(CE_Failure, CPLE_AppDefined,
349 : "Couldn't determine image parameters of FITS file %s (%d).",
350 0 : GetDescription(), status);
351 0 : return CE_Failure;
352 : }
353 :
354 : // Determine data type
355 39 : if (bitpix == BYTE_IMG) {
356 19 : gdalDataType = GDT_Byte;
357 19 : fitsDataType = TBYTE;
358 : }
359 20 : else if (bitpix == SHORT_IMG) {
360 5 : gdalDataType = GDT_Int16;
361 5 : fitsDataType = TSHORT;
362 : }
363 15 : else if (bitpix == LONG_IMG) {
364 5 : gdalDataType = GDT_Int32;
365 5 : fitsDataType = TINT;
366 : }
367 10 : else if (bitpix == FLOAT_IMG) {
368 5 : gdalDataType = GDT_Float32;
369 5 : fitsDataType = TFLOAT;
370 : }
371 5 : else if (bitpix == DOUBLE_IMG) {
372 5 : gdalDataType = GDT_Float64;
373 5 : fitsDataType = TDOUBLE;
374 : }
375 : else {
376 : CPLError(CE_Failure, CPLE_AppDefined,
377 0 : "FITS file %s has unknown data type: %d.", GetDescription(),
378 0 : bitpix);
379 0 : return CE_Failure;
380 : }
381 :
382 : // Determine image dimensions - we assume BSQ ordering
383 39 : if (naxis == 2) {
384 30 : nRasterXSize = naxes[0];
385 30 : nRasterYSize = naxes[1];
386 30 : nBands = 1;
387 : }
388 9 : else if (naxis == 3) {
389 9 : nRasterXSize = naxes[0];
390 9 : nRasterYSize = naxes[1];
391 9 : nBands = naxes[2];
392 : }
393 : else {
394 : CPLError(CE_Failure, CPLE_AppDefined,
395 : "FITS file %s does not have 2 or 3 dimensions.",
396 0 : GetDescription());
397 0 : return CE_Failure;
398 : }
399 :
400 : // Create the bands
401 194 : for (int i = 0; i < nBands; ++i)
402 58 : SetBand(i+1, new FITSRasterBand(this, i+1));
403 :
404 : // Read header information from file and use it to set metadata
405 : // This process understands the CONTINUE standard for long strings.
406 : // We don't bother to capture header names that duplicate information
407 : // already captured elsewhere (e.g. image dimensions and type)
408 39 : int keyNum = 1;
409 : char key[100];
410 : char value[100];
411 39 : bool endReached = false;
412 436 : do {
413 436 : fits_read_keyn(hFITS, keyNum, key, value, 0, &status);
414 436 : if (status) {
415 : CPLError(CE_Failure, CPLE_AppDefined,
416 : "Error while reading key %d from FITS file %s (%d)",
417 0 : keyNum, GetDescription(), status);
418 0 : return CE_Failure;
419 : }
420 : // What we do next depends upon the key and the value
421 436 : if (strcmp(key, "END") == 0) {
422 39 : endReached = true;
423 : }
424 397 : else if (isIgnorableFITSHeader(key)) {
425 : // Ignore it!
426 : }
427 : else { // Going to store something, but check for long strings etc
428 : // Strip off leading and trailing quote if present
429 6 : char* newValue = value;
430 6 : if (value[0] == '\'' && value[strlen(value) - 1] == '\'') {
431 6 : newValue = value + 1;
432 6 : value[strlen(value) - 1] = '\0';
433 : }
434 : // Check for long string
435 6 : if (strrchr(newValue, '&') == newValue + strlen(newValue) - 1) {
436 : // Value string ends in "&", so use long string conventions
437 0 : char* longString = 0;
438 0 : fits_read_key_longstr(hFITS, key, &longString, 0, &status);
439 : // Note that read_key_longst already strips quotes
440 0 : if (status) {
441 : CPLError(CE_Failure, CPLE_AppDefined,
442 : "Error while reading long string for key %s from "
443 0 : "FITS file %s (%d)", key, GetDescription(), status);
444 0 : return CE_Failure;
445 : }
446 0 : SetMetadataItem(key, longString);
447 0 : free(longString);
448 : }
449 : else { // Normal keyword
450 6 : SetMetadataItem(key, newValue);
451 : }
452 : }
453 436 : ++keyNum;
454 : } while (!endReached);
455 :
456 39 : return CE_None;
457 : }
458 :
459 :
460 :
461 : /************************************************************************/
462 : /* Open() */
463 : /************************************************************************/
464 :
465 9643 : GDALDataset* FITSDataset::Open(GDALOpenInfo* poOpenInfo) {
466 :
467 9643 : const char* fitsID = "SIMPLE = T"; // Spaces important!
468 9643 : size_t fitsIDLen = strlen(fitsID); // Should be 30 chars long
469 :
470 9643 : if ((size_t)poOpenInfo->nHeaderBytes < fitsIDLen)
471 8661 : return NULL;
472 982 : if (memcmp(poOpenInfo->pabyHeader, fitsID, fitsIDLen))
473 968 : return NULL;
474 :
475 : // Get access mode and attempt to open the file
476 14 : int status = 0;
477 14 : fitsfile* hFITS = 0;
478 14 : if (poOpenInfo->eAccess == GA_ReadOnly)
479 13 : fits_open_file(&hFITS, poOpenInfo->pszFilename, READONLY, &status);
480 : else
481 1 : fits_open_file(&hFITS, poOpenInfo->pszFilename, READWRITE, &status);
482 14 : if (status) {
483 : CPLError(CE_Failure, CPLE_AppDefined,
484 : "Error while opening FITS file %s (%d).\n",
485 0 : poOpenInfo->pszFilename, status);
486 0 : fits_close_file(hFITS, &status);
487 0 : return NULL;
488 : }
489 :
490 : // Create a FITSDataset object and initialize it from the FITS handle
491 14 : FITSDataset* dataset = new FITSDataset();
492 14 : dataset->eAccess = poOpenInfo->eAccess;
493 :
494 : // Set up the description and initialize the dataset
495 14 : dataset->SetDescription(poOpenInfo->pszFilename);
496 14 : if (dataset->Init(hFITS, true) != CE_None) {
497 0 : delete dataset;
498 0 : return NULL;
499 : }
500 : else
501 : {
502 : /* -------------------------------------------------------------------- */
503 : /* Initialize any PAM information. */
504 : /* -------------------------------------------------------------------- */
505 14 : dataset->SetDescription( poOpenInfo->pszFilename );
506 14 : dataset->TryLoadXML();
507 :
508 14 : return dataset;
509 : }
510 : }
511 :
512 :
513 : /************************************************************************/
514 : /* Create() */
515 : /* */
516 : /* Create a new FITS file. */
517 : /************************************************************************/
518 :
519 37 : GDALDataset *FITSDataset::Create(const char* pszFilename,
520 : int nXSize, int nYSize,
521 : int nBands, GDALDataType eType,
522 : char** papszParmList) {
523 :
524 :
525 : FITSDataset* dataset;
526 : fitsfile* hFITS;
527 37 : int status = 0;
528 :
529 : // No creation options are defined. The BSCALE/BZERO options were
530 : // removed on 2002-07-02 by Simon Perkins because they introduced
531 : // excessive complications and didn't really fit into the GDAL
532 : // paradigm.
533 :
534 : // Create the file - to force creation, we prepend the name with '!'
535 37 : char* extFilename = new char[strlen(pszFilename) + 10]; // 10 for margin!
536 37 : sprintf(extFilename, "!%s", pszFilename);
537 37 : fits_create_file(&hFITS, extFilename, &status);
538 37 : delete[] extFilename;
539 37 : if (status) {
540 : CPLError(CE_Failure, CPLE_AppDefined,
541 0 : "Couldn't create FITS file %s (%d).\n", pszFilename, status);
542 0 : return NULL;
543 : }
544 :
545 : // Determine FITS type of image
546 : int bitpix;
547 37 : if (eType == GDT_Byte)
548 13 : bitpix = BYTE_IMG;
549 24 : else if (eType == GDT_Int16)
550 3 : bitpix = SHORT_IMG;
551 21 : else if (eType == GDT_Int32)
552 3 : bitpix = LONG_IMG;
553 18 : else if (eType == GDT_Float32)
554 3 : bitpix = FLOAT_IMG;
555 15 : else if (eType == GDT_Float64)
556 3 : bitpix = DOUBLE_IMG;
557 : else {
558 : CPLError(CE_Failure, CPLE_AppDefined,
559 12 : "GDALDataType (%d) unsupported for FITS", eType);
560 12 : fits_close_file(hFITS, &status);
561 12 : return NULL;
562 : }
563 :
564 : // Now create an image of appropriate size and type
565 25 : long naxes[3] = {nXSize, nYSize, nBands};
566 25 : int naxis = (nBands == 1) ? 2 : 3;
567 25 : fits_create_img(hFITS, bitpix, naxis, naxes, &status);
568 :
569 : // Check the status
570 25 : if (status) {
571 : CPLError(CE_Failure, CPLE_AppDefined,
572 : "Couldn't create image within FITS file %s (%d).",
573 0 : pszFilename, status);
574 0 : fits_close_file(hFITS, &status);
575 0 : return NULL;
576 : }
577 :
578 25 : dataset = new FITSDataset();
579 25 : dataset->nRasterXSize = nXSize;
580 25 : dataset->nRasterYSize = nYSize;
581 25 : dataset->eAccess = GA_Update;
582 25 : dataset->SetDescription(pszFilename);
583 :
584 : // Init recalculates a lot of stuff we already know, but...
585 25 : if (dataset->Init(hFITS, false) != CE_None) {
586 0 : delete dataset;
587 0 : return NULL;
588 : }
589 : else {
590 25 : return dataset;
591 : }
592 : }
593 :
594 :
595 : /************************************************************************/
596 : /* GDALRegister_FITS() */
597 : /************************************************************************/
598 :
599 338 : void GDALRegister_FITS() {
600 :
601 : GDALDriver* poDriver;
602 :
603 338 : if( GDALGetDriverByName( "FITS" ) == NULL) {
604 336 : poDriver = new GDALDriver();
605 :
606 336 : poDriver->SetDescription( "FITS" );
607 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
608 336 : "Flexible Image Transport System" );
609 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
610 336 : "frmt_various.html#FITS" );
611 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
612 336 : "Byte Int16 Int32 Float32 Float64" );
613 :
614 336 : poDriver->pfnOpen = FITSDataset::Open;
615 336 : poDriver->pfnCreate = FITSDataset::Create;
616 336 : poDriver->pfnCreateCopy = NULL;
617 :
618 336 : GetGDALDriverManager()->RegisterDriver(poDriver);
619 : }
620 338 : }
|
