PyFU is a Pyraf/Python package for mosaicking (co-adding) datacubes that have
mutual pointing and/or spectral offsets. It can also convert linear wavelength
binning to log lambda for kinematic analysis. The input files are required to
be in multi-extension FITS format with one "SCI" image extension each and the
spectral co-ordinate along the third array axis, as with Gemini data. The cubes
must already have linear FITS World co-ordinate systems.

This is not tested as extensively as the official Gemini packages but has been
well used with output from gfcube (and more recently nfcube) in the Gemini IRAF
package.

If you find this work useful, please consider mentioning it in the
acknowledgements or data reduction section of your publication (eg. Version
0.12 of the PyFU 3D data manipulation package was contributed by James Turner
through the Gemini Data Reduction User Forum.").


Requirements
------------

The following dependencies are required to use pyfu:
    NumPy, SciPy, AstroPy >=3.1

The easiest way to get these is to install "pyfu" from Gemini's public conda
channel and conda-forge, as described (for DRAGONS and Gemini IRAF) at:
https://www.gemini.edu/observing/phase-iii/reducing-data/dragons-data-reduction-software


Installation from source  (prefer conda installation instead, as above)
------------------------

Change to the parent directory in which you want to put PyFU's source directory
and unpack the contents of the tar file, eg. to put it in ~/pyfu-0.12/ do:

  cd
  tar zxf pyfu-0.12.tar.gz

PyRAF

  To use the PyRAF interface, set it up in your IRAF loginuser.cl like so
  (then skip to "Usage"):

    reset pyfu = "/where/you/put/it/pyfu-0.12/pyfu/"  # need trailing slash!
    task pyfu.pkg = pyfu$pyfu.cl

  You can also install it into Python as detailed below and adapt the above
  path accordingly (pointing to pythonN.N/site-packages/pyfu/).

Python (including use with NDMapper)

  Configure the Python environment in which you want to install PyFU, eg. 

    For Anaconda:

      conda create -n pyfu_env iraf pyraf   (or similar, IRAF optional)
      conda activate pyfu_env

    For Ureka (very old):

      ur_setup -n pyfu_variant

  Use "which python" to check that you are finding python in the expected
  place, then install the package in the usual way:

    cd pyfu-0.12/
    python setup.py install

  Then do "import pyfu" in Python. This (and the GMOS IRAF package) can
  also be used via the lib.gmos interface in my NDMapper Python package.
  For a bit more information, see the doc strings and/or description of the
  analagous PyRAF interface below.


Usage with PyRAF
----------------

The PyRAF interface is simple to use. Type "pyfu" to load the package.

The mosaicking procedure has 2 steps:

1. Determine any spatial offsets and update the WCS of each input cube
   accordingly, to put them all on the same absolute World co-ordinates, eg.:

     pyfalign cube1,cube2,cube3,cube4 method="correlate"

     - This modifies the FITS header WCS values in place, so make a copy if
       you need a record of the original values.

     - Your pointings must overlap, of course, with some discernable spatial
       feature(s) in common, for offsets to be measurable. There is no support
       for bootstrapping offsets from multiple sources with known relative
       co-ordinates (though you could probably use method="centroid" and add
       an adjustment for the separation to CRVAL1/2 afterwards).

     - Options:

       - To measure the offsets, pyfalign sums each cube over wavelength to
         produce an image. The option method="centroid" then simply finds the
         brightest peak and takes a centroid over its half-maximum region while
         method="correlate" (suitable for more nebulous regions) uses AstroPy's
         convolve_fft to cross-correlate each collapsed image with the first
         one at tenth-of-a-pixel resolution. These algorithms are a bit rough
         and ready but are good enough for registering well-sampled cubes as
         long as they don't pick up spurious sources...

       - The llimit option can be used to exclude wavelength planes towards
         the blue end (which can be quite noisy for GMOS) from the summation.
         An hlimit option may be added later, as needed.

     - You don't have to use pyfalign at all. If your WCS values already take
       pointing differences into account accurately (unlikely), you can just
       proceed to the second step, otherwise you can use any method available
       to figure out the offsets between cubes (eg. IRAF imexam) and update
       the CRVAL/CRPIX header keywords manually such that a given feature has
       the same World co-ordinates in every cube (see Greisen et al., 2002,
       for more information on the WCS conventions). Using pyfalign is much
       easier, however, as long as it finds the right peak(s), which the
       second step can help you confirm.

     - Any spectral offsets between cubes are assumed to be reflected
       accurately in the input headers to begin with, as determined by your
       wavelength calibration.

2. Resample all the cubes onto the same grid as the first one (extended as
   necessary to contain the full mosaic) and "co-add" (actually average) them:

     pyfmosaic cube1,cube2,cube3,cube4 cube_avg

     - Any data quality arrays in the input cubes are used to exclude blank
       pixels from the average (with non-rectangular edges in mind), but DQ
       is currently not propagated to the output.

     - Options:

       - Don't change the posangle parameter for now (see notes).

       - If you specify separate=yes, each input cube will be resampled onto
         the same output grid but written to a separate SCI extension in the
         output file instead of averaging the cubes. This helps verify how well
         the registration worked at step 1; once you're happy, repeat with
         separate=no. It also potentially allows co-addition with rejection
         using an external program such as imcombine, but currently pyfmosaic
         doesn't write out a mask to help track the number of cubes overlapping
         at each output pixel.

       - If you specify var=yes, a variance cube will be generated from VAR
         extensions in the input cubes, as long as they are all available. The
         input variance arrays are resampled in the same way as their
         corresponding science extensions, preserving overall S/N
         characteristics, without accounting for covariance. In the unlikely
         event that the cubes you combine have significantly different
         co-ordinate scales, the calculation will currently be incorrect
         because it does not account for scaling to match the first cube (the
         science arrays themselves should be handled properly but I haven't
         really tested that; I think both should be OK with mutual rotations).

Logarithmic wavelength rebinning is currently performed as a separate step,
like so:

  pyflogbin incube outcube

  - This can take 5-10 minutes to run, so please be patient.

  - You should do any mosaicking beforehand, as pyfmosaic needs linear WCS.
    Also, you should not run pyflogbin more than once on the same data, as I
    haven't really checked that it does the right thing with log input.

  - Options:

    - With flux=yes, pyflogbin will scale its output values according to the
      bin size as a function of wavelength, conserving total flux rather than
      flux density. This is probably not that useful since you will have flux
      calibrated your data (eg. in units per Angstrom) beforehand and variance
      propagation is a better way to track S/N, but the option there if you
      need it... Don't use this unless you're confident you know better.

    - With var=yes, a variance cube will be generated from the input VAR
      extension, as described for pyfmosaic. In the absence of covariance
      propagation, the intention is to conserve overall S/N rather than the
      S/N in a single pixel, as the former would otherwise become quite
      awkward to track accurately during subsequent analysis.

  - The output WCS follows the convention described in Greisen et al. (2005),
    equation 5. This is NOT interpreted correctly by core IRAF tasks (except
    close to the reference pixel) but is easy for the user to evaluate and
    preserves the reference increment and range from the input.


Other notes & limitations
-------------------------

o The expected data format is somewhat Gemini-centric (see introduction).

o Specifying a rotation for the output cube (WRT the input) currently produces
  incorrect output, due to some regression that I haven't had time to fix; you
  should therefore leave pyfmosaic's posangle parameter set to the default.

o This version doesn't do any bad pixel rejection. You should clean up your
  data by the time you generate the input cubes or specify separate=yes and
  then use some other program to combine the output extensions with rejection.
  I did briefly work on a pyfmosaic version with a median combine option
  (which of course uses more memory and is non-optimal in terms of S/N) but
  any further effort will be better invested in AstroPy tools (see Plans).

o Pyfalign doesn't work reliably in the presence of uncorrected cosmic rays.

o I'll probably build a log-rebinning option directly into pyfmosaic to avoid
  the extra step but haven't got around to it yet.

o The spline interpolation used here (from scipy.ndimage) is suitable for
  well-sampled data. If your data are undersampled (eg. NIFS), look at
  gemcube/nifcube in the Gemini IRAF package (you may have to do some work to
  construct the WCS appropriately).


Plans
-----

This code is built on a >1 decade-old experimental data access class/module
that was my first significant venture into Python. The latter includes some
useful code for manipulating simple WCS but it's also a bit incomplete/quirky
and lacking in documentation. Today, it would make more sense to start with
Gemini's astrodata class, on top of AstroPy's nddata and gwcs (or possibly
SunPy's ndcube). Eventually the functionality provided here should be rebuilt
on top of those things. I therefore don't intend to support my libraries for
developing other applications but I'm open to minor improvements in
functionality (my available time being very limited) in the meantime.


Changes
-------

0.12

o Slice with a tuple rather than a list, which is deprecated in recent numpy
  versions.


0.11

o Replace old stsci.imagestats dependency with astropy.stats & numpy, to
  support Python 3.8+.

o Issue warnings rather than failing when (re-)writing non-compliant headers.

o If no wavelength axis CTYPE is found in the headers, assume it's the
  outermost (last FITS) axis by default, rather than failing with an obscure
  indexing error. 


0.10

o Add pyfalign option to specify a lower pixel limit when integrating over
  wavelength to produce an image for centring (from Kathleen Labrie).

o Python 3 compatibility changes. Previous versions only work with Python 2.


0.9

o Restructure slightly & add a setup.py, for direct use as a Python module
  (in particular by NDMapper), as an alternative to the PyRAF interface.

o Use astropy.io.fits instead of pyfits.

o Minor fixes to support API changes in NumPy 1.9 and 1.10. Mask division
  now accounts for fractional edge pixels.


0.8.1

o Don't try & fail to calculate variance when propagation disabled.

o Avoid 1-pixel grid copy mismatch due to repeated rounding.

o Fix recently-broken separate+ option, due to re-used variable name.


Acknowledgements
----------------

Supported by the Gemini Observatory, which is operated by the Association of
Universities for Research in Astronomy, Inc., on behalf of the international
Gemini partnership of Argentina, Brazil, Canada, Chile, and the United
States of America.

This research made use of Astropy, a community-developed core Python package
for Astronomy (Astropy Collaboration, 2013). See http://www.astropy.org.

PyRAF is a product of the Space Telescope Science Institute, which is operated
by AURA for NASA.

The SciPy community (http://scipy.org) seems to have no standard
acknowledgement but provides the underlying numerical Python code on which this
is built. In particular, Stefan van der Walt helpfully resolved an ndimage
edge-artifact bug in 2007 that had been contaminating pyfmosaic results. 


James, March 2025.