In astrophysics and in solar physics we have FITS format as a de-facto standard for 2D images. Of course there are libraries for all the mainstream languages & viewers.
In particular for C and Fortran there exist CFITSIO library (link). A while ago I’ve managed to get read/write access to FITS-files in my programs, but I admit that I don’t understand why it is working.
So, if I don’t get blind, there are no any interfaces defined and probably implicit interfaces were implied in my make (I’ve used default gfortran and ifort ~2019). Now when I’m trying to switch to fpm, I get an error about it, so had to specify implicit-external=true
So my questions are:
If anybody has experience with FITS?
If there exist a more modern approach/library?
Am I missing something in CFITSIO? They provide examples of compiling&linking and a simple program in a f77 style.
If not, is it a good idea (is it straightforward?) to write an interface and generics? For example to read data there is a set of functions like FTGPV[BIJKED] for different data types ( I - short int; J - int; etc)
I have no experience with this type of files and I have not checked the library as yet, but if the interfacing method is an old-style one (i.e. predating Fortran 2003 with ISO_C_BINDING etc.) then indeed the interfaces can be completely implicit and the burden is on you as a programmer to get everything right.
I am working myself on the interface to an extensive C library and I prefer to provide generic interfaces to routines that manipulate the various basic data types. It is a matter of taste, I guess, but I find it more enjoyable to be able to use such generic interfaces instead of having to wonder about the specific names.
As a research astronomer we used FITS files extensively and they are still in use widely, especially in radio astronomy and x-ray astronomy. As far as I know the CFITSIO library is still the only one that is C- and Fortran-callable, but as you note its interfaces were designed in the Fortran77 era. The library was written originally in Fortran77 at NASA Goddard Space Flight Institute by Bill Pence and his team. One problem was that there is, in Fortran77, no standard-conforming way of accessing a binary stream. It was in fact possible to access them direct-access files with a record length of 2880 bytes, and this works well, but specifying a record length in bytes was, before the IOLENGTH descriptor arrived, not easy to do in a standard-conforming way. Thus it was hard to write a library for FITS file access without using non-standard Fortran77.
Anyhow, many years ago a direction came from on high in NASA that they should no longer use Fortran77 so they re-wrote the whole thing in C with a Fortran77-callable interface on top. I tried to persuade them that Fortran90/95 (and later) standards were superior but GSFC never accepted that. So what exists now is a library written in C with a thin Fortran77 style interface.
The problem with providing a sensible modern interface layer is that a few procedures in CFITSIO have one or two arguments where the actual argument can be either a scalar or an array. I think there are work-arounds for this but it is messy. I constructed Fortran90 interfaces for some of the CFITSIO procedures that I needed for my own work, but decided that it would be too much effort to do this for all of them. It might even be easier to write a new library in modern Fortran for accessing FITS files, accessing them via access=“stream”. But as far as I know nobody has done that yet.
I agree with this with one caveat. You have to be extra careful if you need exact byte-for-byte compatability and reproducability with existing files written by CFITSIO.
If that’s the case you should consider writting F2003 interop interfaces for just the routines you need. The issue with the CFITIO F77 interfaces is they appear to be generated using the old cfortran.h macro package. There can be subtle differences depending which operating system and C compiler you use to build the C library. I encountered this 20 plus years ago when I wrote the first version of the C-interop interfaces for netCDF that are now part of the netCDF Fortran package. Based on my experiences with netCDF and later my forSISL interfaces to the SISL NURBS library, Its not very complicated to write the interfaces but its a LOT of work (and I mean a LOT).
Is the issue that the CFITSIO library itself will return a different binary depending on how you compile it or which platform you are on?
If so, then there isn’t much we can do, the binary format itself is not stable across platforms or library versions/compilation options.
However if the binary itself is stable, then I recommend writing our own in pure Fortran and not bother with brittle interfaces and build systems. Yes, it has to be done to be 100% compatible with existing writers/readers. But that’s the same issue for any other binary format out there. Modern Fortran allows precise cross-platform binary writing and reading.
I think the CFITSIO C side binary is stable. The issue is with the cfortran.h macro package that generates the Fortran callable C wrappers. When I did the first netCDF F2003 C-interop interfaces around 2006, I was running the cfortran.h macros through a preprocessor and looking at the generated C code for the wrappers for different Fortran compilers and operating systems. The biggest issue was with the different name mangling (one underscore, two underscores, ALL CAPS) etc along with what was a long vs what was an int etc. Again this was 20 years ago so things have probably changed. My comment about binary compatability was based on my memory (now a vague one) of thinking I got the interfaces correct but seeing slight differences in files generated using my interfaces and the cfortran.h interfaces when I did compare on the resulting files. It took me a while to get them to give the same results but in all honesty I don’t remember what I did. Note that my issues with trying to use cfortran.h on a Cray X1 were the reason I did the C-interop interfaces in the first place. On the X1 the Fortran compiler was the classic Cray vector compiler where both Reals and Integers were both 64 bits by default. However, the C compiler was ANSI C so the default int was 32 bits. This gave gfortran.h fits.
Also, my comments were based on an ignorance of the FITS file format which (after a little research) appears to be a mixture of ASCII and binary (if thats what the user chooses). There is a Fortran 77 package that doesn’t appear to be a wrapper around CFITSIO but I think it has limited capability. Again, I agree that a new from the ground up pure Fortran (or at least as much as possible - CFITSIO appears to support a couple of compression algorithms which would require interfaces to 3rd party libraries) is the optimum solution.
I agree with @certik .
Fortran access=stream provides a way to read/write any file data or record structure.
Stream access is positioned by “file storage units”, which is byte addressable on all Fortran compilers I am now using, eliminating the byte addressing problems of some past Fortran compilers.
With stream=access, I have been able to reproduce Fortran unformatted records for a variety of Fortran compilers, where most formats have horrendous header structures for large records. I have also tested my own simpler record formats.
If you know the file format for FITS IO, I am sure you could clearly create this with stream access.
Yes, it would be nice to have (write) a modern FITS-IO library, but in my case I’m afraid I have neither motivation nor free time to achieve more-or-less mature state. In addition, for my purposes I need very limited set of FITS features: for Fortran-related part I use FITS mainly for storing the results of my optical simulations; otherwise in solar physics we use mostly python (nowadays) and IDL to store&process solar images.
BTW, is the question of various architectures overestimated these days? In my (professional) life I’ve never seen anything significantly different from x86 or arm. All the clusters (all the 4) I’ve seen were Intel-based and running linux, with gfortran and ifort being main compilers. I can’t even imagine one would use something else.
The FITS is relatively simple and standardized, with the header that describes the body quite in detail. In particular 2D image can be stored as 2-4 byte integers or 4-8 byte reals (in big-endian). Thus swapping of bytes is the first “problem”. Do we have a standardized way for that? It is more straightforward for integers, but probably not for reals.
My second question is the following: I’ve asked ChatGPT few times, all the times the answers were completely different (he probably is not completely sure what he is suggesting). But he was suggesting OO-like way of working with the file:
Is there any advantage of using OO-like approach? In IDL I do it like data=readfits("best_solar_image.fits",header) and it works. Of course, for the case of huge files one may need to get the information about the data before reading the actual data; which is normally not the case of solar physics, because we are limited by the size of a CCD (normally up to 4k x 4k). Is there anything else?
There used to be a lot of variations in the binary representation of reals, but nowadays it is indeed merely little-endian (the majority of machines in use, I guess) and big-endian. It is quite possible to swap the bytes using standard Fortran, as that is the only thing that changes. Not the exponent or the mantissa.
If your compiler supports integer(int8), it could be as simple as:
integer(kind=int8) :: bytes(8)
real(kind=real64) :: x, y
x = 1.0d0
bytes = transfer( x, bytes)
bytes = bytes(8:1:-1)
y = transfer( bytes, y )
(As I have no big-endian data sets at hand, I have checked the code by writing x and y to a binary file )
function le_to_be(x) result(y)
use, intrinsic :: iso_fortran_env, only: int64, real64
real(kind=real64) :: x, y
integer(int64) :: i, j
i = transfer(x, i)
j = 0
do k = 0, 7
call mvbits(i, 8*k, 8, j, 8*(7-k))
end do
y = transfer(j, y)
end function
With the Intel ifx compiler it compiles to a bswap (byte swap) instruction:
That style of interface is attractive, but you probably want to use it to e.g. read many FITS files and maybe write some concurrently, so it probably needs a file handle or unit number or something as an additional argument of each call.
That is an impressive reduction of a lot of fortran code down to a single machine instruction. I’m guessing that the array-syntax-based code does something similar, but that code is closer to the target instruction to begin with. Do other compilers achieve a similar reduction down to a single instruction for these codes? The array syntax code can be shortened to just two transfer operations.
bytes = transfer( x, bytes )
y = transfer( bytes(8:1:-1), y )
Is it possible to reflux this down into a single statement?
The problem with providing a sensible modern interface layer is that a few procedures in CFITSIO have one or two arguments where the actual argument can be either a scalar or an array. I think there are work-arounds for this but it is messy.
If it is just an argument or two or if not a lot of calls are made to the routine it is not that bad to make a generic interface or to use an EQUIVANCE between a scalar and a single-value array to avoid having to create and call the procedures without an interface definition but I have had to work with old libraries where there were sometimes a dozen parameters that could be scalars, vectors, or matrices. It made doing a generic interface very cumbersome because of all the permutations so as an interrim I made GitHub - urbanjost/M_flatten: return a rank one array pointer to a scalar or array of any rank · GitHub. Passing anything to a procedure and then treating it in the procedure as a vector was so common I find it (and treating arrays as “typeless storage” where you store integers, arrays and strings in the same array) the hardest thing to deal with when porting old code and really wanting everything to be in a module (which of course means it has to conform strictly to the rules regarding parameter shapes). Does anyone else feel something like the “flatten” procedure should be available as a standard fortran feature? (I am not counting using the ISO_C_BINDING interface and C pointers as standard Fortran in this case)? Note I use M_flatten.f90 as an interrim solution while migrating code, not a final solution, although sometimes I never get around to that final solution. I am not sure how fpt and spagh deal with such code where parameters have different shapes. Maybe they have a better methodology for dealing with the issue.
It is a double edged sword. The type-kind-rank (TKR) argument matching catches many programmer errors, so relaxing that TKR matching would mean that those types of errors could not be detected at compile time. That is why many legacy codes with explicit interfaces have this error to begin with.
I think we talked about this a little in relation to this thread. Spacing of the results of Fortran intrinsic random_number - #36 by RonShepard The “solution” I ended up adopting there was to use assumed rank argument association, which then requires select rank blocks internally. If there is only one such argument, that isn’t so bad, but if there are many such arguments, and all possible rank matching is required for each of them, then the complexity can grow in an exponential way.
There is a special case that is easy to work around. If the actual argument is scalar x, and the dummy argument is intent(in) :: x(:), then you can change the call statement actual argument to [x]. This rank 1 expression matches the TKR dummy argument. But this doesn’t work with other intent, so it is not really a general solution.
Compiler Explorer can be used to check that out very easily, put the code and add a “Compiler” Window, you can add as many as you want and see side-by-side the assembly translation:
Just a few bits of information for anyone interested: The FITS Standard is covered in web pages at fits.gsfc.nasa.gov and a list of software libraries, include cFITSIO, is provided on that site at FITS Libraries Guide
This has a table listing libraries in many languages but as far as I know cFITSIO is the only one providing a really comprehensive set of procedures for reading/writing FITS files. I just counted up those in my (now slightly out-of-date) cFITSIO manual and there are 247 procedures. That total is much larger if one counts variants for the different data types. That suggests it’s quite a big job either to provide (a) a modern Fortran module with wrappers for each of these cFITSIO routines, or (b) to write a new module entirely in Fortran to replace it. It is obvious that cFITSIO was designed in a era when memory was expensive so it made sense to read a file in chunks, that is also why files are blocked in data units of 2880 bytes. That was a design feature being the lowest common multiple of the word-length of all the machines in use at the time it was written and was short enough to work on all known magnetic tape drives. My inclination is to encourage (b), that is a complete replacement, since a lot of features of modern Fortran could be used to simplify the interface considerably, something that would not be as easy with just a set of thin wrappers around cFITSIO.
Maybe a small group of those who have so far contributed here could start thinking about the design principles for such a replacement? One simplification would be to start by supporting just binary images and binary tables. In principle I’d be happy to be involved but I shall be away for most of the next three weeks.
)