#
# wcs.py -- calculations based on world coordinate system.
#
# This is open-source software licensed under a BSD license.
# Please see the file LICENSE.txt for details.
#
"""This module handles calculations based on world coordinate system."""
import math
from collections import OrderedDict
import numpy as np
from ginga.misc import Bunch
__all__ = ['hmsToDeg', 'dmsToDeg', 'decTimeToDeg', 'degToHms', 'degToDms',
'arcsecToDeg', 'hmsStrToDeg', 'dmsStrToDeg',
'trans_coeff', 'eqToEq2000',
'get_xy_rotation_and_scale', 'get_rotation_and_scale',
'get_relative_orientation', 'simple_wcs', 'deg2fmt', 'dispos',
'deltaStarsRaDecDeg1', 'deltaStarsRaDecDeg2', 'get_starsep_RaDecDeg',
'add_offset_radec', 'get_RaDecOffsets', 'lon_to_deg', 'lat_to_deg',
'ra_deg_to_str', 'dec_deg_to_str',
]
class WCSError(Exception):
pass
degPerHMSHour = 15.0 # 360/24
degPerHMSMin = 0.25 # 360.0/24/60
degPerHMSSec = 1.0 / 240.0 # 360.0/24/60/60
degPerDmsMin = 1.0 / 60.0
degPerDmsSec = 1.0 / 3600.0
HMSHourPerDeg = 1.0 / 15.0
HMSMinPerDeg = 4.0
HMSSecPerDeg = 240.0
[docs]
def hmsToDeg(h, m, s):
"""Convert RA hours, minutes, seconds into an angle in degrees."""
return h * degPerHMSHour + m * degPerHMSMin + s * degPerHMSSec
[docs]
def dmsToDeg(sign, deg, min, sec):
"""Convert dec sign, degrees, minutes, seconds into a signed angle in
degrees."""
return sign * (deg + min * degPerDmsMin + sec * degPerDmsSec)
[docs]
def decTimeToDeg(sign_sym, deg, min, sec):
"""Convert dec sign, degrees, minutes, seconds into a signed angle in
degrees.
``sign_sym`` may represent negative as either '-' or numeric -1."""
if sign_sym == -1 or sign_sym == '-':
sign = -1
else:
sign = 1
return dmsToDeg(sign, deg, min, sec)
[docs]
def degToHms(ra):
"""Converts the ra (in degrees) to HMS three tuple.
H and M are in integer and the S part is in float.
"""
assert (ra >= 0.0), WCSError("RA (%f) is negative" % (ra))
assert ra < 360.0, WCSError("RA (%f) > 360.0" % (ra))
rah = ra / degPerHMSHour
ramin = (ra % degPerHMSHour) * HMSMinPerDeg
rasec = (ra % degPerHMSMin) * HMSSecPerDeg
return (int(rah), int(ramin), rasec)
[docs]
def degToDms(dec, isLatitude=True):
"""Convert the dec, in degrees, to an (sign,D,M,S) tuple.
D and M are integer, and sign and S are float.
"""
if isLatitude:
assert dec <= 90, WCSError("DEC (%f) > 90.0" % (dec))
assert dec >= -90, WCSError("DEC (%f) < -90.0" % (dec))
if dec < 0.0:
sign = -1.0
else:
sign = 1.0
dec = dec * sign
mnt, sec = divmod(dec * 3600, 60)
deg, mnt = divmod(mnt, 60)
return (int(sign), int(deg), int(mnt), sec)
[docs]
def arcsecToDeg(arcsec):
"""Convert numeric arcseconds (aka DMS seconds) to degrees of arc.
"""
return arcsec * degPerDmsSec
[docs]
def hmsStrToDeg(ra):
"""Convert a string representation of RA into a float in degrees."""
hour, min, sec = ra.split(':')
ra_deg = hmsToDeg(int(hour), int(min), float(sec))
return ra_deg
[docs]
def dmsStrToDeg(dec):
"""Convert a string representation of DEC into a float in degrees."""
sign_deg, min, sec = dec.split(':')
sign = sign_deg[0:1]
if sign not in ('+', '-'):
sign = '+'
deg = sign_deg
else:
deg = sign_deg[1:]
dec_deg = decTimeToDeg(sign, int(deg), int(min), float(sec))
return dec_deg
[docs]
def ra_deg_to_str(ra_deg, precision=3, format='%02d:%02d:%02d.%03d'):
if ra_deg > 360.0:
ra_deg = math.fmod(ra_deg, 360.0)
ra_hour, ra_min, ra_sec = degToHms(ra_deg)
frac_sec, ra_sec = math.modf(ra_sec)
frac_sec = int(frac_sec * 10 ** precision)
return format % (ra_hour, ra_min, ra_sec, frac_sec)
[docs]
def dec_deg_to_str(dec_deg, precision=2, format='%s%02d:%02d:%02d.%02d'):
sign, dec_degree, dec_min, dec_sec = degToDms(dec_deg)
if sign > 0:
sign_sym = '+'
else:
sign_sym = '-'
frac_sec, dec_sec = math.modf(dec_sec)
frac_s = int(frac_sec * 10 ** precision)
return format % (sign_sym, int(dec_degree), int(dec_min), dec_sec,
frac_sec)
[docs]
def trans_coeff(eq, x, y, z):
"""This function is provided by MOKA2 Development Team (1996.xx.xx)
and used in SOSS system."""
tt = (eq - 2000.0) / 100.0
zeta = 2306.2181 * tt + 0.30188 * tt * tt + 0.017998 * tt * tt * tt
zetto = 2306.2181 * tt + 1.09468 * tt * tt + 0.018203 * tt * tt * tt
theta = 2004.3109 * tt - 0.42665 * tt * tt - 0.041833 * tt * tt * tt
zeta = math.radians(zeta) / 3600.0
zetto = math.radians(zetto) / 3600.0
theta = math.radians(theta) / 3600.0
p11 = (math.cos(zeta) * math.cos(theta) * math.cos(zetto) -
math.sin(zeta) * math.sin(zetto))
p12 = (-math.sin(zeta) * math.cos(theta) * math.cos(zetto) -
math.cos(zeta) * math.sin(zetto))
p13 = -math.sin(theta) * math.cos(zetto)
p21 = (math.cos(zeta) * math.cos(theta) * math.sin(zetto) +
math.sin(zeta) * math.cos(zetto))
p22 = (-math.sin(zeta) * math.cos(theta) * math.sin(zetto) +
math.cos(zeta) * math.cos(zetto))
p23 = -math.sin(theta) * math.sin(zetto)
p31 = math.cos(zeta) * math.sin(theta)
p32 = -math.sin(zeta) * math.sin(theta)
p33 = math.cos(theta)
return (p11, p12, p13, p21, p22, p23, p31, p32, p33)
[docs]
def eqToEq2000(ra_deg, dec_deg, eq):
"""Convert Eq to Eq 2000."""
ra_rad = math.radians(ra_deg)
dec_rad = math.radians(dec_deg)
x = math.cos(dec_rad) * math.cos(ra_rad)
y = math.cos(dec_rad) * math.sin(ra_rad)
z = math.sin(dec_rad)
p11, p12, p13, p21, p22, p23, p31, p32, p33 = trans_coeff(eq, x, y, z)
x0 = p11 * x + p21 * y + p31 * z
y0 = p12 * x + p22 * y + p32 * z
z0 = p13 * x + p23 * y + p33 * z
new_dec = math.asin(z0)
if x0 == 0.0:
new_ra = math.pi / 2.0
else:
new_ra = math.atan(y0 / x0)
if ((y0 * math.cos(new_dec) > 0.0 and x0 * math.cos(new_dec) <= 0.0) or
(y0 * math.cos(new_dec) <= 0.0 and x0 * math.cos(new_dec) < 0.0)):
new_ra += math.pi
elif new_ra < 0.0:
new_ra += 2.0 * math.pi
#new_ra = new_ra * 12.0 * 3600.0 / math.pi
new_ra_deg = new_ra * 12.0 / math.pi * 15.0
#new_dec = new_dec * 180.0 * 3600.0 / math.pi
new_dec_deg = new_dec * 180.0 / math.pi
return (new_ra_deg, new_dec_deg)
[docs]
def get_xy_rotation_and_scale(header):
"""
CREDIT: See IDL code at
https://idlastro.gsfc.nasa.gov/ftp/pro/astrom/getrot.pro
"""
def calc_from_cd(cd1_1, cd1_2, cd2_1, cd2_2):
# TODO: Check if first coordinate in CTYPE is latitude
# if (ctype EQ 'DEC-') or (strmid(ctype, 1) EQ 'LAT') then $
# cd = reverse(cd,1)
det = cd1_1 * cd2_2 - cd1_2 * cd2_1
if det < 0:
sgn = -1
else:
sgn = 1
## if det > 0:
## raise ValueError("Astrometry is for a right-handed coordinate system")
if (cd2_1 == 0.0) or (cd1_2 == 0.0):
# Unrotated coordinates?
xrot = 0.0
yrot = 0.0
cdelt1 = cd1_1
cdelt2 = cd2_2
else:
xrot = math.atan2(sgn * cd1_2, sgn * cd1_1)
yrot = math.atan2(-cd2_1, cd2_2)
cdelt1 = sgn * math.sqrt(cd1_1**2 + cd1_2**2)
cdelt2 = math.sqrt(cd1_1**2 + cd2_1**2)
return xrot, yrot, cdelt1, cdelt2
def calc_from_crota():
try:
crota1 = float(header['CROTA1'])
xrot = crota1
except KeyError:
xrot = None
try:
crota2 = float(header['CROTA2'])
yrot = crota2
except KeyError:
yrot = 0.0
if xrot is None:
xrot = yrot
cdelt1 = float(header.get('CDELT1', 1.0))
cdelt2 = float(header.get('CDELT2', 1.0))
return xrot, yrot, cdelt1, cdelt2
# 1st, check for presence of PC matrix
try:
pc1_1 = header['PC1_1']
pc1_2 = header['PC1_2']
pc2_1 = header['PC2_1']
pc2_2 = header['PC2_2']
cdelt1 = float(header['CDELT1'])
cdelt2 = float(header['CDELT2'])
cd1_1, cd1_2 = cdelt1 * pc1_1, cdelt1 * pc1_2
cd2_1, cd2_2 = cdelt2 * pc2_1, cdelt2 * pc2_2
xrot, yrot, cdelt1, cdelt2 = calc_from_cd(cd1_1, cd1_2,
cd2_1, cd2_2)
except KeyError:
# 2nd, check for presence of CD matrix
try:
cd1_1 = header['CD1_1']
cd1_2 = header['CD1_2']
cd2_1 = header['CD2_1']
cd2_2 = header['CD2_2']
xrot, yrot, cdelt1, cdelt2 = calc_from_cd(cd1_1, cd1_2,
cd2_1, cd2_2)
except KeyError:
# 3rd, check for presence of CROTA keyword
# (or default is north=up)
xrot, yrot, cdelt1, cdelt2 = calc_from_crota()
xrot, yrot = math.degrees(xrot), math.degrees(yrot)
return ((xrot, yrot), (cdelt1, cdelt2))
[docs]
def get_rotation_and_scale(header, skew_threshold=0.001):
"""Calculate rotation and CDELT."""
((xrot, yrot),
(cdelt1, cdelt2)) = get_xy_rotation_and_scale(header)
if math.fabs(xrot) - math.fabs(yrot) > skew_threshold:
raise ValueError("Skew detected: xrot=%.4f yrot=%.4f" % (
xrot, yrot))
rot = yrot
lonpole = float(header.get('LONPOLE', 180.0))
if lonpole != 180.0:
rot += 180.0 - lonpole
return (rot, cdelt1, cdelt2)
[docs]
def get_relative_orientation(image, ref_image):
"""Computes the relative orientation and scale of an image to a reference
image.
Parameters
----------
image
AstroImage based object
ref_image
AstroImage based object
Returns
-------
result
Bunch object containing the relative scale in x and y
and the relative rotation in degrees.
"""
# Get reference image rotation and scale
header = ref_image.get_header()
((xrot_ref, yrot_ref),
(cdelt1_ref, cdelt2_ref)) = get_xy_rotation_and_scale(header)
scale_x, scale_y = math.fabs(cdelt1_ref), math.fabs(cdelt2_ref)
# Get rotation and scale of image
header = image.get_header()
((xrot, yrot),
(cdelt1, cdelt2)) = get_xy_rotation_and_scale(header)
# Determine relative scale of this image to the reference
rscale_x = math.fabs(cdelt1) / scale_x
rscale_y = math.fabs(cdelt2) / scale_y
# Figure out rotation relative to our orientation
rrot_dx, rrot_dy = xrot - xrot_ref, yrot - yrot_ref
# Choose Y rotation as default
rrot_deg = rrot_dy
res = Bunch.Bunch(rscale_x=rscale_x, rscale_y=rscale_y,
rrot_deg=rrot_deg)
return res
[docs]
def simple_wcs(px_x, px_y, ra_deg, dec_deg, px_scale_deg_px, rot_deg,
cdbase=[1, 1]):
"""Calculate a set of WCS keywords for a 2D simple instrument FITS
file with a 'standard' RA/DEC pixel projection.
Parameters
----------
px_x
(ZERO-based) reference pixel of field in X
(usually center of field)
px_y
(ZERO-based) reference pixel of field in Y
(usually center of field)
ra_deg
RA (in deg) for the reference point
dec_deg
DEC (in deg) for the reference point
px_scale_deg_px
Pixel scale (deg/pixel); can be a tuple for different x,y scales
rot_deg
Rotation angle of the field (in deg)
cdbase
CD base
Returns
-------
res : dict
Ordered dictionary object containing WCS headers.
"""
# center of the projection
crpix = (px_x + 1, px_y + 1) # pixel position (WCS is 1 based)
crval = (ra_deg, dec_deg) # RA, Dec (degrees)
# image scale in deg/pix
cdelt = np.array(cdbase) * np.array(px_scale_deg_px)
# Create rotation matrix for position angle of north (radians E of N)
rot_rad = np.radians(rot_deg)
cpa = np.cos(rot_rad)
spa = np.sin(rot_rad)
# a) clockwise rotation
pc = np.array([[cpa, spa], [-spa, cpa]])
# b) counter clockwise
#pc = np.array([[cpa, -spa], [spa, cpa]])
#cd = pc * cdelt
res = OrderedDict((('CRVAL1', crval[0]),
('CRVAL2', crval[1]),
('CRPIX1', crpix[0]),
('CRPIX2', crpix[1]),
('CUNIT1', 'deg'),
('CUNIT2', 'deg'),
('CTYPE1', 'RA---TAN'),
('CTYPE2', 'DEC--TAN'),
('RADESYS', 'FK5'),
# Either PC + CDELT or CD should appear
# PC + CDELT seems to be the preferred approach
# according to the Calabretta papers
('CDELT1', cdelt[0]),
('CDELT2', cdelt[1]),
('PC1_1', pc[0, 0]),
('PC1_2', pc[0, 1]),
('PC2_1', pc[1, 0]),
('PC2_2', pc[1, 1])
## ('CD1_1' , cd[0, 0]),
## ('CD1_2' , cd[0, 1]),
## ('CD2_1' , cd[1, 0]),
## ('CD2_2' , cd[1, 1]),
))
return res
[docs]
def deg2fmt(ra_deg, dec_deg, format):
"""Format coordinates."""
rhr, rmn, rsec = degToHms(ra_deg)
dsgn, ddeg, dmn, dsec = degToDms(dec_deg)
if format == 'hms':
return rhr, rmn, rsec, dsgn, ddeg, dmn, dsec
elif format == 'str':
frac_sec, rsec = math.modf(rsec)
frac_sec = int(frac_sec * 1000)
ra_txt = '%d:%02d:%02d.%03d' % (rhr, rmn, rsec, frac_sec)
if dsgn < 0:
dsgn = '-'
else:
dsgn = '+'
frac_sec, dsec = math.modf(dsec)
frac_s = int(frac_sec * 100)
dec_txt = '%s%d:%02d:%02d.%02d' % (dsgn, ddeg, dmn, dsec, frac_sec)
return ra_txt, dec_txt
[docs]
def dispos(dra0, decd0, dra, decd):
"""Compute distance and position angle solving a spherical
triangle (no approximations).
Source/credit: Skycat
Author: A.P. Martinez
Parameters
----------
dra0 : float
Center RA in decimal degrees.
decd0 : float
Center DEC in decimal degrees.
dra : float
Point RA in decimal degrees.
decd : float
Point DEC in decimal degrees.
Returns
-------
phi : float
Phi in degrees (East of North).
dist : float
Distance in arcmin.
"""
radian = 180.0 / math.pi
# coord transformed in radians
alf = dra / radian
alf0 = dra0 / radian
del_ = decd / radian
del0 = decd0 / radian
sd0 = math.sin(del0)
sd = math.sin(del_)
cd0 = math.cos(del0)
cd = math.cos(del_)
cosda = math.cos(alf - alf0)
cosd = sd0 * sd + cd0 * cd * cosda
dist = math.acos(cosd)
phi = 0.0
if dist > 0.0000004:
sind = math.sin(dist)
cospa = (sd * cd0 - cd * sd0 * cosda) / sind
#if cospa > 1.0:
# cospa=1.0
if math.fabs(cospa) > 1.0:
# 2005-06-02: fix from awicenec@eso.org
cospa = cospa / math.fabs(cospa)
sinpa = cd * math.sin(alf - alf0) / sind
phi = math.acos(cospa) * radian
if sinpa < 0.0:
phi = 360.0 - phi
dist *= radian
dist *= 60.0
if decd0 == 90.0:
phi = 180.0
if decd0 == -90.0:
phi = 0.0
return (phi, dist)
[docs]
def deltaStarsRaDecDeg1(ra1_deg, dec1_deg, ra2_deg, dec2_deg):
"""Spherical triangulation."""
phi, dist = dispos(ra1_deg, dec1_deg, ra2_deg, dec2_deg)
return arcsecToDeg(dist * 60.0)
[docs]
def deltaStarsRaDecDeg2(ra1_deg, dec1_deg, ra2_deg, dec2_deg):
ra1_rad = math.radians(ra1_deg)
dec1_rad = math.radians(dec1_deg)
ra2_rad = math.radians(ra2_deg)
dec2_rad = math.radians(dec2_deg)
sep_rad = math.acos(math.cos(90.0 - dec1_rad) * math.cos(90.0 - dec2_rad) +
math.sin(90.0 - dec1_rad) * math.sin(90.0 - dec2_rad) *
math.cos(ra1_rad - ra2_rad))
res = math.degrees(sep_rad)
return res
# Use spherical triangulation
deltaStarsRaDecDeg = deltaStarsRaDecDeg1
"""Use spherical triangulation."""
[docs]
def get_starsep_RaDecDeg(ra1_deg, dec1_deg, ra2_deg, dec2_deg):
"""Calculate separation."""
sep = deltaStarsRaDecDeg(ra1_deg, dec1_deg, ra2_deg, dec2_deg)
sgn, deg, mn, sec = degToDms(sep)
if deg != 0:
txt = '%02d:%02d:%06.3f' % (deg, mn, sec)
else:
txt = '%02d:%06.3f' % (mn, sec)
return txt
[docs]
def add_offset_radec(ra_deg, dec_deg, delta_deg_ra, delta_deg_dec):
"""
Algorithm to compute RA/Dec from RA/Dec base position plus tangent
plane offsets.
"""
# To radians
x = math.radians(delta_deg_ra)
y = math.radians(delta_deg_dec)
raz = math.radians(ra_deg)
decz = math.radians(dec_deg)
sdecz = math.sin(decz)
cdecz = math.cos(decz)
d = cdecz - y * sdecz
ra2 = math.atan2(x, d) + raz
# Normalize ra into the range 0 to 2*pi
twopi = math.pi * 2
ra2 = math.fmod(ra2, twopi)
if ra2 < 0.0:
ra2 += twopi
dec2 = math.atan2(sdecz + y * cdecz, math.sqrt(x * x + d * d))
# back to degrees
ra2_deg = math.degrees(ra2)
dec2_deg = math.degrees(dec2)
return (ra2_deg, dec2_deg)
[docs]
def get_RaDecOffsets(ra1_deg, dec1_deg, ra2_deg, dec2_deg):
"""Calculate offset."""
delta_ra_deg = ra1_deg - ra2_deg
adj = math.cos(math.radians(dec2_deg))
if delta_ra_deg > 180.0:
delta_ra_deg = (delta_ra_deg - 360.0) * adj
elif delta_ra_deg < -180.0:
delta_ra_deg = (delta_ra_deg + 360.0) * adj
else:
delta_ra_deg *= adj
delta_dec_deg = dec1_deg - dec2_deg
return (delta_ra_deg, delta_dec_deg)
def calc_midpoint_radec(ra1, dec1, ra2, dec2):
delta_ra = ra2 - ra1
bx = np.cos(dec2) * np.cos(delta_ra)
by = np.cos(dec2) * np.sin(delta_ra)
dec_term1 = np.sin(dec1) + np.sin(dec2)
dec_term2 = np.sqrt(np.power((np.cos(dec1) + bx), 2) +
np.power(by, 2))
dec_mid = np.arctan2(dec_term1, dec_term2)
ra_mid = ra1 + np.arctan2(by, np.cos(dec1) + bx)
return (ra_mid, dec_mid)
[docs]
def lon_to_deg(lon):
"""Convert longitude to degrees."""
if isinstance(lon, str) and (':' in lon):
# TODO: handle other coordinate systems
lon_deg = hmsStrToDeg(lon)
else:
lon_deg = float(lon)
return lon_deg
[docs]
def lat_to_deg(lat):
"""Convert latitude to degrees."""
if isinstance(lat, str) and (':' in lat):
# TODO: handle other coordinate systems
lat_deg = dmsStrToDeg(lat)
else:
lat_deg = float(lat)
return lat_deg
def get_ruler_distances(image, p1, p2):
"""Get the distance calculated between two points. A Bunch of
results is returned, containing pixel values and distance values
if the image contains a valid WCS.
"""
x1, y1 = p1[:2]
x2, y2 = p2[:2]
dx, dy = x2 - x1, y2 - y1
res = Bunch.Bunch(x1=x1, y1=y1, x2=x2, y2=y2,
theta=np.arctan2(y2 - y1, x2 - x1),
dx_pix=dx, dy_pix=dy,
dh_pix=np.sqrt(dx**2 + dy**2),
ra_org=None, dec_org=None,
ra_dst=None, dec_dst=None,
ra_heel=None, dec_heel=None,
dx_deg=None, dy_deg=None, dh_deg=None)
if image is not None and hasattr(image, 'wcs') and image.wcs is not None:
# Calculate RA and DEC for the three points
try:
# origination point
ra_org, dec_org = image.pixtoradec(x1, y1)
res.ra_org, res.dec_org = ra_org, dec_org
# destination point
ra_dst, dec_dst = image.pixtoradec(x2, y2)
res.ra_dst, res.dec_dst = ra_dst, dec_dst
# "heel" point making a right triangle
ra_heel, dec_heel = image.pixtoradec(x2, y1)
res.ra_heel, res.dec_heel = ra_heel, dec_heel
res.dh_deg = deltaStarsRaDecDeg(ra_org, dec_org,
ra_dst, dec_dst)
res.dx_deg = deltaStarsRaDecDeg(ra_org, dec_org,
ra_heel, dec_heel)
res.dy_deg = deltaStarsRaDecDeg(ra_heel, dec_heel,
ra_dst, dec_dst)
except Exception as e:
pass
return res
def get_starsep_XY(image, x1, y1, x2, y2):
# source point
ra_org, dec_org = image.pixtoradec(x1, y1)
# destination point
ra_dst, dec_dst = image.pixtoradec(x2, y2)
return get_starsep_RaDecDeg(ra_org, dec_org, ra_dst, dec_dst)
def calc_radius_xy(image, x, y, radius_deg):
"""Calculate a radius (in pixels) from the point (x, y) to a circle
defined by radius in degrees.
"""
# calculate ra/dec of x,y pixel
ra_deg, dec_deg = image.pixtoradec(x, y)
# Calculate position 1 degree from the given one
# NOTE: this needs to add in DEC, not RA
ra2_deg, dec2_deg = add_offset_radec(ra_deg, dec_deg,
0.0, 1.0)
# Calculate the length of this segment--it is pixels/deg
x2, y2 = image.radectopix(ra2_deg, dec2_deg)
px_per_deg_e = math.sqrt(math.fabs(x2 - x)**2 + math.fabs(y2 - y)**2)
# calculate radius based on desired radius_deg
radius_px = px_per_deg_e * radius_deg
return radius_px
def calc_radius_deg2pix(image, ra_deg, dec_deg, delta_deg,
equinox=None):
x, y = image.radectopix(ra_deg, dec_deg, equinox=equinox)
return calc_radius_xy(image, x, y, delta_deg)
def add_offset_xy(image, x, y, delta_deg_x, delta_deg_y):
# calculate ra/dec of x,y pixel
ra_deg, dec_deg = image.pixtoradec(x, y)
# add offsets
ra2_deg, dec2_deg = add_offset_radec(ra_deg, dec_deg,
delta_deg_x, delta_deg_y)
# then back to new pixel coords
x2, y2 = image.radectopix(ra2_deg, dec2_deg)
return (x2, y2)
def calc_radius_center(image, delta_deg):
return calc_radius_xy(image,
float(image.width / 2.0),
float(image.height / 2.0),
delta_deg)
def calc_compass(image, x, y, len_deg_e, len_deg_n):
# Get east and north coordinates
xe, ye = add_offset_xy(image, x, y, len_deg_e, 0.0)
xn, yn = add_offset_xy(image, x, y, 0.0, len_deg_n)
return (x, y, xn, yn, xe, ye)
def calc_compass_radius(image, x, y, radius_px):
# Define an angular length we can use to determine the pixel scale
# along east and north - we want to make sure we use a small length so
# that the offset points still fall inside the image in case they have
# a bounding box set.
delta = 0.1 / 3600 # 0.1 arcsec
xe, ye = add_offset_xy(image, x, y, delta, 0.0)
xn, yn = add_offset_xy(image, x, y, 0.0, delta)
# now calculate the length in pixels of those arcs
# (planar geometry is good enough here)
px_per_delta_e = math.sqrt(math.fabs(ye - y)**2 + math.fabs(xe - x)**2)
px_per_delta_n = math.sqrt(math.fabs(yn - y)**2 + math.fabs(xn - x)**2)
# now calculate the arm length in degrees for each arm
# (this produces same-length arms)
len_deg_e = radius_px / px_per_delta_e * delta
len_deg_n = radius_px / px_per_delta_n * delta
return calc_compass(image, x, y, len_deg_e, len_deg_n)
def calc_compass_center(image):
# calculate center of data
x = float(image.width) / 2.0
y = float(image.height) / 2.0
# radius we want the arms to be (approx 1/4 the smallest dimension)
radius_px = float(min(image.width, image.height)) / 4.0
return calc_compass_radius(image, x, y, radius_px)