github ↗
docs menu
on this page

Module Ymir.Transform

Maps between planes and the sky.

A transform is a list of stages, each one family's parameters and a sense, forward or inverse. Its endpoint types say what it maps: Transform.plane points, Frame.icrs Direction.t directions. Composing appends stages and inverting reverses them; no stage is ever fused, so what a file stated is what a writer prints.

A FITS celestial image stored [NAXIS2; NAXIS1] reads as

ocaml
axes [| 1; 0 |] ~origin:1
>> shift crpix >> sip (a, b) >> linear cd
>> tpv pv >> celestial Tan Frame.icrs ~pv ~native ~crval ~lonpole ~latpole

with Transform.sip and Transform.tpv present when the header has them. Transform.axes holds raven's pixel convention (0-based centres in the data tensor's axis order) and FITS's (1-based, NAXIS1 first) in one stage.

Units. Each planar family reads its input in the unit it states: Transform.axes with an origin, Transform.sip, Transform.linear and Transform.scale in Unit.one, Transform.shift in its offset's unit, Transform.tpv and Transform.celestial in an angle. A list in the wrong order raises Invalid_argument where a stage reads its input.

Domains. A stage is a bijection on its domain, the connected region about its reference point where its map does not fold. TAN's projection, from directions to the plane, needs native θ > 0; ARC's deprojection needs a distance from the reference point of at most 180°; ZPN's ends where its polynomial first turns. A distortion's domain is where its Jacobian's determinant is positive and, inverted, where Newton's method converges. Transform.apply at a finite point outside a domain raises Invalid_argument through Nx.check, naming the stage, the point and its distance or residual, as in "Transform.apply: point 12 is 93.1 deg from the TAN reference (110.8375, -73.4537) deg, outside the projection's domain; Transform.covers gives the mask". NaN maps to NaN. Transform.covers gives the mask for callers that expect points outside.

Precision. Planes are float64: every stage holds its parameters and computes in float64, so a header's numbers are held as the file wrote them. Data at another dtype sit on float64 geometry (Grid, Region).

Batch axes of points and parameters broadcast.

#type ('a, 'b, 'd) t

The type for maps from points of type 'a to points of type 'b.

#type 'd plane = (float, Nx.float64_elt, 'd) Nx.t Quantity.t

The type for points [...; n] in a plane. Pixel coordinates are in Unit.one, intermediate and tangent-plane coordinates in an angle.

#type code =
  1. #| Azp
    (*

    Zenithal perspective from distance μ, onto a plane tilted by γ: [|μ; γ|], PVi_1, PVi_2, defaults 0, 0.

    *)
  2. #| Szp
    (*

    Slant zenithal perspective from distance μ toward (φc, θc): [|μ; φc; θc|], PVi_1 to PVi_3, defaults 0, 0, 90.

    *)
  3. #| Tan
    (*

    Gnomonic: great circles map to lines.

    *)
  4. #| Stg
    (*

    Stereographic: circles map to circles.

    *)
  5. #| Sin
    (*

    Slant orthographic, along (ξ, η): [|ξ; η|], PVi_1, PVi_2, defaults 0, 0.

    *)
  6. #| Arc
    (*

    Zenithal equidistant: distance from the reference is preserved.

    *)
  7. #| Zpn
    (*

    Zenithal polynomial R = Σ Pₘ (90° − θ)ᵐ, θ in radians: 30 coefficients P₀ to P₂₉, PVi_0 to PVi_29, default 0.

    *)
  8. #| Zea
    (*

    Zenithal equal area.

    *)
  9. #| Air
    (*

    Airy's minimum-error zenithal, for a boundary at θb: [|θb|], PVi_1, default 90.

    *)
  10. #| Cyp
    (*

    Cylindrical perspective from distance μ onto a cylinder of radius λ: [|μ; λ|], PVi_1, PVi_2, defaults 1, 1.

    *)
  11. #| Cea
    (*

    Cylindrical equal area, λ the square of the cosine of the standard parallel: [|λ|], PVi_1, default 1.

    *)
  12. #| Car
    (*

    Plate carrée: x is longitude and y latitude.

    *)
  13. #| Mer
    (*

    Mercator's: conformal.

    *)

The type for projections, by their FITS codes. Zenithal projections put the native pole at the plane's origin, cylindrical ones the native point (0°, 0°). Each code's parameters pv are listed with the FITS keywords PVi_m of the latitude axis i that state them, and their defaults; angles are in degrees.

Constructors

#val id : ('a, 'a, 'd) t

id maps every point to itself.

#val axes : int array -> origin:int -> ('d plane, 'd plane, 'd) t

axes p ~origin maps x to x.(p.(k)) + origin in component k. With origin = 0 it permutes plane points of any unit; otherwise they are in Unit.one.

Raises Invalid_argument if p is not a permutation of 0, ..., n - 1.

#val shift : 'd plane -> ('d plane, 'd plane, 'd) t

shift r maps x to x - r: the plane whose origin is r. Points are read in r's unit.

Raises Invalid_argument if r is a scalar.

#val linear : 'd plane -> ('d plane, 'd plane, 'd) t

linear m maps x, in Unit.one, to m · x in m's unit: CD, or PC in Unit.one. m is [...; n; n].

Raises Invalid_argument if m is not square on its last two axes.

#val scale : 'd plane -> ('d plane, 'd plane, 'd) t

scale d maps x, in Unit.one, to x.(k) · d.(k) in component k, in d's unit: CDELT.

Raises Invalid_argument if d is a scalar.

#val sip : ?seed:((float, Nx.float64_elt, 'd) Nx.t * (float, Nx.float64_elt, 'd) Nx.t) -> ((float, Nx.float64_elt, 'd) Nx.t * (float, Nx.float64_elt, 'd) Nx.t) -> ('d plane, 'd plane, 'd) t

sip ~seed (a, b) maps pixel offsets (u, v) from CRPIX, in Unit.one, to (u + f, v + g) with SIP's polynomials f = Σ a.(p).(q) uᵖ v^q and g = Σ b.(p).(q) uᵖ v^q; a and b are [...; n; n], n the order plus one. Its inverse solves for each point with Newton's method, starting from seed, the file's AP and BP of the same form, applied as (u + ap, v + bp), or from the point itself.

Raises Invalid_argument if a matrix is not square on its last two axes.

#val tpv : ?stated:int array -> (float, Nx.float64_elt, 'd) Nx.t -> ('d plane, 'd plane, 'd) t

tpv ~stated pv maps intermediate coordinates (ξ, η), in an angle, to (Σ pv.(0).(k) t_k (ξ, η), Σ pv.(1).(k) t_k (η, ξ)) in degrees, over TPV's forty terms t_k: 1, x, y, r, then the monomials of each degree from 2 to 7 in decreasing powers of x, with r³, r⁵ and r⁷ after degrees 3, 5 and 7, r = √(x² + y²). pv is [...; 2; 40], the longitude axis's PVi_k then the latitude axis's, with the terms the file left out filled: 1 for PV1_1 and PV2_1, 0 otherwise. stated lists the indices of pv flattened on its last two axes that the file gave, all by default. Its inverse solves for each point with Newton's method from the point itself.

Raises Invalid_argument if pv is not [...; 2; 40] or stated not ascending indices below 80.

#val celestial : ?stated:int array -> code -> 'f Frame.t -> pv:(float, Nx.float64_elt, 'd) Nx.t -> native:'d plane -> crval:'d plane -> lonpole:'d plane -> latpole:'d plane -> ('d plane, ('f, 'd) Direction.t, 'd) t

celestial code f ~pv ~native ~crval ~lonpole ~latpole is the projection code with parameters pv, [...; m], followed by the rotation taking the native point native (φ₀, θ₀) to crval (lon, lat) in f, with the celestial pole at native longitude lonpole; latpole picks between two poles where FITS defines two. Angles are [...; 2] for native and crval and scalars for the poles, in any angle unit, as the file gives them. Plane points are angles. stated lists the indices of pv the file gave, all by default.

The plane's origin is the projection's own reference point, (0°, 90°) for a zenithal code and (0°, 0°) for a cylindrical one, whatever native is. The celestial pole is solved as FITS WCS Paper II states, from CRVAL, (φ₀, θ₀) and LONPOLE; where two poles solve, the one nearer latpole is taken.

Where pv does not define the projection or no pole solves, the stage raises Invalid_argument through Nx.check where it is applied, naming the parameters.

Raises Invalid_argument if pv's last axis is not m, if stated is not ascending indices below m, or if an angle is not [...; 2] or not in an angle unit.

#val rotation : ('a Frame.fixed as 'f) Frame.t -> ('b Frame.fixed as 'g) Frame.t -> (('f, 'd) Direction.t, ('g, 'd) Direction.t, 'd) t

rotation f g maps directions in f to the same directions in g, as Direction.rotate does. Its inverse is rotation g f.

#val about : ('f, 'd) Direction.t -> (('f, 'd) Direction.t, 'd plane, 'd) t

about c maps directions to angular offsets about c, in radians, x east and y north: the inverse of celestial Arc at c with LONPOLE 180°. The offset's norm is Direction.separation from c, and its bearing from +y toward +x is Direction.position_angle. A disc of radius ρ about the origin is exactly the cap of angular radius ρ about c. At a pole the meridian is longitude 0, as Direction.lon takes it.

#val gnomonic : ('f, 'd) Direction.t -> (('f, 'd) Direction.t, 'd plane, 'd) t

gnomonic c is the same with Tan: great circles map to lines.

Operations

#val (>>) : ('a, 'b, 'd) t -> ('b, 'c, 'd) t -> ('a, 'c, 'd) t

t >> u applies t, then u.

#val inverse : ('a, 'b, 'd) t -> ('b, 'a, 'd) t

inverse t reverses t's stages and flips each one's sense. inverse (inverse t) has t's stages.

#val apply : ('a, 'b, 'd) t -> 'a -> 'b

apply t x maps x through each stage. Directions are read as the rays they span and returned as unit vectors.

Raises Invalid_argument where a point is outside a stage's domain, or where a stage reads points of the wrong unit or size.

#val covers : ('a, 'b, 'd) t -> 'a -> 'd Nx.bool_t option

covers t x is Some m, m where apply t x is defined: false outside a stage's domain and at NaN, of x's batch shape. It is None for a transform with no stage, which is defined at every point. It does not raise at any point.

#val pp : Stdlib.Format.formatter -> ('a, 'b, 'd) t -> unit

pp formats a transform's stages.

#val ptree : unit -> (('a, 'b, 'd) t, 'd) Nx.Ptree.t

ptree () is the structure of transforms: the number of stages, then each stage at its index, reporting its family, its sense and its static data (permutation, origin, code, frame, stated terms) and walking its parameters. No float is static.