docs menu

Module Ymir.Wcs

FITS world coordinates.

A header's celestial description reads as the stages

ocaml

axes [| 1; 0 |] ~origin:1 (* tensor order, 0-based -> FITS, 1-based *)
>> shift crpix
>> sip (a, b) (* a SIP header *)
>> linear cd (* a CD header *)
(* or: linear pc >> scale cdelt, a PC, CROTA2 or CDELT header *)
>> axes [| 1; 0 |] ~origin:0 (* when latitude is axis 1 *)
>> tpv pv (* a TPV header *)
>> celestial code frame ~pv ~native ~crval ~lonpole ~latpole

CRPIX, CD, PC, CDELT, SIP's coefficients and the PV terms are held as the header writes them, CRVAL in CUNIT (degrees when absent), and LONPOLE, LATPOLE and the native reference point (φ₀, θ₀) = (PVi_1, PVi_2) on the longitude axis i in degrees. A keyword FITS leaves out takes its default: CRPIX and CRVAL 0, CDELT 1, PC the identity, a missing CD element or SIP coefficient 0, a projection's parameters those Transform.code lists, (φ₀, θ₀) = (0°, 90°) for a zenithal projection and (0°, 0°) for a cylindrical one, LATPOLE 90°, and LONPOLE 180° + φ₀ when CRVAL's latitude is below θ₀, φ₀ otherwise. Where θ₀ = 0, CRVAL's latitude is 0 and LONPOLE is 90° from φ₀, LATPOLE alone fixes the celestial pole's latitude and has no default: a header without it is an Error, which astropy (WCSLIB) reads with the pole at 90°. CROTA2 reads as the PC matrix it implies. PVi_3 and PVi_4 on the longitude axis are LONPOLE's and LATPOLE's other spellings outside TPV; a header giving both spellings of one with different values is an Error.

Distortions. SIP reads from A_ORDER, B_ORDER and A_p_q, B_p_q, with AP_ORDER, BP_ORDER, AP_p_q and BP_p_q as the inverse's seed, when the CTYPEs end in -SIP or the primary description has A_ORDER. TPV reads from PVi_0 to PVi_39 on both axes when the CTYPEs say TPV, or say TAN and the latitude axis has PV terms, as SCAMP writes; PV1_1 and PV2_1 default to 1, the other terms to 0.

Frames. The CTYPE prefix chooses the system and RADESYS (or RADECSYS) and EQUINOX (or EPOCH) qualify an equatorial or ecliptic one: absent both, ICRS; EQUINOX alone, FK4 below 1984.0 and FK5 from it; EQUINOX 2000.0 under FK5 when absent; EQUINOX ignored under ICRS.

RA/DEC under ICRS is Frame.icrs, under FK5 at 2000.0 Frame.fk5_j2000.
GLON/GLAT is Frame.galactic, SLON/SLAT Frame.supergalactic.
ELON/ELAT under ICRS at 2000.0 is Frame.ecliptic_j2000.

Other systems and equinoxes are an Error naming what they need.

Scope. The thirteen projections of Transform.code, SIP and TPV are read. Other projections, distortion lookup tables, the fiducial offset PVi_0, PVi_m on the longitude axis beyond m = 1, ..., 4 outside TPV, and a third WCS axis are an Error naming the keyword.

#type keywords = {

#float : string -> (float option, string) Stdlib.result;
#int : string -> (int option, string) Stdlib.result;
#text : string -> (string option, string) Stdlib.result;

}

The type for a header's keywords, as decoded values. kw.float k is the real under keyword k: Ok None when k is absent, an Error when k does not decode as a real or its cards disagree. int is the same for integers, and text for strings, without their quotes. These errors pass through read unchanged. With ymir.fits:

ocaml

let keywords h =
  {
    Wcs.float = (fun k -> Fits.Header.find Fits.Value.float k h);
    int = (fun k -> Fits.Header.find Fits.Value.int k h);
    text = (fun k -> Fits.Header.find Fits.Value.string k h);
  }

#type edit =

#| Float of string * float
(*
Set a keyword to a real.
*)
#| Int of string * int
(*
Set a keyword to an integer.
*)
#| Text of string * string
(*
Set a keyword to a string.
*)
#| Remove of string
(*
Remove a keyword.
*)

The type for keyword edits. With ymir.fits:

ocaml

let edit h = function
  | Wcs.Float (k, x) -> Fits.Header.set Fits.Value.float k x h
  | Int (k, n) -> Fits.Header.set Fits.Value.int k n h
  | Text (k, s) -> Fits.Header.set Fits.Value.string k s h
  | Remove k -> Fits.Header.remove k h

val read : 
  ?alt:char ->
  at:'d Nx.Placement.t ->
  'f Frame.t ->
  keywords ->
  (('d Transform.plane, ('f, 'd) Direction.t, 'd) Transform.t, string)
    Stdlib.result

read ~alt ~at f kw is the map from 0-based pixel indices, in tensor axis order, to directions in f, from kw's primary description or its alternate alt ('A' to 'Z'), its parameters at at. A celestial CUNIT is deg, arcmin, arcsec, mas or rad. It is an Error if kw's frame is not f, as in "SCI: FK5 at equinox 2000.0; the caller expects icrs. Read with Frame.fk5_j2000.", if a keyword's value is outside what WCS allows, if LATPOLE is absent where it alone fixes the pole (see above), or as kw is. Its own errors start with the header's EXTNAME.

Raises Invalid_argument if alt is not ' ' or 'A' to 'Z'.

val write : 
  ?alt:char ->
  ?window:(int * int) array ->
  ('d Transform.plane, ('f, 'd) Direction.t, 'd) Transform.t ->
  keywords ->
  (edit list, string) Stdlib.result

write ~alt ~window t kw is the edits, in order, that make kw's description alt spell t, with window, (start, stop) per tensor axis, taking each CRPIX less its axis's start. A keyword whose value equals t's, bit for bit, gets no edit, nor does one kw lacks whose value is its default and, for a PV term, that the stage does not list as stated. Keywords t no longer spells (CD for a PC header, CROTA, SIP) are removed. A pole keeps the spelling kw gives it. write (read kw) kw is Ok [] for every header read reads, except a CROTA2 header, whose edits spell the PC matrix it reads as, and SCAMP's TAN with PV terms, whose edits spell TPV; their stages are equal.

It is an Error naming the stage if t is not a list read builds, up to an absent PC or a scale stage on its own, or if t's leaves are batched or not finite.

Raises Invalid_argument if alt is not ' ' or 'A' to 'Z', or window does not have two ranges.