PURIFY
Next-generation radio interferometric imaging
Classes | Functions
purify::pfitsio Namespace Reference

Classes

struct  header_params
 

Functions

void write_key (fitsfile *fptr, const std::string &key, const std::string &value, const std::string &comment, int *status)
 write key to fits file header More...
 
void write_key (fitsfile *fptr, const std::string &key, const char *value, const std::string &comment, int *status)
 
void write_history (fitsfile *fptr, const std::string &context, const std::string &history, int *status)
 write history to fits file More...
 
void write2d (const Image< t_real > &eigen_image, const pfitsio::header_params &header, const bool &overwrite)
 Write image to fits file. More...
 
void write2d (const Image< t_real > &image, const std::string &fits_name, const std::string &pix_units="Jy/Beam", const bool &overwrite=true)
 Write image to fits file. More...
 
void write3d (const std::vector< Image< t_real >> &image, const pfitsio::header_params &header, const bool &overwrite=true)
 Write cube to fits file using header information. More...
 
void write3d (const std::vector< Image< t_real >> &image, const std::string &fits_name, const std::string &pix_units="Jy/Beam", const bool &overwrite=true)
 Write cube to fits file. More...
 
std::vector< Image< t_complex > > read3d (const std::string &fits_name)
 Read cube from fits file. More...
 
Image< t_complex > read2d (const std::string &fits_name)
 Read image from fits file. More...
 
template<class T >
std::enable_if< std::is_scalar< T >::value, void >::type write_history (fitsfile *fptr, const std::string &context, const T &history, int *status)
 
template<class T >
std::enable_if< std::is_scalar< T >::value, void >::type write_key (fitsfile *fptr, const std::string &key, const T &value, const std::string &comment, int *status)
 
template<typename T >
T read_key (fitsfile *fptr, const std::string &key, int *status)
 read fits key and return as tuple More...
 
template<typename DERIVED >
void write2d (const Eigen::EigenBase< DERIVED > &input, int nx, int ny, const std::string &fits_name, const std::string &pix_units="Jy/Beam", const bool &overwrite=true)
 
template<typename DERIVED >
void write2d (const Eigen::EigenBase< DERIVED > &input, int nx, int ny, const pfitsio::header_params &header, const bool overwrite=true)
 
template<typename DERIVED >
void write3d (const std::vector< Eigen::EigenBase< DERIVED >> &input, int nx, int ny, const std::string &fits_name, const std::string &pix_units="Jy/Beam", const bool &overwrite=true)
 
template<typename DERIVED >
void write3d (const std::vector< Eigen::EigenBase< DERIVED >> &input, int nx, int ny, const pfitsio::header_params &header, const bool &overwrite=true)
 
template<typename T >
std::enable_if< std::is_same< t_real, typename T::Scalar >::value, void >::type write3d (const Eigen::EigenBase< T > &image, const t_int rows, const t_int cols, const t_int chans, const std::string &fits_name, const std::string &pix_units="Jy/Beam", const bool &overwrite=true)
 
template<typename T >
std::enable_if< std::is_same< t_real, typename T::Scalar >::value, void >::type write3d (const Eigen::EigenBase< T > &image, const t_int rows, const t_int cols, const t_int chans, const pfitsio::header_params &header, const bool &overwrite)
 write 3d fits cube with header More...
 
template<typename T , typename = std::enable_if_t<std::is_same_v<double, typename T::Scalar>>>
void read3d (const std::string &fits_name, Eigen::EigenBase< T > &output, int &rows, int &cols, int &channels, int &pols)
 
utilities::vis_params read_uvfits (const std::vector< std::string > &names, const bool flag=true, const stokes pol=stokes::I)
 Read uvfits files from name of vector. More...
 
utilities::vis_params read_uvfits (const std::string &vis_name2, const utilities::vis_params &u1, const bool flag=true, const stokes pol=stokes::I)
 Reads in combine visiblities from uvfits files. More...
 
utilities::vis_params read_uvfits (const std::string &filename, const bool flag=true, const stokes pol=stokes::I)
 Read uvfits file. More...
 
utilities::vis_params filter_and_combine (const utilities::vis_params &input, const utilities::vis_params &input2, const Vector< t_complex > &stokes_transform, const std::function< bool(t_real, t_real, t_real, t_complex, t_complex, t_real, t_real, t_real, t_complex, t_complex)> &filter=[](const t_real, const t_real, const t_real, const t_complex vis1, const t_complex weight1, const t_real, const t_real, const t_real, const t_complex vis2, const t_complex weight2) { return(weight1.real() > 0.) and(weight2.real() > 0.) and(std::abs(vis1) > 1e-20) and(std::abs(vis2) > 1e-20) and(!std::isnan(vis1.real()) and !std::isnan(vis1.imag())) and(!std::isnan(vis2.real()) and !std::isnan(vis2.imag()));})
 Remove visibilities with zero weighting. More...
 
utilities::vis_params read_polarisation_with_flagging (const Vector< t_real > &data, const Matrix< t_real > &coords, const Vector< t_real > &frequencies, const t_uint pol_index1, const t_uint pol_index2, const t_uint pols, const t_uint baselines, const t_uint channels, const Vector< t_complex > stokes_transform, const std::function< bool(t_complex, t_complex, t_complex, t_complex)> &filter)
 read polarisation with flaggging More...
 
utilities::vis_params read_polarisation (const Vector< t_real > &data, const Matrix< t_real > &coords, const Vector< t_real > &frequencies, const t_uint pol_index1, const t_uint pols, const t_uint baselines, const t_uint channels)
 read polarisation data from uvfits data More...
 
Vector< t_real > read_uvfits_freq (fitsfile *fptr, int *status, const int &col)
 read frequencies for each channel More...
 
void read_uvfits_freq (fitsfile *fptr, int *status, Vector< t_real > &output, const int &col)
 
Vector< t_real > read_uvfits_data (fitsfile *fptr, int *status, const std::vector< int > &naxis, const int &baselines)
 read data from uvfits file More...
 
void read_uvfits_data (fitsfile *fptr, int *status, const std::vector< int > &naxis, const int &baselines, Vector< t_real > &output)
 
Matrix< t_real > read_uvfits_coords (fitsfile *fptr, int *status, const int &groups, const int &pcount)
 read coordinates from uvfits file More...
 
t_real read_value_from_data (const Vector< t_real > &data, const t_uint col, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
 read value from data More...
 
t_complex read_vis_from_data (const Vector< t_real > &data, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
 return visibility for given baseline More...
 
t_complex read_weight_from_data (const Vector< t_real > &data, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
 return weight for given baseline More...
 
void read_uvfits_coords (fitsfile *fptr, int *status, const int &pcount, const int &groups, Matrix< t_real > &output)
 
void read_fits_keys (fitsfile *fptr, int *status)
 Read uvfits keys out. More...
 

Function Documentation

◆ filter_and_combine()

utilities::vis_params purify::pfitsio::filter_and_combine ( const utilities::vis_params &  input,
const utilities::vis_params &  input2,
const Vector< t_complex > &  stokes_transform,
const std::function< bool(t_real, t_real, t_real, t_complex, t_complex, t_real, t_real, t_real, t_complex, t_complex)> &  filter 
)

Remove visibilities with zero weighting.

Definition at line 175 of file uvfits.cc.

179  {
180  t_uint size = 0;
181  if (stokes_transform.size() != 2)
182  throw std::runtime_error("stokes transform is not the right size (!= 2).");
183  if (stokes_transform.isZero()) throw std::runtime_error("stokes transform not valid.");
184  if (input.size() != input2.size())
185  throw std::runtime_error("Cannot apply filter to each data set, sizes do not match.");
186  for (int i = 0; i < input.size(); i++) {
187  if (filter(input.u(i), input.v(i), input.w(i), input.vis(i), input.weights(i), input2.u(i),
188  input2.v(i), input2.w(i), input2.vis(i), input2.weights(i))) {
189  if (std::abs(input.u(i) - input2.u(i)) > 1e-6)
190  throw std::runtime_error("baselines don't match for filter.");
191  if (std::abs(input.v(i) - input2.v(i)) > 1e-6)
192  throw std::runtime_error("baselines don't match for filter.");
193  if (std::abs(input.w(i) - input2.w(i)) > 1e-6)
194  throw std::runtime_error("baselines don't match for filter.");
195  size++;
196  }
197  }
198 
199  PURIFY_LOW_LOG("Applying flags: Keeping {} out of {} data points.", size, input.size());
200  utilities::vis_params output(Vector<t_real>::Zero(size), Vector<t_real>::Zero(size),
201  Vector<t_real>::Zero(size), Vector<t_complex>::Zero(size),
202  Vector<t_complex>::Zero(size), input.units, input.ra, input.dec,
203  input.average_frequency);
204  output.time = Vector<t_real>::Zero(size);
205  output.baseline = Vector<t_uint>::Zero(size);
206  output.frequencies = input.frequencies;
207  t_uint count = 0;
208  for (t_uint i = 0; i < input.size(); i++)
209  if (filter(input.u(i), input.v(i), input.w(i), input.vis(i), input.weights(i), input2.u(i),
210  input2.v(i), input2.w(i), input2.vis(i), input2.weights(i))) {
211  output.u(count) = input.u(i);
212  output.v(count) = input.v(i);
213  output.w(count) = input.w(i);
214  output.vis(count) = input.vis(i) * stokes_transform(0) + input2.vis(i) * stokes_transform(1);
215  output.weights(count) = 1. / std::sqrt(1. / input.weights(i) * stokes_transform(0) +
216  1. / input2.weights(i) * stokes_transform(1));
217  output.baseline(count) = input.baseline(i);
218  output.time(count) = input.time(i);
219  count++;
220  }
221 
222  return output;
223 }
PURIFY_LOW_LOG
#define PURIFY_LOW_LOG(...)
Low priority message.
Definition: logging.h:207

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::baseline, purify::utilities::vis_params::dec, purify::utilities::vis_params::frequencies, PURIFY_LOW_LOG, purify::utilities::vis_params::ra, purify::utilities::vis_params::size(), purify::utilities::vis_params::time, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

◆ read2d()

Image< t_complex > purify::pfitsio::read2d ( const std::string &  fits_name)

Read image from fits file.

Definition at line 109 of file pfitsio.cc.

109  {
110  /*
111  Reads in an image from a fits file and returns the image.
112 
113  fits_name:: name of fits file
114  */
115 
116  const std::vector<Image<t_complex>> images = read3d(fits_name);
117  return images.at(0);
118 }
purify::pfitsio::read3d
std::vector< Image< t_complex > > read3d(const std::string &fits_name)
Read cube from fits file.
Definition: pfitsio.cc:95

References read3d().

Referenced by getInputData(), main(), and TEST_CASE().

◆ read3d() [1/2]

std::vector< Image< t_complex > > purify::pfitsio::read3d ( const std::string &  fits_name)

Read cube from fits file.

Definition at line 95 of file pfitsio.cc.

95  {
96  std::vector<Image<t_complex>> eigen_images;
97  Vector<double> image;
98  int rows, cols, channels, pols = 1;
99  read3d<Vector<double>>(fits_name, image, rows, cols, channels, pols);
100  for (int i = 0; i < channels; i++) {
101  Vector<t_complex> eigen_image = Vector<t_complex>::Zero(rows * cols);
102  eigen_image.real() = image.segment(i * rows * cols, rows * cols);
103  eigen_images.push_back(Image<t_complex>::Map(eigen_image.data(), rows, cols));
104  }
105  return eigen_images;
106 }

Referenced by getInputData(), and read2d().

◆ read3d() [2/2]

template<typename T , typename = std::enable_if_t<std::is_same_v<double, typename T::Scalar>>>
void purify::pfitsio::read3d ( const std::string &  fits_name,
Eigen::EigenBase< T > &  output,
int &  rows,
int &  cols,
int &  channels,
int &  pols 
)

Definition at line 291 of file pfitsio.h.

292  {
293  /*
294  Reads in a cube from a fits file and returns the vector of images.
295  fits_name:: name of fits file
296  */
297 
298  fitsfile *fptr;
299  int status = 0;
300  PURIFY_LOW_LOG("Reading fits file {}", fits_name);
301  if (fits_open_image(&fptr, fits_name.c_str(), READONLY, &status))
302  throw std::runtime_error("Error opening image " + fits_name);
303  const int naxis = pfitsio::read_key<int>(fptr, "NAXIS", &status);
304  if (naxis < 1) throw std::runtime_error("Image contains zero axes.");
305  rows = pfitsio::read_key<int>(fptr, "NAXIS1", &status);
306  cols = (naxis > 1) ? pfitsio::read_key<int>(fptr, "NAXIS2", &status) : 1;
307  channels = (naxis > 2) ? pfitsio::read_key<int>(fptr, "NAXIS3", &status) : 1;
308  pols = (naxis > 3) ? pfitsio::read_key<int>(fptr, "NAXIS4", &status) : 1;
309  PURIFY_LOW_LOG("Axes {}", naxis);
310  std::vector<long> fpixel(naxis, 1);
311  PURIFY_LOW_LOG("Dimensions {}x{}x{}x{}", rows, cols, channels, pols);
312  if (pols > 1) throw std::runtime_error("Too many polarisations when reading " + fits_name);
313  t_real nulval = 0;
314  int anynul = 0;
315  output.derived() = Vector<typename T::Scalar>::Zero(rows * cols * channels * pols);
316  if (fits_read_pix(fptr, TDOUBLE, fpixel.data(), static_cast<long>(output.size()), &nulval,
317  output.derived().data(), &anynul, &status))
318  throw std::runtime_error("Error reading data from image " + fits_name);
319  if (anynul) PURIFY_LOW_LOG("There are bad values when reading " + fits_name);
320  if (fits_close_file(fptr, &status))
321  throw std::runtime_error("Problem closing fits file: " + fits_name);
322 }

References PURIFY_LOW_LOG.

◆ read_fits_keys()

void purify::pfitsio::read_fits_keys ( fitsfile *  fptr,
int *  status 
)

Read uvfits keys out.

Definition at line 405 of file uvfits.cc.

405  {
406  std::shared_ptr<char> card =
407  std::shared_ptr<char>(new char[FLEN_CARD], [](char *ptr) { delete[] ptr; });
408  int nkeys;
409  fits_get_hdrspace(fptr, &nkeys, NULL, status);
410 
411  for (int ii = 1; ii <= nkeys; ii++) {
412  fits_read_record(fptr, ii, card.get(), status); /* read keyword */
413  std::printf("%s\n", card.get());
414  }
415  std::printf("END\n\n"); /* terminate listing with END */
416 }

◆ read_key()

template<typename T >
T purify::pfitsio::read_key ( fitsfile *  fptr,
const std::string &  key,
int *  status 
)

read fits key and return as tuple

Definition at line 122 of file pfitsio.h.

122  {
123  T value;
124  char comment[FLEN_COMMENT];
125  if constexpr (std::is_same_v<T, std::string>) {
126  if (fits_read_key(fptr, TSTRING, key.data(), value.data(), comment, status)) {
127  throw std::runtime_error("Error reading value from key " + key);
128  }
129  } else {
130  int datatype = 0;
131  if (std::is_same<T, double>::value)
132  datatype = TDOUBLE;
133  else if (std::is_same_v<T, float>)
134  datatype = TFLOAT;
135  else if (std::is_same_v<T, int>)
136  datatype = TINT;
137  else if (std::is_same_v<T, bool>)
138  datatype = TLOGICAL;
139  else {
140  throw std::runtime_error("Key type of value not supported by PURIFY template for " + key);
141  }
142 
143  if (fits_read_key(fptr, datatype, key.data(), &value, comment, status)) {
144  throw std::runtime_error("Error reading value from key " + key);
145  }
146  }
147  return value;
148 }

◆ read_polarisation()

utilities::vis_params purify::pfitsio::read_polarisation ( const Vector< t_real > &  data,
const Matrix< t_real > &  coords,
const Vector< t_real > &  frequencies,
const t_uint  pol_index1,
const t_uint  pols,
const t_uint  baselines,
const t_uint  channels 
)

read polarisation data from uvfits data

Definition at line 284 of file uvfits.cc.

287  {
288  const t_uint stride = channels * baselines;
289  utilities::vis_params uv_data;
290  uv_data.u = Vector<t_real>::Zero(stride);
291  uv_data.v = Vector<t_real>::Zero(stride);
292  uv_data.w = Vector<t_real>::Zero(stride);
293  uv_data.time = Vector<t_real>::Zero(stride);
294  uv_data.baseline = Vector<t_uint>::Zero(stride);
295  uv_data.vis = Vector<t_complex>::Zero(stride);
296  uv_data.weights = Vector<t_complex>::Zero(stride);
297  if (pol_index1 >= pols) throw std::runtime_error("Polarisation index out of bounds.");
298  // reading in data
299  t_uint count = 0;
300  for (t_uint c = 0; c < channels; c++)
301  for (t_uint b = 0; b < baselines; b++) {
302  uv_data.vis(count) = read_vis_from_data(data, pol_index1, pols, c, channels, b, baselines);
303  uv_data.weights(count) =
304  read_weight_from_data(data, pol_index1, pols, c, channels, b, baselines);
305  count++;
306  }
307  uv_data.average_frequency = frequencies.array().mean();
308  for (t_uint i = 0; i < frequencies.size(); i++) {
309  uv_data.u.segment(i * baselines, baselines) = coords.row(0) * frequencies(i);
310  uv_data.v.segment(i * baselines, baselines) = -coords.row(1) * frequencies(i);
311  uv_data.w.segment(i * baselines, baselines) = coords.row(2) * frequencies(i);
312  uv_data.baseline.segment(i * baselines, baselines) = coords.row(3).cast<t_uint>();
313  uv_data.time.segment(i * baselines, baselines) = coords.row(4);
314  }
315  uv_data.frequencies = frequencies;
316  uv_data.units = utilities::vis_units::lambda;
317  return uv_data;
318 }
purify::constant::c
const t_real c
speed of light in vacuum
Definition: types.h:72
purify::pfitsio::read_weight_from_data
t_complex read_weight_from_data(const Vector< t_real > &data, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
return weight for given baseline
Definition: uvfits.cc:388
purify::pfitsio::read_vis_from_data
t_complex read_vis_from_data(const Vector< t_real > &data, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
return visibility for given baseline
Definition: uvfits.cc:381

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::baseline, purify::constant::c, purify::utilities::vis_params::frequencies, purify::utilities::lambda, read_vis_from_data(), read_weight_from_data(), purify::utilities::vis_params::time, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

◆ read_polarisation_with_flagging()

utilities::vis_params purify::pfitsio::read_polarisation_with_flagging ( const Vector< t_real > &  data,
const Matrix< t_real > &  coords,
const Vector< t_real > &  frequencies,
const t_uint  pol_index1,
const t_uint  pol_index2,
const t_uint  pols,
const t_uint  baselines,
const t_uint  channels,
const Vector< t_complex >  stokes_transform,
const std::function< bool(t_complex, t_complex, t_complex, t_complex)> &  filter 
)

read polarisation with flaggging

Definition at line 225 of file uvfits.cc.

229  {
230  t_uint count = 0;
231  for (t_uint c = 0; c < channels; c++)
232  for (t_uint b = 0; b < baselines; b++) {
233  t_complex const weight1 =
234  read_weight_from_data(data, pol_index1, pols, c, channels, b, baselines);
235  t_complex const weight2 =
236  read_weight_from_data(data, pol_index2, pols, c, channels, b, baselines);
237  t_complex const vis1 = read_vis_from_data(data, pol_index1, pols, c, channels, b, baselines);
238  t_complex const vis2 = read_vis_from_data(data, pol_index2, pols, c, channels, b, baselines);
239  if (filter(vis1, weight1, vis2, weight2)) count++;
240  }
241  const t_uint stride = count;
242  utilities::vis_params uv_data;
243  uv_data.u = Vector<t_real>::Zero(stride);
244  uv_data.v = Vector<t_real>::Zero(stride);
245  uv_data.w = Vector<t_real>::Zero(stride);
246  uv_data.time = Vector<t_real>::Zero(stride);
247  uv_data.baseline = Vector<t_uint>::Zero(stride);
248  uv_data.vis = Vector<t_complex>::Zero(stride);
249  uv_data.weights = Vector<t_complex>::Zero(stride);
250  uv_data.average_frequency = frequencies.array().mean();
251  if (pol_index1 >= pols) throw std::runtime_error("Polarisation index out of bounds.");
252  if (pol_index2 >= pols) throw std::runtime_error("Polarisation index out of bounds.");
253  PURIFY_LOW_LOG("Applying flags: Keeping {} out of {} data points.", stride, channels * baselines);
254  // reading in data
255  count = 0;
256  for (t_uint c = 0; c < channels; c++)
257  for (t_uint b = 0; b < baselines; b++) {
258  t_complex const weight1 =
259  read_weight_from_data(data, pol_index1, pols, c, channels, b, baselines);
260  t_complex const weight2 =
261  read_weight_from_data(data, pol_index2, pols, c, channels, b, baselines);
262  t_complex const vis1 = read_vis_from_data(data, pol_index1, pols, c, channels, b, baselines);
263  t_complex const vis2 = read_vis_from_data(data, pol_index2, pols, c, channels, b, baselines);
264  // checking flags of both polarisations
265  if (filter(vis1, weight1, vis2, weight2)) {
266  uv_data.vis(count) = uv_data.vis(count) =
267  vis1 * stokes_transform(0) + vis2 * stokes_transform(1);
268  uv_data.weights(count) =
269  1. / std::sqrt(1. / weight1 * stokes_transform(0) + 1. / weight2 * stokes_transform(1));
270  uv_data.u(count) = coords(0, b) * frequencies(c);
271  uv_data.v(count) = -coords(1, b) * frequencies(c);
272  uv_data.w(count) = coords(2, b) * frequencies(c);
273  uv_data.baseline(count) = static_cast<t_uint>(coords(3, b));
274  uv_data.time(count) = coords(4, b);
275  count++;
276  }
277  }
278  assert(count == stride);
279  uv_data.frequencies = frequencies;
280  uv_data.units = utilities::vis_units::lambda;
281  PURIFY_MEDIUM_LOG("All Data Read!");
282  return uv_data;
283 }
PURIFY_MEDIUM_LOG
#define PURIFY_MEDIUM_LOG(...)
Medium priority message.
Definition: logging.h:205

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::baseline, purify::constant::c, purify::utilities::vis_params::frequencies, purify::utilities::lambda, PURIFY_LOW_LOG, PURIFY_MEDIUM_LOG, read_vis_from_data(), read_weight_from_data(), purify::utilities::vis_params::time, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by read_uvfits().

◆ read_uvfits() [1/3]

utilities::vis_params purify::pfitsio::read_uvfits ( const std::string &  filename,
const bool  flag,
const stokes  pol 
)

Read uvfits file.

Definition at line 49 of file uvfits.cc.

49  {
50  fitsfile *fptr;
51  int status = 0;
52  int hdupos;
53  int baselines;
54  int naxes;
55  int pcount;
56  t_real ra = 0;
57  t_real dec = 0;
58  std::shared_ptr<char> comment =
59  std::shared_ptr<char>(new char[FLEN_CARD], [](char *ptr) { delete[] ptr; });
60  PURIFY_MEDIUM_LOG("Reading uvfits {}", filename);
61  if (fits_open_file(&fptr, filename.c_str(), READONLY, &status))
62  throw std::runtime_error("Could not open file " + filename);
63 
64  int hdutype;
65  if (fits_movabs_hdu(fptr, 1, &hdutype, &status)) throw std::runtime_error("Error changing HDU.");
66  if (hdutype != IMAGE_HDU)
67  throw std::runtime_error("HDU 1 not expected type " + std::to_string(hdutype));
68 
69  fits_read_key(fptr, TDOUBLE, "CRVAL5", &ra, comment.get(), &status);
70  fits_read_key(fptr, TDOUBLE, "CRVAL6", &dec, comment.get(), &status);
71  fits_read_key(fptr, TINT, "GCOUNT", &baselines, comment.get(), &status);
72  fits_read_key(fptr, TINT, "NAXIS", &naxes, comment.get(), &status);
73  fits_read_key(fptr, TINT, "PCOUNT", &pcount, comment.get(), &status);
74  if (naxes == 0) throw std::runtime_error("No axes in header... ");
75  if (pcount == 0) throw std::runtime_error("No uvw or time coordinates in header... ");
76  t_uint total = 1;
77  std::vector<int> naxis(naxes, 0);
78  // filling up axis sizes
79  for (int i = 1; i < naxes; i++) {
80  int n;
81  std::string key = "NAXIS" + std::to_string(i + 1);
82  fits_read_key(fptr, TINT, key.c_str(), &n, comment.get(), &status);
83  naxis[i] = n;
84  PURIFY_LOW_LOG("NAXIS {} Size: {}", i, n);
85  total *= n;
86  }
87  const t_uint channels = naxis.at(3);
88  const t_uint pols = naxis.at(2);
89  const t_uint ifs = naxis.at(4);
90  const t_uint pointings_num = naxis.at(5);
91  PURIFY_LOW_LOG("RA: {} deg, Dec.: {} deg", ra, dec);
92  PURIFY_MEDIUM_LOG("Pointings: {}", pointings_num);
93  PURIFY_MEDIUM_LOG("IFs: {}", ifs);
94  PURIFY_MEDIUM_LOG("Baselines: {}", baselines);
95  PURIFY_MEDIUM_LOG("Channels: {}", channels);
96  PURIFY_MEDIUM_LOG("Polarisations: {}", pols);
97  PURIFY_LOW_LOG(
98  "Purify currently assumes uvfits files are in XX YY XY YX format. Also will only read Stokes "
99  "I for now.");
100  PURIFY_MEDIUM_LOG("Total data per baseline: {}", total);
101  if (pointings_num > 1) throw std::runtime_error("More than one pointing is not supported.");
102  if (ifs > 1) throw std::runtime_error("More than one IF is not supported.");
103  PURIFY_LOW_LOG("Reading Data...");
104 
105  const Matrix<t_real> coords = read_uvfits_coords(fptr, &status, baselines, pcount);
106 
107  // (real, imag, weight, pol, frequency) is the order
108  const Vector<t_real> data = read_uvfits_data(fptr, &status, naxis, baselines);
109  const Vector<t_real> frequencies = read_uvfits_freq(fptr, &status, 4);
110  if (frequencies.size() != channels)
111  throw std::runtime_error("Number of frequencies doesn't match number of channels. " +
112  std::to_string(frequencies.size()) + "!=" + std::to_string(channels));
113  int pol_index1;
114  int pol_index2;
115  Vector<t_complex> stokes_transform = Vector<t_complex>::Zero(2);
116  switch (pol) {
117  case stokes::I:
118  pol_index1 = 0;
119  pol_index2 = 1;
120  stokes_transform(0) = 1. / 2;
121  stokes_transform(1) = 1. / 2;
122  break;
123  case stokes::Q:
124  throw std::runtime_error("Polarisation not supported for reading uvfits.");
125  pol_index1 = 0;
126  pol_index2 = 1;
127  stokes_transform(0) = 1. / 2;
128  stokes_transform(1) = -1. / 2;
129  break;
130  case stokes::U:
131  throw std::runtime_error("Polarisation not supported for reading uvfits.");
132  pol_index1 = 0;
133  pol_index2 = 1;
134  stokes_transform(0) = 1. / 2 * t_complex(0, -1);
135  stokes_transform(1) = 1. / 2 * t_complex(0, -1);
136  break;
137  case stokes::V:
138  throw std::runtime_error("Polarisation not supported for reading uvfits.");
139  pol_index1 = 0;
140  pol_index2 = 1;
141  stokes_transform(0) = 1. / 2 * t_complex(0, -1);
142  stokes_transform(1) = -1. / 2 * t_complex(0, -1);
143  break;
144  default:
145  throw std::runtime_error("Polarisation not supported for reading uvfits.");
146  pol_index1 = 0;
147  pol_index2 = 1;
148  stokes_transform(0) = 1. / 2;
149  stokes_transform(1) = 1. / 2;
150  break;
151  }
152  fits_close_file(fptr, &status);
153  if (status) { /* print any error messages */
154  fits_report_error(stderr, status);
155  throw std::runtime_error("Error reading fits file...");
156  }
157  auto uv_data = read_polarisation_with_flagging(
158  data, coords, frequencies, pol_index1, pol_index2, pols, baselines, channels,
159  stokes_transform,
160  [flag](const t_complex vis1, const t_complex weight1, const t_complex vis2,
161  const t_complex weight2) {
162  if (flag)
163  return (weight1.real() > 0.) and (weight2.real() > 0.) and (std::abs(vis1) > 1e-20) and
164  (std::abs(vis2) > 1e-20) and
165  (!std::isnan(vis1.real()) and !std::isnan(vis1.imag())) and
166  (!std::isnan(vis2.real()) and !std::isnan(vis2.imag()));
167  else
168  return true;
169  });
170  uv_data.ra = ra;
171  uv_data.dec = dec;
172  return uv_data;
173 }
purify::pfitsio::read_uvfits_coords
void read_uvfits_coords(fitsfile *fptr, int *status, const int &pcount, const int &groups, Matrix< t_real > &output)
Definition: uvfits.cc:394
purify::pfitsio::read_uvfits_data
void read_uvfits_data(fitsfile *fptr, int *status, const std::vector< int > &naxis, const int &baselines, Vector< t_real > &output)
Definition: uvfits.cc:353
purify::pfitsio::read_polarisation_with_flagging
utilities::vis_params read_polarisation_with_flagging(const Vector< t_real > &data, const Matrix< t_real > &coords, const Vector< t_real > &frequencies, const t_uint pol_index1, const t_uint pol_index2, const t_uint pols, const t_uint baselines, const t_uint channels, const Vector< t_complex > stokes_transform, const std::function< bool(t_complex, t_complex, t_complex, t_complex)> &filter)
read polarisation with flaggging
Definition: uvfits.cc:225
purify::pfitsio::read_uvfits_freq
void read_uvfits_freq(fitsfile *fptr, int *status, Vector< t_real > &output, const int &col)
Definition: uvfits.cc:326
purify::utilities::vis_params::ra
t_real ra
Definition: uvw_utilities.h:25

References purify::I, PURIFY_LOW_LOG, PURIFY_MEDIUM_LOG, purify::Q, purify::utilities::vis_params::ra, read_polarisation_with_flagging(), read_uvfits_coords(), read_uvfits_data(), read_uvfits_freq(), purify::U, and purify::V.

◆ read_uvfits() [2/3]

utilities::vis_params purify::pfitsio::read_uvfits ( const std::string &  vis_name2,
const utilities::vis_params &  uv1,
const bool  flag,
const stokes  pol 
)

Reads in combine visiblities from uvfits files.

Definition at line 20 of file uvfits.cc.

21  {
22  utilities::vis_params uv;
23  const bool w_term = not uv1.w.isZero(0);
24  const auto uv2 = read_uvfits(vis_name2, flag, pol);
25  if (std::abs(uv1.ra - uv2.ra) > 1e-6)
26  throw std::runtime_error(vis_name2 + ": wrong RA in pointing.");
27  if (std::abs(uv1.dec - uv2.dec) > 1e-6)
28  throw std::runtime_error(vis_name2 + ": wrong DEC in pointing.");
29  uv.ra = uv1.ra;
30  uv.dec = uv1.dec;
31  uv.u = Vector<t_real>::Zero(uv1.size() + uv2.size());
32  uv.v = Vector<t_real>::Zero(uv1.size() + uv2.size());
33  uv.w = Vector<t_real>::Zero(uv1.size() + uv2.size());
34  uv.vis = Vector<t_complex>::Zero(uv1.size() + uv2.size());
35  uv.weights = Vector<t_complex>::Zero(uv1.size() + uv2.size());
36  uv.u.segment(0, uv1.size()) = uv1.u;
37  uv.v.segment(0, uv1.size()) = uv1.v;
38  uv.w.segment(0, uv1.size()) = uv1.w;
39  uv.vis.segment(0, uv1.size()) = uv1.vis;
40  uv.weights.segment(0, uv1.size()) = uv1.weights;
41  uv.u.segment(uv1.size(), uv2.size()) = uv2.u;
42  uv.v.segment(uv1.size(), uv2.size()) = uv2.v;
43  uv.w.segment(uv1.size(), uv2.size()) = uv2.w;
44  uv.vis.segment(uv1.size(), uv2.size()) = uv2.vis;
45  uv.weights.segment(uv1.size(), uv2.size()) = uv2.weights;
46  return uv;
47 }
purify::pfitsio::read_uvfits
utilities::vis_params read_uvfits(const std::string &filename, const bool flag, const stokes pol)
Read uvfits file.
Definition: uvfits.cc:49

References purify::utilities::vis_params::dec, purify::utilities::vis_params::ra, read_uvfits(), purify::utilities::vis_params::size(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

◆ read_uvfits() [3/3]

utilities::vis_params purify::pfitsio::read_uvfits ( const std::vector< std::string > &  names,
const bool  flag,
const stokes  pol 
)

Read uvfits files from name of vector.

Definition at line 12 of file uvfits.cc.

13  {
14  utilities::vis_params output = read_uvfits(names.at(0), flag, pol);
15  if (names.size() == 1) return output;
16  for (int i = 1; i < names.size(); i++) output = read_uvfits(names.at(i), output);
17  return output;
18 }

Referenced by main(), purify::read_measurements::read_measurements(), read_uvfits(), and TEST_CASE().

◆ read_uvfits_coords() [1/2]

Matrix< t_real > purify::pfitsio::read_uvfits_coords ( fitsfile *  fptr,
int *  status,
const int &  groups,
const int &  pcount 
)

read coordinates from uvfits file

Definition at line 367 of file uvfits.cc.

368  {
369  Matrix<t_real> output;
370  read_uvfits_coords(fptr, status, pcount, groups, output);
371  return output;
372 }

Referenced by read_uvfits().

◆ read_uvfits_coords() [2/2]

void purify::pfitsio::read_uvfits_coords ( fitsfile *  fptr,
int *  status,
const int &  pcount,
const int &  groups,
Matrix< t_real > &  output 
)

Definition at line 394 of file uvfits.cc.

395  {
396  output = Matrix<t_real>::Zero(pcount, groups);
397  int nulval = 0;
398  int anynul;
399  // reading in parameters per baseline
400  for (int i = 0; i < groups; i++)
401  fits_read_col(fptr, TDOUBLE, 1, 1 + i, 1, pcount, &nulval, output.col(i).data(), &anynul,
402  status);
403 }

◆ read_uvfits_data() [1/2]

Vector< t_real > purify::pfitsio::read_uvfits_data ( fitsfile *  fptr,
int *  status,
const std::vector< int > &  naxis,
const int &  baselines 
)

read data from uvfits file

Definition at line 346 of file uvfits.cc.

347  {
348  Vector<t_real> output;
349  read_uvfits_data(fptr, status, naxis, baselines, output);
350  return output;
351 }

Referenced by read_uvfits().

◆ read_uvfits_data() [2/2]

void purify::pfitsio::read_uvfits_data ( fitsfile *  fptr,
int *  status,
const std::vector< int > &  naxis,
const int &  baselines,
Vector< t_real > &  output 
)

Definition at line 353 of file uvfits.cc.

354  {
355  long nelements = 1;
356  for (int i = 2; i < naxis.size(); i++) {
357  nelements *= static_cast<long>(naxis.at(i));
358  }
359  if (nelements == 0) throw std::runtime_error("Zero number of elements.");
360  output = Vector<t_real>::Zero(naxis.at(1) * nelements * baselines);
361  int nulval = 0;
362  int anynul = 0;
363  fits_read_col(fptr, TDOUBLE, 2, 1, 1, static_cast<long>(output.size()), &nulval, output.data(),
364  &anynul, status);
365 }

◆ read_uvfits_freq() [1/2]

Vector< t_real > purify::pfitsio::read_uvfits_freq ( fitsfile *  fptr,
int *  status,
const int &  col 
)

read frequencies for each channel

Definition at line 320 of file uvfits.cc.

320  {
321  Vector<t_real> output;
322  read_uvfits_freq(fptr, status, output, col);
323  return output;
324 }

Referenced by read_uvfits().

◆ read_uvfits_freq() [2/2]

void purify::pfitsio::read_uvfits_freq ( fitsfile *  fptr,
int *  status,
Vector< t_real > &  output,
const int &  col 
)

Definition at line 326 of file uvfits.cc.

326  {
327  t_real cfreq;
328  t_real dfreq;
329  int crpix;
330  int nfreq;
331  std::shared_ptr<char> comment =
332  std::shared_ptr<char>(new char[FLEN_CARD], [](char *ptr) { delete[] ptr; });
333  std::string key = "CRVAL" + std::to_string(col);
334  fits_read_key(fptr, TDOUBLE, key.c_str(), &cfreq, comment.get(), status);
335  key = "CDELT" + std::to_string(col);
336  fits_read_key(fptr, TDOUBLE, key.c_str(), &dfreq, comment.get(), status);
337  key = "CRPIX" + std::to_string(col);
338  fits_read_key(fptr, TINT, key.c_str(), &crpix, comment.get(), status);
339  key = "NAXIS" + std::to_string(col);
340  fits_read_key(fptr, TINT, key.c_str(), &nfreq, comment.get(), status);
341  if (nfreq == 0) throw std::runtime_error("Wrong number of channels read from header.");
342  output = Vector<t_real>::Zero(nfreq);
343  for (int i = 0; i < output.size(); i++) output(i) = (i - nfreq * 0.5) * dfreq + cfreq;
344 }

◆ read_value_from_data()

t_real purify::pfitsio::read_value_from_data ( const Vector< t_real > &  data,
const t_uint  col,
const t_uint  pol,
const t_uint  pols,
const t_uint  chan,
const t_uint  chans,
const t_uint  baseline,
const t_uint  baselines 
)

read value from data

Definition at line 374 of file uvfits.cc.

376  {
377  const t_uint index = col + 3 * (pol + pols * (chan + chans * baseline));
378  return data(index);
379 }

Referenced by read_vis_from_data(), and read_weight_from_data().

◆ read_vis_from_data()

t_complex purify::pfitsio::read_vis_from_data ( const Vector< t_real > &  data,
const t_uint  pol,
const t_uint  pols,
const t_uint  chan,
const t_uint  chans,
const t_uint  baseline,
const t_uint  baselines 
)

return visibility for given baseline

Definition at line 381 of file uvfits.cc.

383  {
384  return t_complex(read_value_from_data(data, 0, pol, pols, chan, chans, baseline, baselines),
385  read_value_from_data(data, 1, pol, pols, chan, chans, baseline, baselines));
386 }
purify::pfitsio::read_value_from_data
t_real read_value_from_data(const Vector< t_real > &data, const t_uint col, const t_uint pol, const t_uint pols, const t_uint chan, const t_uint chans, const t_uint baseline, const t_uint baselines)
read value from data
Definition: uvfits.cc:374

References read_value_from_data().

Referenced by read_polarisation(), and read_polarisation_with_flagging().

◆ read_weight_from_data()

t_complex purify::pfitsio::read_weight_from_data ( const Vector< t_real > &  data,
const t_uint  pol,
const t_uint  pols,
const t_uint  chan,
const t_uint  chans,
const t_uint  baseline,
const t_uint  baselines 
)

return weight for given baseline

Definition at line 388 of file uvfits.cc.

390  {
391  return t_complex(read_value_from_data(data, 2, pol, pols, chan, chans, baseline, baselines), 0);
392 }

References read_value_from_data().

Referenced by read_polarisation(), and read_polarisation_with_flagging().

◆ write2d() [1/4]

template<typename DERIVED >
void purify::pfitsio::write2d ( const Eigen::EigenBase< DERIVED > &  input,
int  nx,
int  ny,
const pfitsio::header_params &  header,
const bool  overwrite = true 
)

Definition at line 165 of file pfitsio.h.

166  {
167  Image<t_real> const data = input.derived().real().template cast<t_real>();
168  write2d(Image<t_real>::Map(data.data(), nx, ny), header, overwrite);
169 }
purify::pfitsio::write2d
void write2d(const Eigen::EigenBase< DERIVED > &input, int nx, int ny, const pfitsio::header_params &header, const bool overwrite=true)
Definition: pfitsio.h:165

References write2d().

◆ write2d() [2/4]

template<typename DERIVED >
void purify::pfitsio::write2d ( const Eigen::EigenBase< DERIVED > &  input,
int  nx,
int  ny,
const std::string &  fits_name,
const std::string &  pix_units = "Jy/Beam",
const bool &  overwrite = true 
)

Definition at line 159 of file pfitsio.h.

160  {
161  Image<t_real> const data = input.derived().real().template cast<t_real>();
162  write2d(Image<t_real>::Map(data.data(), nx, ny), fits_name, pix_units, overwrite);
163 }

References write2d().

◆ write2d() [3/4]

void purify::pfitsio::write2d ( const Image< t_real > &  eigen_image,
const pfitsio::header_params &  header,
const bool &  overwrite 
)

Write image to fits file.

Write image to fits file using header information.

Definition at line 30 of file pfitsio.cc.

31  {
32  /*
33  Writes an image to a fits file.
34 
35  image:: image data, a 2d Image.
36  header:: structure containing header information
37  overwrite:: if true, overwrites old fits file with same name
38 
39  */
40 
41  write3d(std::vector<Image<t_real>>(1, eigen_image), header, overwrite);
42 }
purify::pfitsio::write3d
void write3d(const std::vector< Image< t_real >> &eigen_images, const std::string &fits_name, const std::string &pix_units, const bool &overwrite)
Write cube to fits file.
Definition: pfitsio.cc:77

References write3d().

Referenced by purify::factory::add_updater(), main(), padmm(), padmm_factory(), saveDirtyImage(), saveMeasurementEigenVector(), savePSF(), TEST_CASE(), and write2d().

◆ write2d() [4/4]

void purify::pfitsio::write2d ( const Image< t_real > &  eigen_image,
const std::string &  fits_name,
const std::string &  pix_units,
const bool &  overwrite 
)

Write image to fits file.

Definition at line 44 of file pfitsio.cc.

45  {
46  /*
47  Writes an image to a fits file.
48 
49  image:: image data, a 2d Image.
50  fits_name:: string containing the file name of the fits file.
51  pix_units:: units of flux
52  ra:: centre pixel coordinate in ra
53  dec:: centre pixel coordinate in dec
54 
55  */
56  pfitsio::header_params header;
57  header.fits_name = fits_name;
58  header.pix_units = pix_units;
59 
60  write2d(eigen_image, header, overwrite);
61 }
purify::pfitsio::write2d
void write2d(const Image< t_real > &eigen_image, const std::string &fits_name, const std::string &pix_units, const bool &overwrite)
Write image to fits file.
Definition: pfitsio.cc:44

References purify::pfitsio::header_params::fits_name, purify::pfitsio::header_params::pix_units, and write2d().

◆ write3d() [1/6]

template<typename T >
std::enable_if<std::is_same<t_real, typename T::Scalar>::value, void>::type purify::pfitsio::write3d ( const Eigen::EigenBase< T > &  image,
const t_int  rows,
const t_int  cols,
const t_int  chans,
const pfitsio::header_params &  header,
const bool &  overwrite 
)

write 3d fits cube with header

Definition at line 214 of file pfitsio.h.

216  {
217  /*
218  Writes an image to a fits file.
219  image:: image data, a 2d Image.
220  header:: structure containing header information
221  overwrite:: if true, overwrites old fits file with same name
222  */
223  if (image.size() != rows * cols * chans)
224  throw std::runtime_error("image or cube size does not match dimensions.");
225 
226  PURIFY_LOW_LOG("Writing fits file {}", header.fits_name);
227  if (overwrite == true) {
228  remove(header.fits_name.c_str());
229  }
230  fitsfile *fptr;
231  int status = 0;
232  std::vector<long> naxes = {static_cast<long>(rows), static_cast<long>(cols),
233  static_cast<long>(chans), 1};
234  std::vector<long> fpixel = {1, 1, 1, 1};
235 
236  if (fits_create_file(&fptr, header.fits_name.c_str(), &status))
237  throw std::runtime_error("Problem creating fits file:" + header.fits_name);
238  if (fits_create_img(fptr, DOUBLE_IMG, static_cast<t_int>(naxes.size()), naxes.data(), &status))
239  throw std::runtime_error("Problem creating HDU for image in fits file:" + header.fits_name);
240  if (fits_write_pix(fptr, TDOUBLE, fpixel.data(), static_cast<long>(image.size()),
241  const_cast<t_real *>(image.derived().data()), &status))
242  throw std::runtime_error("Problem writing image in fits file:" + header.fits_name);
243 
244 #define PURIFY_MACRO(KEY, VALUE, COMMENT) write_key(fptr, KEY, VALUE, COMMENT, &status);
245 
246  // Writing to fits header
247  PURIFY_MACRO("BUNIT", header.pix_units, ""); // units
248  PURIFY_MACRO("NAXIS", static_cast<t_int>(naxes.size()), ""); // number of axes
249  PURIFY_MACRO("NAXIS1", static_cast<t_int>(naxes.at(0)), "");
250  PURIFY_MACRO("NAXIS2", static_cast<t_int>(naxes.at(1)), "");
251  PURIFY_MACRO("NAXIS3", static_cast<t_int>(naxes.at(2)), "");
252  PURIFY_MACRO("NAXIS4", static_cast<t_int>(naxes.at(3)), "");
253  PURIFY_MACRO("CRPIX1", static_cast<t_int>(std::floor(naxes.at(0) / 2)) + 1, "");
254  PURIFY_MACRO("CRPIX2", static_cast<t_int>(std::floor(naxes.at(1) / 2)) + 1, "");
255  PURIFY_MACRO("CRPIX3", 1, "");
256  PURIFY_MACRO("CRPIX4", 1, "");
257  PURIFY_MACRO("CRVAL1", static_cast<t_real>(header.ra * 180. / purify::constant::pi), "");
258  PURIFY_MACRO("CRVAL2", static_cast<t_real>(header.dec * 180. / purify::constant::pi), "");
259  PURIFY_MACRO("CRVAL3", static_cast<t_real>(header.mean_frequency * 1.), "");
260  PURIFY_MACRO("CRVAL4", static_cast<t_real>(1), "");
261  PURIFY_MACRO("CTYPE1", "RA--SIN", "");
262  PURIFY_MACRO("CTYPE2", "DEC--SIN", "");
263  PURIFY_MACRO("CTYPE3", "FREQ", "");
264  PURIFY_MACRO("CTYPE4", "STOKES", "");
265  PURIFY_MACRO("CDELT1", static_cast<t_real>(-header.cell_x / 3600.), "");
266  PURIFY_MACRO("CDELT2", static_cast<t_real>(header.cell_y / 3600.), "");
267  PURIFY_MACRO("CDELT3", static_cast<t_real>(header.channel_width * 1.), "");
268  PURIFY_MACRO("CDELT4", 1., "");
269  PURIFY_MACRO("BTYPE", "intensity", "");
270  PURIFY_MACRO("EQUINOX", 2000, "");
271 
272 #undef PURIFY_MACRO
273  write_history(fptr, "SOFTWARE", "Purify", &status);
274  write_history(fptr, "PURIFY-NITERS", header.niters, &status);
275  if (header.hasconverged) {
276  write_history(fptr, "PURIFY-CONVERGED", "True", &status);
277  } else {
278  write_history(fptr, "PURIFY-CONVERGED", "False", &status);
279  }
280  write_history(fptr, "PURIFY-RELATIVEVARIATION", header.relative_variation, &status);
281  write_history(fptr, "PURIFY-RESIDUALCONVERGENCE", header.residual_convergence, &status);
282  write_history(fptr, "PURIFY-EPSILON", header.epsilon, &status);
283  if (fits_write_date(fptr, &status))
284  throw std::runtime_error("Problem writing date in fits file:" + header.fits_name);
285 
286  if (fits_close_file(fptr, &status))
287  throw std::runtime_error("Problem closing fits file:" + header.fits_name);
288 }
purify::constant::pi
const t_real pi
mathematical constant
Definition: types.h:70
purify::pfitsio::write_history
std::enable_if< std::is_scalar< T >::value, void >::type write_history(fitsfile *fptr, const std::string &context, const T &history, int *status)
Definition: pfitsio.h:88
PURIFY_MACRO
#define PURIFY_MACRO(VAR, H2)
Definition: pfitsio.h:62

References purify::pfitsio::header_params::cell_x, purify::pfitsio::header_params::cell_y, purify::pfitsio::header_params::channel_width, purify::pfitsio::header_params::dec, purify::pfitsio::header_params::epsilon, purify::pfitsio::header_params::fits_name, purify::pfitsio::header_params::hasconverged, purify::pfitsio::header_params::mean_frequency, purify::pfitsio::header_params::niters, purify::constant::pi, purify::pfitsio::header_params::pix_units, PURIFY_LOW_LOG, PURIFY_MACRO, purify::pfitsio::header_params::ra, purify::pfitsio::header_params::relative_variation, purify::pfitsio::header_params::residual_convergence, and write_history().

◆ write3d() [2/6]

template<typename T >
std::enable_if<std::is_same<t_real, typename T::Scalar>::value, void>::type purify::pfitsio::write3d ( const Eigen::EigenBase< T > &  image,
const t_int  rows,
const t_int  cols,
const t_int  chans,
const std::string &  fits_name,
const std::string &  pix_units = "Jy/Beam",
const bool &  overwrite = true 
)

Definition at line 202 of file pfitsio.h.

205  {
206  pfitsio::header_params header;
207  header.fits_name = fits_name;
208  header.pix_units = pix_units;
209  write3d<T>(image, rows, cols, chans, header, overwrite);
210 }

References purify::pfitsio::header_params::fits_name, and purify::pfitsio::header_params::pix_units.

◆ write3d() [3/6]

template<typename DERIVED >
void purify::pfitsio::write3d ( const std::vector< Eigen::EigenBase< DERIVED >> &  input,
int  nx,
int  ny,
const pfitsio::header_params &  header,
const bool &  overwrite = true 
)

Definition at line 191 of file pfitsio.h.

192  {
193  std::vector<Image<t_real>> images;
194  for (int i = 0; i < input.size(); i++) {
195  Image<t_real> const data = input.derived().real().template cast<t_real>();
196  images.push_back(Image<t_real>::Map(data.data(), nx, ny));
197  }
198  write3d(images, header, overwrite);
199 }
purify::pfitsio::write3d
std::enable_if< std::is_same< t_real, typename T::Scalar >::value, void >::type write3d(const Eigen::EigenBase< T > &image, const t_int rows, const t_int cols, const t_int chans, const pfitsio::header_params &header, const bool &overwrite)
write 3d fits cube with header
Definition: pfitsio.h:214

References write3d().

◆ write3d() [4/6]

template<typename DERIVED >
void purify::pfitsio::write3d ( const std::vector< Eigen::EigenBase< DERIVED >> &  input,
int  nx,
int  ny,
const std::string &  fits_name,
const std::string &  pix_units = "Jy/Beam",
const bool &  overwrite = true 
)

Definition at line 180 of file pfitsio.h.

182  {
183  std::vector<Image<t_real>> images;
184  for (int i = 0; i < input.size(); i++) {
185  Image<t_real> const data = input.derived().real().template cast<t_real>();
186  images.push_back(Image<t_real>::Map(data.data(), nx, ny));
187  }
188  write3d(images, fits_name, pix_units, overwrite);
189 }

References write3d().

◆ write3d() [5/6]

void purify::pfitsio::write3d ( const std::vector< Image< t_real >> &  eigen_images,
const pfitsio::header_params &  header,
const bool &  overwrite 
)

Write cube to fits file using header information.

Definition at line 63 of file pfitsio.cc.

64  {
65  std::vector<long> naxes = {static_cast<long>(eigen_images.at(0).rows()),
66  static_cast<long>(eigen_images.at(0).cols()),
67  static_cast<long>(eigen_images.size()), 1};
68  std::vector<long> fpixel = {1, 1, 1, 1};
69 
70  Vector<t_real> image = Vector<t_real>::Zero(naxes.at(0) * naxes.at(1) * naxes.at(2));
71  for (int i = 0; i < eigen_images.size(); i++)
72  image.segment(i * naxes.at(0) * naxes.at(1), naxes.at(0) * naxes.at(1)) =
73  Vector<t_real>::Map(eigen_images.at(i).data(), naxes.at(0) * naxes.at(1));
74  write3d<Vector<t_real>>(image, naxes.at(0), naxes.at(1), naxes.at(2), header, overwrite);
75 }

Referenced by write2d(), and write3d().

◆ write3d() [6/6]

void purify::pfitsio::write3d ( const std::vector< Image< t_real >> &  eigen_images,
const std::string &  fits_name,
const std::string &  pix_units,
const bool &  overwrite 
)

Write cube to fits file.

Definition at line 77 of file pfitsio.cc.

78  {
79  /*
80  Writes a vector of images to a fits file.
81  image:: image data, a 3d Image.
82  fits_name:: string containing the file name of the fits file.
83  pix_units:: units of flux
84  ra:: centre pixel coordinate in ra
85  dec:: centre pixel coordinate in dec
86 */
87  pfitsio::header_params header;
88  header.fits_name = fits_name;
89  header.pix_units = pix_units;
90 
91  write3d(eigen_images, header, overwrite);
92 }

References purify::pfitsio::header_params::fits_name, purify::pfitsio::header_params::pix_units, and write3d().

◆ write_history() [1/2]

void purify::pfitsio::write_history ( fitsfile *  fptr,
const std::string &  context,
const std::string &  history,
int *  status 
)

write history to fits file

Definition at line 22 of file pfitsio.cc.

23  {
24  const std::string entry = context + ": " + history;
25  if (fits_write_history(fptr, const_cast<char *>(entry.c_str()), status))
26  throw std::runtime_error("Problem writing comments in fits file");
27 }

Referenced by write3d(), and write_history().

◆ write_history() [2/2]

template<class T >
std::enable_if<std::is_scalar<T>::value, void>::type purify::pfitsio::write_history ( fitsfile *  fptr,
const std::string &  context,
const T &  history,
int *  status 
)

Definition at line 88 of file pfitsio.h.

89  {
90  T value = history;
91  write_history(fptr, context, std::to_string(value), status);
92 }

References write_history().

◆ write_key() [1/3]

void purify::pfitsio::write_key ( fitsfile *  fptr,
const std::string &  key,
const char *  value,
const std::string &  comment,
int *  status 
)

Definition at line 17 of file pfitsio.cc.

18  {
19  write_key(fptr, key, std::string(value), comment, status);
20 }
purify::pfitsio::write_key
void write_key(fitsfile *fptr, const std::string &key, const char *value, const std::string &comment, int *status)
Definition: pfitsio.cc:17

References write_key().

◆ write_key() [2/3]

void purify::pfitsio::write_key ( fitsfile *  fptr,
const std::string &  key,
const std::string &  value,
const std::string &  comment,
int *  status 
)

write key to fits file header

Definition at line 8 of file pfitsio.cc.

9  {
10  std::string entry = value;
11  if (fits_write_key(fptr, TSTRING, const_cast<char *>(key.c_str()),
12  reinterpret_cast<void *>(const_cast<char *>(entry.c_str())),
13  const_cast<char *>(comment.c_str()), status))
14  throw std::runtime_error("Problem writing key in fits file: " + key);
15 }

Referenced by write_key().

◆ write_key() [3/3]

template<class T >
std::enable_if<std::is_scalar<T>::value, void>::type purify::pfitsio::write_key ( fitsfile *  fptr,
const std::string &  key,
const T &  value,
const std::string &  comment,
int *  status 
)

Definition at line 95 of file pfitsio.h.

99  {
100  int datatype = 0;
101  if (std::is_same<T, double>::value)
102  datatype = TDOUBLE;
103  else if (std::is_same<T, float>::value)
104  datatype = TFLOAT;
105  else if (std::is_same<T, int>::value)
106  datatype = TINT;
107  else if (std::is_same<T, t_int>::value)
108  datatype = TINT;
109  else if (std::is_same<T, bool>::value)
110  datatype = TLOGICAL;
111  else
112  throw std::runtime_error("Key type of value not supported by PURIFY template for " + key);
113 
114  T entry = value;
115  if (fits_write_key(fptr, datatype, const_cast<char *>(key.c_str()),
116  reinterpret_cast<void *>(&entry), const_cast<char *>(comment.c_str()), status))
117  throw std::runtime_error("Problem writing key in fits file: " + key);
118 }
Generated by doxygen 1.9.1