[fitsbits] Adaptive mesh refinement

Folks:

This is a modified version of an email I sent earlier to the fitswcs mailing
list. I thought it probably made sense to open up the discussion to the wider
FITS community.

I was asked how to store data from model calculations using adaptive mesh
refinement in FITS files. To my mind, there are two possibilities, pixel lists
or the -TAB projection. Personally, I prefer the -TAB projection, for reasons
listed below, but I think that one sentence in the definition makes it
impracticle for this application. Below I present my understanding of each of
the options. Please let me know if my understanding is incorrect or incomplete.

Scenario 1: The -TAB projection.

The advantages of the -TAB projection are that it explicitly ties together the
data array with the coordinate information, and that the coordinates are stored
directly without needing any additional computations.

Suppose that we have a model with 1000 points in it, each of which has
coordinates x(i), y(i), z(i). My belief is that I could store this using the
following keywords. (Some keywords, like CUNITn, have been left out for brevity.)

NAXIS = 1 /One data axis
NAXIS1 = 1000
WCSAXES = 3 /Three coordinate axes

CTYPE1 = 'HEQX-TAB' /Heliocentric Earth Equatorial coordinates
CRPIX1 = 1
CRVAL1 = 1
CDELT1 = 1
CTYPE2 = 'HEQY-TAB'
CRPIX2 = 1
CRVAL2 = 1
CDELT2 = 1
CTYPE3 = 'HEQZ-TAB'
CRPIX3 = 1
CRVAL3 = 1
CDELT3 = 1

PC1_1 = 1 /Use PC matrix to tie coordinates together.
PC1_2 = 0
PC1_3 = 0
PC2_1 = 1
PC2_2 = 1
PC2_3 = 0
PC3_1 = 1
PC3_2 = 0
PC3_3 = 1

PS1_0 = 'WCS-TAB'
PV1_1 = 1
PV1_2 = 1
PS1_1 = 'COORDS'
PV1_3 = 1

PS2_0 = 'WCS-TAB'
PV2_1 = 1
PV2_2 = 1
PS2_1 = 'COORDS'
PV2_3 = 2

PS3_0 = 'WCS-TAB'
PV3_1 = 1
PV3_2 = 1
PS3_1 = 'COORDS'
PV3_3 = 3

The coordinate array in the binary table should have the dimensions
(3,1000,1,1). However, the documentation disallows degenerate axes, so we'd
have to substitute (3,1000,2,2), of which only 25% is actually used. This is
highly inefficient, and makes the scheme impracticle. Is this a correct
assessment? Is anything wrong about the above header?


Scenario 2: Pixel lists

The most straightforward way to use pixel lists would seem to be to set
TCRPXn=0, TCRVLn=0, TCDLTn=1, so that the coordinates can be stored directly in
the table. One would then have a binary table header such as the following

NAXIS = 2 /Binary table
NAXIS1 = 16 /Number of bytes per row
NAXIS2 = 1000 /Number of rows
PCOUNT = 0 /Random parameter count
GCOUNT = 1 /Group count
TFIELDS = 4 /Number of columns

TFORM1 = '1E ' /Real*4 (floating point)
TTYPE1 = 'DENSITY ' /Label for column 1

TFORM2 = '1E ' /Real*4 (floating point)
TCTYP2 = 'HEQX ' /
TCRPX2 = 0 /
TCRVL2 = 0.00000 /
TCDLT2 = 1.00000 /

TFORM3 = '1E ' /Real*4 (floating point)
TCTYP3 = 'HEQY ' /
TCRPX3 = 0 /
TCRVL3 = 0.00000 /
TCDLT3 = 1.00000 /

TFORM4 = '1E ' /Real*4 (floating point)
TCTYP4 = 'HEQZ ' /
TCRPX4 = 0 /
TCRVL4 = 0.00000 /
TCDLT4 = 1.00000 /

Some care would have to be exercised when spherical coordinates are used. To
keep the philosophy of storing the coordinates directly, I think that one would
have to treat the data as being in the plate carree (-CAR) projection.

Thank you,

Bill Thompson



--
William Thompson
NASA Goddard Space Flight Center
Code 612.1
Greenbelt, MD 20771
USA

301-286-2040