Package wsatools :: Module Astrometry
[hide private]

Module Astrometry

source code

Astrometric calculations. Routines for positional calculations in the celestial sphere derived from FORTRAN SLALIB of P.T. Wallace. Functions to convert between pixel (X, Y) co-ordinates in the detector plane and sky (RA, Dec) co-ordinates in the celestial sphere, use our PyWCS wrapper to the WCS routines of Starlink's AST library.

Pairing module replacement functions:

   import wsatools.pairing as pairing
   pairing.multiPair(inputList, matchRadius, workPath)
   pairing.pairOff(set0, set1, pFile0, pFile1, pairCriterion)

becomes:

   import wsatools.Astrometry as astro
   astro.pairObjects(inputList, matchRadius, workPath)
   astro.pairObjects([set0, set1], pairCriterion, pFile0, oneWay=True)

Author: R.S. Collins, N.C. Hambly

Organization: WFAU, IfA, University of Edinburgh

Requires: Starlink-AST, Numpy, PyWCS

Classes [hide private]
  Pixel
Pixel position data.
  SkyPos
Sky position data.
Functions [hide private]
float
angDiff(angle1, angle2, minVal=-180, cycleLength=360)
Calculates the angular difference between the given two angles accounting for zero crossings.
source code
float
angularSep(pos1, pos2)
Angle between two points on a sphere.
source code
float
angularSepArr(v1, v2, isEqu=False)
Angle between two arrays of points on a sphere or plane.
source code
tuple(float, float, float, float)
calcGroupCentre(cartPosList, isWeighted=False)
Takes list of positions in radians and returns centre in radians.
source code
list
calcHSpaces(coords, initPos, nCoordInHSDef)
Calculates hSpaces for a set of coordinates
source code
long
calcHtmIndex(ra, dec, level=20)
Calculate an HTM index for the given position at the given level.
source code
dict(int: numpy.array)
calcMosaicProductCorners(requiredMosaics)
Takes a set of required mosaics and calculates a set of positions for the corners.
source code
tuple(float, float)
convCartesianToGalactic(skyVector)
Converts a J2000 (FK5) Cartesian vector to a Galactic coordinate.
source code
tuple(float, float)
convCartesianToSky(skyVector)
Converts Cartesian vector to a spherical coordinate.
source code
numpy.array(numpy.array(float, float, float))
convPixelToCartSky(pixels, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')
Converts from pixel (X,Y) co-ordinates in the detector plane to Cartesian sky co-ordinates in the celestial sphere, for a given detector's observation (e.g.
source code
list(SkyPos)
convPixelToSky(pixels, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')
Converts from pixel (X,Y) co-ordinates in the detector plane to sky (RA, Dec) co-ordinates in the celestial sphere, for a given detector's observation (e.g.
source code
list(Pixel)
convSkyToPixel(positions, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')
Converts from sky (RA, Dec) co-ordinates in the celestial sphere to pixel (X,Y) co-ordinates in the detector plane, for a given detector's observation (e.g.
source code
numpy.array(float, float, float) or numpy.array(numpy.array(float, float, float))
convSkyToCartesian(coords)
Converts spherical coordinate(s) to Cartesian vector(s).
source code
tuple(float, float)
convSkyToGalactic(coords)
Converts B1950 (FK4) sky coordinate to Galactic coordinate.
source code
float, float
convSkyToTangent(pos, axis)
Projection of spherical coordinates onto the tangent plane, also known as the gnomic projection.
source code
float, float
convTangentToSky(pos, axis)
Transform tangent plane coordinates back to spherical given the gnomic projection (standard) and tangent point coordinates.
source code
list(SkyPos)
findCorners(multiframeID, extNum, db)
Calculates the sky (RA, Dec) co-ordinates of the given detector's corners.
source code
 
findTileCentre(offsetCoords, offsetPos, posAngle)
Takes coordinates of pawprint, offsetPos and posAngle.
source code
float
getMedAngle(angles, expRange, minVal=-180, cycleLength=360)
Returns: Median angle, taking into account cyclical nature of measurement.
source code
 
getOffSetsFromPos(offPos)
Get position within tile.
source code
 
getOffSetPos(offX, offY)
Get position within tile.
source code
list(tuple)
matchFrames(frames, raIndex, centrePos, tolerance, extNum=None, extIndex=None)
Returns a list of frames matching an input position.
source code
numpy.array(int)
_matchPairs(sourcesM, sourcesS, decsS, tolerance, cartSep)
Matches given lists of detections and returns a pointing array.
source code
float
normalAngle(radians)
Normalize an angle into the range 0 to 2xPI radians.
source code
 
pairObjects(inputList, matchRadius, outPathName, oneWay=False)
Pairs object positions from given binary files, saving the resulting pointings in binary files.
source code
tuple(float)
_queryWCS(db, multiframeID, extNum)
Looks up WCS values in the database for given observation.
source code
 
convertCoordEso(raDeg, decDeg, outString=False)
Returns correct format of tileRa and tileDec
source code
Variables [hide private]
  __package__ = 'wsatools'
Function Details [hide private]

angDiff(angle1, angle2, minVal=-180, cycleLength=360)

source code 

Calculates the angular difference between the given two angles accounting for zero crossings.

Parameters:
  • angle1 (float) - First angle (default is degrees).
  • angle2 (float) - Second angle (default is degrees).
  • minVal (float) - Value of the angle at the "zero crossing".
  • cycleLength (float) - Difference between minimum and maximum values, i.e. separation of values at "zero crossing".
Returns: float
Angular difference between two angles.

angularSep(pos1, pos2)

source code 

Angle between two points on a sphere.

Parameters:
  • pos1 (tuple(float, float)) - Co-ordinate on celestial sphere in radians of the first point, (RA, Dec), though can be any spherical co-ordinate.
  • pos2 (tuple(float, float)) - Co-ordinate on celestial sphere in radians of the second point, (RA, Dec), though can be any spherical co-ordinate.
Returns: float
Angle in radians (always positive) between the two points.

Note: Based on the SLALIB algorithms dsep() and dsepv() by P.T. Wallace.

angularSepArr(v1, v2, isEqu=False)

source code 

Angle between two arrays of points on a sphere or plane.

Parameters:
  • v1 (numpy.array) - Co-ordinates in radians of the first points, e.g. (RA, Dec).
  • v2 (numpy.array) - Co-ordinates in radians of the second points, e.g. (RA, Dec).
  • isEqu (bool) - True, if equatorial (spherical/celestial) co-ordinates.
Returns: float
Angle in radians (always positive) between the two points.

Note: Based on the SLALIB algorithms dsep() and dsepv() by P.T. Wallace.

calcGroupCentre(cartPosList, isWeighted=False)

source code 

Takes list of positions in radians and returns centre in radians.

Parameters:
  • cartPosList (numpy.array) - Array of RA, Dec co-ordinate positions in radians.
  • isWeighted (bool) - If True, calculate using weights.
Returns: tuple(float, float, float, float)
Centre of list of positions in RA, Dec co-ordinates with errors (radians).

calcHSpaces(coords, initPos, nCoordInHSDef)

source code 

Calculates hSpaces for a set of coordinates

Parameters:
  • coords (list(tuple)) - List of Cartesian coordinates on the sky
  • initPos (list(int)) - Initial positions of each new hspace in coords
  • nCoordInHSDef (int) - Number of coordinates used to define an H-space
Returns: list
list of values, 4 for each half-space. The first 3 are the Cartesian coordinates of the normal vector to the plane and the 4th is the offset of the plane from the origin.

calcHtmIndex(ra, dec, level=20)

source code 

Calculate an HTM index for the given position at the given level.

>>> calcHtmIndex(ra=200.981, dec=-1.06456, level=20)
11030377854233L
Parameters:
  • ra (float) - Right ascension in degrees.
  • dec (float) - Declination in degrees.
  • level (int) - HTM level to calculate down to.
Returns: long
HTM index of given position.

calcMosaicProductCorners(requiredMosaics)

source code 

Takes a set of required mosaics and calculates a set of positions for the corners.

Parameters:
  • requiredMosaics (list(namedtuple)) - List of namedtuples of RequiredMosaics attributes
Returns: dict(int: numpy.array)
Dictionary of 2D numpy arrays of Cartesian coordinates for the corners of the mosaic of the given productID

convCartesianToGalactic(skyVector)

source code 

Converts a J2000 (FK5) Cartesian vector to a Galactic coordinate.

Parameters:
  • skyVector (numpy.array(float, float, float)) - 3-vector of Cartesian coordinates (FK5 J2000).
Returns: tuple(float, float)
Galactic coordinate (l, b) in degrees.

convCartesianToSky(skyVector)

source code 

Converts Cartesian vector to a spherical coordinate. The spherical coordinate is longitude (+ve anticlockwise looking from the positive latitude pole) and latitude. The Cartesian coordinate is right handed, with the x-axis at zero longitude and latitude, and the z-axis at the positive latitude pole.

Parameters:
  • skyVector (tuple(float, float, float)) - 3-vector of right-handed Cartesian coordinates.
Returns: tuple(float, float)
Coordinate on celestial sphere in degrees (RA, Dec), though can be any spherical coordinate.

convPixelToCartSky(pixels, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')

source code 

Converts from pixel (X,Y) co-ordinates in the detector plane to Cartesian sky co-ordinates in the celestial sphere, for a given detector's observation (e.g. multiframeID, extNum). If multiframeID and extNum supplied it looks up WCS values in the database for this observation, otherwise WCS must be supplied. Uses PyWCS to perform the conversion.

Parameters:
  • pixels (list(tuple(float, float))) - List of pixel (X,Y) co-ordinates.
  • multiframeID (int) - MultiframeID of the observation.
  • extNum (int) - Detector extension number.
  • db (DbSession) - Optional database connection to obtain WCS co-ordinates.
  • wcs (tuple(float)) - World Co-ordinate System for detector.
  • projection (str) - Projection, either ZPN or TAN.
Returns: numpy.array(numpy.array(float, float, float))
Array of 3-vectors of right-handed Cartesian sky coordinates ordered by and corresponding to the list of pixel co-ordinates.

To Do: Would be more efficient if PyWCS performed multiple co-ordinate conversions for a given WCS value set, using ctypes.

convPixelToSky(pixels, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')

source code 

Converts from pixel (X,Y) co-ordinates in the detector plane to sky (RA, Dec) co-ordinates in the celestial sphere, for a given detector's observation (e.g. multiframeID, extNum). If multiframeID and extNum supplied it looks up WCS values in the database for this observation, otherwise WCS must be supplied. Uses PyWCS to perform the conversion.

Parameters:
  • pixels (list(tuple(float, float))) - List of pixel (X,Y) co-ordinates.
  • multiframeID (int) - MultiframeID of the observation.
  • extNum (int) - Detector extension number.
  • db (DbSession) - Optional database connection to obtain WCS co-ordinates.
  • wcs (tuple(float)) - World Co-ordinate System for detector.
  • projection (str) - Projection, either ZPN or TAN.
Returns: list(SkyPos)
Sky (RA, Dec) co-ordinates in degrees ordered by and corresponding to the list of pixel co-ordinates.

To Do: Would be more efficient if PyWCS performed multiple co-ordinate conversions for a given WCS value set, using ctypes.

convSkyToPixel(positions, multiframeID=None, extNum=None, db=None, wcs=None, projection='ZPN')

source code 

Converts from sky (RA, Dec) co-ordinates in the celestial sphere to pixel (X,Y) co-ordinates in the detector plane, for a given detector's observation (e.g. multiframeID, extNum). If multiframeID and extNum supplied it looks up WCS values in the database for this observation, otherwise WCS must be supplied. Uses PyWCS to perform the conversion.

Parameters:
  • positions (list(tuple(float, float))) - List of sky (RA, Dec) co-ordinates in degrees.
  • multiframeID (int) - MultiframeID of the observation.
  • extNum (int) - Detector extension number.
  • db (DbSession) - Optional database connection to obtain WCS co-ordinates.
  • wcs (tuple(float)) - World Co-ordinate System for detector.
  • projection (str) - Projection, either ZPN or TAN.
Returns: list(Pixel)
Pixel (X, Y) co-ordinates ordered by and corresponding to the list of sky position co-ordinates.

To Do: Would be more efficient if PyWCS performed multiple co-ordinate conversions for a given WCS value set, using ctypes.

convSkyToCartesian(coords)

source code 

Converts spherical coordinate(s) to Cartesian vector(s). The spherical coordinates are longitude (+ve anticlockwise looking from the positive latitude pole) and latitude. The Cartesian coordinates are right handed, with the x-axis at zero longitude and latitude, and the z-axis at the positive latitude pole.

Parameters:
  • coords (tuple(float, float) or numpy.array(numpy.array(float, float))) - Co-ordinate(s) on celestial sphere in radians (RA, Dec), though can be any spherical co-ordinate.
Returns: numpy.array(float, float, float) or numpy.array(numpy.array(float, float, float))
(Array of) 3-vector(s) of right-handed Cartesian coordinates.

Note: Based on the SLALIB algorithm dcs2c() by P.T. Wallace.

convSkyToGalactic(coords)

source code 

Converts B1950 (FK4) sky coordinate to Galactic coordinate.

Parameters:
  • coords (tuple(float, float)) - Co-ordinate on celestial sphere in radians (RA, Dec).
Returns: tuple(float, float)
Galactic coordinate (l, b) in degrees.

Note: Based on Practical Astronomy With Your Calculator. Section 29 p43.

convSkyToTangent(pos, axis)

source code 

Projection of spherical coordinates onto the tangent plane, also known as the gnomic projection. Function translated from the SLALIB FORTRAN source.

Parameters:
  • pos (tuple(float, float)) - Co-ordinate on celestial sphere in radians to be projected (RA, Dec), though can be any spherical co-ordinate.
  • axis (tuple(float, float)) - Co-ordinate on celestial sphere in radians of the tangent point (RA, Dec), though can be any spherical co-ordinate.
Returns: float, float
Tangent plane co-ordinates (Xi, Eta) in radians.

Note: Based on the SLALIB algorithm ds2tp() by P.T. Wallace.

convTangentToSky(pos, axis)

source code 

Transform tangent plane coordinates back to spherical given the gnomic projection (standard) and tangent point coordinates. Function translated from the SLALIB Fortran source.

Parameters:
  • pos (tuple(float, float)) - Co-ordinate on tangent plane in radians (Xi, Eta).
  • axis (tuple(float, float)) - Co-ordinate on celestial sphere in radians of the tangent point (RA, Dec), though can be any spherical co-ordinate.
Returns: float, float
Co-ordinate on celestial sphere in radians (RA, Dec), though can be any spherical co-ordinate.

Note: Based on the SLALIB algorithm dtp2s() by P.T. Wallace.

findCorners(multiframeID, extNum, db)

source code 

Calculates the sky (RA, Dec) co-ordinates of the given detector's corners.

Parameters:
  • multiframeID (int) - MultiframeID of the observation.
  • extNum (int) - Detector extension number.
  • db (DbSession) - Database connection to use for metadata queries.
Returns: list(SkyPos)
Sky (RA, Dec) co-ordinates in degrees of each corner.

getMedAngle(angles, expRange, minVal=-180, cycleLength=360)

source code 
Returns: float
Median angle, taking into account cyclical nature of measurement.

matchFrames(frames, raIndex, centrePos, tolerance, extNum=None, extIndex=None)

source code 

Returns a list of frames matching an input position. Takes a position, a positional tolerance and a list of frames (or a frame-like list) and return the subset of those matching that position within the specified tolerance.

Method: identifies the first frame that could be a match based on Dec in the (ordered) input list, then cycles through the list checking all possible matches up to the most distant Dec possible. Note that if the RA and Dec come from a frame that is contained within the given list, that frame will appear as a matched frame in the returned list in addition to any other matches. The optional arguments are provided so that if frame matching required that only common extension numbers are matched, and that information is provided in those input arguments, then the function will only match those.

Parameters:
  • frames (list(tuple)) - List of frames to match.
  • raIndex (int) - Index in frame tuple of the RA co-ordinate (Dec is assumed to be index+1).
  • centrePos (tuple(float, float)) - Co-ordinate of central position (RA, Dec) in degrees.
  • tolerance (float) - Positional matching tolerance for frames in a set (degrees).
  • extNum (int) - Optionally specify a particular image extension to be matched.
  • extIndex (int) - If matching by image extension then provide the index of the frame tuple that contains the extension number.
Returns: list(tuple)
A subset of frames that match the given position.

_matchPairs(sourcesM, sourcesS, decsS, tolerance, cartSep)

source code 

Matches given lists of detections and returns a pointing array.

Parameters:
  • sourcesM (list(tuple(float, (int, float, float, float)))) - List of sorted master detection declinations, record number and Cartesian co-ordinates.
  • sourcesS (list(tuple(float, (int, float, float, float)))) - List of sorted secondary detection declinations, record number and Cartesian co-ordinates.
  • decsS (list(float)) - List of sorted secondary detection declinations.
  • tolerance (float) - Pairing criterion in arcseconds with additional 10%.
  • cartSep (float) - Pairing criterion in Cartesian distance.
Returns: numpy.array(int)
Pointing array of detection indices.

Note: Using xrange here instead of range gives major performance benefits.

To Do:
  • Need to try a fully vectorised calculation in NumPy for performance.
  • Should be able to combine all three sub for-loops in this function into one - only doing Cartesian separation calculations when in the declination tolerance range. Also try slicing the list, but need to be certain perf. is good, also can possibly do this as a find type func.

normalAngle(radians)

source code 

Normalize an angle into the range 0 to 2xPI radians.

Parameters:
  • radians (float) - The angle to normalize (in radians).
Returns: float
The angle normalized to the range 0 to 2xPI.

pairObjects(inputList, matchRadius, outPathName, oneWay=False)

source code 

Pairs object positions from given binary files, saving the resulting pointings in binary files.

Parameters:
  • inputList - List of detection data sets, one for each passband, as a tuple of (filePath, (raIndex, raSize, raUnit), (decIndex, decSize, decUnit), recordLength)
  • inputList - list(tuple)
  • matchRadius (float) - Pairing separation criterion in arcseconds.
  • outPathName (str) - Path to write the pointing files too. File name must be included if oneWay is True.
  • oneWay (bool) - If True, don't prepare a reverse pointing file.
To Do:
  • Make completely schema driven for both data types and units, basically the caller gets raIndex/raSize/raUnit from schema file and also passes raType, e.g. 'd' in this case.
  • Can be made even more Pythonic and faster if we include the merge table in with this function, avoiding the need for index look-ups and reading/writing to binary files... assuming memory req. are OK.

Note: NumPy is useful and fast for the pp array, but slow for the Cartesian vectors since they are stored as a list. If the calculations were fully vectorised then this wouldn't be an issue.

_queryWCS(db, multiframeID, extNum)

source code 

Looks up WCS values in the database for given observation.

Parameters:
  • db (DbSession) - Database connection to obtain WCS co-ordinates.
  • multiframeID (int) - MultiframeID of the observation.
  • extNum (int) - Detector extension number.
Returns: tuple(float)
WCS co-ordinates tuple.