CoordinateUtil.h

Go to the documentation of this file.

//# CoordinateUtils.h: static functions dealing with coordinates
//# Copyright (C) 1997,1998,1999,2000,2001,2002,2004
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//# $Id: CoordinateUtil.h 20652 2009-07-06 05:04:32Z Malte.Marquarding $

#ifndef COORDINATES_COORDINATEUTIL_H
#define COORDINATES_COORDINATEUTIL_H


#include <casa/aips.h>
#include <casa/Arrays/Vector.h>
#include <measures/Measures/Stokes.h>
#include <coordinates/Coordinates/Coordinate.h>

#include <measures/Measures/MDirection.h>        //# For enums
#include <measures/Measures/MFrequency.h>
#include <measures/Measures/MeasConvert.h>

namespace casa { //# NAMESPACE CASA - BEGIN

class CoordinateSystem;
class DirectionCoordinate;
class ObsInfo;
class String;
class LogIO;
class MEpoch;
class MPosition;
class Unit;


// <summary>Functions for creating default CoordinateSystems</summary>
// <use visibility=export>

// <reviewed reviewer="" date="" tests="" demos="">
// </reviewed>

// <prerequisite>
//   <li> <linkto class="CoordinateSystem">CoordinateSystem</linkto>
// </prerequisite>
//
// <etymology> 
// CoordinateUtils follows the AIPS++ naming convention for static functions
// that are associated with a class.
// </etymology>
//
// <synopsis>
// This file contains declarations for static functions that manipulate
// coordinate systems. It currently contains functions for:
// <ul>
// <li> Adding default axes to a CoordinateSystem
// <li> Creating a default CoordinateSystem
// <li> Finding specified axes in a CoordinateSystem
// </ul>
// 
// The functions for adding default axes to a CoordinateSystem can add
// either a RA/DEC pair of axes, a Polarisation Axis, or a Spectral Axis to
// a user supplied coordinate system. The default values for these functions
// are:
// <ul>
// <li> <src>addDirAxes</src> this adds a DirectionCoordinate with a
// reference pixel of (0,0) corresponding to an RA/DEC of (0,0) in a
// J2000 reference frame. The pixel increment is 1 arc-minute.
// <li> <src>addIQUVAxis</src> this adds a polarization axis with four
// elements corresponding to the Stokes (I,Q,U,V) components.
// <li> <src>addIAxis</src> this adds a polarization axis with one
// element corresponding to the Stokes I component only
// <li> <src>addFreqAxis</src> this adds a spectral axis with a reference
// frequency of 1.415GHz on channel 0. The channel bandwidth (pixel
// increment) is 1kHz, and the reference frame is the kinematical Local Standard of
// rest (<linkto class="MFrequency">MFrequency</linkto>::LSRK). 
// </ul>
//
// The <src>defaultCoords</src> functions, create from scratch a
// CoordinateSystem using the above described <src>addXXXAxis</src>
// functions to add the required number of dimensions to the
// CoordinateSystem. Only 2, 3 or 4 dimensional coordinate systems can be
// constructed using these functions. The coordinate systems always have
// RA/Dec axes. Three dimensional Systems add a spectral axis and
// four-dimensional systems add an IQUV  polarization axis. An exception
// (AipsError) is thrown if <src>defaultCoords(uInt)</src> is called with a
// parameter that is not 2, 3, or 4. 
//
// The <src>defaultCoordsXX</src> functions return the coordinate system by
// value (which involves a copy of the CoordinateSystem) and hence are not
// as effcient as the <src>addXXXAxis</src> functions.
//
// If the default axes provided by these functions are not quite what is
// required it is possible to use member functions of the 
// <linkto class="CoordinateSystem">CoordinateSystem</linkto>
// and <linkto class="Coordinate">Coordinate</linkto> classes 
// (<linkto class="DirectionCoordinate">DirectionCoordinate</linkto>,
// <linkto class="StokesCoordinate">StokesCoordinate</linkto>,
// <linkto class="SpectralCoordinate">SpectralCoordinate</linkto> etc.)
// to tweak the appropriate parameters of the specified axis.
//
// Now we turn to the functions for finding axes in a CoordinateSystem. With
// a CoordinateSystem object it is not required that the first Coordinate
// axis in the the CoordinateSystem map to the first pixel axis in an
// image. Hence it is necessary to determine which pixel axis corresponds to a
// specified Coordinate and this can be done using these functions. Some
// coordinate types, in particular DirectionCoordinate, usually map to more
// than one pixel axis (DirectionsCoordinates are inherently two-dimensional).
// 
// This group contains declarations for static functions that search
// CoordinateSystem's for a coordinate of the specified type. It returns the
// pixel axis (zero relative) of the specified coordinate type. If the supplied
// Coordinate system does not contain the specified coordinate type the
// returned value is function specific (but usually -1). If the supplied   
// CoordinateSystem contains two or more of the specified coordinateType then
// an exception (AipsError) is thrown.
//
// Finally functions are provided for removing lists of pixel/world axes
// from a CoordinateSystem.
// This process is made a little awkward by the fact that when you
// remove one axis, all the rest shuffle down one, so it is
// provided here.  Generally, one only needs to remove one axis
// (in which case you should use the CoordinateSystem::removeWorldAxis and
// CoordinateSystem::removcePixelAxis functions), but on occaision,
// the multiple need is there.
// </synopsis>
//
// <example>
// I use these functions when creating test images. 
// <srcblock>
// PagedImage(IPosition(4,256,256,4,32), CoordinateUtil::defaultCoords4D(),
//            String("test.image"));
// </srcblock>
// </example>
//
// <example>
// Functions are needed to handle images without specifying a canonical
// coordinate order. For example suppose we want to find the spectral aixs
// of a PagedImage object.
//     
// <srcblock>
//   const Int spectralAxis = CoordinateUtil::findSpectralAxis(image.coordinates());
//   cout << "The spectral axis is of shape " << image.shape()(spectralAxis) << endl;
// </srcblock>
// </example>
//
// <example>
// Here we remove the first and last world axes, and their associated
// pixel axes from a 3D CoordinateSystem.  The reference values and
// reference pixels are used for the replacement values.
//     
// <srcblock>
//   CoordinateSystem cSys = CoordinateUtil::defaultCoords3D();
//   Vector<Int> worldAxes(2);
//   worldAxes(0) = 0; worldAxes(1) = cSys.nWorldAxes()-1;
//   Vector<Double> worldRep;
//   Bool ok = CoordinateUtil::removeAxes(cSys, worldRep, worldAxes, True);
//   cout << "For world axes used " << worldRep << " for replacement" << endl;
// </srcblock>
// </example>
//
//
// <motivation>
// I got fed up writing small functions to create and find coordinates when writing
// test programs involving Images and ComponentModels.
// </motivation>
//
// <thrown>
//    <li> AipsError
// </thrown>
//
// <todo asof="1997/01/23">
//   <li> This code does all I want at the moment
// </todo>

//  <linkfrom anchor=defaultAxes classes="CoordinateSystem">
//      Static functions for creating <here>default</here> coordinate systems
//  </linkfrom>

class CoordinateUtil
{
public: 

// Add a RA/DEC pair of direction axes (ie. a DirectionCoordinate) to the
// user supplied CoordinateSystem. See the synopsis above for the current
// default values.
static void addDirAxes(CoordinateSystem& coords);

// Add a Stokes I,Q,U,V axis to the user supplied CoordinateSystem.
static void addIQUVAxis(CoordinateSystem& coords);

// Add a Stokes I (only) axis to the user supplied CoordinateSystem.
static void addIAxis(CoordinateSystem& coords);

// Add a Stokes axis of length 1 to 4 selected from I,Q,U,V
// E.g. if shape=2 you get IQ.   Returns False if shape
// is not in the range 1 to 4
static Bool addStokesAxis(CoordinateSystem& coords, uInt shape);

// Add Linear axes.  The LinearCoordinate can have > 1 axes (like
// the DirectionCoordinate has 2).  The number of axes is given
// by the length of the names argument.   If you supply a shape,
// it will be used to set the reference pixel to 1/2 the shape.
// If the shape does not have the same number of elements as
// the names variable, the reference pixel will be 0
static void addLinearAxes (CoordinateSystem & coords,
                           const Vector<String>& names,
                           const IPosition& shape);

// Add a spectral axis to the user supplied CoordinateSystem. See the
// synopsis above for the current default values.
static void addFreqAxis(CoordinateSystem& coords);

// Return a 2-dimensional coordinate system with RA/DEC axes only. 
static CoordinateSystem defaultCoords2D();

// Return a 3-dimensional coordinate system with RA/DEC axes and a spectral axis.
static CoordinateSystem defaultCoords3D();

// Return a 4-dimensional coordinate system with RA/DEC axes, an IQUV
// polarisation axis  and a spectral axis.
static CoordinateSystem defaultCoords4D();

// Calls one of the above three functions depending of the arguement. An
// AipsError is thrown if dims is not 2, 3, or 4.
static CoordinateSystem defaultCoords(uInt dims);

// If doLinear=False, Tries to make a standard RA/DEC/Stokes/Frequency CoordinateSystem
// depending upon the shape.   The shape for the Stokes axis
// must be <= 4.   If axis 2 can't be Stokes it will be a Spectral
// axis instead.  AFter the standard types, the rest (if any)
// of the CoordinateSystem consists of LinearCoordinates.
// If doLinear=True, then you just get a linear coordinate system
static CoordinateSystem makeCoordinateSystem(const IPosition& shape,
                                             Bool doLinear=False);

//
// Find which pixel axis in the CoordinateSystem corresponds to the
// SpectralCoordinate. If there is no SpectralCoordinate in the coordinate
// system then return -1.
static Int findSpectralAxis(const CoordinateSystem & coords);

// Find the SpectralCoordinate in the CoordinateSystem, and then
// return the most general description of where it is.  
// If there is no SpectralCoordinate in the CoordinateSystem then return 
// -1 for coordinate.  If the world or pixel axis has been removed,
// return -1 for that value. 
static void findSpectralAxis(Int& pixelAxis, Int& worldAxis, Int& coordinate,
                             const CoordinateSystem & coords);

// Find which pixel axes correspond to the DirectionCoordinate in the 
// supplied coordinate system and return this as a Vector. If there is no 
// DirectionCoordinate in the CoordinateSystem then return a Vector of zero 
// length. Normally the returned Vector will have a length of two.  
// However, if the pixel axis has been removed, then the resultant
// vector will take the value -1 for that axis.
static Vector<Int> findDirectionAxes(const CoordinateSystem & coords);

// Find which pixel axes correspond to the DirectionCoordinate in the supplied coordinate
// system and return the most general description of where it is. If there is 
// no DirectionCoordinate then coordinate is returned with value -1.
// Values of -1 in the returned vectors indicate an axis has been removed.
static void findDirectionAxes(Vector<Int>& pixelAxes, Vector<Int>& worldAxes,
                              Int& coordinate, const CoordinateSystem & coords);

// Find which pixel axis is the polarisation axis in the supplied
// CoordinateSystem and return this. If there is no StokesCoordinate in the
// CoordinateSystem return a negative number. The actual polarisations on the
// returned pixel axis are returned in the whichPols Vector. Each element of
// this Vector is a Stokes::StokesTypes enumerator and the length of the Vector
// is the same as the length of the polarisation axis. If there is no
// polarisation axis the whichPols returns a unit length Vector containing
// Stokes::I
static Int findStokesAxis(Vector<Stokes::StokesTypes>& whichPols, 
                          const CoordinateSystem& coords);

// Find the StokesCoordinate in the CoordinateSystem, and then
// return the most general description of where it is.  
// If there is no StokesCoordinate in the CoordinateSystem then return 
// -1 for coordinate.  If the world or pixel axis has been removed,
// return -1 for that value. 
static void findStokesAxis(Int& pixelAxis, Int& worldAxis, Int& coordinate,
                           const CoordinateSystem & coords);

// Find Coordinate type for this pixel or world axis
// <group>
static Coordinate::Type findPixelAxis (const CoordinateSystem& cSys, Int axis);
static Coordinate::Type findWorldAxis (const CoordinateSystem& cSys, Int axis);
// </group>

// Remove a list of world axes and their associated
// pixel axes from a <src>CoordinateSystem</src>. The list of world
// axes to be removed is derived from a list giving either axes to remove, 
// or axes to keep (controlled by whether <src>remove</src> 
// is <src>True</src> or <src>False</src>.  The replacement values (see functions 
// <src>CoordinateSystem::removeWorldAxis</src>) for the world axes
// can be given.  For the associated pixel axes, the pixel replacement
// coordinate is found by converting the world coordinate 
// to a pixel coordinate. If the length of the replacement value 
// vector is not the number of world axes to be removed then
// the reference values will be used (e.g. use zero length
// vectors).
static Bool removeAxes(CoordinateSystem& cSys,
                       Vector<Double>& worldReplacement,
                       const Vector<Int>& worldAxes,
                       const Bool remove);

// Remove a list of pixel axes but not their associated
// world axes from a <src>CoordinateSystem</src>. 
// The list of pixel axes to be removed is derived from a 
// list giving either axes to remove, 
// or axes to keep (controlled by whether <src>remove</src> 
// is <src>True</src> or <src>False</src>.  The replacement values (see functions 
// <src>CoordinateSystem::removePixelAxis</src>) for the pixel axes
// can be given.  If the length of the replacement value 
// vector is not the number of pixel axes to be removed then
// the reference pixel will be used (e.g. use zero length
// vectors).
static Bool removePixelAxes(CoordinateSystem& cSys,
                            Vector<Double>& pixelReplacement,
                            const Vector<Int>& pixelAxes,
                            const Bool remove);

// Physically (nont just virtually) drop coordinates from the CoordinateSystem
// if all axes are fully removed. For coordinates with axes partially removed
// (world/pixel) preserve that removal state in the output CS.  No effort
// is made to deal in any way with transposed systems.
static Bool dropRemovedAxes (CoordinateSystem& cSysOut, const CoordinateSystem& cSysIn);

// Setup Measures conversion machine for MDirections.
// Returns True if the machine was needed and set.  Returns False
// if the machine was not needed and not set.  
  static Bool makeDirectionMachine(LogIO& os, MDirection::Convert& machine,
                                   const DirectionCoordinate& dirCoordTo,
                                   const DirectionCoordinate& dirCoordFrom,
                                   const ObsInfo& obsTo,
                                   const ObsInfo& obsFrom);

// Setup Measures conversion machines for MFrequencies. 
// Returns False if a trial conversion failed, else returns True.
// There must be both a Direction and a Spectral
// Coordinate in the CoordinateSystem when making the Frequency machine,
// else an exception occurs.
   static Bool makeFrequencyMachine(LogIO& os, MFrequency::Convert& machine, 
                                    Int coordinateTo, Int coordinateFrom, 
                                    const CoordinateSystem& coordsTo, 
                                    const CoordinateSystem& coordsFrom,
                                    const Unit& unit=Unit(String("Hz")));

// Setup Measures conversion machines for MFrequencies.
// Returns False if a trial conversion failed, else returns True.
   static Bool makeFrequencyMachine(LogIO& os, MFrequency::Convert& machine,
                                    MFrequency::Types typeTo, MFrequency::Types typeFrom,
                                    const MDirection& dirTo, const MDirection& dirFrom,
                                    const MEpoch& epochTo, const MEpoch& epochFrom,
                                    const MPosition& posTo, const MPosition& posFrom,
                                    const Unit& unit=Unit(String("Hz")));

// Find the Sky in the CoordinateSystem. Assumes only one DirectionCoordinate.
// <src>pixelAxes</src> and <src>worldAxes</src>  say where
// in the CS the DirectionCoordinate axes are (long then lat).
// Returns False and an error message if it can't find the sky.
   static Bool findSky(String& errorMessage, Int& dirCoord, Vector<Int>& pixelAxes,
                       Vector<Int>& worldAxes, const CoordinateSystem& cSys);
//

// Does the CoordinateSystem hold just the sky ?  Exception if not.
// Returns True if CS pixel axis 0 is the longitude and 1 latitude  
// else returns False
   static Bool isSky (LogIO& os, const CoordinateSystem& cSys);

// Do the specified axes hold the sky ?  Returns False if no DirectionCoordinate
// or if only one axis of the DirectionCoordinate is held or the specified
// pixel axes don't pertain to the DirectionCoordinate.  
   static Bool holdsSky (Bool& holdsOneSkyAxis, const CoordinateSystem& cSys, 
                         Vector<Int> pixelAxes);


// Find the Stokes for the specified pixel. If there is no Stokes in the
// CoordinateSystem, returns Stokes::I
   static Stokes::StokesTypes findSingleStokes (LogIO& os, const CoordinateSystem& cSys,
                                                uInt pixel=0);

// Set the world axis units in the CS to 'deg' for Direction. For Spectral
// set the velocity handling to use 'km/s' units.  Other coordinates
// are not touched.
   static void setNiceAxisLabelUnits(CoordinateSystem& cSys);

// Set world axis units for specific Coordinate.  Returnd False if fails to set units
// with error in cSys.errorMessage().  
   static Bool setCoordinateUnits (CoordinateSystem& cSys, const Vector<String>& units,
                                   uInt which);

// Set a unit for all unremoved world axes in the DirectionCoordinate in the
// CS.  Returns False if fails to set unit with error in cSys.  If no DC
// returns True
   static Bool setDirectionUnit (CoordinateSystem& cSys, const String& unit, Int which=-1);

// Set Direction conversion layer of DirectionCoordinate in CoordinateSystem
// so that pixel<->world go to the specified direction system (a valid
// MDirection::Types string).  Returns False with error if direction
// system invalid.  If no DirectionCoordinate returns True
   static Bool setDirectionConversion (String& errorMsg, CoordinateSystem& cSys,
                                      const String directionSystem);

// Set spectral state of SpectralCoordinate in CoordinateSystem.
// Unit must be consistent with Hz or m/s and the doppler a valid MDoppler string.
// For no change, leave either String empty.
// Returns False if invalid inputs (and CS not changed) and an error message. 
   static Bool setSpectralState (String& errorMsg, CoordinateSystem& cSys, 
                                 const String& unit, const String& doppler);


// Set velocity state of SpectralCoordinate in CoordinateSystem.
// Unit must be consistent m/s and the doppler a valid MDoppler string.
// For no change, leave either String empty.
// Returns False if invalid inputs (and CS not changed) and an error message. 
   static Bool setVelocityState (String& errorMsg, CoordinateSystem& cSys, 
                                 const String& unit, const String& doppler);

// Set Spectral conversion layer of SpectralCoordinate in CoordinateSystem
// so that pixel<->world go to the specified frequency system (a valid
// MFrequency::Types string).  Returns False if frequency system invalid
// or if no DirectionCoordinate or if cant get Date/Epoch
   static Bool setSpectralConversion (String& errorMsg, CoordinateSystem& cSys,
                                      const String frequencySystem);

// Set default format unit and doppler velocity state of SpectralCoordinate in CoordinateSystem.
// Unit can be consistent with Hz or m/s
// Returns False if invalid inputs (and CS not changed) and an error message. 
   static Bool setSpectralFormatting (String& errorMsg, CoordinateSystem& cSys, 
                                      const String& unit, const String& doppler);

// Convert an absolute pixel coordinate to world and format with 
// default Coordinate formatting
// <group>
   static String formatCoordinate(const IPosition& pixel, CoordinateSystem& cSys);
   static String formatCoordinate(const Vector<Double>& pixel, CoordinateSystem& cSys);
// </group>

// Generate axis label String from coordinate. Specify coordinate axis,
// whether world or pixel labels required, whether absolute or
// relative.   For spectral coordinates, doVel says if you want to 
// use the velocity information contained in it to generate the label
   static String axisLabel (const Coordinate& coord, uInt axisInCoordinate=0,
                            Bool doWorld=True, Bool doAbs=True, Bool doVel=False);

  // <group name=Coordinate comparison>
  // Check how the coordinates of this and that compare.
  // The return value tells how they compare.
  // <br>-1: left is subset
  // <br>0: equal 
  // <br>1: left is superset
  // <br>9: invalid (mismatch)
  static Int compareCoordinates (const CoordinateSystem& thisCsys,
                                 const CoordinateSystem& thatCsys);

  // Convert the world axes map given in worldAxes to a pixel axes map.
  static Vector<Int> toPixelAxes (const CoordinateSystem& thisCsys,
                                  const CoordinateSystem& thatCsys,
                                  const Vector<Int>& worldAxes);

  // Check if the axes in the pixel axes map are in ascending order.
  static Bool checkOrder (const Vector<Int>& pixelAxes);

  // Find the new and stretch axes when comparing the old and new
  // coordinates and shapes (helper for ExtendImage).
  static Bool findExtendAxes (IPosition& newAxes,
                              IPosition& stretchAxes,
                              const IPosition& newShape,
                              const IPosition& oldShape,
                              const CoordinateSystem& newCsys,
                              const CoordinateSystem& oldCsys);
  // </group>

  // Fix up Cylindrical parameters in any DirectionCoordinate for when the longitude 
  // is outside of [-180,180] range.  If it returns False, it failed and an error 
  // message is returned as well.  This function should be called on any
  // CS made from an imported image like FITS
  static Bool cylindricalFix (CoordinateSystem& cSys, String& errorMessage, const IPosition& shape);

  // Apply the binning factors to the CS and create a new one reflecting the binning
  // You can optionally throw an exception if factors is non-unit for any Stokes axis
  static CoordinateSystem makeBinnedCoordinateSystem (const IPosition& factors,
                                                      const CoordinateSystem& cSysIn,
                                                      Bool failOnStokes=False);
private:
  // Sets pos to the position found for tel in the database, or
  // raises an exception + error message.
  static void findObservatoryOrRaiseException(LogIO& os, MPosition& pos,
                                              const String& tel);
};


} //# NAMESPACE CASA - END

#endif