#![allow(dead_code)]
#![allow(missing_docs)]
use crate::error::{IoError, Result};
use byteorder::{BigEndian, ReadBytesExt, WriteBytesExt};
use scirs2_core::ndarray::Array2;
use std::collections::HashMap;
use std::fs::File;
use std::io::{BufWriter, Read, Seek, SeekFrom, Write};
use std::path::Path;
pub struct FitsFile {
file_path: String,
hdus: Vec<HDU>,
}
#[derive(Debug, Clone)]
pub struct HDU {
pub hdu_type: HDUType,
pub header: FitsHeader,
pub data_offset: u64,
pub data_size: usize,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum HDUType {
Primary,
Image,
AsciiTable,
BinaryTable,
}
#[derive(Debug, Clone)]
pub struct FitsHeader {
cards: Vec<HeaderCard>,
card_map: HashMap<String, usize>,
}
#[derive(Debug, Clone)]
pub struct HeaderCard {
pub keyword: String,
pub value: CardValue,
pub comment: Option<String>,
}
#[derive(Debug, Clone, PartialEq)]
pub enum CardValue {
Boolean(bool),
Integer(i64),
Float(f64),
String(String),
Complex(f64, f64),
None,
}
impl FitsHeader {
pub fn new() -> Self {
Self {
cards: Vec::new(),
card_map: HashMap::new(),
}
}
pub fn add_card(&mut self, card: HeaderCard) {
let index = self.cards.len();
self.card_map.insert(card.keyword.clone(), index);
self.cards.push(card);
}
pub fn get_card(&self, keyword: &str) -> Option<&HeaderCard> {
self.card_map.get(keyword).map(|&idx| &self.cards[idx])
}
pub fn get_bool(&self, keyword: &str) -> Result<bool> {
match self.get_card(keyword) {
Some(card) => match &card.value {
CardValue::Boolean(b) => Ok(*b),
_ => Err(IoError::ParseError(format!(
"Keyword {keyword} is not a boolean"
))),
},
None => Err(IoError::ParseError(format!("Keyword {keyword} not found"))),
}
}
pub fn get_i64(&self, keyword: &str) -> Result<i64> {
match self.get_card(keyword) {
Some(card) => match &card.value {
CardValue::Integer(i) => Ok(*i),
_ => Err(IoError::ParseError(format!(
"Keyword {keyword} is not an integer"
))),
},
None => Err(IoError::ParseError(format!("Keyword {keyword} not found"))),
}
}
pub fn get_f64(&self, keyword: &str) -> Result<f64> {
match self.get_card(keyword) {
Some(card) => match &card.value {
CardValue::Float(f) => Ok(*f),
CardValue::Integer(i) => Ok(*i as f64),
_ => Err(IoError::ParseError(format!(
"Keyword {keyword} is not a number"
))),
},
None => Err(IoError::ParseError(format!("Keyword {keyword} not found"))),
}
}
pub fn get_string(&self, keyword: &str) -> Result<String> {
match self.get_card(keyword) {
Some(card) => match &card.value {
CardValue::String(s) => Ok(s.clone()),
_ => Err(IoError::ParseError(format!(
"Keyword {keyword} is not a string"
))),
},
None => Err(IoError::ParseError(format!("Keyword {keyword} not found"))),
}
}
pub fn cards(&self) -> &[HeaderCard] {
&self.cards
}
pub fn has_keyword(&self, keyword: &str) -> bool {
self.card_map.contains_key(keyword)
}
}
impl Default for FitsHeader {
fn default() -> Self {
Self::new()
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FitsDataType {
UInt8,
Int16,
Int32,
Int64,
Float32,
Float64,
}
impl FitsDataType {
pub fn byte_size(&self) -> usize {
match self {
FitsDataType::UInt8 => 1,
FitsDataType::Int16 => 2,
FitsDataType::Int32 => 4,
FitsDataType::Int64 => 8,
FitsDataType::Float32 => 4,
FitsDataType::Float64 => 8,
}
}
pub fn bitpix(&self) -> i32 {
match self {
FitsDataType::UInt8 => 8,
FitsDataType::Int16 => 16,
FitsDataType::Int32 => 32,
FitsDataType::Int64 => 64,
FitsDataType::Float32 => -32,
FitsDataType::Float64 => -64,
}
}
pub fn frombitpix(bitpix: i32) -> Result<Self> {
match bitpix {
8 => Ok(FitsDataType::UInt8),
16 => Ok(FitsDataType::Int16),
32 => Ok(FitsDataType::Int32),
64 => Ok(FitsDataType::Int64),
-32 => Ok(FitsDataType::Float32),
-64 => Ok(FitsDataType::Float64),
_ => Err(IoError::ParseError(format!(
"Invalid BITPIX value: {bitpix}"
))),
}
}
}
impl FitsFile {
pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
let file_path = path.as_ref().to_string_lossy().to_string();
let mut file =
File::open(path.as_ref()).map_err(|_e| IoError::FileNotFound(file_path.clone()))?;
let mut hdus = Vec::new();
let mut offset = 0u64;
loop {
file.seek(SeekFrom::Start(offset))
.map_err(|e| IoError::ParseError(format!("Failed to seek: {e}")))?;
let header = Self::readheader(&mut file)?;
let hdu_type = if hdus.is_empty() {
HDUType::Primary
} else {
match header.get_string("XTENSION").ok().as_deref() {
Some("IMAGE") => HDUType::Image,
Some("TABLE") => HDUType::AsciiTable,
Some("BINTABLE") => HDUType::BinaryTable,
_ => HDUType::Image,
}
};
let data_size = Self::calculate_data_size(&header)?;
let header_blocks = ((header.cards.len() + 35) / 36) as u64; let data_offset = offset + header_blocks * 2880;
hdus.push(HDU {
hdu_type,
header,
data_offset,
data_size,
});
let data_blocks = ((data_size + 2879) / 2880) as u64;
offset = data_offset + data_blocks * 2880;
if hdus
.last()
.expect("Operation failed")
.header
.cards
.iter()
.any(|c| c.keyword == "END")
{
if file.seek(SeekFrom::Start(offset)).is_err() {
break;
}
let mut test_buf = [0u8; 80];
match file.read_exact(&mut test_buf) {
Ok(_) => {
let test_str = String::from_utf8_lossy(&test_buf);
if test_str.trim().is_empty() {
break;
}
}
Err(_) => break,
}
}
}
Ok(Self { file_path, hdus })
}
fn readheader<R: Read>(reader: &mut R) -> Result<FitsHeader> {
let mut header = FitsHeader::new();
let mut card_buf = [0u8; 80];
loop {
reader
.read_exact(&mut card_buf)
.map_err(|e| IoError::ParseError(format!("Failed to read header card: {e}")))?;
let cardstr = String::from_utf8_lossy(&card_buf);
if let Some(card) = Self::parseheader_card(&cardstr) {
if card.keyword == "END" {
break;
}
header.add_card(card);
}
}
Ok(header)
}
fn parseheader_card(cardstr: &str) -> Option<HeaderCard> {
if cardstr.len() < 8 {
return None;
}
let keyword = cardstr[0..8].trim().to_string();
if keyword.is_empty() || keyword == "COMMENT" || keyword == "HISTORY" {
return Some(HeaderCard {
keyword,
value: CardValue::None,
comment: Some(cardstr[8..].trim().to_string()),
});
}
if keyword == "END" {
return Some(HeaderCard {
keyword,
value: CardValue::None,
comment: None,
});
}
if cardstr.len() > 9 && &cardstr[8..9] == "=" {
let value_comment = &cardstr[10..];
if let Some(slash_pos) = value_comment.find('/') {
let valuestr = value_comment[..slash_pos].trim();
let comment = value_comment[slash_pos + 1..].trim().to_string();
let value = Self::parse_card_value(valuestr);
Some(HeaderCard {
keyword,
value,
comment: Some(comment),
})
} else {
let value = Self::parse_card_value(value_comment.trim());
Some(HeaderCard {
keyword,
value,
comment: None,
})
}
} else {
None
}
}
fn parse_card_value(valuestr: &str) -> CardValue {
if valuestr == "T" {
return CardValue::Boolean(true);
}
if valuestr == "F" {
return CardValue::Boolean(false);
}
if valuestr.starts_with('\'') && valuestr.ends_with('\'') {
let s = valuestr[1..valuestr.len() - 1].trim().to_string();
return CardValue::String(s);
}
if let Ok(i) = valuestr.parse::<i64>() {
return CardValue::Integer(i);
}
if let Ok(f) = valuestr.parse::<f64>() {
return CardValue::Float(f);
}
CardValue::String(valuestr.to_string())
}
fn calculate_data_size(header: &FitsHeader) -> Result<usize> {
if let Ok(bitpix) = header.get_i64("BITPIX") {
let mut size = (bitpix.abs() / 8) as usize;
let naxis = header.get_i64("NAXIS").unwrap_or(0) as usize;
for i in 1..=naxis {
let axis_key = format!("NAXIS{i}");
if let Ok(axis_size) = header.get_i64(&axis_key) {
size *= axis_size as usize;
}
}
return Ok(size);
}
if let Ok(naxis2) = header.get_i64("NAXIS2") {
if let Ok(naxis1) = header.get_i64("NAXIS1") {
return Ok((naxis1 * naxis2) as usize);
}
}
Ok(0)
}
pub fn primaryheader(&self) -> &FitsHeader {
&self.hdus[0].header
}
pub fn hdu_count(&self) -> usize {
self.hdus.len()
}
pub fn gethdu(&self, index: usize) -> Result<&HDU> {
self.hdus
.get(index)
.ok_or_else(|| IoError::ParseError(format!("HDU index {index} out of range")))
}
pub fn read_image<T: FitsNumeric>(&self) -> Result<Array2<T>> {
self.readhdu_image(0)
}
pub fn readhdu_image<T: FitsNumeric>(&self, hduindex: usize) -> Result<Array2<T>> {
let hdu = self.gethdu(hduindex)?;
if hdu.hdu_type != HDUType::Primary && hdu.hdu_type != HDUType::Image {
return Err(IoError::ParseError(format!(
"HDU {hduindex} is not an image"
)));
}
let bitpix = hdu.header.get_i64("BITPIX")?;
let naxis = hdu.header.get_i64("NAXIS")?;
if naxis != 2 {
return Err(IoError::ParseError(format!(
"Expected 2D image, got {naxis}D"
)));
}
let naxis1 = hdu.header.get_i64("NAXIS1")? as usize;
let naxis2 = hdu.header.get_i64("NAXIS2")? as usize;
let mut file = File::open(&self.file_path)
.map_err(|_e| IoError::FileNotFound(self.file_path.clone()))?;
file.seek(SeekFrom::Start(hdu.data_offset))
.map_err(|e| IoError::ParseError(format!("Failed to seek to data: {e}")))?;
let datatype = FitsDataType::frombitpix(bitpix as i32)?;
let mut values = Vec::with_capacity(naxis1 * naxis2);
for _ in 0..(naxis1 * naxis2) {
let value = T::read_fits(&mut file, datatype)?;
values.push(value);
}
let array = Array2::from_shape_vec((naxis2, naxis1), values)
.map_err(|e| IoError::ParseError(format!("Failed to create array: {e}")))?;
Ok(array.t().to_owned())
}
pub fn image_dimensions(&self, hduindex: usize) -> Result<Vec<usize>> {
let hdu = self.gethdu(hduindex)?;
let naxis = hdu.header.get_i64("NAXIS")? as usize;
let mut dims = Vec::with_capacity(naxis);
for i in 1..=naxis {
let axis_key = format!("NAXIS{i}");
let size = hdu.header.get_i64(&axis_key)? as usize;
dims.push(size);
}
Ok(dims)
}
}
pub trait FitsNumeric: Default + Clone {
fn read_fits<R: Read>(reader: &mut R, datatype: FitsDataType) -> Result<Self>;
fn write_fits<W: Write>(&self, writer: &mut W, datatype: FitsDataType) -> Result<()>;
}
impl FitsNumeric for f32 {
fn read_fits<R: Read>(reader: &mut R, datatype: FitsDataType) -> Result<Self> {
match datatype {
FitsDataType::Float32 => reader
.read_f32::<BigEndian>()
.map_err(|e| IoError::ParseError(format!("Failed to read f32: {e}"))),
FitsDataType::Float64 => reader
.read_f64::<BigEndian>()
.map(|v| v as f32)
.map_err(|e| IoError::ParseError(format!("Failed to read f64: {e}"))),
FitsDataType::Int16 => reader
.read_i16::<BigEndian>()
.map(|v| v as f32)
.map_err(|e| IoError::ParseError(format!("Failed to read i16: {e}"))),
FitsDataType::Int32 => reader
.read_i32::<BigEndian>()
.map(|v| v as f32)
.map_err(|e| IoError::ParseError(format!("Failed to read i32: {e}"))),
_ => Err(IoError::ParseError(format!(
"Unsupported conversion from {datatype:?} to f32"
))),
}
}
fn write_fits<W: Write>(&self, writer: &mut W, datatype: FitsDataType) -> Result<()> {
match datatype {
FitsDataType::Float32 => writer
.write_f32::<BigEndian>(*self)
.map_err(|e| IoError::FileError(format!("Failed to write f32: {e}"))),
_ => Err(IoError::FileError(format!(
"Unsupported conversion from f32 to {datatype:?}"
))),
}
}
}
impl FitsNumeric for f64 {
fn read_fits<R: Read>(reader: &mut R, datatype: FitsDataType) -> Result<Self> {
match datatype {
FitsDataType::Float64 => reader
.read_f64::<BigEndian>()
.map_err(|e| IoError::ParseError(format!("Failed to read f64: {e}"))),
FitsDataType::Float32 => reader
.read_f32::<BigEndian>()
.map(|v| v as f64)
.map_err(|e| IoError::ParseError(format!("Failed to read f32: {e}"))),
FitsDataType::Int32 => reader
.read_i32::<BigEndian>()
.map(|v| v as f64)
.map_err(|e| IoError::ParseError(format!("Failed to read i32: {e}"))),
_ => Err(IoError::ParseError(format!(
"Unsupported conversion from {datatype:?} to f64"
))),
}
}
fn write_fits<W: Write>(&self, writer: &mut W, datatype: FitsDataType) -> Result<()> {
match datatype {
FitsDataType::Float64 => writer
.write_f64::<BigEndian>(*self)
.map_err(|e| IoError::FileError(format!("Failed to write f64: {e}"))),
_ => Err(IoError::FileError(format!(
"Unsupported conversion from f64 to {datatype:?}"
))),
}
}
}
pub struct FitsWriter {
writer: BufWriter<File>,
currenthdu: usize,
}
impl FitsWriter {
pub fn create<P: AsRef<Path>>(path: P) -> Result<Self> {
let file = File::create(path.as_ref())
.map_err(|e| IoError::FileError(format!("Failed to create file: {e}")))?;
Ok(Self {
writer: BufWriter::new(file),
currenthdu: 0,
})
}
pub fn write_image_2d<T: FitsNumeric>(
&mut self,
data: &Array2<T>,
datatype: FitsDataType,
) -> Result<()> {
let mut header = FitsHeader::new();
header.add_card(HeaderCard {
keyword: "SIMPLE".to_string(),
value: CardValue::Boolean(true),
comment: Some("Standard FITS format".to_string()),
});
header.add_card(HeaderCard {
keyword: "BITPIX".to_string(),
value: CardValue::Integer(datatype.bitpix() as i64),
comment: Some("Number of bits per pixel".to_string()),
});
header.add_card(HeaderCard {
keyword: "NAXIS".to_string(),
value: CardValue::Integer(2),
comment: Some("Number of axes".to_string()),
});
let (rows, cols) = data.dim();
header.add_card(HeaderCard {
keyword: "NAXIS1".to_string(),
value: CardValue::Integer(cols as i64),
comment: Some("Length of axis 1".to_string()),
});
header.add_card(HeaderCard {
keyword: "NAXIS2".to_string(),
value: CardValue::Integer(rows as i64),
comment: Some("Length of axis 2".to_string()),
});
self.writeheader(&header)?;
for col in 0..cols {
for row in 0..rows {
data[[row, col]].write_fits(&mut self.writer, datatype)?;
}
}
let data_bytes = rows * cols * datatype.byte_size();
let padding = (2880 - (data_bytes % 2880)) % 2880;
if padding > 0 {
let pad_bytes = vec![0u8; padding];
self.writer
.write_all(&pad_bytes)
.map_err(|e| IoError::FileError(format!("Failed to write padding: {e}")))?;
}
self.currenthdu += 1;
Ok(())
}
fn writeheader(&mut self, header: &FitsHeader) -> Result<()> {
for card in &header.cards {
self.writeheader_card(card)?;
}
self.writeheader_card(&HeaderCard {
keyword: "END".to_string(),
value: CardValue::None,
comment: None,
})?;
let cards_written = header.cards.len() + 1; let cards_per_block = 36;
let remaining = cards_per_block - (cards_written % cards_per_block);
if remaining < cards_per_block {
let blank_card = vec![b' '; 80];
for _ in 0..remaining {
self.writer
.write_all(&blank_card)
.map_err(|e| IoError::FileError(format!("Failed to write blank card: {e}")))?;
}
}
Ok(())
}
fn writeheader_card(&mut self, card: &HeaderCard) -> Result<()> {
let mut cardstr = String::with_capacity(80);
cardstr.push_str(&format!("{:<8}", card.keyword));
match &card.value {
CardValue::Boolean(b) => {
cardstr.push_str(&format!("= {:>20}", if *b { "T" } else { "F" }));
}
CardValue::Integer(i) => {
cardstr.push_str(&format!("= {i:>20}"));
}
CardValue::Float(f) => {
cardstr.push_str(&format!("= {f:>20.10E}"));
}
CardValue::String(s) => {
cardstr.push_str(&format!("= '{s:<18}'"));
}
CardValue::None => {
}
_ => {}
}
if let Some(comment) = &card.comment {
if cardstr.len() < 31 {
cardstr.push_str(&" ".repeat(31 - cardstr.len()));
}
cardstr.push_str(" / ");
cardstr.push_str(comment);
}
match cardstr.len().cmp(&80) {
std::cmp::Ordering::Less => cardstr.push_str(&" ".repeat(80 - cardstr.len())),
std::cmp::Ordering::Greater => cardstr.truncate(80),
std::cmp::Ordering::Equal => {}
}
self.writer
.write_all(cardstr.as_bytes())
.map_err(|e| IoError::FileError(format!("Failed to write header card: {e}")))
}
pub fn flush(&mut self) -> Result<()> {
self.writer
.flush()
.map_err(|e| IoError::FileError(format!("Failed to flush: {e}")))
}
pub fn close(mut self) -> Result<()> {
self.flush()
}
}
pub struct VOTable {
pub metadata: HashMap<String, String>,
pub columns: Vec<VOTableColumn>,
pub data: Vec<Vec<VOTableValue>>,
}
#[derive(Debug, Clone)]
pub struct VOTableColumn {
pub name: String,
pub datatype: String,
pub arraysize: Option<String>,
pub unit: Option<String>,
pub description: Option<String>,
pub ucd: Option<String>,
}
#[derive(Debug, Clone, PartialEq)]
pub enum VOTableValue {
Boolean(bool),
Integer(i64),
Float(f64),
Double(f64),
String(String),
Array(Vec<VOTableValue>),
Null,
}
impl VOTable {
pub fn new() -> Self {
Self {
metadata: HashMap::new(),
columns: Vec::new(),
data: Vec::new(),
}
}
pub fn add_column(&mut self, column: VOTableColumn) {
self.columns.push(column);
}
pub fn add_row(&mut self, row: Vec<VOTableValue>) -> Result<()> {
if row.len() != self.columns.len() {
return Err(IoError::FileError(format!(
"Row has {} values but table has {} columns",
row.len(),
self.columns.len()
)));
}
self.data.push(row);
Ok(())
}
pub fn get_column(&self, name: &str) -> Option<usize> {
self.columns.iter().position(|c| c.name == name)
}
pub fn get_column_data(&self, columnindex: usize) -> Result<Vec<&VOTableValue>> {
if columnindex >= self.columns.len() {
return Err(IoError::ParseError(format!(
"Column _index {columnindex} out of range"
)));
}
Ok(self.data.iter().map(|row| &row[columnindex]).collect())
}
pub fn read<P: AsRef<Path>>(path: P) -> Result<Self> {
Ok(Self::new())
}
pub fn write<P: AsRef<Path>>(&self, path: P) -> Result<()> {
let file = File::create(path.as_ref())
.map_err(|e| IoError::FileError(format!("Failed to create file: {e}")))?;
let mut writer = BufWriter::new(file);
writeln!(writer, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>")
.map_err(|e| IoError::FileError(format!("Failed to write XML header: {e}")))?;
writeln!(
writer,
"<VOTABLE version=\"1.4\" xmlns=\"http://www.ivoa.net/xml/VOTable/v1.3\">"
)
.map_err(|e| IoError::FileError(format!("Failed to write VOTABLE tag: {e}")))?;
writeln!(writer, " <RESOURCE>")
.map_err(|e| IoError::FileError(format!("Failed to write RESOURCE tag: {e}")))?;
writeln!(writer, " <TABLE>")
.map_err(|e| IoError::FileError(format!("Failed to write TABLE tag: {e}")))?;
for column in &self.columns {
write!(
writer,
" <FIELD name=\"{}\" datatype=\"{}\"",
column.name, column.datatype
)
.map_err(|e| IoError::FileError(format!("Failed to write FIELD: {e}")))?;
if let Some(unit) = &column.unit {
write!(writer, " unit=\"{unit}\"")
.map_err(|e| IoError::FileError(format!("Failed to write unit: {e}")))?;
}
writeln!(writer, "/>")
.map_err(|e| IoError::FileError(format!("Failed to write FIELD close: {e}")))?;
}
writeln!(writer, " <DATA>")
.map_err(|e| IoError::FileError(format!("Failed to write DATA tag: {e}")))?;
writeln!(writer, " <TABLEDATA>")
.map_err(|e| IoError::FileError(format!("Failed to write TABLEDATA tag: {e}")))?;
for row in &self.data {
write!(writer, " <TR>")
.map_err(|e| IoError::FileError(format!("Failed to write TR: {e}")))?;
for value in row {
match value {
VOTableValue::String(s) => write!(writer, "<TD>{s}</TD>"),
VOTableValue::Integer(i) => write!(writer, "<TD>{i}</TD>"),
VOTableValue::Float(f) | VOTableValue::Double(f) => {
write!(writer, "<TD>{f}</TD>")
}
VOTableValue::Boolean(b) => {
write!(writer, "<TD>{}</TD>", if *b { "true" } else { "false" })
}
VOTableValue::Null => write!(writer, "<TD/>"),
VOTableValue::Array(_) => write!(writer, "<TD>[]</TD>"), }
.map_err(|e| IoError::FileError(format!("Failed to write TD: {e}")))?;
}
writeln!(writer, "</TR>")
.map_err(|e| IoError::FileError(format!("Failed to write TR close: {e}")))?;
}
writeln!(writer, " </TABLEDATA>")
.map_err(|e| IoError::FileError(format!("Failed to close TABLEDATA: {e}")))?;
writeln!(writer, " </DATA>")
.map_err(|e| IoError::FileError(format!("Failed to close DATA: {e}")))?;
writeln!(writer, " </TABLE>")
.map_err(|e| IoError::FileError(format!("Failed to close TABLE: {e}")))?;
writeln!(writer, " </RESOURCE>")
.map_err(|e| IoError::FileError(format!("Failed to close RESOURCE: {e}")))?;
writeln!(writer, "</VOTABLE>")
.map_err(|e| IoError::FileError(format!("Failed to close VOTABLE: {e}")))?;
writer
.flush()
.map_err(|e| IoError::FileError(format!("Failed to flush: {e}")))
}
}
impl Default for VOTable {
fn default() -> Self {
Self::new()
}
}
#[derive(Debug, Clone)]
pub struct GeoTransform {
pub ref_lon: f64,
pub ref_lat: f64,
pub lon_scale: f64,
pub latscale: f64,
}
impl GeoTransform {
pub fn new(ref_lon: f64, ref_lat: f64, lon_scale: f64, latscale: f64) -> Self {
Self {
ref_lon,
ref_lat,
lon_scale,
latscale,
}
}
pub fn pixel_to_geo(&self, px: f64, py: f64) -> (f64, f64) {
let lon = self.ref_lon + px * self.lon_scale;
let lat = self.ref_lat + py * self.latscale;
(lon, lat)
}
pub fn geo_to_pixel(&self, lon: f64, lat: f64) -> (f64, f64) {
let px = (lon - self.ref_lon) / self.lon_scale;
let py = (lat - self.ref_lat) / self.latscale;
(px, py)
}
pub fn apply_wcs(&self, wcs: &WCSTransform) -> GeoTransform {
Self {
ref_lon: wcs.crval1,
ref_lat: wcs.crval2,
lon_scale: wcs.cdelt1,
latscale: wcs.cdelt2,
}
}
}
#[derive(Debug, Clone)]
pub struct WCSTransform {
pub crval1: f64,
pub crval2: f64,
pub crpix1: f64,
pub crpix2: f64,
pub cdelt1: f64,
pub cdelt2: f64,
pub cd_matrix: Option<[[f64; 2]; 2]>,
pub ctype1: String,
pub ctype2: String,
}
impl WCSTransform {
pub fn from_fitsheader(header: &FitsHeader) -> Result<Self> {
Ok(Self {
crval1: header.get_f64("CRVAL1").unwrap_or(0.0),
crval2: header.get_f64("CRVAL2").unwrap_or(0.0),
crpix1: header.get_f64("CRPIX1").unwrap_or(1.0),
crpix2: header.get_f64("CRPIX2").unwrap_or(1.0),
cdelt1: header.get_f64("CDELT1").unwrap_or(1.0),
cdelt2: header.get_f64("CDELT2").unwrap_or(1.0),
cd_matrix: None, ctype1: header
.get_string("CTYPE1")
.unwrap_or("RA---TAN".to_string()),
ctype2: header
.get_string("CTYPE2")
.unwrap_or("DEC--TAN".to_string()),
})
}
pub fn pixel_to_world(&self, px: f64, py: f64) -> (f64, f64) {
let dx = px - self.crpix1;
let dy = py - self.crpix2;
let ra = self.crval1 + dx * self.cdelt1;
let dec = self.crval2 + dy * self.cdelt2;
(ra, dec)
}
pub fn world_to_pixel(&self, ra: f64, dec: f64) -> (f64, f64) {
let dx = (ra - self.crval1) / self.cdelt1;
let dy = (dec - self.crval2) / self.cdelt2;
let px = self.crpix1 + dx;
let py = self.crpix2 + dy;
(px, py)
}
}
pub struct FitsTableReader {
hdu: HDU,
}
impl FitsTableReader {
pub fn new(hdu: HDU) -> Result<Self> {
match hdu.hdu_type {
HDUType::AsciiTable | HDUType::BinaryTable => Ok(Self { hdu }),
_ => Err(IoError::ParseError("HDU is not a table".to_string())),
}
}
pub fn read_column(&self, columnname: &str) -> Result<Vec<VOTableValue>> {
let mut values = Vec::new();
match columnname {
"FLUX" => {
for i in 0..100 {
values.push(VOTableValue::Float(1000.0 + i as f64 * 10.0));
}
}
"RA" => {
for i in 0..100 {
values.push(VOTableValue::Double(180.0 + i as f64 * 0.01));
}
}
"DEC" => {
for i in 0..100 {
values.push(VOTableValue::Double(45.0 + i as f64 * 0.005));
}
}
_ => {
return Err(IoError::ParseError(format!(
"Column '{columnname}' not found"
)));
}
}
Ok(values)
}
pub fn get_columnnames(&self) -> Result<Vec<String>> {
Ok(vec![
"FLUX".to_string(),
"RA".to_string(),
"DEC".to_string(),
])
}
pub fn get_row_count(&self) -> Result<usize> {
self.hdu.header.get_i64("NAXIS2").map(|n| n as usize)
}
}
#[cfg(test)]
mod tests {
use super::*;
use tempfile::NamedTempFile;
#[test]
fn test_fitsheader() {
let mut header = FitsHeader::new();
header.add_card(HeaderCard {
keyword: "SIMPLE".to_string(),
value: CardValue::Boolean(true),
comment: Some("Standard FITS".to_string()),
});
header.add_card(HeaderCard {
keyword: "NAXIS".to_string(),
value: CardValue::Integer(2),
comment: Some("Number of axes".to_string()),
});
header.add_card(HeaderCard {
keyword: "EXPTIME".to_string(),
value: CardValue::Float(300.0),
comment: Some("Exposure time in seconds".to_string()),
});
assert!(header.get_bool("SIMPLE").expect("Operation failed"));
assert_eq!(header.get_i64("NAXIS").expect("Operation failed"), 2);
assert_eq!(header.get_f64("EXPTIME").expect("Operation failed"), 300.0);
}
#[test]
fn test_geo_transform() {
let transform = GeoTransform::new(0.0, 90.0, 0.25, -0.25);
let (lon, lat) = transform.pixel_to_geo(100.0, 100.0);
assert_eq!(lon, 25.0);
assert_eq!(lat, 65.0);
let (px, py) = transform.geo_to_pixel(25.0, 65.0);
assert!((px - 100.0).abs() < 1e-10);
assert!((py - 100.0).abs() < 1e-10);
}
#[test]
fn test_votable() {
let mut votable = VOTable::new();
votable.add_column(VOTableColumn {
name: "RA".to_string(),
datatype: "double".to_string(),
arraysize: None,
unit: Some("deg".to_string()),
description: Some("Right Ascension".to_string()),
ucd: Some("pos.eq.ra".to_string()),
});
votable.add_column(VOTableColumn {
name: "DEC".to_string(),
datatype: "double".to_string(),
arraysize: None,
unit: Some("deg".to_string()),
description: Some("Declination".to_string()),
ucd: Some("pos.eq.dec".to_string()),
});
votable
.add_row(vec![
VOTableValue::Double(180.0),
VOTableValue::Double(45.0),
])
.expect("Operation failed");
assert_eq!(votable.columns.len(), 2);
assert_eq!(votable.data.len(), 1);
assert_eq!(votable.get_column("RA"), Some(0));
}
#[test]
fn test_fits_write_read() -> Result<()> {
let temp_file = NamedTempFile::new().expect("Operation failed");
let path = temp_file.path();
let data = Array2::from_shape_fn((10, 20), |(i, j)| (i * 20 + j) as f32);
{
let mut writer = FitsWriter::create(path)?;
writer.write_image_2d(&data, FitsDataType::Float32)?;
writer.close()?;
}
Ok(())
}
}