wcs/
lib.rs

1extern crate mapproj;
2#[macro_use]
3extern crate quick_error;
4
5#[doc = include_str!("../readme.md")]
6pub mod error;
7
8use coo_system::CooSystem;
9use error::Error;
10pub mod coo_system;
11pub mod params;
12mod projection;
13mod sip;
14mod utils;
15
16use crate::projection::WCSCanonicalProjection;
17pub use params::WCSParams;
18
19// Imports
20use mapproj::{
21    conic::{cod::Cod, coe::Coe, coo::Coo, cop::Cop},
22    cylindrical::{car::Car, cea::Cea, cyp::Cyp, mer::Mer},
23    hybrid::hpx::Hpx,
24    img2celestial::Img2Celestial,
25    img2proj::{WcsImgXY2ProjXY, WcsWithSipImgXY2ProjXY},
26    pseudocyl::{ait::Ait, mol::Mol, par::Par, sfl::Sfl},
27    zenithal::{
28        air::Air, arc::Arc, azp::Azp, ncp::Ncp, sin::Sin, stg::Stg, szp::Szp, tan::Tan, zea::Zea,
29        zpn::Zpn,
30    },
31    XYZ,
32};
33
34use paste::paste;
35/// macro
36macro_rules! create_specific_proj {
37    ( $proj_name:ident, $params:expr, $ctype1:expr, $naxis1:expr, $naxis2:expr, $crpix1:expr, $crpix2:expr, $img2proj:expr ) => {{
38        let (positional_angle, proj) = $proj_name::parse_proj(&$params)?;
39
40        let is_sip_found = &$ctype1[($ctype1.len() - 3)..] == "SIP";
41        if is_sip_found {
42            let sip = sip::parse_sip($params, $naxis1, $naxis2, $crpix1, $crpix2)?;
43            let img2proj = WcsWithSipImgXY2ProjXY::new($img2proj, sip);
44
45            paste! {
46                Ok((WCSCelestialProj::[ <$proj_name Sip> ](Img2Celestial::new(img2proj, proj)), positional_angle))
47            }
48        } else {
49            Ok((
50                WCSCelestialProj::$proj_name(Img2Celestial::new($img2proj, proj)),
51                positional_angle,
52            ))
53        }
54    }};
55}
56
57/// Structure alias coming from mapproj defining
58/// image space pixel coordinates
59pub type ImgXY = mapproj::ImgXY;
60/// Structure alias coming from mapproj defining
61/// longitude and latitude expressed in degrees
62pub type LonLat = mapproj::LonLat;
63
64#[derive(Debug)]
65pub struct WCS {
66    /* Metadata keywords */
67    /// Size of the image in pixels in its i-th dimension
68    naxisi: Box<[i64]>,
69    /// Field of view of the image along NAXIS1
70    fov1: f64,
71    /// Field of view of the image along NAXIS2
72    fov2: f64,
73    /// Main sub structure defining the projection
74    proj: WCSProj,
75}
76
77/// Main object structure descripting a WCS object
78/// Once created, the user can proceed two operation on it
79/// * The projection of a (lon, lat) tuple onto the image space.
80///   Results are given in pixels
81/// * The unprojection of a (x, y) tuple given in pixel coordinates onto the sphere.
82///   Results are given as a (lon, lat) tuple expressed in degrees
83impl WCS {
84    pub fn new(params: &WCSParams) -> Result<Self, Error> {
85        // Check for ZNAXISi first (case of tiled compressed image stored in bin tables)
86        let naxisi = match params.naxis {
87            2 => {
88                let naxis1 = params
89                    .znaxis1
90                    .or(params.naxis1)
91                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS1"))?;
92
93                let naxis2 = params
94                    .znaxis2
95                    .or(params.naxis2)
96                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS2"))?;
97
98                Ok(vec![naxis1, naxis2].into_boxed_slice())
99            }
100            3 => {
101                let naxis1 = params
102                    .znaxis1
103                    .or(params.naxis1)
104                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS1"))?;
105                let naxis2 = params
106                    .znaxis2
107                    .or(params.naxis2)
108                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS2"))?;
109                let naxis3 = params
110                    .znaxis3
111                    .or(params.naxis3)
112                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS3"))?;
113
114                Ok(vec![naxis1, naxis2, naxis3].into_boxed_slice())
115            }
116            4 => {
117                let naxis1 = params
118                    .znaxis1
119                    .or(params.naxis1)
120                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS1"))?;
121                let naxis2 = params
122                    .znaxis2
123                    .or(params.naxis2)
124                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS2"))?;
125                let naxis3 = params
126                    .znaxis3
127                    .or(params.naxis3)
128                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS3"))?;
129                let naxis4 = params
130                    .znaxis4
131                    .or(params.naxis4)
132                    .ok_or(Error::MandatoryWCSKeywordsMissing("NAXIS4"))?;
133
134                Ok(vec![naxis1, naxis2, naxis3, naxis4].into_boxed_slice())
135            }
136            _ => Err(Error::NotSupportedNaxis(params.naxis)),
137        }?;
138
139        // At least NAXIS >= 2
140        let proj = WCSProj::new(naxisi[0], naxisi[1], params)?;
141
142        let fov1 = proj.s_lon * (naxisi[0] as f64);
143        let fov2 = proj.s_lat * (naxisi[1] as f64);
144
145        Ok(WCS {
146            naxisi,
147            fov1,
148            fov2,
149            proj,
150        })
151    }
152
153    /// Returns the dimensions of the image given by the NAXIS1 x NAXIS2 keyword
154    pub fn img_dimensions(&self) -> &[i64] {
155        &self.naxisi[..]
156    }
157
158    pub fn field_of_view(&self) -> (f64, f64) {
159        (self.fov1, self.fov2)
160    }
161
162    /// Project a (lon, lat) in ICRS
163    ///
164    /// The result is given a (X, Y) tuple expressed in pixel coordinates.
165    ///
166    /// # Arguments
167    ///
168    /// * `lonlat`: the 3D sphere vertex expressed as a (lon, lat) tuple given in degrees
169    pub fn proj(&self, lonlat: &LonLat) -> Option<ImgXY> {
170        self.proj.proj_lonlat(lonlat)
171    }
172
173    /// Unproject a (X, Y) point to get a position on the sky in ICRS system
174    ///
175    /// # Arguments
176    ///
177    /// * `img_pos`: the image space point expressed as a (X, Y) tuple given en pixels
178    pub fn unproj(&self, img_pos: &ImgXY) -> Option<LonLat> {
179        self.proj.unproj_lonlat(img_pos)
180    }
181
182    /// Get the coordinate system frame
183    pub fn coo_system(&self) -> &CooSystem {
184        self.proj.coo_system()
185    }
186}
187
188use std::ops::Deref;
189impl Deref for WCS {
190    type Target = WCSProj;
191
192    fn deref(&self) -> &Self::Target {
193        &self.proj
194    }
195}
196
197pub struct WCSProj {
198    /// The right part of the CTYPE keyword
199    /// The projection type
200    proj: WCSCelestialProj,
201    /// The left part of the CTYPE keyword
202    /// The coordinate system
203    coo_system: CooSystem,
204    pos_angle: f64,
205    s_lon: f64, // scale in degrees along the longitude axis
206    s_lat: f64, // scale in degrees along the latitude axis
207}
208
209impl std::fmt::Debug for WCSProj {
210    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
211        f.debug_struct("WCS")
212            .field("coosys", &self.coo_system)
213            .field("euler angle", &self.pos_angle)
214            .finish()
215    }
216}
217
218/// Main enum structure descripting a WCS object
219/// Once created, the user can proceed two operation on it
220/// * The projection of a (lon, lat) tuple onto the image space.
221///   Results are given in pixels
222/// * The unprojection of a (x, y) tuple given in pixel coordinates onto the sphere.
223///   Results are given as a (lon, lat) tuple expressed in degrees
224pub enum WCSCelestialProj {
225    // Zenithal
226    Azp(Img2Celestial<Azp, WcsImgXY2ProjXY>),
227    Szp(Img2Celestial<Szp, WcsImgXY2ProjXY>),
228    Tan(Img2Celestial<Tan, WcsImgXY2ProjXY>),
229    Stg(Img2Celestial<Stg, WcsImgXY2ProjXY>),
230    Sin(Img2Celestial<Sin, WcsImgXY2ProjXY>),
231    Arc(Img2Celestial<Arc, WcsImgXY2ProjXY>),
232    Zpn(Img2Celestial<Zpn, WcsImgXY2ProjXY>),
233    Zea(Img2Celestial<Zea, WcsImgXY2ProjXY>),
234    Air(Img2Celestial<Air, WcsImgXY2ProjXY>),
235    Ncp(Img2Celestial<Ncp, WcsImgXY2ProjXY>),
236    // Cylindrical
237    Cyp(Img2Celestial<Cyp, WcsImgXY2ProjXY>),
238    Cea(Img2Celestial<Cea, WcsImgXY2ProjXY>),
239    Car(Img2Celestial<Car, WcsImgXY2ProjXY>),
240    Mer(Img2Celestial<Mer, WcsImgXY2ProjXY>),
241    // Pseudo-Cylindrical
242    Sfl(Img2Celestial<Sfl, WcsImgXY2ProjXY>),
243    Par(Img2Celestial<Par, WcsImgXY2ProjXY>),
244    Mol(Img2Celestial<Mol, WcsImgXY2ProjXY>),
245    Ait(Img2Celestial<Ait, WcsImgXY2ProjXY>),
246    // Conic
247    Cop(Img2Celestial<Cop, WcsImgXY2ProjXY>),
248    Cod(Img2Celestial<Cod, WcsImgXY2ProjXY>),
249    Coe(Img2Celestial<Coe, WcsImgXY2ProjXY>),
250    Coo(Img2Celestial<Coo, WcsImgXY2ProjXY>),
251    // Hybrid
252    Hpx(Img2Celestial<Hpx, WcsImgXY2ProjXY>),
253
254    // SIP handling
255    // Zenithal
256    AzpSip(Img2Celestial<Azp, WcsWithSipImgXY2ProjXY>),
257    SzpSip(Img2Celestial<Szp, WcsWithSipImgXY2ProjXY>),
258    TanSip(Img2Celestial<Tan, WcsWithSipImgXY2ProjXY>),
259    StgSip(Img2Celestial<Stg, WcsWithSipImgXY2ProjXY>),
260    SinSip(Img2Celestial<Sin, WcsWithSipImgXY2ProjXY>),
261    ArcSip(Img2Celestial<Arc, WcsWithSipImgXY2ProjXY>),
262    ZpnSip(Img2Celestial<Zpn, WcsWithSipImgXY2ProjXY>),
263    ZeaSip(Img2Celestial<Zea, WcsWithSipImgXY2ProjXY>),
264    AirSip(Img2Celestial<Air, WcsWithSipImgXY2ProjXY>),
265    NcpSip(Img2Celestial<Ncp, WcsWithSipImgXY2ProjXY>),
266    // Cylindrical
267    CypSip(Img2Celestial<Cyp, WcsWithSipImgXY2ProjXY>),
268    CeaSip(Img2Celestial<Cea, WcsWithSipImgXY2ProjXY>),
269    CarSip(Img2Celestial<Car, WcsWithSipImgXY2ProjXY>),
270    MerSip(Img2Celestial<Mer, WcsWithSipImgXY2ProjXY>),
271    // Pseudo-Cylindrical
272    SflSip(Img2Celestial<Sfl, WcsWithSipImgXY2ProjXY>),
273    ParSip(Img2Celestial<Par, WcsWithSipImgXY2ProjXY>),
274    MolSip(Img2Celestial<Mol, WcsWithSipImgXY2ProjXY>),
275    AitSip(Img2Celestial<Ait, WcsWithSipImgXY2ProjXY>),
276    // Conic
277    CopSip(Img2Celestial<Cop, WcsWithSipImgXY2ProjXY>),
278    CodSip(Img2Celestial<Cod, WcsWithSipImgXY2ProjXY>),
279    CoeSip(Img2Celestial<Coe, WcsWithSipImgXY2ProjXY>),
280    CooSip(Img2Celestial<Coo, WcsWithSipImgXY2ProjXY>),
281    // Hybrid
282    HpxSip(Img2Celestial<Hpx, WcsWithSipImgXY2ProjXY>),
283}
284
285fn parse_pc_matrix(params: &WCSParams) -> Option<(f64, f64, f64, f64)> {
286    let pc11 = params.pc1_1;
287    let pc12 = params.pc1_2;
288    let pc21 = params.pc2_1;
289    let pc22 = params.pc2_2;
290
291    let pc_matrix_found = match (&pc11, &pc12, &pc21, &pc22) {
292        (None, None, None, None) => false,
293        // The CD1_1 keyword has been found
294        // We are in a case where the CDij are given
295        _ => true,
296    };
297
298    if pc_matrix_found {
299        let pc11 = pc11.unwrap_or(1.0);
300        let pc12 = pc12.unwrap_or(0.0);
301        let pc21 = pc21.unwrap_or(0.0);
302        let pc22 = pc22.unwrap_or(1.0);
303
304        Some((pc11, pc12, pc21, pc22))
305    } else {
306        None
307    }
308}
309
310fn parse_cd_matrix(params: &WCSParams) -> Option<(f64, f64, f64, f64)> {
311    let cd11 = params.cd1_1;
312    let cd12 = params.cd1_2;
313    let cd21 = params.cd2_1;
314    let cd22 = params.cd2_2;
315
316    let cd_matrix_found = match (&cd11, &cd12, &cd21, &cd22) {
317        (None, None, None, None) => false,
318        // The CD1_1 keyword has been found
319        // We are in a case where the CDij are given
320        _ => true,
321    };
322
323    if cd_matrix_found {
324        let cd11 = cd11.unwrap_or(1.0);
325        let cd12 = cd12.unwrap_or(0.0);
326        let cd21 = cd21.unwrap_or(0.0);
327        let cd22 = cd22.unwrap_or(1.0);
328
329        Some((cd11, cd12, cd21, cd22))
330    } else {
331        None
332    }
333}
334
335impl WCSProj {
336    /// Create a WCS from a specific fits header parsed with fitsrs
337    /// # Param
338    /// * `naxis1` - Size of the image in its first dimension (in pixels)
339    /// * `naxis2` - Size of the image in its second dimension (in pixels)
340    /// * `params` - Header unit coming from fitsrs.
341    ///   This contains all the cards of one HDU.
342    pub fn new(naxis1: i64, naxis2: i64, params: &WCSParams) -> Result<Self, Error> {
343        // 1. Identify the image <-> intermediate projection
344        // a. Linear transformation matrix cases:
345        // - CRPIXi + CDij
346        // - CRPIXi + CDELTi + CROTA2
347        // - CRPIXi + CDELTi + PCij
348        let crpix1 = params.crpix1.unwrap_or(0.0);
349        let crpix2 = params.crpix2.unwrap_or(0.0);
350
351        // Choice of the wcs order:
352        // 1 - Priority to define the projection is given to CD
353        // 2 - Then, to the couple PC + CDELT
354        // 3 - Finally to the old CROTA + CDELT convention
355        let (img2proj, s_lon, s_lat) =
356            if let Some((cd11, cd12, cd21, cd22)) = parse_cd_matrix(params) {
357                // CDij case
358                (
359                    WcsImgXY2ProjXY::from_cd(crpix1, crpix2, cd11, cd12, cd21, cd22),
360                    cd11.abs(),
361                    cd22.abs(),
362                )
363            } else {
364                // Search for CDELTi
365                let cdelt1 = params.cdelt1.unwrap_or(1.0);
366                let cdelt2 = params.cdelt2.unwrap_or(1.0);
367
368                if let Some((pc11, pc12, pc21, pc22)) = parse_pc_matrix(params) {
369                    // CDELTi + PCij case
370                    (
371                        WcsImgXY2ProjXY::from_pc(
372                            crpix1, crpix2, pc11, pc12, pc21, pc22, cdelt1, cdelt2,
373                        ),
374                        (cdelt1 * pc11).abs(),
375                        (cdelt2 * pc22).abs(),
376                    )
377                } else {
378                    // CDELTi + CROTA2 case
379                    let crota2 = params.crota2.unwrap_or(0.0);
380                    let cosc = crota2.to_radians().cos();
381
382                    (
383                        WcsImgXY2ProjXY::from_cr(crpix1, crpix2, crota2, cdelt1, cdelt2),
384                        (cdelt1 * cosc).abs(),
385                        (cdelt2 * cosc).abs(),
386                    )
387                }
388            };
389
390        // 2. Identify the projection type
391        let ctype1 = &params.ctype1;
392        let proj_name = &ctype1[5..=7];
393
394        let (proj, pos_angle) = match proj_name.as_bytes() {
395            // Zenithal
396            b"AZP" => {
397                create_specific_proj!(Azp, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
398            }
399            b"SZP" => {
400                create_specific_proj!(Szp, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
401            }
402            b"TAN" => {
403                create_specific_proj!(Tan, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
404            }
405            b"STG" => {
406                create_specific_proj!(Stg, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
407            }
408            b"SIN" => {
409                create_specific_proj!(Sin, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
410            }
411            b"ARC" => {
412                create_specific_proj!(Arc, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
413            }
414            b"ZPN" => {
415                create_specific_proj!(Zpn, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
416            }
417            b"ZEA" => {
418                create_specific_proj!(Zea, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
419            }
420            b"AIR" => {
421                create_specific_proj!(Air, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
422            }
423            b"NCP" => {
424                create_specific_proj!(Ncp, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
425            }
426            // Cylindrical
427            b"CYP" => {
428                create_specific_proj!(Cyp, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
429            }
430            b"CEA" => {
431                create_specific_proj!(Cea, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
432            }
433            b"CAR" => {
434                create_specific_proj!(Car, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
435            }
436            b"MER" => {
437                create_specific_proj!(Mer, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
438            }
439            // Pseudo-cylindrical
440            b"SFL" => {
441                create_specific_proj!(Sfl, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
442            }
443            b"PAR" => {
444                create_specific_proj!(Par, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
445            }
446            b"MOL" => {
447                create_specific_proj!(Mol, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
448            }
449            b"AIT" => {
450                create_specific_proj!(Ait, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
451            }
452            // Conic
453            b"COP" => {
454                create_specific_proj!(Cop, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
455            }
456            b"COD" => {
457                create_specific_proj!(Cod, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
458            }
459            b"COE" => {
460                create_specific_proj!(Coe, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
461            }
462            b"COO" => {
463                create_specific_proj!(Coo, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
464            }
465            // HEALPix
466            b"HPX" => {
467                create_specific_proj!(Hpx, params, ctype1, naxis1, naxis2, crpix1, crpix2, img2proj)
468            }
469            _ => Err(Error::NotImplementedProjection(proj_name.to_string())),
470        }?;
471
472        let coo_system = CooSystem::parse(&params)?;
473
474        Ok(WCSProj {
475            proj,
476            coo_system,
477            pos_angle,
478            s_lon,
479            s_lat,
480        })
481    }
482
483    /// Project a (lon, lat) given in ICRS frame to get its corresponding location on the image
484    ///
485    /// The result is given a (X, Y) tuple expressed in pixel coordinates.
486    ///
487    /// # Arguments
488    ///
489    /// * `lonlat`: a coo expressed as (lon, lat) tuple given in degrees and in ICRS system
490    pub fn proj_lonlat(&self, lonlat: &LonLat) -> Option<ImgXY> {
491        let lonlat = &self.coo_system.from_icrs(lonlat.clone());
492
493        match &self.proj {
494            // Zenithal
495            WCSCelestialProj::Azp(wcs) => wcs.lonlat2img(lonlat),
496            WCSCelestialProj::Szp(wcs) => wcs.lonlat2img(lonlat),
497            WCSCelestialProj::Tan(wcs) => wcs.lonlat2img(lonlat),
498            WCSCelestialProj::Stg(wcs) => wcs.lonlat2img(lonlat),
499            WCSCelestialProj::Sin(wcs) => wcs.lonlat2img(lonlat),
500            WCSCelestialProj::Arc(wcs) => wcs.lonlat2img(lonlat),
501            WCSCelestialProj::Zpn(wcs) => wcs.lonlat2img(lonlat),
502            WCSCelestialProj::Zea(wcs) => wcs.lonlat2img(lonlat),
503            WCSCelestialProj::Air(wcs) => wcs.lonlat2img(lonlat),
504            WCSCelestialProj::Ncp(wcs) => wcs.lonlat2img(lonlat),
505            // Pseudo-cyl
506            WCSCelestialProj::Cyp(wcs) => wcs.lonlat2img(lonlat),
507            WCSCelestialProj::Cea(wcs) => wcs.lonlat2img(lonlat),
508            WCSCelestialProj::Car(wcs) => wcs.lonlat2img(lonlat),
509            WCSCelestialProj::Mer(wcs) => wcs.lonlat2img(lonlat),
510            // Cylindrical
511            WCSCelestialProj::Sfl(wcs) => wcs.lonlat2img(lonlat),
512            WCSCelestialProj::Par(wcs) => wcs.lonlat2img(lonlat),
513            WCSCelestialProj::Mol(wcs) => wcs.lonlat2img(lonlat),
514            WCSCelestialProj::Ait(wcs) => wcs.lonlat2img(lonlat),
515            // Conic
516            WCSCelestialProj::Cop(wcs) => wcs.lonlat2img(lonlat),
517            WCSCelestialProj::Cod(wcs) => wcs.lonlat2img(lonlat),
518            WCSCelestialProj::Coe(wcs) => wcs.lonlat2img(lonlat),
519            WCSCelestialProj::Coo(wcs) => wcs.lonlat2img(lonlat),
520            // Hybrid
521            WCSCelestialProj::Hpx(wcs) => wcs.lonlat2img(lonlat),
522
523            /* Sip variants */
524            // Zenithal
525            WCSCelestialProj::AzpSip(wcs) => wcs.lonlat2img(lonlat),
526            WCSCelestialProj::SzpSip(wcs) => wcs.lonlat2img(lonlat),
527            WCSCelestialProj::TanSip(wcs) => wcs.lonlat2img(lonlat),
528            WCSCelestialProj::StgSip(wcs) => wcs.lonlat2img(lonlat),
529            WCSCelestialProj::SinSip(wcs) => wcs.lonlat2img(lonlat),
530            WCSCelestialProj::ArcSip(wcs) => wcs.lonlat2img(lonlat),
531            WCSCelestialProj::ZpnSip(wcs) => wcs.lonlat2img(lonlat),
532            WCSCelestialProj::ZeaSip(wcs) => wcs.lonlat2img(lonlat),
533            WCSCelestialProj::AirSip(wcs) => wcs.lonlat2img(lonlat),
534            WCSCelestialProj::NcpSip(wcs) => wcs.lonlat2img(lonlat),
535            // Pseudo-cyl
536            WCSCelestialProj::CypSip(wcs) => wcs.lonlat2img(lonlat),
537            WCSCelestialProj::CeaSip(wcs) => wcs.lonlat2img(lonlat),
538            WCSCelestialProj::CarSip(wcs) => wcs.lonlat2img(lonlat),
539            WCSCelestialProj::MerSip(wcs) => wcs.lonlat2img(lonlat),
540            // Cylindrical
541            WCSCelestialProj::SflSip(wcs) => wcs.lonlat2img(lonlat),
542            WCSCelestialProj::ParSip(wcs) => wcs.lonlat2img(lonlat),
543            WCSCelestialProj::MolSip(wcs) => wcs.lonlat2img(lonlat),
544            WCSCelestialProj::AitSip(wcs) => wcs.lonlat2img(lonlat),
545            // Conic
546            WCSCelestialProj::CopSip(wcs) => wcs.lonlat2img(lonlat),
547            WCSCelestialProj::CodSip(wcs) => wcs.lonlat2img(lonlat),
548            WCSCelestialProj::CoeSip(wcs) => wcs.lonlat2img(lonlat),
549            WCSCelestialProj::CooSip(wcs) => wcs.lonlat2img(lonlat),
550            // Hybrid
551            WCSCelestialProj::HpxSip(wcs) => wcs.lonlat2img(lonlat),
552        }
553    }
554
555    pub fn proj_xyz(&self, xyz: &(f64, f64, f64)) -> Option<ImgXY> {
556        let xyz = &self.coo_system.from_icrs_xyz(XYZ::new(xyz.0, xyz.1, xyz.2));
557
558        match &self.proj {
559            // Zenithal
560            WCSCelestialProj::Azp(wcs) => wcs.xyz2img(xyz),
561            WCSCelestialProj::Szp(wcs) => wcs.xyz2img(xyz),
562            WCSCelestialProj::Tan(wcs) => wcs.xyz2img(xyz),
563            WCSCelestialProj::Stg(wcs) => wcs.xyz2img(xyz),
564            WCSCelestialProj::Sin(wcs) => wcs.xyz2img(xyz),
565            WCSCelestialProj::Arc(wcs) => wcs.xyz2img(xyz),
566            WCSCelestialProj::Zpn(wcs) => wcs.xyz2img(xyz),
567            WCSCelestialProj::Zea(wcs) => wcs.xyz2img(xyz),
568            WCSCelestialProj::Air(wcs) => wcs.xyz2img(xyz),
569            WCSCelestialProj::Ncp(wcs) => wcs.xyz2img(xyz),
570            // Pseudo-cyl
571            WCSCelestialProj::Cyp(wcs) => wcs.xyz2img(xyz),
572            WCSCelestialProj::Cea(wcs) => wcs.xyz2img(xyz),
573            WCSCelestialProj::Car(wcs) => wcs.xyz2img(xyz),
574            WCSCelestialProj::Mer(wcs) => wcs.xyz2img(xyz),
575            // Cylindrical
576            WCSCelestialProj::Sfl(wcs) => wcs.xyz2img(xyz),
577            WCSCelestialProj::Par(wcs) => wcs.xyz2img(xyz),
578            WCSCelestialProj::Mol(wcs) => wcs.xyz2img(xyz),
579            WCSCelestialProj::Ait(wcs) => wcs.xyz2img(xyz),
580            // Conic
581            WCSCelestialProj::Cop(wcs) => wcs.xyz2img(xyz),
582            WCSCelestialProj::Cod(wcs) => wcs.xyz2img(xyz),
583            WCSCelestialProj::Coe(wcs) => wcs.xyz2img(xyz),
584            WCSCelestialProj::Coo(wcs) => wcs.xyz2img(xyz),
585            // Hybrid
586            WCSCelestialProj::Hpx(wcs) => wcs.xyz2img(xyz),
587
588            /* Sip variants */
589            // Zenithal
590            WCSCelestialProj::AzpSip(wcs) => wcs.xyz2img(xyz),
591            WCSCelestialProj::SzpSip(wcs) => wcs.xyz2img(xyz),
592            WCSCelestialProj::TanSip(wcs) => wcs.xyz2img(xyz),
593            WCSCelestialProj::StgSip(wcs) => wcs.xyz2img(xyz),
594            WCSCelestialProj::SinSip(wcs) => wcs.xyz2img(xyz),
595            WCSCelestialProj::ArcSip(wcs) => wcs.xyz2img(xyz),
596            WCSCelestialProj::ZpnSip(wcs) => wcs.xyz2img(xyz),
597            WCSCelestialProj::ZeaSip(wcs) => wcs.xyz2img(xyz),
598            WCSCelestialProj::AirSip(wcs) => wcs.xyz2img(xyz),
599            WCSCelestialProj::NcpSip(wcs) => wcs.xyz2img(xyz),
600            // Pseudo-cyl
601            WCSCelestialProj::CypSip(wcs) => wcs.xyz2img(xyz),
602            WCSCelestialProj::CeaSip(wcs) => wcs.xyz2img(xyz),
603            WCSCelestialProj::CarSip(wcs) => wcs.xyz2img(xyz),
604            WCSCelestialProj::MerSip(wcs) => wcs.xyz2img(xyz),
605            // Cylindrical
606            WCSCelestialProj::SflSip(wcs) => wcs.xyz2img(xyz),
607            WCSCelestialProj::ParSip(wcs) => wcs.xyz2img(xyz),
608            WCSCelestialProj::MolSip(wcs) => wcs.xyz2img(xyz),
609            WCSCelestialProj::AitSip(wcs) => wcs.xyz2img(xyz),
610            // Conic
611            WCSCelestialProj::CopSip(wcs) => wcs.xyz2img(xyz),
612            WCSCelestialProj::CodSip(wcs) => wcs.xyz2img(xyz),
613            WCSCelestialProj::CoeSip(wcs) => wcs.xyz2img(xyz),
614            WCSCelestialProj::CooSip(wcs) => wcs.xyz2img(xyz),
615            // Hybrid
616            WCSCelestialProj::HpxSip(wcs) => wcs.xyz2img(xyz),
617        }
618    }
619
620    pub fn unproj_xyz(&self, img_pos: &ImgXY) -> Option<XYZ> {
621        let xyz = match &self.proj {
622            // Zenithal
623            WCSCelestialProj::Azp(wcs) => wcs.img2xyz(&img_pos),
624            WCSCelestialProj::Szp(wcs) => wcs.img2xyz(&img_pos),
625            WCSCelestialProj::Tan(wcs) => wcs.img2xyz(&img_pos),
626            WCSCelestialProj::Stg(wcs) => wcs.img2xyz(&img_pos),
627            WCSCelestialProj::Sin(wcs) => wcs.img2xyz(&img_pos),
628            WCSCelestialProj::Arc(wcs) => wcs.img2xyz(&img_pos),
629            WCSCelestialProj::Zpn(wcs) => wcs.img2xyz(&img_pos),
630            WCSCelestialProj::Zea(wcs) => wcs.img2xyz(&img_pos),
631            WCSCelestialProj::Air(wcs) => wcs.img2xyz(&img_pos),
632            WCSCelestialProj::Ncp(wcs) => wcs.img2xyz(&img_pos),
633            // Pseudo-cyl
634            WCSCelestialProj::Cyp(wcs) => wcs.img2xyz(&img_pos),
635            WCSCelestialProj::Cea(wcs) => wcs.img2xyz(&img_pos),
636            WCSCelestialProj::Car(wcs) => wcs.img2xyz(&img_pos),
637            WCSCelestialProj::Mer(wcs) => wcs.img2xyz(&img_pos),
638            // Cylindrical
639            WCSCelestialProj::Sfl(wcs) => wcs.img2xyz(&img_pos),
640            WCSCelestialProj::Par(wcs) => wcs.img2xyz(&img_pos),
641            WCSCelestialProj::Mol(wcs) => wcs.img2xyz(&img_pos),
642            WCSCelestialProj::Ait(wcs) => wcs.img2xyz(&img_pos),
643            // Conic
644            WCSCelestialProj::Cop(wcs) => wcs.img2xyz(&img_pos),
645            WCSCelestialProj::Cod(wcs) => wcs.img2xyz(&img_pos),
646            WCSCelestialProj::Coe(wcs) => wcs.img2xyz(&img_pos),
647            WCSCelestialProj::Coo(wcs) => wcs.img2xyz(&img_pos),
648            // Hybrid
649            WCSCelestialProj::Hpx(wcs) => wcs.img2xyz(&img_pos),
650
651            /* Sip variants */
652            // Zenithal
653            WCSCelestialProj::AzpSip(wcs) => wcs.img2xyz(&img_pos),
654            WCSCelestialProj::SzpSip(wcs) => wcs.img2xyz(&img_pos),
655            WCSCelestialProj::TanSip(wcs) => wcs.img2xyz(&img_pos),
656            WCSCelestialProj::StgSip(wcs) => wcs.img2xyz(&img_pos),
657            WCSCelestialProj::SinSip(wcs) => wcs.img2xyz(&img_pos),
658            WCSCelestialProj::ArcSip(wcs) => wcs.img2xyz(&img_pos),
659            WCSCelestialProj::ZpnSip(wcs) => wcs.img2xyz(&img_pos),
660            WCSCelestialProj::ZeaSip(wcs) => wcs.img2xyz(&img_pos),
661            WCSCelestialProj::AirSip(wcs) => wcs.img2xyz(&img_pos),
662            WCSCelestialProj::NcpSip(wcs) => wcs.img2xyz(&img_pos),
663            // Pseudo-cyl
664            WCSCelestialProj::CypSip(wcs) => wcs.img2xyz(&img_pos),
665            WCSCelestialProj::CeaSip(wcs) => wcs.img2xyz(&img_pos),
666            WCSCelestialProj::CarSip(wcs) => wcs.img2xyz(&img_pos),
667            WCSCelestialProj::MerSip(wcs) => wcs.img2xyz(&img_pos),
668            // Cylindrical
669            WCSCelestialProj::SflSip(wcs) => wcs.img2xyz(&img_pos),
670            WCSCelestialProj::ParSip(wcs) => wcs.img2xyz(&img_pos),
671            WCSCelestialProj::MolSip(wcs) => wcs.img2xyz(&img_pos),
672            WCSCelestialProj::AitSip(wcs) => wcs.img2xyz(&img_pos),
673            // Conic
674            WCSCelestialProj::CopSip(wcs) => wcs.img2xyz(&img_pos),
675            WCSCelestialProj::CodSip(wcs) => wcs.img2xyz(&img_pos),
676            WCSCelestialProj::CoeSip(wcs) => wcs.img2xyz(&img_pos),
677            WCSCelestialProj::CooSip(wcs) => wcs.img2xyz(&img_pos),
678            // Hybrid
679            WCSCelestialProj::HpxSip(wcs) => wcs.img2xyz(&img_pos),
680        };
681        xyz.map(|v| self.coo_system.to_icrs_xyz(v))
682    }
683
684    /// Unproject a (X, Y) point from the image space to get its corresponding location on the sphere
685    ///
686    /// The result is (lon, lat) tuple expressed in degrees in ICRS
687    ///
688    /// # Arguments
689    ///
690    /// * `img_pos`: the image space point expressed as a (X, Y) tuple given en pixels
691    pub fn unproj_lonlat(&self, img_pos: &ImgXY) -> Option<LonLat> {
692        let lonlat = match &self.proj {
693            // Zenithal
694            WCSCelestialProj::Azp(wcs) => wcs.img2lonlat(&img_pos),
695            WCSCelestialProj::Szp(wcs) => wcs.img2lonlat(&img_pos),
696            WCSCelestialProj::Tan(wcs) => wcs.img2lonlat(&img_pos),
697            WCSCelestialProj::Stg(wcs) => wcs.img2lonlat(&img_pos),
698            WCSCelestialProj::Sin(wcs) => wcs.img2lonlat(&img_pos),
699            WCSCelestialProj::Arc(wcs) => wcs.img2lonlat(&img_pos),
700            WCSCelestialProj::Zpn(wcs) => wcs.img2lonlat(&img_pos),
701            WCSCelestialProj::Zea(wcs) => wcs.img2lonlat(&img_pos),
702            WCSCelestialProj::Air(wcs) => wcs.img2lonlat(&img_pos),
703            WCSCelestialProj::Ncp(wcs) => wcs.img2lonlat(&img_pos),
704            // Pseudo-cyl
705            WCSCelestialProj::Cyp(wcs) => wcs.img2lonlat(&img_pos),
706            WCSCelestialProj::Cea(wcs) => wcs.img2lonlat(&img_pos),
707            WCSCelestialProj::Car(wcs) => wcs.img2lonlat(&img_pos),
708            WCSCelestialProj::Mer(wcs) => wcs.img2lonlat(&img_pos),
709            // Cylindrical
710            WCSCelestialProj::Sfl(wcs) => wcs.img2lonlat(&img_pos),
711            WCSCelestialProj::Par(wcs) => wcs.img2lonlat(&img_pos),
712            WCSCelestialProj::Mol(wcs) => wcs.img2lonlat(&img_pos),
713            WCSCelestialProj::Ait(wcs) => wcs.img2lonlat(&img_pos),
714            // Conic
715            WCSCelestialProj::Cop(wcs) => wcs.img2lonlat(&img_pos),
716            WCSCelestialProj::Cod(wcs) => wcs.img2lonlat(&img_pos),
717            WCSCelestialProj::Coe(wcs) => wcs.img2lonlat(&img_pos),
718            WCSCelestialProj::Coo(wcs) => wcs.img2lonlat(&img_pos),
719            // Hybrid
720            WCSCelestialProj::Hpx(wcs) => wcs.img2lonlat(&img_pos),
721
722            /* Sip variants */
723            // Zenithal
724            WCSCelestialProj::AzpSip(wcs) => wcs.img2lonlat(&img_pos),
725            WCSCelestialProj::SzpSip(wcs) => wcs.img2lonlat(&img_pos),
726            WCSCelestialProj::TanSip(wcs) => wcs.img2lonlat(&img_pos),
727            WCSCelestialProj::StgSip(wcs) => wcs.img2lonlat(&img_pos),
728            WCSCelestialProj::SinSip(wcs) => wcs.img2lonlat(&img_pos),
729            WCSCelestialProj::ArcSip(wcs) => wcs.img2lonlat(&img_pos),
730            WCSCelestialProj::ZpnSip(wcs) => wcs.img2lonlat(&img_pos),
731            WCSCelestialProj::ZeaSip(wcs) => wcs.img2lonlat(&img_pos),
732            WCSCelestialProj::AirSip(wcs) => wcs.img2lonlat(&img_pos),
733            WCSCelestialProj::NcpSip(wcs) => wcs.img2lonlat(&img_pos),
734            // Pseudo-cyl
735            WCSCelestialProj::CypSip(wcs) => wcs.img2lonlat(&img_pos),
736            WCSCelestialProj::CeaSip(wcs) => wcs.img2lonlat(&img_pos),
737            WCSCelestialProj::CarSip(wcs) => wcs.img2lonlat(&img_pos),
738            WCSCelestialProj::MerSip(wcs) => wcs.img2lonlat(&img_pos),
739            // Cylindrical
740            WCSCelestialProj::SflSip(wcs) => wcs.img2lonlat(&img_pos),
741            WCSCelestialProj::ParSip(wcs) => wcs.img2lonlat(&img_pos),
742            WCSCelestialProj::MolSip(wcs) => wcs.img2lonlat(&img_pos),
743            WCSCelestialProj::AitSip(wcs) => wcs.img2lonlat(&img_pos),
744            // Conic
745            WCSCelestialProj::CopSip(wcs) => wcs.img2lonlat(&img_pos),
746            WCSCelestialProj::CodSip(wcs) => wcs.img2lonlat(&img_pos),
747            WCSCelestialProj::CoeSip(wcs) => wcs.img2lonlat(&img_pos),
748            WCSCelestialProj::CooSip(wcs) => wcs.img2lonlat(&img_pos),
749            // Hybrid
750            WCSCelestialProj::HpxSip(wcs) => wcs.img2lonlat(&img_pos),
751        };
752
753        lonlat.map(|ll| self.coo_system.to_icrs(ll))
754    }
755
756    /// Getter of the coordinate system
757    pub fn coo_system(&self) -> &CooSystem {
758        &self.coo_system
759    }
760}
761
762#[cfg(test)]
763mod tests {
764    use super::WCS;
765    use crate::Error;
766    use crate::WCSParams;
767
768    use crate::mapproj::Projection;
769
770    use fitsrs::card::Value;
771    use fitsrs::fits::Fits;
772    use fitsrs::hdu::header::{extension::image::Image, Header};
773    use fitsrs::Pixels;
774
775    use glob::glob;
776    use mapproj::{CanonicalProjection, ImgXY, LonLat};
777    use serde::Deserialize;
778    use std::f64::consts::PI;
779    use std::fs::File;
780    use std::io::BufReader;
781
782    use std::convert::TryFrom;
783
784    use crate::Error::MandatoryWCSKeywordsMissing;
785    fn parse_card<'de, T: Deserialize<'de>>(
786        header: &'de Header<Image>,
787        key: &'static str,
788    ) -> Result<T, Error> {
789        header
790            .get_parsed::<T>(key)
791            .map_err(|_| MandatoryWCSKeywordsMissing("Card cannot be parsed"))
792    }
793
794    fn parse_opt_card<'de, T: Deserialize<'de>>(
795        header: &'de Header<Image>,
796        key: &'static str,
797    ) -> Option<T> {
798        header.get_parsed::<T>(key).ok()
799    }
800
801    impl<'a> TryFrom<&'a Header<Image>> for WCS {
802        type Error = Error;
803
804        fn try_from(h: &'a Header<Image>) -> Result<Self, Self::Error> {
805            let params = WCSParams {
806                naxis: parse_card::<i64>(h, "NAXIS").unwrap(),
807                ctype1: parse_card::<String>(h, "CTYPE1").unwrap(),
808
809                naxis1: parse_opt_card::<i64>(h, "NAXIS1"),
810                naxis2: parse_opt_card::<i64>(h, "NAXIS2"),
811
812                znaxis1: parse_opt_card::<i64>(h, "ZNAXIS1"),
813                znaxis2: parse_opt_card::<i64>(h, "ZNAXIS2"),
814                znaxis3: parse_opt_card::<i64>(h, "ZNAXIS3"),
815                znaxis4: parse_opt_card::<i64>(h, "ZNAXIS4"),
816
817                ctype2: parse_opt_card::<String>(h, "CTYPE2"),
818                ctype3: parse_opt_card::<String>(h, "CTYPE3"),
819
820                a_order: parse_opt_card::<i64>(h, "A_ORDER"),
821                b_order: parse_opt_card::<i64>(h, "B_ORDER"),
822                ap_order: parse_opt_card::<i64>(h, "AP_ORDER"),
823                bp_order: parse_opt_card::<i64>(h, "BP_ORDER"),
824                crpix1: parse_opt_card::<f64>(h, "CRPIX1"),
825                crpix2: parse_opt_card::<f64>(h, "CRPIX2"),
826                crpix3: parse_opt_card::<f64>(h, "CRPIX3"),
827                crval1: parse_opt_card::<f64>(h, "CRVAL1"),
828                crval2: parse_opt_card::<f64>(h, "CRVAL2"),
829                crval3: parse_opt_card::<f64>(h, "CRVAL3"),
830                crota1: parse_opt_card::<f64>(h, "CROTA1"),
831                crota2: parse_opt_card::<f64>(h, "CROTA2"),
832                crota3: parse_opt_card::<f64>(h, "CROTA3"),
833                cdelt1: parse_opt_card::<f64>(h, "CDELT1"),
834                cdelt2: parse_opt_card::<f64>(h, "CDELT2"),
835                cdelt3: parse_opt_card::<f64>(h, "CDELT3"),
836                naxis3: parse_opt_card::<i64>(h, "NAXIS3"),
837                naxis4: parse_opt_card::<i64>(h, "NAXIS4"),
838                lonpole: parse_opt_card::<f64>(h, "LONPOLE"),
839                latpole: parse_opt_card::<f64>(h, "LATPOLE"),
840                equinox: parse_opt_card::<f64>(h, "EQUINOX"),
841                epoch: parse_opt_card::<f64>(h, "EPOCH"),
842                radesys: parse_opt_card::<String>(h, "RADESYS"),
843                pv1_0: parse_opt_card::<f64>(h, "PV1_0"),
844                pv1_1: parse_opt_card::<f64>(h, "PV1_1"),
845                pv1_2: parse_opt_card::<f64>(h, "PV1_2"),
846                pv2_0: parse_opt_card::<f64>(h, "PV2_0"),
847                pv2_1: parse_opt_card::<f64>(h, "PV2_1"),
848                pv2_2: parse_opt_card::<f64>(h, "PV2_2"),
849                pv2_3: parse_opt_card::<f64>(h, "PV2_3"),
850                pv2_4: parse_opt_card::<f64>(h, "PV2_4"),
851                pv2_5: parse_opt_card::<f64>(h, "PV2_5"),
852                pv2_6: parse_opt_card::<f64>(h, "PV2_6"),
853                pv2_7: parse_opt_card::<f64>(h, "PV2_7"),
854                pv2_8: parse_opt_card::<f64>(h, "PV2_8"),
855                pv2_9: parse_opt_card::<f64>(h, "PV2_9"),
856                pv2_10: parse_opt_card::<f64>(h, "PV2_10"),
857                pv2_11: parse_opt_card::<f64>(h, "PV2_11"),
858                pv2_12: parse_opt_card::<f64>(h, "PV2_12"),
859                pv2_13: parse_opt_card::<f64>(h, "PV2_13"),
860                pv2_14: parse_opt_card::<f64>(h, "PV2_14"),
861                pv2_15: parse_opt_card::<f64>(h, "PV2_15"),
862                pv2_16: parse_opt_card::<f64>(h, "PV2_16"),
863                pv2_17: parse_opt_card::<f64>(h, "PV2_17"),
864                pv2_18: parse_opt_card::<f64>(h, "PV2_18"),
865                pv2_19: parse_opt_card::<f64>(h, "PV2_19"),
866                pv2_20: parse_opt_card::<f64>(h, "PV2_20"),
867                cd1_1: parse_opt_card::<f64>(h, "CD1_1"),
868                cd1_2: parse_opt_card::<f64>(h, "CD1_2"),
869                cd1_3: parse_opt_card::<f64>(h, "CD1_3"),
870                cd2_1: parse_opt_card::<f64>(h, "CD2_1"),
871                cd2_2: parse_opt_card::<f64>(h, "CD2_2"),
872                cd2_3: parse_opt_card::<f64>(h, "CD2_3"),
873                cd3_1: parse_opt_card::<f64>(h, "CD3_1"),
874                cd3_2: parse_opt_card::<f64>(h, "CD3_2"),
875                cd3_3: parse_opt_card::<f64>(h, "CD3_3"),
876                pc1_1: parse_opt_card::<f64>(h, "PC1_1"),
877                pc1_2: parse_opt_card::<f64>(h, "PC1_2"),
878                pc1_3: parse_opt_card::<f64>(h, "PC1_3"),
879                pc2_1: parse_opt_card::<f64>(h, "PC2_1"),
880                pc2_2: parse_opt_card::<f64>(h, "PC2_2"),
881                pc2_3: parse_opt_card::<f64>(h, "PC2_3"),
882                pc3_1: parse_opt_card::<f64>(h, "PC3_1"),
883                pc3_2: parse_opt_card::<f64>(h, "PC3_2"),
884                pc3_3: parse_opt_card::<f64>(h, "PC3_3"),
885                a_0_0: parse_opt_card::<f64>(h, "A_0_0"),
886                a_1_0: parse_opt_card::<f64>(h, "A_1_0"),
887                a_2_0: parse_opt_card::<f64>(h, "A_2_0"),
888                a_3_0: parse_opt_card::<f64>(h, "A_3_0"),
889                a_4_0: parse_opt_card::<f64>(h, "A_4_0"),
890                a_5_0: parse_opt_card::<f64>(h, "A_5_0"),
891                a_6_0: parse_opt_card::<f64>(h, "A_6_0"),
892                a_0_1: parse_opt_card::<f64>(h, "A_0_1"),
893                a_1_1: parse_opt_card::<f64>(h, "A_1_1"),
894                a_2_1: parse_opt_card::<f64>(h, "A_2_1"),
895                a_3_1: parse_opt_card::<f64>(h, "A_3_1"),
896                a_4_1: parse_opt_card::<f64>(h, "A_4_1"),
897                a_5_1: parse_opt_card::<f64>(h, "A_5_1"),
898                a_0_2: parse_opt_card::<f64>(h, "A_0_2"),
899                a_1_2: parse_opt_card::<f64>(h, "A_1_2"),
900                a_2_2: parse_opt_card::<f64>(h, "A_2_2"),
901                a_3_2: parse_opt_card::<f64>(h, "A_3_2"),
902                a_4_2: parse_opt_card::<f64>(h, "A_4_2"),
903                a_0_3: parse_opt_card::<f64>(h, "A_0_3"),
904                a_1_3: parse_opt_card::<f64>(h, "A_1_3"),
905                a_2_3: parse_opt_card::<f64>(h, "A_2_3"),
906                a_3_3: parse_opt_card::<f64>(h, "A_3_3"),
907                a_0_4: parse_opt_card::<f64>(h, "A_0_4"),
908                a_1_4: parse_opt_card::<f64>(h, "A_1_4"),
909                a_2_4: parse_opt_card::<f64>(h, "A_2_4"),
910                a_0_5: parse_opt_card::<f64>(h, "A_0_5"),
911                a_1_5: parse_opt_card::<f64>(h, "A_1_5"),
912                a_0_6: parse_opt_card::<f64>(h, "A_0_6"),
913                ap_0_0: parse_opt_card::<f64>(h, "AP_0_0"),
914                ap_1_0: parse_opt_card::<f64>(h, "AP_1_0"),
915                ap_2_0: parse_opt_card::<f64>(h, "AP_2_0"),
916                ap_3_0: parse_opt_card::<f64>(h, "AP_3_0"),
917                ap_4_0: parse_opt_card::<f64>(h, "AP_4_0"),
918                ap_5_0: parse_opt_card::<f64>(h, "AP_5_0"),
919                ap_6_0: parse_opt_card::<f64>(h, "AP_6_0"),
920                ap_0_1: parse_opt_card::<f64>(h, "AP_0_1"),
921                ap_1_1: parse_opt_card::<f64>(h, "AP_1_1"),
922                ap_2_1: parse_opt_card::<f64>(h, "AP_2_1"),
923                ap_3_1: parse_opt_card::<f64>(h, "AP_3_1"),
924                ap_4_1: parse_opt_card::<f64>(h, "AP_4_1"),
925                ap_5_1: parse_opt_card::<f64>(h, "AP_5_1"),
926                ap_0_2: parse_opt_card::<f64>(h, "AP_0_2"),
927                ap_1_2: parse_opt_card::<f64>(h, "AP_1_2"),
928                ap_2_2: parse_opt_card::<f64>(h, "AP_2_2"),
929                ap_3_2: parse_opt_card::<f64>(h, "AP_3_2"),
930                ap_4_2: parse_opt_card::<f64>(h, "AP_4_2"),
931                ap_0_3: parse_opt_card::<f64>(h, "AP_0_3"),
932                ap_1_3: parse_opt_card::<f64>(h, "AP_1_3"),
933                ap_2_3: parse_opt_card::<f64>(h, "AP_2_3"),
934                ap_3_3: parse_opt_card::<f64>(h, "AP_3_3"),
935                ap_0_4: parse_opt_card::<f64>(h, "AP_0_4"),
936                ap_1_4: parse_opt_card::<f64>(h, "AP_1_4"),
937                ap_2_4: parse_opt_card::<f64>(h, "AP_2_4"),
938                ap_0_5: parse_opt_card::<f64>(h, "AP_0_5"),
939                ap_1_5: parse_opt_card::<f64>(h, "AP_1_5"),
940                ap_0_6: parse_opt_card::<f64>(h, "AP_0_6"),
941                b_0_0: parse_opt_card::<f64>(h, "B_0_0"),
942                b_1_0: parse_opt_card::<f64>(h, "B_1_0"),
943                b_2_0: parse_opt_card::<f64>(h, "B_2_0"),
944                b_3_0: parse_opt_card::<f64>(h, "B_3_0"),
945                b_4_0: parse_opt_card::<f64>(h, "B_4_0"),
946                b_5_0: parse_opt_card::<f64>(h, "B_5_0"),
947                b_6_0: parse_opt_card::<f64>(h, "B_6_0"),
948                b_0_1: parse_opt_card::<f64>(h, "B_0_1"),
949                b_1_1: parse_opt_card::<f64>(h, "B_1_1"),
950                b_2_1: parse_opt_card::<f64>(h, "B_2_1"),
951                b_3_1: parse_opt_card::<f64>(h, "B_3_1"),
952                b_4_1: parse_opt_card::<f64>(h, "B_4_1"),
953                b_5_1: parse_opt_card::<f64>(h, "B_5_1"),
954                b_0_2: parse_opt_card::<f64>(h, "B_0_2"),
955                b_1_2: parse_opt_card::<f64>(h, "B_1_2"),
956                b_2_2: parse_opt_card::<f64>(h, "B_2_2"),
957                b_3_2: parse_opt_card::<f64>(h, "B_3_2"),
958                b_4_2: parse_opt_card::<f64>(h, "B_4_2"),
959                b_0_3: parse_opt_card::<f64>(h, "B_0_3"),
960                b_1_3: parse_opt_card::<f64>(h, "B_1_3"),
961                b_2_3: parse_opt_card::<f64>(h, "B_2_3"),
962                b_3_3: parse_opt_card::<f64>(h, "B_3_3"),
963                b_0_4: parse_opt_card::<f64>(h, "B_0_4"),
964                b_1_4: parse_opt_card::<f64>(h, "B_1_4"),
965                b_2_4: parse_opt_card::<f64>(h, "B_2_4"),
966                b_0_5: parse_opt_card::<f64>(h, "B_0_5"),
967                b_1_5: parse_opt_card::<f64>(h, "B_1_5"),
968                b_0_6: parse_opt_card::<f64>(h, "B_0_6"),
969                bp_0_0: parse_opt_card::<f64>(h, "BP_0_0"),
970                bp_1_0: parse_opt_card::<f64>(h, "BP_1_0"),
971                bp_2_0: parse_opt_card::<f64>(h, "BP_2_0"),
972                bp_3_0: parse_opt_card::<f64>(h, "BP_3_0"),
973                bp_4_0: parse_opt_card::<f64>(h, "BP_4_0"),
974                bp_5_0: parse_opt_card::<f64>(h, "BP_5_0"),
975                bp_6_0: parse_opt_card::<f64>(h, "BP_6_0"),
976                bp_0_1: parse_opt_card::<f64>(h, "BP_0_1"),
977                bp_1_1: parse_opt_card::<f64>(h, "BP_1_1"),
978                bp_2_1: parse_opt_card::<f64>(h, "BP_2_1"),
979                bp_3_1: parse_opt_card::<f64>(h, "BP_3_1"),
980                bp_4_1: parse_opt_card::<f64>(h, "BP_4_1"),
981                bp_5_1: parse_opt_card::<f64>(h, "BP_5_1"),
982                bp_0_2: parse_opt_card::<f64>(h, "BP_0_2"),
983                bp_1_2: parse_opt_card::<f64>(h, "BP_1_2"),
984                bp_2_2: parse_opt_card::<f64>(h, "BP_2_2"),
985                bp_3_2: parse_opt_card::<f64>(h, "BP_3_2"),
986                bp_4_2: parse_opt_card::<f64>(h, "BP_4_2"),
987                bp_0_3: parse_opt_card::<f64>(h, "BP_0_3"),
988                bp_1_3: parse_opt_card::<f64>(h, "BP_1_3"),
989                bp_2_3: parse_opt_card::<f64>(h, "BP_2_3"),
990                bp_3_3: parse_opt_card::<f64>(h, "BP_3_3"),
991                bp_0_4: parse_opt_card::<f64>(h, "BP_0_4"),
992                bp_1_4: parse_opt_card::<f64>(h, "BP_1_4"),
993                bp_2_4: parse_opt_card::<f64>(h, "BP_2_4"),
994                bp_0_5: parse_opt_card::<f64>(h, "BP_0_5"),
995                bp_1_5: parse_opt_card::<f64>(h, "BP_1_5"),
996                bp_0_6: parse_opt_card::<f64>(h, "BP_0_6"),
997            };
998
999            WCS::new(&params)
1000        }
1001    }
1002
1003    fn wcs_from_fits_header(header: &Header<Image>) -> Result<WCS, Error> {
1004        header.try_into()
1005    }
1006
1007    #[test]
1008    fn test_visualize() {
1009        let f = File::open("examples/panstarrs-rotated-around-orion.fits").unwrap();
1010
1011        let reader = BufReader::new(f);
1012        let mut fits = Fits::from_reader(reader);
1013        let hdu = fits.next().unwrap().unwrap();
1014
1015        match hdu {
1016            HDU::XImage(hdu) | HDU::Primary(hdu) => {
1017                let header = hdu.get_header();
1018
1019                // Parse data
1020                let data = match fits.get_data(&hdu).pixels() {
1021                    Pixels::F32(it) => it.collect::<Vec<_>>(),
1022                    _ => unreachable!(),
1023                };
1024
1025                let wcs = wcs_from_fits_header(&header).unwrap();
1026                reproject_fits_image(mapproj::zenithal::azp::Azp::new(), &wcs, &header, &data);
1027                reproject_fits_image(mapproj::zenithal::szp::Szp::new(), &wcs, &header, &data);
1028                reproject_fits_image(mapproj::zenithal::tan::Tan::new(), &wcs, &header, &data);
1029                reproject_fits_image(mapproj::zenithal::stg::Stg::new(), &wcs, &header, &data);
1030                reproject_fits_image(mapproj::zenithal::sin::Sin::new(), &wcs, &header, &data);
1031                reproject_fits_image(mapproj::zenithal::arc::Arc::new(), &wcs, &header, &data);
1032                reproject_fits_image(mapproj::zenithal::zea::Zea::new(), &wcs, &header, &data);
1033                reproject_fits_image(mapproj::zenithal::air::Air::new(), &wcs, &header, &data);
1034                reproject_fits_image(mapproj::zenithal::ncp::Ncp::new(), &wcs, &header, &data);
1035
1036                reproject_fits_image(mapproj::pseudocyl::mol::Mol::new(), &wcs, &header, &data);
1037                reproject_fits_image(mapproj::pseudocyl::ait::Ait::new(), &wcs, &header, &data);
1038                reproject_fits_image(mapproj::pseudocyl::par::Par::new(), &wcs, &header, &data);
1039                reproject_fits_image(mapproj::pseudocyl::sfl::Sfl::new(), &wcs, &header, &data);
1040
1041                reproject_fits_image(mapproj::cylindrical::cyp::Cyp::new(), &wcs, &header, &data);
1042                reproject_fits_image(mapproj::cylindrical::cea::Cea::new(), &wcs, &header, &data);
1043                reproject_fits_image(mapproj::cylindrical::car::Car::new(), &wcs, &header, &data);
1044                reproject_fits_image(mapproj::cylindrical::mer::Mer::new(), &wcs, &header, &data);
1045
1046                reproject_fits_image(mapproj::conic::cod::Cod::new(), &wcs, &header, &data);
1047                reproject_fits_image(mapproj::conic::cop::Cop::new(), &wcs, &header, &data);
1048                reproject_fits_image(mapproj::conic::coo::Coo::new(), &wcs, &header, &data);
1049                reproject_fits_image(mapproj::conic::coe::Coe::new(), &wcs, &header, &data);
1050
1051                reproject_fits_image(mapproj::hybrid::hpx::Hpx::new(), &wcs, &header, &data);
1052            }
1053            _ => unreachable!(),
1054        }
1055    }
1056
1057    fn reproject_fits_image<'a, T: CanonicalProjection>(
1058        proj: T,
1059        wcs: &WCS,
1060        header: &Header<Image>,
1061        data: &[f32],
1062    ) {
1063        let scale = header.get_parsed::<f32>("BSCALE").unwrap_or(1.0);
1064        let offset = header.get_parsed::<f32>("BZERO").unwrap_or(0.0);
1065
1066        let xtension = header.get_xtension();
1067        let naxis = xtension.get_naxis();
1068        let width = naxis[0];
1069        let height = naxis[1];
1070
1071        //let proj = mapproj::zenithal::sin::Sin::new();
1072        let bounds = proj.bounds();
1073        let x_bounds = bounds.x_bounds().as_ref().unwrap_or(&((-PI)..=PI));
1074        let y_bounds = bounds.y_bounds().as_ref().unwrap_or(&((-PI)..=PI));
1075
1076        let x_off = x_bounds.start();
1077        let x_len = x_bounds.end() - x_bounds.start();
1078
1079        let y_off = y_bounds.start();
1080        let y_len = y_bounds.end() - y_bounds.start();
1081
1082        const WIDTH_IMAGE: usize = 1024;
1083        const HEIGHT_IMAGE: usize = 1024;
1084        // Create a new ImgBuf with width: imgx and height: imgy
1085        let mut imgbuf = image::ImageBuffer::new(WIDTH_IMAGE as u32, HEIGHT_IMAGE as u32);
1086
1087        for y in 0..height {
1088            for x in 0..width {
1089                let grayscale_val = (data[(y * width + x) as usize] * scale + offset) as u8;
1090
1091                let img_xy = ImgXY::new(x as f64, y as f64);
1092                if let Some(lonlat) = wcs.unproj(&img_xy) {
1093                    if let Some(proj_xy) = proj.proj_lonlat(&lonlat) {
1094                        let proj_x = ((proj_xy.x() as f64) - x_off) / x_len; // between 0 and 1
1095                        let proj_y = ((proj_xy.y() as f64) - y_off) / y_len; // between 0 and 1
1096
1097                        if (0.0..1.0).contains(&proj_x) && (0.0..1.0).contains(&proj_y) {
1098                            let ix = (proj_x * (WIDTH_IMAGE as f64)) as usize;
1099                            let iy = (proj_y * (HEIGHT_IMAGE as f64)) as usize;
1100
1101                            let pixel = imgbuf
1102                                .get_pixel_mut(ix as u32, (HEIGHT_IMAGE as u32) - iy as u32 - 1);
1103                            *pixel = image::Rgb([grayscale_val, grayscale_val, grayscale_val]);
1104                        }
1105                    }
1106                }
1107            }
1108        }
1109
1110        let filename = &format!(
1111            "tests/reproj/pans-{}.jpeg",
1112            <T as CanonicalProjection>::WCS_NAME
1113        );
1114        imgbuf.save(filename).unwrap();
1115    }
1116
1117    macro_rules! assert_delta {
1118        ($x:expr, $y:expr, $d:expr) => {
1119            if ($x - $y).abs() > $d {
1120                panic!();
1121            }
1122        };
1123    }
1124
1125    #[test]
1126    fn astropy_comparison() {
1127        use std::fs;
1128
1129        let fits_file_paths = fs::read_dir("./examples")
1130            .unwrap()
1131            .filter_map(|res| res.ok())
1132            // Map the directory entries to paths
1133            .map(|dir_entry| dir_entry.path())
1134            // Filter out all paths with extensions other than `csv`
1135            .filter_map(|path| {
1136                if path.extension().map_or(false, |ext| ext == "fits") {
1137                    Some(path)
1138                } else {
1139                    None
1140                }
1141            });
1142
1143        for path in fits_file_paths {
1144            println!("Test {:?}", path.display());
1145
1146            let f = File::open(path.clone()).unwrap();
1147
1148            let reader = BufReader::new(f);
1149            let mut fits = Fits::from_reader(reader);
1150            let hdu = fits.next().unwrap().unwrap();
1151            match hdu {
1152                HDU::XImage(hdu) | HDU::Primary(hdu) => {
1153                    let header = hdu.get_header();
1154                    let wcs = wcs_from_fits_header(header).unwrap();
1155
1156                    // add the astropy opened results
1157                    let path_astropy = path.with_extension("fits.csv");
1158                    let f = File::open(path_astropy).unwrap();
1159
1160                    let mut rdr = csv::Reader::from_reader(BufReader::new(f));
1161                    for result in rdr.records() {
1162                        let record = result.unwrap();
1163
1164                        let ra: f64 = record[0].parse().unwrap();
1165                        let dec: f64 = record[1].parse().unwrap();
1166                        let x: f64 = record[2].parse().unwrap();
1167                        let y: f64 = record[3].parse().unwrap();
1168
1169                        if ra.is_finite() && dec.is_finite() {
1170                            if let Some(img_xy) = wcs.proj(&LonLat::new(ra, dec)) {
1171                                //dbg!(img_xy.x() - x);
1172                                //dbg!(img_xy.y() - y);
1173
1174                                assert_delta!(img_xy.x(), x, 1e-4);
1175                                assert_delta!(img_xy.y(), y, 1e-4);
1176                            }
1177                        }
1178                    }
1179                }
1180                _ => unreachable!(),
1181            };
1182        }
1183    }
1184
1185    #[test]
1186    fn crval_to_crpix() {
1187        for entry in glob("examples/*.fits").unwrap() {
1188            if let Ok(path) = dbg!(entry) {
1189                let f = File::open(path).unwrap();
1190                let reader = BufReader::new(f);
1191                let mut fits = Fits::from_reader(reader);
1192                let hdu = fits.next().unwrap().unwrap();
1193
1194                match hdu {
1195                    HDU::XImage(hdu) | HDU::Primary(hdu) => {
1196                        let header = hdu.get_header();
1197                        let crval1 = header.get_parsed::<f64>("CRVAL1").unwrap_or(0.0);
1198                        let crval2 = header.get_parsed::<f64>("CRVAL2").unwrap_or(0.0);
1199                        let crpix1 =
1200                            if let Some(Value::Integer { value, .. }) = header.get("CRPIX1") {
1201                                *value as f64
1202                            } else if let Some(Value::Float { value, .. }) = header.get("CRPIX1") {
1203                                *value
1204                            } else {
1205                                0.0
1206                            };
1207
1208                        let crpix2 =
1209                            if let Some(Value::Integer { value, .. }) = header.get("CRPIX2") {
1210                                *value as f64
1211                            } else if let Some(Value::Float { value, .. }) = header.get("CRPIX2") {
1212                                *value
1213                            } else {
1214                                0.0
1215                            };
1216
1217                        let wcs = wcs_from_fits_header(&header).unwrap();
1218
1219                        // crval to crpix
1220                        let proj_px = wcs
1221                            .proj(&LonLat::new(
1222                                dbg!(crval1).to_radians(),
1223                                dbg!(crval2).to_radians(),
1224                            ))
1225                            .unwrap();
1226                        assert_delta!(proj_px.x(), crpix1, 1e-6);
1227                        assert_delta!(proj_px.y(), crpix2, 1e-6);
1228
1229                        // crpix to crval
1230                        let lonlat = wcs.unproj_lonlat(&ImgXY::new(crpix1, crpix2)).unwrap();
1231                        assert_delta!(lonlat.lon(), crval1.to_radians(), 1e-6);
1232                        assert_delta!(lonlat.lat(), crval2.to_radians(), 1e-6);
1233                    }
1234                    _ => unreachable!(),
1235                };
1236            }
1237        }
1238    }
1239
1240    use fitsrs::hdu::HDU;
1241    #[test]
1242    fn open_fits() {
1243        let f = File::open("examples/neowise.fits").unwrap();
1244
1245        let reader = BufReader::new(f);
1246        let mut fits = Fits::from_reader(reader);
1247        let hdu = fits.next().unwrap().unwrap();
1248        match hdu {
1249            HDU::XImage(hdu) | HDU::Primary(hdu) => {
1250                let header = hdu.get_header();
1251                let wcs = wcs_from_fits_header(header).unwrap();
1252                dbg!(wcs.unproj(&ImgXY::new(0.0, 1200.0)));
1253            }
1254            _ => unreachable!(),
1255        }
1256    }
1257}
wcs/lib.rs

wcs/
lib.rs
